2
$\begingroup$

I have a gray scale image of fibres in different orientations. My goal is to mark the area where the fibres have a specific angle and neglect the rest. In the future it should be an automated process with similar images but different width of the areas to mark. That's my input image with the marked area to extract:

Input

At first I removed the vertical white threads which distort the spectrum a lot. I removed them by finding outliers with thresholding and interpolating the "holes" by directional interpolation. The result is here:

Threads Interpolated

Now I did some FFT Filtering and removed all occurring orientations but the -45°, so that it looks as following:

FFT Filtered

Now it looks like an easy task to extract the "non fibre area", but I'm struggling with that. Simple thresholding isn't possible, cause the grayvalues are similar all across the image. I tried with gradient but the result isn't satisfying. If you have any ideas, please let me know.

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ It looks like you need to find the correct algorithm using a heuristically guided stochastic search. In other words -- make wild-ass guesses (the stochastic part), ask yourself if they might work (the heuristic part), try them out, then keep the ones that look best (the search part). Personally, instead of looking for what I don't want and then looking for uniformity, I think I'd look for what I do want, which is lines at a more or less 45 degree angle in the up-right, down-left direction. Filter for that, and then look for lots of contrast. $\endgroup$
    – TimWescott
    Nov 30, 2022 at 17:15

2 Answers 2

Reset to default
2
$\begingroup$

1.- Input image and fix contrast

A=imread('001.jpg');

3 layers RGB to single layer Y:

A1=rgb2gray(A); 

[sz1,sz2]=size(A1)  % sz1:Y  sz2:X
h1=figure(1); imshow(A1); title('input image A')

Choosing best contrast can be done manually with imcontrast or automatically with command imadjust

A2=imadjust(A1,[0/255,11/255]);

figure(2);imshow(A2);title('tune contrast A')

imcontrast(h1)

enter image description here

enter image description here

2.- FFT2(A)

fftA2=fft2(A2);               % image spectrum not centered
fftcA2=fftshift(fftA2);      % image spectrum centered

figure(3); imshow(log(1+abs(fftcA2)),[]); title('centered |FFT(A)|')

enter image description here

A3=imgradientxy(A2,'intermediate');  % sobe (default) | prewitt | central | intermediate

Binarizing

B1=A3;
B1(B1>0)=255;B1(B1<0)=0;
B1=logical(B1);
A3=B1;
figure(4);h2=imshow(A3);title('A3 binarized image')

enter image description here

Centering spectrum

fftA3=fft2(A3);             
fftcA3=fftshift(fftA3);      % centered spectrum

figure(5); imshow(abs(fftcA3),[]); title('centered |FFT(A)|')

enter image description here

surf(|FFT2(A)|)

