Engauge Digitizer  2
 All Classes Files Functions Variables Enumerations Enumerator Friends Pages
PointMatchAlgorithm.cpp
1 #include "ColorFilter.h"
2 #include "DocumentModelPointMatch.h"
3 #include "EngaugeAssert.h"
4 #include <iostream>
5 #include "Logger.h"
6 #include "PointMatchAlgorithm.h"
7 #include <QFile>
8 #include <QImage>
9 #include <qmath.h>
10 #include <QTextStream>
11 
12 using namespace std;
13 
14 #define FOLD2DINDEX(i,j,jmax) ((i)*(jmax)+j)
15 
16 const int PIXEL_OFF = -1; // Negative of PIXEL_ON so two off pixels are just as valid as two on pixels when
17  // multiplied. One off pixel and one on pixel give +1 * -1 = -1 which reduces the correlation
18 const int PIXEL_ON = 1; // Arbitrary value as long as negative of PIXEL_OFF
19 
21  m_isGnuplot (isGnuplot)
22 {
23  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::PointMatchAlgorithm";
24 }
25 
26 void PointMatchAlgorithm::allocateMemory(double** array,
27  fftw_complex** arrayPrime,
28  int width,
29  int height)
30 {
31  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::allocateMemory";
32 
33  *array = new double [width * height];
34  ENGAUGE_CHECK_PTR(*array);
35 
36  *arrayPrime = new fftw_complex [width * height];
37  ENGAUGE_CHECK_PTR(*arrayPrime);
38 }
39 
40 void PointMatchAlgorithm::assembleLocalMaxima(double* convolution,
41  PointMatchList& listCreated,
42  int width,
43  int height,
44  int sampleXCenter,
45  int sampleYCenter)
46 {
47  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::assembleLocalMaxima";
48 
49  // Ignore tiny correlation values near zero by applying this threshold
50  const double SINGLE_PIXEL_CORRELATION = 1.0;
51 
52  for (int i = 0; i < width; i++) {
53  for (int j = 0; j < height; j++) {
54 
55  double convIJ = convolution[FOLD2DINDEX(i, j, height)];
56 
57  // Loop through nearest neighbor points
58  bool isLocalMax = true;
59  for (int iDelta = -1; (iDelta <= 1) && isLocalMax; iDelta++) {
60 
61  int iNeighbor = i + iDelta;
62  if ((0 <= iNeighbor) && (iNeighbor < width)) {
63 
64  for (int jDelta = -1; (jDelta <= 1) && isLocalMax; jDelta++) {
65 
66  int jNeighbor = j + jDelta;
67  if ((0 <= jNeighbor) && (jNeighbor < height)) {
68 
69  double convNeighbor = convolution[FOLD2DINDEX(iNeighbor, jNeighbor, height)];
70  if (convIJ < convNeighbor) {
71 
72  isLocalMax = false;
73 
74  } else if (convIJ == convNeighbor) {
75 
76  // Rare situation. In the event of a tie, the lower row/column wins (an arbitrary convention)
77  if ((jDelta < 0) || (jDelta == 0 && iDelta < 0)) {
78 
79  isLocalMax = false;
80  }
81  }
82  }
83  }
84  }
85  }
86 
87  if (isLocalMax && (convIJ > SINGLE_PIXEL_CORRELATION)) {
88 
89  // Save new local maximum
90  PointMatchTriplet t (i + sampleXCenter,
91  j + sampleYCenter,
92  convolution [FOLD2DINDEX(i, j, height)]);
93 
94  listCreated.append(t);
95  }
96  }
97  }
98 }
99 
100 void PointMatchAlgorithm::computeConvolution(fftw_complex* imagePrime,
101  fftw_complex* samplePrime,
102  int width, int height,
103  double** convolution)
104 {
105  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::computeConvolution";
106 
107  fftw_complex* convolutionPrime;
108 
109  allocateMemory(convolution,
110  &convolutionPrime,
111  width,
112  height);
113 
114  // Perform in-place conjugation of the sample since equation is F-1 {F(f) * F*(g)}
115  conjugateMatrix(width,
116  height,
117  samplePrime);
118 
119  // Perform the convolution in transform space
120  multiplyMatrices(width,
121  height,
122  imagePrime,
123  samplePrime,
124  convolutionPrime);
125 
126  // Backward transform the convolution
127  fftw_plan pConvolution = fftw_plan_dft_c2r_2d(width,
128  height,
129  convolutionPrime,
130  *convolution,
131  FFTW_ESTIMATE);
132 
133  fftw_execute (pConvolution);
134 
135  releasePhaseArray(convolutionPrime);
136 }
137 
138 void PointMatchAlgorithm::conjugateMatrix(int width,
139  int height,
140  fftw_complex* matrix)
141 {
142  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::conjugateMatrix";
143 
144  ENGAUGE_CHECK_PTR(matrix);
145 
146  for (int x = 0; x < width; x++) {
147  for (int y = 0; y < height; y++) {
148 
149  int index = FOLD2DINDEX(x, y, height);
150  matrix [index] [1] = -1.0 * matrix [index] [1];
151  }
152  }
153 }
154 
155 void PointMatchAlgorithm::dumpToGnuplot (double* convolution,
156  int width,
157  int height,
158  const QString &filename) const
159 {
160  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::dumpToGnuplot";
161 
162  cout << "Writing gnuplot file: " << filename.toLatin1().data() << "\n";
163 
164  QFile file (filename);
165  if (file.open (QIODevice::WriteOnly | QIODevice::Text)) {
166 
167  QTextStream str (&file);
168 
169  str << "# Suggested gnuplot commands:" << endl;
170  str << "# set hidden3d" << endl;
171  str << "# splot \"" << filename << "\" u 1:2:3 with pm3d" << endl;
172  str << endl;
173 
174  str << "# I J Convolution" << endl;
175  for (int i = 0; i < width; i++) {
176  for (int j = 0; j < height; j++) {
177 
178  double convIJ = convolution[FOLD2DINDEX(i, j, height)];
179  str << i << " " << j << " " << convIJ << endl;
180  }
181  str << endl; // pm3d likes blank lines between rows
182  }
183  }
184 
185  file.close();
186 }
187 
188 QList<QPoint> PointMatchAlgorithm::findPoints (const QList<PointMatchPixel> &samplePointPixels,
189  const QImage &imageProcessed,
190  const DocumentModelPointMatch &modelPointMatch,
191  const Points &pointsExisting)
192 {
193  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::findPoints"
194  << " samplePointPixels=" << samplePointPixels.count();
195 
196  // Use larger arrays for computations, if necessary, to improve fft performance
197  int originalWidth = imageProcessed.width();
198  int originalHeight = imageProcessed.height();
199  int width = optimizeLengthForFft(originalWidth);
200  int height = optimizeLengthForFft(originalHeight);
201 
202  // The untransformed (unprimed) and transformed (primed) storage arrays can be huge for big pictures, so minimize
203  // the number of allocated arrays at every point in time
204  double *image, *sample, *convolution;
205  fftw_complex *imagePrime, *samplePrime;
206 
207  // Compute convolution=F(-1){F(image)*F(*)(sample)}
208  int sampleXCenter, sampleYCenter, sampleXExtent, sampleYExtent;
209  loadImage(imageProcessed,
210  modelPointMatch,
211  pointsExisting,
212  width,
213  height,
214  &image,
215  &imagePrime);
216  loadSample(samplePointPixels,
217  width,
218  height,
219  &sample,
220  &samplePrime,
221  &sampleXCenter,
222  &sampleYCenter,
223  &sampleXExtent,
224  &sampleYExtent);
225  computeConvolution(imagePrime,
226  samplePrime,
227  width,
228  height,
229  &convolution);
230 
231  if (m_isGnuplot) {
232 
233  dumpToGnuplot(image,
234  width,
235  height,
236  "image.gnuplot");
237  dumpToGnuplot(sample,
238  width,
239  height,
240  "sample.gnuplot");
241  dumpToGnuplot(convolution,
242  width,
243  height,
244  "convolution.gnuplot");
245  }
246 
247  // Assemble local maxima, where each is the maxima centered in a region
248  // having a width of sampleWidth and a height of sampleHeight
249  PointMatchList listCreated;
250  assembleLocalMaxima(convolution,
251  listCreated,
252  width,
253  height,
254  sampleXCenter,
255  sampleYCenter);
256  qSort (listCreated);
257 
258  // Copy sorted match points to output
259  QList<QPoint> pointsCreated;
260  PointMatchList::iterator itr;
261  for (itr = listCreated.begin(); itr != listCreated.end(); itr++) {
262 
263  PointMatchTriplet triplet = *itr;
264  pointsCreated.push_back (triplet.point ());
265 
266  // Current order of maxima would be fine if they never overlapped. However, they often overlap so as each
267  // point is pulled off the list, and its pixels are removed from the image, we might consider updating all
268  // succeeding maxima here if those maximax overlap the just-removed maxima. The maxima list is kept
269  // in descending order according to correlation value
270  }
271 
272  releaseImageArray(image);
273  releasePhaseArray(imagePrime);
274  releaseImageArray(sample);
275  releasePhaseArray(samplePrime);
276  releaseImageArray(convolution);
277 
278  return pointsCreated;
279 }
280 
281 void PointMatchAlgorithm::loadImage(const QImage &imageProcessed,
282  const DocumentModelPointMatch &modelPointMatch,
283  const Points &pointsExisting,
284  int width,
285  int height,
286  double** image,
287  fftw_complex** imagePrime)
288 {
289  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::loadImage";
290 
291  allocateMemory(image,
292  imagePrime,
293  width,
294  height);
295 
296  populateImageArray(imageProcessed,
297  width,
298  height,
299  image);
300 
301  removePixelsNearExistingPoints(*image,
302  width,
303  height,
304  pointsExisting,
305  modelPointMatch.maxPointSize());
306 
307  // Forward transform the image
308  fftw_plan pImage = fftw_plan_dft_r2c_2d(width,
309  height,
310  *image,
311  *imagePrime,
312  FFTW_ESTIMATE);
313  fftw_execute(pImage);
314 }
315 
316 void PointMatchAlgorithm::loadSample(const QList<PointMatchPixel> &samplePointPixels,
317  int width,
318  int height,
319  double** sample,
320  fftw_complex** samplePrime,
321  int* sampleXCenter,
322  int* sampleYCenter,
323  int* sampleXExtent,
324  int* sampleYExtent)
325 {
326  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::loadImage";
327 
328  // Populate 2d sample array with same size (width x height) as image so fft transforms will have same
329  // dimensions, which means their transforms can be multiplied element-to-element
330  allocateMemory(sample,
331  samplePrime,
332  width,
333  height);
334 
335  populateSampleArray(samplePointPixels,
336  width,
337  height,
338  sample,
339  sampleXCenter,
340  sampleYCenter,
341  sampleXExtent,
342  sampleYExtent);
343 
344  // Forward transform the sample
345  fftw_plan pSample = fftw_plan_dft_r2c_2d(width,
346  height,
347  *sample,
348  *samplePrime,
349  FFTW_ESTIMATE);
350  fftw_execute(pSample);
351 }
352 
353 void PointMatchAlgorithm::multiplyMatrices(int width,
354  int height,
355  fftw_complex* in1,
356  fftw_complex* in2,
357  fftw_complex* out)
358 {
359  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::multiplyMatrices";
360 
361  for (int x = 0; x < width; x++) {
362  for (int y = 0; y < height; y++) {
363 
364  int index = FOLD2DINDEX(x, y, height);
365 
366  out [index] [0] = in1 [index] [0] * in2 [index] [0] - in1 [index] [1] * in2 [index] [1];
367  out [index] [1] = in1 [index] [0] * in2 [index] [1] + in1 [index] [1] * in2 [index] [0];
368  }
369  }
370 }
371 
372 int PointMatchAlgorithm::optimizeLengthForFft(int originalLength)
373 {
374  // LOG4CPP_INFO_S is below
375 
376  const int INITIAL_CLOSEST_LENGTH = 0;
377 
378  // Loop through powers, looking for the lowest multiple of 2^a * 3^b * 5^c * 7^d that is at or above the original
379  // length. Since the original length is expected to usually be less than 2000, we use only the smaller primes
380  // (2, 3, 5 and 7) and ignore 11 and 13 even though fftw can benefit from those as well
381  int closestLength = INITIAL_CLOSEST_LENGTH;
382  for (int power2 = 1; power2 < originalLength; power2 *= 2) {
383  for (int power3 = 1; power3 < originalLength; power3 *= 3) {
384  for (int power5 = 1; power5 < originalLength; power5 *= 5) {
385  for (int power7 = 1; power7 < originalLength; power7 *= 7) {
386 
387  int newLength = power2 * power3 * power5 * power7;
388  if (originalLength <= newLength) {
389 
390  if ((closestLength == INITIAL_CLOSEST_LENGTH) ||
391  (newLength < closestLength)) {
392 
393  // This is the best so far, so save it. No special weighting is given to powers of 2, although those
394  // can be more efficient than other
395  closestLength = newLength;
396  }
397  }
398  }
399  }
400  }
401  }
402 
403  if (closestLength == INITIAL_CLOSEST_LENGTH) {
404 
405  // No closest length was found, so just return the original length and expect slow fft performance
406  closestLength = originalLength;
407  }
408 
409  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::optimizeLengthForFft"
410  << " originalLength=" << originalLength
411  << " newLength=" << closestLength;
412 
413  return closestLength;
414 }
415 
416 void PointMatchAlgorithm::populateImageArray(const QImage &imageProcessed,
417  int width,
418  int height,
419  double** image)
420 {
421  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::populateImageArray";
422 
423  // Initialize memory with original image in real component, and imaginary component set to zero
424  ColorFilter colorFilter;
425  for (int x = 0; x < width; x++) {
426  for (int y = 0; y < height; y++) {
427  bool pixelIsOn = colorFilter.pixelFilteredIsOn (imageProcessed,
428  x,
429  y);
430 
431  (*image) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ?
432  PIXEL_ON :
433  PIXEL_OFF);
434  }
435  }
436 }
437 
438 void PointMatchAlgorithm::populateSampleArray(const QList<PointMatchPixel> &samplePointPixels,
439  int width,
440  int height,
441  double** sample,
442  int* sampleXCenter,
443  int* sampleYCenter,
444  int* sampleXExtent,
445  int* sampleYExtent)
446 {
447  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::populateSampleArray";
448 
449  // Compute bounds
450  bool first = true;
451  unsigned int i;
452  int xMin = width, yMin = height, xMax = 0, yMax = 0;
453  for (i = 0; i < (unsigned int) samplePointPixels.size(); i++) {
454 
455  int x = (samplePointPixels.at(i)).xOffset();
456  int y = (samplePointPixels.at(i)).yOffset();
457  if (first || (x < xMin))
458  xMin = x;
459  if (first || (x > xMax))
460  xMax = x;
461  if (first || (y < yMin))
462  yMin = y;
463  if (first || (y > yMax))
464  yMax = y;
465 
466  first = false;
467  }
468 
469  const int border = 1; // #pixels in border on each side
470 
471  xMin -= border;
472  yMin -= border;
473  xMax += border;
474  yMax += border;
475 
476  // Initialize memory with original image in real component, and imaginary component set to zero
477  int x, y;
478  for (x = 0; x < width; x++) {
479  for (y = 0; y < height; y++) {
480  (*sample) [FOLD2DINDEX(x, y, height)] = PIXEL_OFF;
481  }
482  }
483 
484  // We compute the center of mass of the on pixels. This means user does not have to precisely align
485  // the encompassing circle when selecting the sample point, since surrounding off pixels will not
486  // affect the center of mass computed only from on pixels
487  double xSumOn = 0, ySumOn = 0, countOn = 0;
488 
489  for (i = 0; i < (unsigned int) samplePointPixels.size(); i++) {
490 
491  // Place, quite arbitrarily, the sample image up against the top left corner
492  x = (samplePointPixels.at(i)).xOffset() - xMin;
493  y = (samplePointPixels.at(i)).yOffset() - yMin;
494  ENGAUGE_ASSERT((0 < x) && (x < width));
495  ENGAUGE_ASSERT((0 < y) && (y < height));
496 
497  bool pixelIsOn = samplePointPixels.at(i).pixelIsOn();
498 
499  (*sample) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ? PIXEL_ON : PIXEL_OFF);
500 
501  if (pixelIsOn) {
502  xSumOn += x;
503  ySumOn += y;
504  ++countOn;
505  }
506  }
507 
508  // Compute location of center of mass, which will represent the center of the point
509  countOn = qMax (1.0, countOn);
510  *sampleXCenter = (int) (0.5 + xSumOn / countOn);
511  *sampleYCenter = (int) (0.5 + ySumOn / countOn);
512 
513  // Dimensions of portion of array actually used by sample (rest is empty)
514  *sampleXExtent = xMax - xMin + 1;
515  *sampleYExtent = yMax - yMin + 1;
516 }
517 
518 void PointMatchAlgorithm::releaseImageArray(double* array)
519 {
520  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::releaseImageArray";
521 
522  ENGAUGE_CHECK_PTR(array);
523  delete[] array;
524 }
525 
526 void PointMatchAlgorithm::releasePhaseArray(fftw_complex* arrayPrime)
527 {
528  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::releasePhaseArray";
529 
530  ENGAUGE_CHECK_PTR(arrayPrime);
531  delete[] arrayPrime;
532 }
533 
534 void PointMatchAlgorithm::removePixelsNearExistingPoints(double* image,
535  int imageWidth,
536  int imageHeight,
537  const Points &pointsExisting,
538  int pointSeparation)
539 {
540  LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::removePixelsNearExistingPoints";
541 
542  for (int i = 0; i < pointsExisting.size(); i++) {
543 
544  int xPoint = pointsExisting.at(i).posScreen().x();
545  int yPoint = pointsExisting.at(i).posScreen().y();
546 
547  // Loop through rows of pixels
548  int yMin = yPoint - pointSeparation;
549  if (yMin < 0)
550  yMin = 0;
551  int yMax = yPoint + pointSeparation;
552  if (imageHeight < yMax)
553  yMax = imageHeight;
554 
555  for (int y = yMin; y < yMax; y++) {
556 
557  // Pythagorean theorem gives range of x values
558  int radical = pointSeparation * pointSeparation - (y - yPoint) * (y - yPoint);
559  if (0 < radical) {
560 
561  int xMin = (int) (xPoint - qSqrt((double) radical));
562  if (xMin < 0)
563  xMin = 0;
564  int xMax = xPoint + (xPoint - xMin);
565  if (imageWidth < xMax)
566  xMax = imageWidth;
567 
568  // Turn off pixels in this row of pixels
569  for (int x = xMin; x < xMax; x++) {
570 
571  image [FOLD2DINDEX(x, y, imageHeight)] = PIXEL_OFF;
572 
573  }
574  }
575  }
576  }
577 }
double maxPointSize() const
Get method for max point size.
Model for DlgSettingsPointMatch and CmdSettingsPointMatch.
Class for filtering image to remove unimportant information.
Definition: ColorFilter.h:12
QPoint point() const
Return (x,y) coordinates as a point.
Representation of one matched point as produced from the point match algorithm.
bool pixelFilteredIsOn(const QImage &image, int x, int y) const
Return true if specified filtered pixel is on.
QList< QPoint > findPoints(const QList< PointMatchPixel > &samplePointPixels, const QImage &imageProcessed, const DocumentModelPointMatch &modelPointMatch, const Points &pointsExisting)
Find points that match the specified sample point pixels. They are sorted by best-to-worst match...
PointMatchAlgorithm(bool isGnuplot)
Single constructor.