Fourier Image Analysis

Settings and Helper Functions

Avoid error messages for similar symbols (e.g. viewI, viewF, viewC, viewM)

Off[General::spell]; Off[General::spell1];

Re-align images

Rotate2D[array_, 
          offset_] := RotateLeft[#, 
              (Subscript[offset, [[2]]])] & /@ RotateLeft[array, 
            (Subscript[offset, [[1]]])];
recenter[array_] := 
      Rotate2D[array, Floor[Dimensions[array]/2]]

Create {image,fourier} from source image or fourier-transform image

SetOptions[Fourier, 
        FourierParameters -> {1, 1}];
SetOptions[InverseFourier, 
        FourierParameters -> {1, 1}];
ItoIF[imageI_] := {imageI, 
        Chop[Fourier[imageI]]}
FtoIF[imageF_] := {Chop[InverseFourier[imageF]], imageF}

Import color image and convert to {grayscale,fourier}

loadIF[imagePath_String] := 
      ItoIF[recenter[Import[imagePath][[1, 
              1]] /. {{r_, g_, 
                  b_} -> (r + g + b)/(3*255. )}]]

Create {image,fourier} from function[x,y]

createIF[imageFunc_Function, size_] := 
      ItoIF[Table[
N[
            imageFunc[
              Mod[x + size[[1]]/2, 
                  size[[1]]] - 
                size[[1]]/2,  
              Mod[y + size[[2]]/2, 
                  size[[2]]] - 
                size[[2]]/
                  2]], 
{x, 0, 
            size[[1]] - 1}, {y, 0, 
            size[[2]] - 1}]]

Functions for viewing image and Fourier transform data.
First argument = {image,fourier} data to display. Second argument is optional. It gives function of one argument to convert 0-1 value to color; default is built-in function GrayLevel, which gives a simple grayscale image. Third argument is also optional. It gives the display function to use. Most useful non-default value is Identity, which generates the image but does not display it. Use to stack multiple images together using GraphicsArray[] without also displaying the individual images when they are generated.

viewI[im_, cf_:  GrayLevel, 
          disp_:  $DisplayFunction] := 
Show[
          Graphics[
            Raster[recenter[
                im[[1]]], 
{{0, 
                  0}, {Dimensions[
                      im[[1]]][[2]], Dimensions[
                      im[[1]]][[1]]}}, 
{Min[
                  im[[1]]], 
                Max[im[[1]]]}, 
              ColorFunction -> 
                cf], 
AspectRatio -> 
              Automatic], 
DisplayFunction -> 
            disp];
viewF[im_, cf_:  GrayLevel, disp_:  $DisplayFunction] := 
      ListDensityPlot[
        recenter[Abs[im[[2]]]], 
        Mesh -> False, Frame -> False, ColorFunction -> cf, 
        DisplayFunction -> disp]

Display image, transform and inverse(transform(image)) side-by-side

viewIF[im_, icf_:  GrayLevel, fcf_:  GrayLevel, 
        disp_:  $DisplayFunction] := 
Show[
        GraphicsArray[{
viewI[im, icf, 
              Identity], 
viewF[im, fcf, 
              Identity]
}], 
        DisplayFunction -> disp]

View convolution of two {image,fourier} objects; displays both source images, fourier transforms of both, fourier transform of convolution, and result. Both images should be the same size.
Returns {image,fourier} result of convolution

viewC[imageA_, imageB_,  icf_:  GrayLevel, fcf_:  GrayLevel, 
        disp_:  $DisplayFunction] := 
      Block[
{imageC = 
            FtoIF[imageA[[2]]*
                imageB[[2]]]}, 
Show[
          GraphicsArray[{
 (* 
                stack showing image A image/fourier and image B image/
                    fourier *) 
GraphicsArray[{
{viewIF[imageA, icf, fcf, 
                      Identity]}, 
{viewIF[imageB, icf, 
                      fcf, Identity]}
}], 
 (* 
                image C image/
                    fourier *) 
GraphicsArray[
{viewF[imageC, fcf, Identity], 
                  viewI[imageC, icf, 
                    Identity]}
]
}], 
DisplayFunction -> 
            disp]; 
imageC]

