function FILTERED=filtr(MATRIX,HEIGHT,WIDTH,Xc,Yc)

% Erases selected areas of the PSD.
%
% Considering the PSD visualization and its 2 two caracteristics legs,
% you can define 4 areas : 3, 6, 9, 12 o'clock.
%
% filtr allows you to erase 1 or 2 pairs of opposite areas,
% staying remote of a distance gap/2 (that you are invited to input)
% from the legs you want to keep intact.
%
% EXAMPLES:
%
% >> filtr('cell.jpg')
%    where cell.jpg is the cellular pattern image.
% OR
% >> filtr(MATRIX)
%    where MATRIX is the double-type matrix of
%    the cellular pattern image, i.e. MATRIX=double(imread('cell.jpg'));
%
% You can also specify the HEIGHT, the WIDTH, and the center coordinates Xc and Yc
% of your image. If you do so, type:
% >> filtr(...,HEIGHT,WIDTH,Xc,Yc)
% Enter either 1 argument (the image or its matrix), or 5.
%
% >> FILTERED=filtr(...) returns as FILTERED the filtered cellular pattern matrix
% (CAUTION: it's not the PSD matrix!).


if isstr(MATRIX)
   MATRIX=double(imread(MATRIX)); % Turns the image into a matrix
end
% The previous loop allows us to input the image either as a function argument
% (then CELLIMAGE is a string, for instance cell.jpg),
% or the image matrix (then CELLIMAGE is a matrix).

if nargin==1
   FINAL=calcpsd(MATRIX);
   [HEIGHT WIDTH]=size(FINAL);
   [Xc Yc]=center(FINAL);
end

if nargin~=1 & nargin~=5
   error (' Number of expected input arguments: 1 or 5');
end


% The difficulty is to paint black at the right location of the PSD,
% but we can't do it on the actual PSD, since it is then impossible
% to come back to the image matrix (because of .^2, log,... operators).
% So we process the filter on the NONFINAL matrix, whose pixels are
% on the same location as in the PSD (with of course not the same intensity),
% but which allows us to come back to the image afterwards.

m=0;

WIND=MATRIX.*fspecial('gaussian',[HEIGHT WIDTH],0.2*min(HEIGHT,WIDTH)); % Windows the image with a gaussian filter
WIND=grayscale(WIND); % Rescales the entries between 0 and 255.
NONFINAL=fftshift(fft2(WIND));

while m~=2
   
   disp(' ');
   disp(' Please enter the angles of the 2 legs of the PSD :');
   TH1=2*pi/360*input('\n Positive angle (in degrees) = ');
   TH2=2*pi/360*input('\n Negative angle (in degrees) = ');
   
   g=input('\n Surviving gap = ');
   
   g1=g/2/cos(TH1);
   g2=g/2/cos(TH2);
   
   sm=menu('What area would you like to erase?','Around 3 and 9 o''clock', 'Around 6 and 12 o''clock','Around 3,6,9 and 12 o''clock');
   
   MASK=ones(HEIGHT,WIDTH);
   
   I=0;
   h=waitbar(0,'Please wait, filtering in progress...');
   
   for Y=1:WIDTH
      waitbar(I/WIDTH,h)
      I=I+1;
      X1=Xc-(Y-Yc)*tan(TH1);
      X2=Xc-(Y-Yc)*tan(TH2);
      
      for X=1:HEIGHT
         
         if (sm==1 | sm==3) & ((Y<=Yc & X2+g2<X & X<X1-g1) | (Y>=Yc & X1+g1<X & X<X2-g2)), MASK(X,Y)=0; end
         if (sm==2 | sm==3) & ((X<X2-g2 & X<X1-g1) | (X1+g1<X & X2+g2<X)), MASK(X,Y)=0; end
         
      end
   end
   
   close(h)
   
   MASKED=MASK.*NONFINAL;
   FILTERED=real(ifft2(ifftshift(MASKED)));
   LOCALFINAL=grayscale(calcpsd(FILTERED));
   
   figure;colormap(gray);imshow(LOCALFINAL);image(LOCALFINAL);title('Filtr -- Filtered Power Spectral Density');axis on;zoom on
   figure;colormap(gray);imshow(FILTERED);image(FILTERED);title('Filtr -- Filtered Cellular Pattern');axis on;zoom on % Imshow>>Image gives better results than Image alone. I don't know why...
   
   m=menu('What would you like to do now?','Try again','Quit or return to main menu');
   
end

return