figure(5); hs11=surf(10*log10(abs(fftcA3))); title('surf centered |FFT(A)|');hs11.EdgeColor='none'`

enter image description here

check min max mean values of |image spectrum|

min(abs(fftcA3(:)))
max(abs(fftcA3(:)))
mean(abs(fftcA3(:)))

% fftcA32_m=abs(fftcA3); % /mean(abs(fftcA3(:)));  % -min(abs(fftcA3/mean(abs(fftcA3(:)))));
fftcA32_m=abs(fftcA3); % /mean(abs(fftcA3(:)));
fftcA32_a=angle(fftcA3);

check splitting mod angle real imag and then combining keeps image

[fftcA32r,fftcA32i]=pol2cart(fftcA32_a,fftcA32_m);
fftcA32_2=fftcA32r+1j*fftcA32i;
A32_2=ifft2(fftcA32_2);
A32_2=real(A32_2); % ifft2 returns imaginary amounts 1e-15 

figure(6);imshow(A32_2,[])

enter image description here

The portions of spectrum that I removed didn't work

3.- nulling FFT(A)<th1

th1=500;

fftcA33_m=fftcA32_m;
fftcA33_m(fftcA33_m<th1)=0;

figure(7); ax7=gca;
hs11=surf(ax7,fftcA33_m);hs11.EdgeColor='none';title('centered |FFT(A)| small values removed');
hold(ax7,'on')

enter image description here

[fftcA33r,fftcA33i]=pol2cart(fftcA32_a,fftcA33_m);
fftcA33=fftcA33r+1j*fftcA33i;
A33=ifft2(fftcA33);
A33=real(A33);

figure(8);imshow(A33,[]); title('effect nulling |FFT(A)|<th1 []')

enter image description here

Now, although grey one can appreciate that there are 2 zone, one stripes at 45° (area of interest) and the larger zone at 135°.

4.- How to find spectrum peaks

[pks,locs,W,P]=findpeaks(fftcA32_m(:),'Threshold',th1/2);
[ylocs,xlocs]=ind2sub(size(fftcA32_m),locs);
plot3(ax7,xlocs,ylocs,pks,'r*')

enter image description here

xylocs=[xlocs ylocs]

Here I tried to zero single peaks (a small square around) or all peaks in each quadrant, different combinations, but no improvement.

5.- Filtering with pattern samples After trying different samples and filters, the smaller H1 the better

H1=[0 0 1;0 1 0;1 0 0];

B2_1=imfilter(A2,H1); 
figure(9); imshow(B2_1);

enter image description here

B2_12=~imbinarize(B2_1);

figure(10);imshow(B2_12);hold on;title('sought area now available')
[ny,nx]=find(B2_12==1);
plot(nx,ny,'r*')

fitobj1=fit(nx,ny,'poly1') 
p1=fitobj1.p1;p2=fitobj1.p2;
plot([1 sz1],[p1+p2 p1*sz1+p2],'b','LineWidth',1.5)

plot(fitobj1,nx,ny,'g-')

enter image description here

6.- Regression line

x0=[1:1:-sz1]';
y0=(p1*x0+p2);
7.- Calculating point-to-line distances
D=[];
a=-p1;b=1;c=-p2;r=1/(a^2+b^2)^.5;
for k1=1:1:numel(nx)
    D=[D; abs(a*nx(k1)+b*ny(k1)+c)/r];
end

figure(11);hh1=histogram(D(:,1),numel(nx));grid on;title('histogram distances to regression line')

enter image description here

point-to-line distances sample

D([1:20])

8.- Standard deviation

sgm=(sum((D(:,1)-mean(D(:,1))).^2)/(numel(nx))).^.5;  
% same as
std(D(:,1));

D2 1st column is normalized distances

D2=[D/sgm D nx ny];

D3=sortrows(D2,1,'descend');

D3([1:20],:)

95% confidence interval : 2 sigma up 2 sigma down removing all points outside this interval

D3(D3(:,1)>2,:)=[];

now all remaining points fall within +2*sgm -2*sgm distance to the regression line,sample

D3([1:20],:)

Removing outliers in nx ny

nx2=D3(:,3);ny2=D3(:,4);

9.- Fitting line

fitobj2=fit(nx2,ny2,'poly1'); % ,'Startpoint',[1 1],'Exclude',outL)

p12=fitobj2.p1;p22=fitobj2.p2;

figure(12);imshow(A);hold on
% [ny,nx]=find(B2_12==1);
plot(nx,ny,'r*')
hold on
plot(nx2,ny2,'go')

plot([1 sz1],[p1+p2 p1*sz1+p2],'r','LineWidth',1.5)
plot([1 sz1],[p1+p2 p1*sz1+p2]+sgm,'r','LineWidth',3)
plot([1 sz1],[p1+p2 p1*sz1+p2]-sgm,'r','LineWidth',3)
plot([1 sz1],[p12+p22 p12*sz1+p22]+sgm,'b','LineWidth',3) % corrected
pplot([1 sz1],[p12+p22 p12*sz1+p22]-sgm,'b','LineWidth',3)

enter image description here

Improve this answer
$\endgroup$
8
  • $\begingroup$ Wow, thanks so much for this answer! I also tried a lot of fft filtering but it was never convincing. Your approach to find just a couple of values in the area and then fit a line through them is great. Unfortunately, I need to determine the borders of the area very precisely, so the 95% intervall is a good estimation but not the end. Still, from your last image I can go on now, thanks again! $\endgroup$
    – Till
    Dec 9, 2022 at 10:13
  • $\begingroup$ Hi Till, also tried partial spectrum zeroing, some peaks, some quadrants, .. but didn't work, and I agree with you; directly manipulating the spectrum without a clear direction is not easy. The next step is neural network processing. Found a NN based obsolete MATLAB script for SAR (satellite imagery processing) that given a set of samples: trees, rocks, sea, oil spillage, .. the script colours imagery accordingly. But it may take me a while to upgrade it and make it fly. I will post if I fix it. $\endgroup$
    – John Bofarull
    Dec 10, 2022 at 16:48
  • $\begingroup$ Thanks, John! I wanted to avoid using a NN, cause I found it's overkill for this task. But apprantly it's more difficult than I thought. I could also train a network instead of using a predefined one, cause we have hundreds of those images. Let me know if you fix your NN. $\endgroup$
    – Till
    Dec 13, 2022 at 14:06
  • $\begingroup$ did TimWestcott use NN? $\endgroup$
    – John Bofarull
    Dec 13, 2022 at 21:43
  • $\begingroup$ Maybe you didn't see that it was me who answered my own question and not Tim^^ I didn't use any kind of NN and I didn't share any code cause I didn't think that question would raise any further interest. But if you want, I can edit my answer and explain more detailed what exactly did I do. $\endgroup$
    – Till
    Dec 14, 2022 at 7:36
1
$\begingroup$

Thanks, TimWescott! The heurisitc approach was a great idea. I divided the image in rectangles and evaluated the quotient of standard deviations in +45 and -45 direction. That gave me a good division of feature and non-feature area:

enter image description here

Then I just did a thresholding and that was it :)

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ Hi Tim, any chance to read, if ony a bit of it, the code you used to clearly differentiate perpendicular 'wave fronts'? $\endgroup$
    – John Bofarull
    Dec 10, 2022 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.