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