program FourierFilters;

// By Bruce Comber
// Western Australia Police, Perth, WA, Australia

{
First press the FFT button to perform a Fourier Transform on an
image of your choice.

Next, place one or two flags on an image in the frequency
domain for the purpose of producing either a low, high or band
pass filter image. The frequency of the cut off is determined by
one flag for the high and low pass, and two for the band pass.

The low pass procedure generates a sine wave of value 1 at the
centre of the image going to 0 at the flag. The high pass image
is the low pass subtracted from 1.

The band pass is a sine wave from 0 to 0 between the flags.

The filter image is mutiplied by the user's active image to 
produce the filtered result. 

This routine only accomodates square images as it does not detect
radii that may be longer than a side of the image.

Dimensions must be a power of 2 (ie. 64, 128, 256, 512 etc).
}


var
actimg,high,low,band,image;    {activeimage, output from low, high and bandpass procedures}
centre,centx, centy;     {centre of active image}
flags,radius1,radius2;   {placed flags an the radii of each from the centre}
i,j;                     {looping increments}
length;                  {difference between inner and outer radii}
position;                {location of a point on 'length'}


procedure mutilpy_and_inverse(image);
begin
image:=ifft(image*actimg);
show(image);
end;


procedure Transform ;
button btn_F,'Calculate Fourier Transform' ;
var
  new ;
begin
getactiveimage(actimg);
new := FFT( actimg ) ;
show(new);
setdisplaymode( new,dm_Random+dm_Magnitude ) ;
end;


procedure lowpassfilter;
button btn_l,'LOW PASS Fourier filter, place flag to mark radius';
begin
getactiveimage(actimg);
flags:=getflags(actimg);
low:=createimage(typ_single,getxsize(actimg),getysize(actimg));

if OnDesktop(low) then
  Update(low)
else
  Show(low,'Low pass filter');

centx:=getxsize(low)/2;
centy:=getysize(low)/2;

radius1:= integer(sqrt(sqr(centx-flags[0,0])+sqr(centy-flags[1,0])));

for i:=(0-radius1)to radius1 do
begin
    for j:=(0-radius1) to radius1 do
    begin
    if (sqrt(sqr(i)+sqr(j)))<radius1 then
       begin
	       if (i=0) and (j=0) then
    	   low[centx+i,centx+j]:=1
		   else
       	   low[centx+i,centx+j]:=1-(1+(sin(((180*sqrt(sqr(i)+sqr(j))/radius1)-90)*(pi/180))))/2;
       end;
    end;
    end;
update(low);
mutilpy_and_inverse(low);
end;


procedure bandpassfilter;
button btn_b,'BAND PASS Fourier filter, place 2 flags for the radii';
begin
getactiveimage(actimg);
centre:=getxsize(actimg)/2;
flags:=getflags(actimg);

band:=single(actimg-actimg);

if OnDesktop(band) then
  Update(band)
else
  Show(band,'Band pass filter');

radius1:=sqrt(sqr(centre-flags[0,0])+sqr(centre-flags[1,0]));
radius2:=sqrt(sqr(centre-flags[0,1])+sqr(centre-flags[1,1]));

if radius2>radius1 then
begin
i:=radius2;
radius2:=radius1;
radius1:=i;
end;

length:=radius1-radius2;

for i:=(centre-integer(radius1)) to (centre+integer(radius1))do
    begin
    for j:=(centre-integer(radius1)) to (centre+integer(radius1)) do
        begin
        if (sqrt(sqr(centre-i)+sqr(centre-j))<radius1) and (sqrt(sqr(centre-i)+sqr(centre-j))>radius2) then
        begin
             position:= sqrt(sqr(centre-i)+sqr(centre-j)) -radius2;
             band[i,j]:=(sin(((360*position/length)-90)*3.1415/180)+1)/2;
        end;
        end;
    end;
update(band);
mutilpy_and_inverse(band);
end;


procedure highpassfilter;
button btn_H,'HIGH PASS Fourier filter, place flag to mark radius';
begin
getactiveimage(actimg);
flags:=getflags(actimg);
high:=createimage(typ_single,getxsize(actimg),getysize(actimg));

if OnDesktop(high) then
  Update(high)
else
  Show(high,'High pass filter');

centx:=getxsize(high)/2;
centy:=getysize(high)/2;

radius1:= integer(sqrt(sqr(centx-flags[0,0])+sqr(centy-flags[1,0])));

for i:=(0-radius1)to radius1 do
begin
    for j:=(0-radius1) to radius1 do
    begin
    if (sqrt(sqr(i)+sqr(j)))<radius1 then
       begin
       	   if (i=0) and (j=0) then
    	   high[centx+i,centx+j]:=1
		   else
	   	   high[centx+i,centx+j]:=1-(1+(sin(((180*sqrt(sqr(i)+sqr(j))/radius1)-90)*(pi/180))))/2;
       end;
    end;
    end;
high:=1-high;
update(high);
mutilpy_and_inverse(high);
end;



begin
end