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