2
$\begingroup$

The code below:

  1. reads in a number of png images, where the background is dark and objects of interest (here 3) are bright dots

  2. determines in each image the coordinates and other information about the objects

  3. tracks the objects by using Nearest

Single image (reduced size, cropped) with 3 tiny bright objects (at upper corners and lower center):

enter image description here

Pixels of a single object:

enter image description here

The code below can only be used if the number of objects in each image is the same.

As output the objects' coordinates are printed before and after tracking.

9 test images (zipped) are located here: http://bit.ly/28J56Pd

The notebook file which contains the code can be downloaded here: http://goo.gl/WhS7eP

My question is: Where can the code be improved?

Especially: How can I combine all ComponentMeasurement calls and assign the results to the set of used arrays (now in single command lines)?

ChoiceDialog[{FileNameSetter[Dynamic[imageDir], "Directory"], Dynamic[imageDir]}];

SetDirectory[imageDir];
dirBase = FileNameTake[imageDir];

(*read in file names and determine number of files*)

fNames = FileNames["*.png"];
numFiles = Length[fNames];

(*find number of objects in first image*)

i = 1;
image = Import[fNames[[i]]];
t = FindThreshold[image, Method -> "Entropy"];

binImage = DeleteSmallComponents[Binarize[image, t], 8];
newImage = ImageMultiply[image, binImage];

numObjects = Length@ComponentMeasurements[newImage, "IntensityCentroid"][[All, 2]];

(*define arrays*)

coordinates = Array[0 &, {numFiles, numObjects, 2}];
threshold = Array[0 &, {numFiles, numObjects}];
equivalentDiskRadius = Array[0 &, {numFiles, numObjects}];
totalIntensity = Array[0 &, {numFiles, numObjects}];
meanIntensity = Array[0 &, {numFiles, numObjects}];
standardDeviationIntensity = Array[0 &, {numFiles, numObjects}];
boundingBoxArea = Array[0 &, {numFiles, numObjects}];
numPixels = Array[0 &, {numFiles, numObjects}];

(*find coordinates and other relevant data of each object in first frame*)

i = 1;
coordinates[[i, All]] = 
  ComponentMeasurements[newImage, "IntensityCentroid"][[All, 2]];
threshold[[i, All]] = t;
equivalentDiskRadius[[i, All]] = 
  ComponentMeasurements[newImage, "EquivalentDiskRadius"][[All, 2]];
totalIntensity[[i, All]] = 
  ComponentMeasurements[newImage, "TotalIntensity"][[All, 2]];
meanIntensity[[i, All]] = 
  ComponentMeasurements[newImage, "MeanIntensity"][[All, 2]];
standardDeviationIntensity[[i]] = 
  ComponentMeasurements[newImage, "StandardDeviationIntensity"][[All, 
   2]];
boundingBoxArea[[i, All]] = 
  ComponentMeasurements[newImage, "BoundingBoxArea"][[All, 2]];
numPixels[[i, All]] = 
  ComponentMeasurements[newImage, "Count"][[All, 2]];

Print[i, " ", coordinates[[i, All]], "\n"];

(*loop over all images*) 

