image FFT

FFTnd.bat

@:=' -- $RCSfile: FFTnd.bat,v $
@start c:\strawberry\perl\bin\perl.exe -x -S %0 %* & goto endofperl '; # vim: filetype=perl nowrap
#!perl

# note by doing an FFT I go from spacial domain to frequency domain...
#  i.e. bridging the physical world to the spiritual one <3
#
my ($source) = q$Source: /my/perl/script/fftnd.pl $ =~ /: (\S+)/;

eval { our $pwd = Win32::GetCwd() }; # current directory

#$ENV{PATH} = '\usr\tools\netpbm-10.27\bin;\usr\tools\jpeg-6b-4\bin;\windows';

use PDL;
use PDL::FFT qw(:Func);
use PDL::ImageND;             # provides kernctr()
use PDL::NiceSlice; # MATLAB style of arrays ...
use YAML::Syck qw(DumpFile);
use Brewed::PERMA qw(basen);

# -------------------------------------
# an attempted to do the FFT of love !
my $file = shift || 'love.jpeg'; # 526x526 image :)
my $outdir = 'results'; mkdir $outdir unless -d $outdir;


my $setdir = $pwd;
my @content = ();
if (-d $file) {
  $setdir = $file;
local *D;
opendir D,$setdir; @content = grep /.*\.(?:jpg|png|gif)/, readdir(D); closedir D;
} else {
 @content = ($file);
}

my $workdir = $ENV{TEMP}.'\workbox';
unlink $workdir; system sprintf'junction "%s" "%s"',$workdir,$setdir;

my $score = ();
my $nbfiles = $#content;
# ----------------------------------- [[
foreach my $filename ( @content[0 .. $nbfiles] ) {
my $file = $workdir.'\\'.$filename;
my ($fpath,$basen,$ext) = &basen($file);
$basen =~ s/ /_/g; # remove spaces !
# -------------------------------------
next if (-e "$outdir\\${basen}_filtered.jpeg");

# grab the piddle for the image
my $photo = rpic($file);
my ($depth, $xsize, $ysize) = dims($photo);

my $image = float($photo);
wpic($image,'image.jpeg');
my $spectrum = zeroes(3,$xsize,$ysize);
my $phase = zeroes(3,$xsize,$ysize);


# filter :
$green = $image(1,,;-); # green plane
$tmp = rvals($green)<50;        # Radially-symmetric filter function
$filter = kernctr $tmp, $tmp; # Shift origin to 0,0
#&displayk('filter.pgm',$filter);


$PDL::debug = 1;
$PDL::verbose = 1;

for my $plane (0 .. $depth-1) {
# create both real and imaginary plane ...
my $real = $image($plane,,;-);  # the ";-" means leave out dimensions of order 1 
my $imag = $real * 0;

fftnd $real, $imag;      # Perform Fourier transform

if (1) {
wpic(kernctr($imag,$imag), "imag${plane}.jpg");
$phase($plane,,) .= $imag(*,,);
$spectrum($plane,,) .= $real(*,,);
}

if (0) { use PDL::Image2D;
$rot = $real->rot2d(45,0,1);
# need to resize a/o crop
}

($real2,$imag2) = cmul $real,$imag , $filter, 0;

#&displayk("spectrum${plane}.jpeg",$real2);
#&displayk("phase${plane}.jpeg",$imag2);
ifftnd $real2, $imag2;      # Perform Inverse Fourier transform

$image($plane,,) .= $real2(*,,); # put this data back in image ( the "*" adds a dummy dimension, default order 1)

}
&displayk("phase.jpg",$phase);
&displayk("spectrum.jpg",$spectrum);
#link 'phase.jpg',"${basen}_phase.jpeg";
#link 'spectrum.jpg',"${basen}_freq.jpeg";

my $max = max($image);
printf "max: %g\n",$max;
$image->where( $image<0.0 ) .= 0.0;
$image = $image*255/$max;
wpic($image, "filtered.jpeg");


system 'PicasaPhotoViewer.exe filtered.jpeg';
rename 'phase.jpg', "$outdir\\${basen}_phase.jpg";
rename 'spectrum.jpg', "$outdir\\${basen}_spectrum.jpg";
rename 'filtered.jpeg', "$outdir\\${basen}_filtered.jpeg";

}
# ----------------------------------- ]]


sub displayk {
  my ($file,$piddle) = @_;
  my $centered = kernctr $piddle, $piddle; # Shift origin to 0,0
  my $scaled = byte($centered*255/max($centered));
  wpic($scaled, $file);
  system $file;
  return $?;
}
sub display {
  my ($file,$piddle) = @_;
  my $scaled = byte($piddle*255/max($piddle));
  wpic($scaled, $file);
  system $file;
  return $?;
}

1; # vim: ts=2 et noai
__END__
:endofperl

I was not able to get PDL::Graphics::PGPLOT to work on my machine so I use the image viewer form Picasa to display the result (work arround !)

therefore this code needs a few dependencies,