Message Boards Message Boards

Elliptical Fourier Descriptors

GROUPS:

Elliptical Fourier Descriptors for Biological Pattern Recognition

Elliptical Fourier Descriptors is a mathematical technique which utilizes a sum of ellipses with various harmonics to fit contours and silhouettes in an image. Broadly speaking the technique is widely used in Computer Vision and visual Pattern Matching. The first mode is by default an ellipse. Even the second or higher harmonics are all all ellipses, however for many practical purposes it is the sum of ellipses that fits the pattern. An exception to this might be when classifying different shapes from each other.

An example might be, say, we are interested to find the contour that matches the object. enter image description here

Given the image above (courtesy of PalaeoMath) a natural evolution of pattern identification will be as follows: evolution of ellipses with higher harmonics fitting the silhouette

courtesy of PalaeoMath)

Lately I have been experimenting with aggregates (clusters) of cells that naturally evolve, grow and undergo a change in morphology during their course of development. An attempt to classify the different aggregates is important. It may enable us to understand the mechanical and chemical interactions that determine cell behaviour and fate as a group.

I tried searching for a code in Mathematica (my language of choice) with a robust implementation of Elliptical Fourier Descriptors. Failing to find one I wrote one. The code is a translation from the works of "FEATURE EXTRACTION METHODS FOR CHARACTER RECOGNITION--A SURVEY - OIVIND DUE TRIER et al, Pattern Recognition Vol 29" and also I would like to give credit to Henrik Blidh for a nice Python implementation.

Options[EllipticalFourier] = {normalize -> True, 
   sizeInvariance -> True, order -> 10, ptsToplot -> 300, 
   locus -> {0., 0.}};