Do[

  image = Import[fNames[[i]]];
  t = FindThreshold[image, Method -> "Entropy"];
  binImage = 
   DeleteSmallComponents[Binarize[image, t], 
    8];(*8 has to be manually set*);
  newImage = ImageMultiply[image, binImage];

  coordinates[[i, All]] = 
   ComponentMeasurements[newImage, "IntensityCentroid"][[All, 2]];
  threshold[[i, All]] = t;
  equivalentDiskRadius[[i, All]] = 
   ComponentMeasurements[newImage, "EquivalentDiskRadius"][[All, 2]];
  totalIntensity[[i, All]] = 
   ComponentMeasurements[newImage, "TotalIntensity"][[All, 2]];
  meanIntensity[[i, All]] = 
   ComponentMeasurements[newImage, "MeanIntensity"][[All, 2]];
  standardDeviationIntensity[[i]] = 
   ComponentMeasurements[newImage, "StandardDeviationIntensity"][[All,
     2]];
  boundingBoxArea[[i, All]] = 
   ComponentMeasurements[newImage, "BoundingBoxArea"][[All, 2]];
  numPixels[[i, All]] = 
   ComponentMeasurements[newImage, "Count"][[All, 2]];

  Print[i, " ", coordinates[[i, All]]];

  (*sort coordinate order according to nearest neighbors, so that objects are tracked*)

  (* This is wrong: order = Flatten@Nearest[coordinates[[i - 1, All]] -> Automatic][coordinates[[i, All]]];*)
  order = Flatten@Nearest[coordinates[[i, All]] -> Automatic][coordinates[[i - 1, All]]];
  coordinates[[i, All]] = coordinates[[i, #]] & /@ order;
  threshold[[i, All]] = threshold[[i, #]] & /@ order;
  equivalentDiskRadius[[i, All]] = 
   equivalentDiskRadius[[i, #]] & /@ order;
  totalIntensity[[i, All]] = totalIntensity[[i, #]] & /@ order;
  meanIntensity[[i, All]] = meanIntensity[[i, #]] & /@ order;
  standardDeviationIntensity[[i]] = 
   standardDeviationIntensity[[i, #]] & /@ order;
  boundingBoxArea[[i, All]] = boundingBoxArea[[i, #]] & /@ order;
  numPixels[[i, All]] = numPixels[[i, #]] & /@ order;

  Print[i, " ", coordinates[[i, All]], "\n"],

  {i, 2, numFiles}
  ];

1 {{600.924,1321.05},{1282.57,1319.71},{927.97,726.014}}


2 {{601.472,1320.6},{1282.83,1320.66},{928.494,726.171}}

2 {{601.472,1320.6},{1282.83,1320.66},{928.494,726.171}}


3 {{601.599,1320.62},{1283.29,1321.66},{928.725,726.025}}

3 {{601.599,1320.62},{1283.29,1321.66},{928.725,726.025}}


4 {{1284.07,1322.6},{601.375,1320.79},{929.205,724.945}}

4 {{601.375,1320.79},{1284.07,1322.6},{929.205,724.945}}


5 {{1284.84,1323.69},{600.729,1321.26},{929.304,723.538}}

5 {{600.729,1321.26},{1284.84,1323.69},{929.304,723.538}}


6 {{1285.88,1324.76},{599.797,1321.9},{929.635,721.943}}

6 {{599.797,1321.9},{1285.88,1324.76},{929.635,721.943}}


7 {{1287.23,1325.46},{598.597,1322.92},{930.251,720.146}}

7 {{598.597,1322.92},{1287.23,1325.46},{930.251,720.146}}


8 {{1288.29,1325.99},{596.995,1324.19},{930.652,718.686}}

8 {{596.995,1324.19},{1288.29,1325.99},{930.652,718.686}}


9 {{595.348,1325.49},{1289.63,1326.32},{931.156,717.095}}

9 {{595.348,1325.49},{1289.63,1326.32},{931.156,717.095}}

Plot of the tracks:

ListPlot[# & /@ Transpose[coordinates], 
  AspectRatio -> Automatic, PlotRange -> Automatic, 
  PlotStyle -> AbsolutePointSize[10], Frame -> True, 
  AspectRatio -> 1, FrameLabel -> {{"y (m)", ""}, {"x (m)", ""}}, 
  BaseStyle -> {FontWeight -> "Bold", FontSize -> 20, 
    FontFamily -> "Calibri"}, ImageSize -> 600]

enter image description here

Improve this question
$\endgroup$
2
  • $\begingroup$ goo.gl/Niy2B6 is not working with "It was found to be violating our Terms of Service" $\endgroup$
    – Sumit
    Commented Jun 20, 2016 at 11:23
  • $\begingroup$ now it should work ... bit.ly/28J56Pd $\endgroup$
    – mrz
    Commented Jun 20, 2016 at 12:45

2 Answers 2

Reset to default
2
$\begingroup$

Where can the code be improved?

First of all: Try to avoid loops (Do, For, While) where possible - they're usually much more verbose than functional alternatives.

Second: Read the documentation. ComponentMeasurements can take a list of measurements, and will then return a list of values for each component.

Third: Use functions to encapsulate units of code.

For example, this function takes an image filename, finds components in the image and returns a list of measurements for each component:

Clear[getComponents]
getComponents[file_] := Module[{image, t, binImage},
  image = Import[file];
  t = FindThreshold[image, Method -> "Entropy"];
  binImage = DeleteSmallComponents[Binarize[image, t], 8];
  Sow[t];
  ComponentMeasurements[{MorphologicalComponents[binImage], image},
    {"IntensityCentroid", "EquivalentDiskRadius", "TotalIntensity", 
     "MeanIntensity", "StandardDeviationIntensity", "BoundingBoxArea",
      "Count"}][[All, 2]]
  ]

You can then use:

files = FileNames["some folder\\*.png"];

getComponents /@ files

to execute this function for all images in that folder. This will return a multidimensional array, with the indices [frame index, component index, measurement index].

You can use transpose to get the measurement index in front, and assign it to a list of symbols:

{intensityCentroid, equivalentDiskRadius, totalIntensity, meanIntensity, 
 standardDeviationIntensity, boundingBoxArea, count} = 
 Transpose[getComponents /@ files, {2, 3, 1}]

now intensityCentroid contains a 2d-array of intensity centroids (first index for the frame, second index for the component), equivalentDiskRadius contains a 2d-array of radii and so on.

You'll notice that we're missing the threshold - that's not a measurement, and it doesn't fit into our neat array structure. So I returned it using Sow, and I can get the values that were sowed in each call using Reap, so this line:

{{intensityCentroid, equivalentDiskRadius, totalIntensity, 
    meanIntensity, standardDeviationIntensity, boundingBoxArea, 
    count}, {thresholds}} = 
  Reap[Transpose[getComponents /@ files, {2, 3, 1}]];

I didn't reorder the components using Nearest above, to keep the code simple. It's not hard to add that, though:

Clear[getComponents, reorderComponents]
reorderComponents[components_, expectedPositions_] := Module[{indices},
  indices = (Nearest[expectedPositions -> Automatic] /@ 
      components[[All, 1]])[[All, 1]];
  components[[indices]]]

trackingCenters = {};
getComponents[file_] := Module[{image, t, binImage, components},
  image = Import[file];
  t = FindThreshold[image, Method -> "Entropy"];
  binImage = DeleteSmallComponents[Binarize[image, t], 8];
  Sow[t];
  components = 
   ComponentMeasurements[{MorphologicalComponents[binImage], image},
     {"IntensityCentroid", "EquivalentDiskRadius", "TotalIntensity", 
      "MeanIntensity", "StandardDeviationIntensity", 
      "BoundingBoxArea", "Count"}][[All, 2]];

  If[trackingCenters === {}, trackingCenters = components[[All, 1]]];
  reorderComponents[components, trackingCenters]
  ]
Improve this answer
$\endgroup$
3
  • $\begingroup$ This is great ... I will test you code earliest this evening or tomorrow. I am recording images with a 4Megapixels (1000Hz) Photron camera (CMOS). In the experiment (a radio frequency plasma discharge) I illuminate micron sized spherical particles with a sheet of laser light and observe the reflected light under 90 degrees.) As you see with the upper code I can only track particles if the number per image is constant. Important will be to develop a tracking code if the number varies (single particle appear or disappear) - as e.g. trackpy (Python) does. Do you have experience in this field? $\endgroup$
    – mrz
    Commented Jun 21, 2016 at 10:45
  • $\begingroup$ Perfect, both parts: object detection alone and in the second part combined with tracking work very well. $\endgroup$
    – mrz
    Commented Jun 24, 2016 at 10:52
  • $\begingroup$ I found some problems and would like to ask you for help. Here are some images tinyurl.com/jm4uuhy where 3 objects are moving. Somehow reorderComponents is not giving the correct assignments. With my upper code (see correction: order = Flatten@Nearest[coordinates[[i, All]] -> Automatic][coordinates[[i-1, All]]]; ) I plotted the trajectories (subfolder analysis). I used DeleteSmallComponents[Binarize[image, t], 1];. $\endgroup$
    – mrz
    Commented Oct 19, 2016 at 8:05
1
$\begingroup$ 
img = Import["https://i.sstatic.net/bkNTj.png"];
pos = ComponentMeasurements[img// Binarize,
        {"Centroid", "EquivalentDiskRadius"}][[All, 2, 1]];
Length[pos]
Show[img, Graphics[{Red, Circle[#, 10]} & /@ pos]]

3

enter image description here

You might be interested in Count flowers in an image and Count Elements in Image.

For your complete folder,

SetDirectory[NotebookDirectory[]];
Table[img[n] = Import["image_" <> ToString[n] <> ".png"];
      pos[n] = ComponentMeasurements[img[n] // Binarize,
           {"Centroid", "EquivalentDiskRadius"}][[All, 2, 1]];
      Show[img[n], Graphics[{Red, Circle[#, 50]} & /@ pos[n]], ImageSize -> 200]
,{n, 9}]
Improve this answer
$\endgroup$
2
  • $\begingroup$ This code is very interesting and helpful ... a memory problem could be that my image sets consist of about 3000 images, each 4 Megapixel. $\endgroup$
    – mrz
    Commented Jun 20, 2016 at 14:50
  • $\begingroup$ then I would suggest using img instead of img[n]. You don't need to store all the images. $\endgroup$
    – Sumit
    Commented Jun 20, 2016 at 15:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.