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