Extract diagonal area from image

Question

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:

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:

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

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.

Answer 1

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)

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)|')

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')

Centering spectrum

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

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

surf(|FFT2(A)|)

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

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,[])

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')

[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 []')

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*')

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);

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-')

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')

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)

Answer 2

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:

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