using Image Correlation for Particle Image Velocimetry

Question

EDIT: can this approach be extended to make a robust Particle Image Velocimetry implementation in Mathematica? As of now, there is no package that does PIV in Mathematica (to my knowledge) and there are several in Matlab (which use Fourier Transform to do cross-correlation between images). A crude implementation of PIV can be found @ https://github.com/alihashmiii/simple-piv

Note: posted the same question on Wolfram Community (along with the notebook): http://community.wolfram.com/groups/-/m/t/1124830?p_p_auth=ERlfgQS9

I have been trying to implement a code for determining flow-field using Particle Image Velocimetry.

In this technique a user can take two images. Using small windows from the first image (which act as kernels) and search windows from the second image one can determine the cross-correlation which simply tells where the small window moves within a given search window. This process can be repeated between the second and the third image and so on.

A clear detail can be found in the second paragraph: http://www.physics.emory.edu/faculty/weeks//idl/piv.html

I have two images here (posting as a gif, you can save this and import it in mathematica as a list of two images):

I use the following code to generate the flow-field.

windowsize = 32; (* select window size *)
imgDim = ImageDimensions[images[[1]]]; (* dimensions for the images *)
imgone = ImageCrop[images[[1]], imgDim - (2*windowsize)]; (* removing
border from first image: we dont want to create windows at the borders *)

firstimgsplits = ImagePartition[imgone, windowsize]; 
(* breaking the first image into small windows *)
searchwindows = ImagePartition[images[[2]], windowsize*3, {windowsize, windowsize}];
(* breaking the second image into search windows *)

{dim1, dim2} = Dimensions@searchwindows;
H = Last@ImageDimensions[imgone];

(* get midpoints of the windows of the first frame *)
midptsFirst = Flatten[Table[{i windowsize + windowsize/2, 
 j (windowsize) + windowsize/2}, {i, 1, dim1}, {j, 1, dim2}], 1];

(* pts in the second image where correlation is max *)
correlationPts = Table[MorphologicalComponents[ImagePad[
   ImageAdjust@ImageCorrelate[searchwindows[[i + 1, j + 1]], 
     firstimgsplits[[i + 1, j + 1]], NormalizedSquaredEuclideanDistance, 
     PerformanceGoal -> "Quality"], {{j*windowsize, H - windowsize (j + 1)},
{H - windowsize (i + 1), windowsize i}}, White]]~Position~0,
{i, 0, dim1 - 1}, {j, 0, dim2 - 1}]~Flatten~2;

now when i create a flow-field from the displacement of points (red pts in the first image and cyan pts in second image) I can see that something is not right. My eyes tell me that the particles have move in a direction different from the ones found using ImageCorrelate

This should be rather straightforward for Mathematica. I do not know what is wrong in this simple piece of code. I will appreciate some help here.

ListAnimate@{Show[images[[1]], Graphics[{Red, Point@midptsFirst}]], 
Show[images[[2]], Graphics[{Cyan, PointSize[Medium], Point@correlationPts,
{Pink, Arrowheads[Small], MapThread[Arrow[{#1, #2}] &, {midptsFirst, correlationPts}]}}]]}

Answer 1

Full Code on github: https://github.com/alihashmiii/simple-piv/blob/master/flowtrack.m

PIV[image1_?ImageQ,image2_?ImageQ,win_Integer,pivmethod_]:=Module[{windowsize=win,
imgDim=ImageDimensions[image1],img1NoBorder,interrogateWin,searchWins,
dim,midPtsImg1,correlationPts,h,f,width,img1,img2,imgdata,height},

{img1,img2} = imagePreprocess[image1,image2]; (* image preprocessing *)

{width,height}=imgDim;
img1NoBorder=ImageCrop[img1,imgDim-(2*windowsize)]; (*removing border from the first image*)
interrogateWin=img1NoBorder~ImagePartition~windowsize; 
searchWins=ImagePartition[img2,3*windowsize,{windowsize,windowsize}];
dim=Dimensions[searchWins];
{f,h}= ImageDimensions[img1NoBorder];
midPtsImg1=Transpose[Table[Abs[{iwindowsize+windowsize/2,j(windowsize)+windowsize/2}-
{0,height}],{i,Last@dim},{j,First@dim}]]~Flatten~1;
correlationPts=Table[
imgdata=ImageData@interrogateWin[[i+1,j+1]];
If[And@@Map[First[imgdata]==#&,imgdata,{1}],
Indeterminate,
Composition[{#[[2]],height-#[[1]]}&,
FirstPosition[0], MorphologicalComponents]@ImagePad[ImageAdjust@ImageCorrelate[
searchWins[[i+1,j+1]],interrogateWin[[i+1,j+1]],NormalizedSquaredEuclideanDistance,
PerformanceGoal->"Quality"],
{{j windowsize,f-windowsize(j+1)},{h-windowsize(i+1),i windowsize}},White]],
{i,0,First[dim]-1},{j,0,Last[dim]-1}];

correlationPts=Cases[correlationPts,{__Integer}|Indeterminate,Infinity];
{midPtsImg1,correlationPts}=DeleteCases[Thread[{midPtsImg1,correlationPts}],
{_,Indeterminate}]\[Transpose];
Graphics[{Arrowheads[0.01],MapThread[Arrow[{##}]&,{midPtsImg1,correlationPts}]},
Frame->True]
];

This module is used to generate the following results:

Answer 2

Edit: based on conversation with Sander Huisman

This method is regarded as PTV and can work well for simple flows and sparse cases. For complex flows and dense cases PIV is preferred approach (see the question and the other answer)

Old

as mentioned by Henrik Schachner (http://community.wolfram.com/groups/-/m/t/1126003) we can use FindGeometricTransformto compute the flow-field with ease. The characteristic shape of the flow-field obtained using this approach is the same as the one obtained with ImageCorrelate method (that I used earlier). Nevertheless, this approach stands out to be the faster way of determining flows

images = Import["https://i.sstatic.net/TimgF.gif"];
imgDim = ImageDimensions[First@images];
windowsize = 32;
imgCorrD = First@imgDim;
img = ImageCrop[#, imgDim - (2*windowsize)] & /@ images; 
gtrf = Last[FindGeometricTransform @@ img];
rpts = RandomInteger[{windowsize,windowsize+First@ImageDimensions[First@img]}, {1000, 2}];
Graphics[{Arrowheads[.01], Darker@Cyan,  Arrow /@ Transpose[{Map[Abs[# - {0, 
imgCorrD}] &, #] &@rpts,
 Map[Abs[# - {0, imgCorrD}] &, #] &@gtrf[rpts]}]}, ImageSize -> 600]