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,
35 int width,
36 int height)
37{
38 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::allocateMemory";
39
40 *array = new double [unsigned (width * height)];
42
43 *arrayPrime = new fftw_complex [unsigned (width * height)];
45}
46
47void PointMatchAlgorithm::assembleLocalMaxima(double* convolution,
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)
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 &&
96
97 // Save new local maximum
99 j,
101
102 listCreated.append(t);
103 }
104 }
105 }
106}
107
108void PointMatchAlgorithm::computeConvolution(fftw_complex* imagePrime,
110 int width, int height,
111 double** convolution,
112 int sampleXCenter,
113 int sampleYCenter)
114{
115 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::computeConvolution";
116
118
119 allocateMemory(convolution,
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,
128
129 // Perform the convolution in transform space
130 multiplyMatrices(width,
131 height,
135
136 // Backward transform the convolution
138 height,
142
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)];
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,
171{
172 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::conjugateMatrix";
173
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
195 if (file.open (QIODevice::WriteOnly | QIODevice::Text)) {
196
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
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;
236
237 // Compute convolution=F(-1){F(image)*F(*)(sample)}
239 loadImage(imageProcessed,
240 modelPointMatch,
242 width,
243 height,
244 &image,
245 &imagePrime);
246 loadSample(samplePointPixels,
247 width,
248 height,
249 &sample,
255 computeConvolution(imagePrime,
257 width,
258 height,
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
282 assembleLocalMaxima(convolution,
284 width,
285 height);
286 std::sort (listCreated.begin(),
287 listCreated.end());
288
289 // Copy sorted match points to output
291 PointMatchList::iterator itr;
292 for (itr = listCreated.begin(); itr != listCreated.end(); itr++) {
293
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,
319{
320 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::loadImage";
321
322 allocateMemory(image,
324 width,
325 height);
326
327 populateImageArray(imageProcessed,
328 width,
329 height,
330 image);
331
332 removePixelsNearExistingPoints(*image,
333 width,
334 height,
336 qFloor (modelPointMatch.maxPointSize()));
337
338 // Forward transform the image
340 height,
341 *image,
342 *imagePrime,
345}
346
347void PointMatchAlgorithm::loadSample(const QList<PointMatchPixel> &samplePointPixels,
348 int width,
349 int height,
350 double** sample,
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,
363 width,
364 height);
365
366 populateSampleArray(samplePointPixels,
367 width,
368 height,
369 sample,
374
375 // Forward transform the sample
377 height,
378 *sample,
382}
383
384void PointMatchAlgorithm::multiplyMatrices(int width,
385 int height,
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
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
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
427 }
428 }
429 }
430 }
431 }
432 }
433
435
436 // No closest length was found, so just return the original length and expect slow fft performance
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
456 for (int x = 0; x < width; x++) {
457 for (int y = 0; y < height; y++) {
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
554 delete[] array;
555}
556
557void PointMatchAlgorithm::releasePhaseArray(fftw_complex* arrayPrime)
558{
559 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::releasePhaseArray";
560
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)
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)
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}
const int INNER_RADIUS_MIN
#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
Class for filtering image to remove unimportant information.
Definition ColorFilter.h:21
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.
#define LOG4CPP_INFO_S(logger)
Definition convenience.h:18
const QString GNUPLOT_FILE_MESSAGE