EllipticalFourier[img_, OptionsPattern[]] := 
 Module[{dXY, dt, t, T, \[Phi], contourpositions, contour, xi, 
   A0, \[Delta], C0, dCval, coeffL},
  contourpositions = 
   PixelValuePositions[MorphologicalPerimeter[img], 1];
  contour = 
   Part[contourpositions, Last@FindShortestTour[contourpositions]];
  dXY = Differences[contour];
  dt = (x \[Function] N@Sqrt[x])[Total /@ (dXY^2)];
  t = Prepend[Accumulate[dt], 0];
  T = Last@t;
  \[Phi] = 2 \[Pi] t/T;

  normalizeEFD[coeff_] := With[
    {\[Theta] = 
      0.5 ArcTan[(coeff[[1, 1]]^2 - coeff[[1, 2]]^2 + 
          coeff[[1, 3]]^2 - coeff[[1, 4]]^2),
        2 (coeff[[1, 1]] coeff[[1, 2]] + 
           coeff[[1, 3]] coeff[[1, 4]])]},
    Block[{coeffList = coeff, \[Psi], \[Psi]r, a, c, \[Theta]r},
     \[Theta]r[
       i_Integer] := {{Cos[i \[Theta]], -Sin[i \[Theta]]}, {Sin[
         i \[Theta]], Cos[i \[Theta]]}};
     coeffList = Partition[#, 2, 2] & /@ coeffList;
     coeffList = MapIndexed[(#1.\[Theta]r[First@#2]) &, coeffList];
     a = First@Flatten@First@coeffList;
     c = (Flatten@First@coeffList)[[3]];
     \[Psi] = ArcTan[a, c];
     \[Psi]r = {{Cos[\[Psi]], Sin[\[Psi]]}, {-Sin[\[Psi]], 
        Cos[\[Psi]]}};
     coeffList = MapIndexed[Flatten[\[Psi]r.#] &, coeffList];
     If[OptionValue[sizeInvariance], a = First@First@coeffList; 
      coeffList /= a, coeffList]]
    ];

  ellipticalCoefficients[order_] := Module[{coeffs, func},
    func[list : {{__} ..}, n_Integer] := 
     Block[{listtemp = list, const = T/(2 (n^2) ( \[Pi]^2)), 
       phiN = \[Phi] *n,
       dCosPhiN, dSinPhiN, an, bn, cn, dn},
      dCosPhiN = Cos[Rest@phiN] - Cos[Most@phiN];
      dSinPhiN = Sin[Rest@phiN] - Sin[Most@phiN];
      an = const*Total@(dXY[[All, 1]]/dt *dCosPhiN);
      bn = const*Total@(dXY[[All, 1]]/dt*dSinPhiN);
      cn = const*Total@(dXY[[All, 2]]/dt*dCosPhiN);
      dn = const*Total@(dXY[[All, 2]]/dt*dSinPhiN);
      listtemp[[n]] = {an, bn, cn, dn};
      listtemp
      ];
    coeffs = ConstantArray[0, {order, 4}];
    coeffL = 
     If[OptionValue@normalize, normalizeEFD[#], #] &@
      Fold[func, coeffs, Range@order]
    ];

  dcCoefficients := (
    xi =  Accumulate[Part[dXY, All, 1]] - (dXY[[All, 1]]/dt)*Rest@t;
    A0 = (1/T)*
      Total@(dXY[[All, 1]]/(2 dt) * Differences[t^2] +  xi*dt);
    \[Delta] = 
     Accumulate[Part[dXY, All, 2]] - (dXY[[All, 2]]/dt)*Rest@t;
    C0 = (1/T)*
      Total@(dXY[[All, 2]]/(2 dt) * Differences[t^2] +  \[Delta]*dt);
    {First@First@contour + A0, Last@First@contour + C0}
    );

  plotEFD[coeff_] := 
   Last@Reap@Block[{n = Length@coeff, nhalf, nrow = 2, ti,
       xt, yt, func, \[Theta], i = OptionValue[ptsToplot]},
      nhalf = Ceiling[n/2];
      ti = Array[# &, i, {0., 1}];
      xt = ConstantArray[OptionValue[locus][[1]], i];
      yt = ConstantArray[OptionValue[locus][[2]], i];
      Do[
       xt += 
        coeff[[j, 1]] Cos[2 (j) \[Pi] ti] + 
         coeff[[j, 2]] Sin[2 (j) \[Pi] ti];
       yt += 
        coeff[[j, 3]] Cos[2 (j) \[Pi] ti] + 
         coeff[[j, 4]] Sin[2 (j) \[Pi] ti];
       Sow@
        Graphics[{{Thick, XYZColor[0, 0, 1, 0.8], 
           Line@contour}, {Thick, XYZColor[1, 0, 0, 0.6], 
           Line@Thread[{xt, yt}]}, 
          Inset[Style[#, Bold, Darker@Green, FontSize -> 12] &@
            Text["Harmonic: " <> ToString@j]]}];
       , {j, 1, Length@coeff}];
      ];

  ellipticalCoefficients[OptionValue[order]];
  plotEFD[coeffL]
  ]

The code is initialized as follows:

img = Import@"C:\\Users\\Ali Hashmi\\Desktop\\f7masks\\72.tif";
results = 
  EllipticalFourier[img, normalize -> False, order -> 6, 
   ptsToplot -> 200, locus -> dcCoefficients];
GraphicsGrid[ArrayReshape[results, {2, 3}], ImageSize -> Large]

A couple of tests were performed and the results can be seen below:

enter image description here

It is quite evident that the red curve tends to converge toward the blue curve (cat) as the harmonics get higher.

I also happened to test them on my aggregates and the results look good.enter image description here

Next I will try to use the shape descriptors to quantitatively analyze the cell clusters

Attachments:
POSTED BY: Ali Hashmi
Answer
1 month ago

enter image description here - Congratulations! This post is now a Staff Pick! Thank you for your wonderful contributions. Please, keep them coming!

POSTED BY: Moderation Team
Answer
1 month ago

Thank you very much !

POSTED BY: Ali Hashmi
Answer
1 month ago

one can use the notebook below which is only a slight modification. It removes the need to use Clear on dcCoefficients each time which was defined in the above post as a global variable.

Attachments:
POSTED BY: Ali Hashmi
Answer
15 days ago

Group Abstract Group Abstract