View image-space multiiplication of two {image,fourier} objects; displays both source images, fourier transforms of both, result and fourier transform of result. Both images should be the same size.
Returns resulting {image,fourier} of multiplication

viewM[imageA_, imageB_,  icf_:  GrayLevel, fcf_:  GrayLevel, 
        disp_:  $DisplayFunction] := 
      Block[
{imageC = 
            ItoIF[imageA[[1]]*
                imageB[[1]]]}, 
Show[
          GraphicsArray[{
 (* 
                stack showing image A image/fourier and image B image/
                    fourier *) 
GraphicsArray[{
{viewIF[imageA, icf, fcf, 
                      Identity]}, 
{viewIF[imageB, icf, 
                      fcf, Identity]}
}], 
 (* 
                image C image/
                    fourier *) 
GraphicsArray[
{viewIF[imageC, icf, fcf, 
                    Identity]}
]
}], 
DisplayFunction -> 
            disp]; 
imageC]

Basic transform/inverse transform

viewIF[
        createIF[Cos[#1*π/16 + #2*π/8] &, {64, 64}]];

\(\(image = 
        loadIF["/Users/olano/Class/demos/fourier/mathematica/maggy.jpg"];\)\)
viewIF[image];

square = 
        createIF[If[Abs[#1] < 2.5 && Abs[#2] < 2.5, 1, 0] &, 
          Dimensions[
            image[[1]]]];
viewIF[square];

viewC[image, square];

viewM[image, square];

circle = 
        createIF[If[#(1^2) + #(2^2) < (2.5^2), 1, 0] &, 
          Dimensions[
            image[[1]]]];
viewIF[circle];

viewC[image, circle];

gaussian = 
        createIF[
          N[(1/(5 (Sqrt[ 2  π]))) 
                Exp[-((#(1^2) + #(2^2))/(2*(5^2)))]] &, 
          Dimensions[
            image[[1]]]];
viewIF[gaussian];

viewC[image, gaussian];

dg = 
        createIF[Derivative[0, 
              1][(1/(5 (Sqrt[ 2  π]))) 
                Exp[-((#(1^2) + #(2^2))/(2*(5^2)))] &], 
          Dimensions[
            image[[1]]]];
viewIF[dg];

viewC[image, dg];

sample = 
        createIF[
          If[Mod[#1, 5] == 2  &&  Mod[#2, 5] == 2, 1, 0] &, 
          Dimensions[
            image[[1]]]];
viewIF[sample];

sampled = viewM[image, sample];

viewC[sampled, square];

viewC[sampled, gaussian];

Pattern Matching

mask = 
        createIF[N[If[#(1^2) + #(2^2) < (7^2), 1, 0]] &, 
          Dimensions[
            image[[1]]]];
viewI[{( .5 + mask[[1]])*
            Rotate2D[
              image[[1]], {11, 
                7}]}];

pattern = 
        ItoIF[mask[[1]]*
            Rotate2D[
              image[[1]], {11, 
                7}]];
crossCorrelation = viewC[image, Conjugate[pattern]];

patternSq = 
        Total[pattern[[1]]^2, 
          2];
imageSq = 
        FtoIF[ItoIF[
                image[[1]]^2][[2]] * 
            mask[[2]]];
match = 
        imageSq[[1]] - 
          2*crossCorrelation[[1]] + 
          patternSq;
Show[
        GraphicsArray[{GraphicsArray[
{{viewI[image, 
                    GrayLevel, Identity], 
                  viewI[imageSq, GrayLevel, 
                    Identity]}, 
{viewI[pattern, GrayLevel, 
                    Identity], 
                  viewI[crossCorrelation, GrayLevel, 
                    Identity]}}], 
viewI[{match}, 
              If[# <  .1, Hue[1, 1, 1], Hue[1, 0, #]] &, 
              Identity]}]];