Engauge Digitizer 2
Loading...
Searching...
No Matches
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
19using namespace std;
20
21#define FOLD2DINDEX(i,j,jmax) ((i)*(jmax)+j)
22
23const 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
25const 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
33void 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
47void 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
108void 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
168void 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
185void 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
218QList<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 std::sort (listCreated.begin(),
287 listCreated.end());
288
289 // Copy sorted match points to output
290 QList<QPoint> pointsCreated;
291 PointMatchList::iterator itr;
292 for (itr = listCreated.begin(); itr != listCreated.end(); itr++) {
293
294 PointMatchTriplet triplet = *itr;
295 pointsCreated.push_back (triplet.point ());
296
297 // Current order of maxima would be fine if they never overlapped. However, they often overlap so as each
298 // point is pulled off the list, and its pixels are removed from the image, we might consider updating all
299 // succeeding maxima here if those maximax overlap the just-removed maxima. The maxima list is kept
300 // in descending order according to correlation value
301 }
302
303 releaseImageArray(image);
304 releasePhaseArray(imagePrime);
305 releaseImageArray(sample);
306 releasePhaseArray(samplePrime);
307 releaseImageArray(convolution);
308
309 return pointsCreated;
310}
311
312void PointMatchAlgorithm::loadImage(const QImage &imageProcessed,
313 const DocumentModelPointMatch &modelPointMatch,
314 const Points &pointsExisting,
315 int width,
316 int height,
317 double** image,
318 fftw_complex** imagePrime)
319{
320 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::loadImage";
321
322 allocateMemory(image,
323 imagePrime,
324 width,
325 height);
326
327 populateImageArray(imageProcessed,
328 width,
329 height,
330 image);
331
332 removePixelsNearExistingPoints(*image,
333 width,
334 height,
335 pointsExisting,
336 qFloor (modelPointMatch.maxPointSize()));
337
338 // Forward transform the image
339 fftw_plan pImage = fftw_plan_dft_r2c_2d(width,
340 height,
341 *image,
342 *imagePrime,
343 FFTW_ESTIMATE);
344 fftw_execute(pImage);
345}
346
347void PointMatchAlgorithm::loadSample(const QList<PointMatchPixel> &samplePointPixels,
348 int width,
349 int height,
350 double** sample,
351 fftw_complex** samplePrime,
352 int* sampleXCenter,
353 int* sampleYCenter,
354 int* sampleXExtent,
355 int* sampleYExtent)
356{
357 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::loadImage";
358
359 // Populate 2d sample array with same size (width x height) as image so fft transforms will have same
360 // dimensions, which means their transforms can be multiplied element-to-element
361 allocateMemory(sample,
362 samplePrime,
363 width,
364 height);
365
366 populateSampleArray(samplePointPixels,
367 width,
368 height,
369 sample,
370 sampleXCenter,
371 sampleYCenter,
372 sampleXExtent,
373 sampleYExtent);
374
375 // Forward transform the sample
376 fftw_plan pSample = fftw_plan_dft_r2c_2d(width,
377 height,
378 *sample,
379 *samplePrime,
380 FFTW_ESTIMATE);
381 fftw_execute(pSample);
382}
383
384void PointMatchAlgorithm::multiplyMatrices(int width,
385 int height,
386 fftw_complex* in1,
387 fftw_complex* in2,
388 fftw_complex* out)
389{
390 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::multiplyMatrices";
391
392 for (int x = 0; x < width; x++) {
393 for (int y = 0; y < height; y++) {
394
395 int index = FOLD2DINDEX(x, y, height);
396
397 out [index] [0] = in1 [index] [0] * in2 [index] [0] - in1 [index] [1] * in2 [index] [1];
398 out [index] [1] = in1 [index] [0] * in2 [index] [1] + in1 [index] [1] * in2 [index] [0];
399 }
400 }
401}
402
403int PointMatchAlgorithm::optimizeLengthForFft(int originalLength)
404{
405 // LOG4CPP_INFO_S is below
406
407 const int INITIAL_CLOSEST_LENGTH = 0;
408
409 // Loop through powers, looking for the lowest multiple of 2^a * 3^b * 5^c * 7^d that is at or above the original
410 // length. Since the original length is expected to usually be less than 2000, we use only the smaller primes
411 // (2, 3, 5 and 7) and ignore 11 and 13 even though fftw can benefit from those as well
412 int closestLength = INITIAL_CLOSEST_LENGTH;
413 for (int power2 = 1; power2 < originalLength; power2 *= 2) {
414 for (int power3 = 1; power3 < originalLength; power3 *= 3) {
415 for (int power5 = 1; power5 < originalLength; power5 *= 5) {
416 for (int power7 = 1; power7 < originalLength; power7 *= 7) {
417
418 int newLength = power2 * power3 * power5 * power7;
419 if (originalLength <= newLength) {
420
421 if ((closestLength == INITIAL_CLOSEST_LENGTH) ||
422 (newLength < closestLength)) {
423
424 // This is the best so far, so save it. No special weighting is given to powers of 2, although those
425 // can be more efficient than other
426 closestLength = newLength;
427 }
428 }
429 }
430 }
431 }
432 }
433
434 if (closestLength == INITIAL_CLOSEST_LENGTH) {
435
436 // No closest length was found, so just return the original length and expect slow fft performance
437 closestLength = originalLength;
438 }
439
440 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::optimizeLengthForFft"
441 << " originalLength=" << originalLength
442 << " newLength=" << closestLength;
443
444 return closestLength;
445}
446
447void PointMatchAlgorithm::populateImageArray(const QImage &imageProcessed,
448 int width,
449 int height,
450 double** image)
451{
452 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::populateImageArray";
453
454 // Initialize memory with original image in real component, and imaginary component set to zero
455 ColorFilter colorFilter;
456 for (int x = 0; x < width; x++) {
457 for (int y = 0; y < height; y++) {
458 bool pixelIsOn = colorFilter.pixelFilteredIsOn (imageProcessed,
459 x,
460 y);
461
462 (*image) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ?
463 PIXEL_ON :
464 PIXEL_OFF);
465 }
466 }
467}
468
469void PointMatchAlgorithm::populateSampleArray(const QList<PointMatchPixel> &samplePointPixels,
470 int width,
471 int height,
472 double** sample,
473 int* sampleXCenter,
474 int* sampleYCenter,
475 int* sampleXExtent,
476 int* sampleYExtent)
477{
478 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::populateSampleArray";
479
480 // Compute bounds
481 bool first = true;
482 unsigned int i;
483 int xMin = width, yMin = height, xMax = 0, yMax = 0;
484 for (i = 0; i < static_cast<unsigned int> (samplePointPixels.size()); i++) {
485
486 int x = samplePointPixels.at(signed (i)).xOffset();
487 int y = samplePointPixels.at(signed (i)).yOffset();
488 if (first || (x < xMin))
489 xMin = x;
490 if (first || (x > xMax))
491 xMax = x;
492 if (first || (y < yMin))
493 yMin = y;
494 if (first || (y > yMax))
495 yMax = y;
496
497 first = false;
498 }
499
500 const int border = 1; // #pixels in border on each side
501
502 xMin -= border;
503 yMin -= border;
504 xMax += border;
505 yMax += border;
506
507 // Initialize memory with original image in real component, and imaginary component set to zero
508 int x, y;
509 for (x = 0; x < width; x++) {
510 for (y = 0; y < height; y++) {
511 (*sample) [FOLD2DINDEX(x, y, height)] = PIXEL_OFF;
512 }
513 }
514
515 // We compute the center of mass of the on pixels. This means user does not have to precisely align
516 // the encompassing circle when selecting the sample point, since surrounding off pixels will not
517 // affect the center of mass computed only from on pixels
518 double xSumOn = 0, ySumOn = 0, countOn = 0;
519
520 for (i = 0; i < static_cast<unsigned int> (samplePointPixels.size()); i++) {
521
522 // Place, quite arbitrarily, the sample image up against the top left corner
523 x = (samplePointPixels.at(signed (i))).xOffset() - xMin;
524 y = (samplePointPixels.at(signed (i))).yOffset() - yMin;
525 ENGAUGE_ASSERT((0 < x) && (x < width));
526 ENGAUGE_ASSERT((0 < y) && (y < height));
527
528 bool pixelIsOn = samplePointPixels.at(signed (i)).pixelIsOn();
529
530 (*sample) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ? PIXEL_ON : PIXEL_OFF);
531
532 if (pixelIsOn) {
533 xSumOn += x;
534 ySumOn += y;
535 ++countOn;
536 }
537 }
538
539 // Compute location of center of mass, which will represent the center of the point
540 countOn = qMax (1.0, countOn);
541 *sampleXCenter = qFloor (0.5 + xSumOn / countOn);
542 *sampleYCenter = qFloor (0.5 + ySumOn / countOn);
543
544 // Dimensions of portion of array actually used by sample (rest is empty)
545 *sampleXExtent = xMax - xMin + 1;
546 *sampleYExtent = yMax - yMin + 1;
547}
548
549void PointMatchAlgorithm::releaseImageArray(double* array)
550{
551 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::releaseImageArray";
552
553 ENGAUGE_CHECK_PTR(array);
554 delete[] array;
555}
556
557void PointMatchAlgorithm::releasePhaseArray(fftw_complex* arrayPrime)
558{
559 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::releasePhaseArray";
560
561 ENGAUGE_CHECK_PTR(arrayPrime);
562 delete[] arrayPrime;
563}
564
565void PointMatchAlgorithm::removePixelsNearExistingPoints(double* image,
566 int imageWidth,
567 int imageHeight,
568 const Points &pointsExisting,
569 int pointSeparation)
570{
571 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::removePixelsNearExistingPoints";
572
573 for (int i = 0; i < pointsExisting.size(); i++) {
574
575 int xPoint = qFloor (pointsExisting.at(signed (i)).posScreen().x());
576 int yPoint = qFloor (pointsExisting.at(signed (i)).posScreen().y());
577
578 // Loop through rows of pixels
579 int yMin = yPoint - pointSeparation;
580 if (yMin < 0)
581 yMin = 0;
582 int yMax = yPoint + pointSeparation;
583 if (imageHeight < yMax)
584 yMax = imageHeight;
585
586 for (int y = yMin; y < yMax; y++) {
587
588 // Pythagorean theorem gives range of x values
589 int radical = pointSeparation * pointSeparation - (y - yPoint) * (y - yPoint);
590 if (0 < radical) {
591
592 int xMin = qFloor (xPoint - qSqrt(double (radical)));
593 if (xMin < 0)
594 xMin = 0;
595 int xMax = xPoint + (xPoint - xMin);
596 if (imageWidth < xMax)
597 xMax = imageWidth;
598
599 // Turn off pixels in this row of pixels
600 for (int x = xMin; x < xMax; x++) {
601
602 image [FOLD2DINDEX(x, y, imageHeight)] = PIXEL_OFF;
603
604 }
605 }
606 }
607 }
608}
#define ENGAUGE_ASSERT(cond)
Drop in replacement for Q_ASSERT.
#define ENGAUGE_CHECK_PTR(ptr)
Drop in replacement for Q_CHECK_PTR.
log4cpp::Category * mainCat
Definition Logger.cpp:14
const int PIXEL_ON
#define FOLD2DINDEX(i, j, jmax)
const int PIXEL_OFF
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.
Model for DlgSettingsPointMatch and CmdSettingsPointMatch.
double maxPointSize() const
Get method for max point size.
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.
Representation of one matched point as produced from the point match algorithm.
QPoint point() const
Return (x,y) coordinates as a point.
#define LOG4CPP_INFO_S(logger)
Definition convenience.h:18
const QString GNUPLOT_FILE_MESSAGE