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);
287
288 // Copy sorted match points to output
290 PointMatchList::iterator itr;
291 for (itr = listCreated.begin(); itr != listCreated.end(); itr++) {
292
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
311void PointMatchAlgorithm::loadImage(const QImage &imageProcessed,
312 const DocumentModelPointMatch &modelPointMatch,
313 const Points &pointsExisting,
314 int width,
315 int height,
316 double** image,
318{
319 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::loadImage";
320
321 allocateMemory(image,
323 width,
324 height);
325
326 populateImageArray(imageProcessed,
327 width,
328 height,
329 image);
330
331 removePixelsNearExistingPoints(*image,
332 width,
333 height,
335 qFloor (modelPointMatch.maxPointSize()));
336
337 // Forward transform the image
339 height,
340 *image,
341 *imagePrime,
344}
345
346void PointMatchAlgorithm::loadSample(const QList<PointMatchPixel> &samplePointPixels,
347 int width,
348 int height,
349 double** sample,
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,
362 width,
363 height);
364
365 populateSampleArray(samplePointPixels,
366 width,
367 height,
368 sample,
373
374 // Forward transform the sample
376 height,
377 *sample,
381}
382
383void PointMatchAlgorithm::multiplyMatrices(int width,
384 int height,
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
402int 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
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
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
426 }
427 }
428 }
429 }
430 }
431 }
432
434
435 // No closest length was found, so just return the original length and expect slow fft performance
437 }
438
439 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::optimizeLengthForFft"
440 << " originalLength=" << originalLength
441 << " newLength=" << closestLength;
442
443 return closestLength;
444}
445
446void 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
455 for (int x = 0; x < width; x++) {
456 for (int y = 0; y < height; y++) {
458 x,
459 y);
460
461 (*image) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ?
462 PIXEL_ON :
463 PIXEL_OFF);
464 }
465 }
466}
467
468void 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
548void PointMatchAlgorithm::releaseImageArray(double* array)
549{
550 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::releaseImageArray";
551
553 delete[] array;
554}
555
556void PointMatchAlgorithm::releasePhaseArray(fftw_complex* arrayPrime)
557{
558 LOG4CPP_INFO_S ((*mainCat)) << "PointMatchAlgorithm::releasePhaseArray";
559
561 delete[] arrayPrime;
562}
563
564void 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)
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)
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}
const int INNER_RADIUS_MIN
#define ENGAUGE_ASSERT(cond)
Drop in replacement for Q_ASSERT if defined(QT_NO_DEBUG) && !defined(QT_FORCE_ASSERTS) define ENGAUGE...
#define ENGAUGE_CHECK_PTR(ptr)
#endif
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