7 #include "ColorFilter.h"
8 #include "DocumentModelPointMatch.h"
9 #include "EngaugeAssert.h"
12 #include "PointMatchAlgorithm.h"
16 #include <QTextStream>
20 #define FOLD2DINDEX(i,j,jmax) ((i)*(jmax)+j)
22 const int PIXEL_OFF = -1;
24 const int PIXEL_ON = 1;
27 m_isGnuplot (isGnuplot)
29 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::PointMatchAlgorithm";
32 void PointMatchAlgorithm::allocateMemory(
double** array,
33 fftw_complex** arrayPrime,
37 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::allocateMemory";
39 *array =
new double [width * height];
40 ENGAUGE_CHECK_PTR(*array);
42 *arrayPrime =
new fftw_complex [width * height];
43 ENGAUGE_CHECK_PTR(*arrayPrime);
46 void PointMatchAlgorithm::assembleLocalMaxima(
double* convolution,
47 PointMatchList& listCreated,
51 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::assembleLocalMaxima";
54 const double SINGLE_PIXEL_CORRELATION = 1.0;
56 for (
int i = 0; i < width; i++) {
57 for (
int j = 0; j < height; j++) {
60 double convIJ = log10 (convolution[FOLD2DINDEX(i, j, height)]);
63 bool isLocalMax =
true;
64 for (
int iDelta = -1; (iDelta <= 1) && isLocalMax; iDelta++) {
66 int iNeighbor = i + iDelta;
67 if ((0 <= iNeighbor) && (iNeighbor < width)) {
69 for (
int jDelta = -1; (jDelta <= 1) && isLocalMax; jDelta++) {
71 int jNeighbor = j + jDelta;
72 if ((0 <= jNeighbor) && (jNeighbor < height)) {
75 double convNeighbor = log10 (convolution[FOLD2DINDEX(iNeighbor, jNeighbor, height)]);
76 if (convIJ < convNeighbor) {
80 }
else if (convIJ == convNeighbor) {
83 if ((jDelta < 0) || (jDelta == 0 && iDelta < 0)) {
94 (convIJ > SINGLE_PIXEL_CORRELATION) ) {
99 convolution [FOLD2DINDEX(i, j, height)]);
101 listCreated.append(t);
107 void PointMatchAlgorithm::computeConvolution(fftw_complex* imagePrime,
108 fftw_complex* samplePrime,
109 int width,
int height,
110 double** convolution,
114 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::computeConvolution";
116 fftw_complex* convolutionPrime;
118 allocateMemory(convolution,
124 conjugateMatrix(width,
129 multiplyMatrices(width,
136 fftw_plan pConvolution = fftw_plan_dft_c2r_2d(width,
142 fftw_execute (pConvolution);
144 releasePhaseArray(convolutionPrime);
148 double *temp =
new double [width * height];
149 ENGAUGE_CHECK_PTR(temp);
151 for (
int i = 0; i < width; i++) {
152 for (
int j = 0; j < height; j++) {
153 temp [FOLD2DINDEX(i, j, height)] = (*convolution) [FOLD2DINDEX(i, j, height)];
156 for (
int iFrom = 0; iFrom < width; iFrom++) {
157 for (
int jFrom = 0; jFrom < height; jFrom++) {
159 int iTo = (iFrom + sampleXCenter) % width;
160 int jTo = (jFrom + sampleYCenter) % height;
161 (*convolution) [FOLD2DINDEX(iTo, jTo, height)] = temp [FOLD2DINDEX(iFrom, jFrom, height)];
167 void PointMatchAlgorithm::conjugateMatrix(
int width,
169 fftw_complex* matrix)
171 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::conjugateMatrix";
173 ENGAUGE_CHECK_PTR(matrix);
175 for (
int x = 0; x < width; x++) {
176 for (
int y = 0; y < height; y++) {
178 int index = FOLD2DINDEX(x, y, height);
179 matrix [index] [1] = -1.0 * matrix [index] [1];
184 void PointMatchAlgorithm::dumpToGnuplot (
double* convolution,
187 const QString &filename)
const
189 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::dumpToGnuplot";
191 cout <<
"Writing gnuplot file: " << filename.toLatin1().data() <<
"\n";
193 QFile file (filename);
194 if (file.open (QIODevice::WriteOnly | QIODevice::Text)) {
196 QTextStream str (&file);
198 str <<
"# Suggested gnuplot commands:" << endl;
199 str <<
"# set hidden3d" << endl;
200 str <<
"# splot \"" << filename <<
"\" u 1:2:3 with pm3d" << endl;
203 str <<
"# I J Convolution" << endl;
204 for (
int i = 0; i < width; i++) {
205 for (
int j = 0; j < height; j++) {
207 double convIJ = convolution[FOLD2DINDEX(i, j, height)];
208 str << i <<
" " << j <<
" " << convIJ << endl;
218 const QImage &imageProcessed,
220 const Points &pointsExisting)
222 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::findPoints"
223 <<
" samplePointPixels=" << samplePointPixels.count();
226 int originalWidth = imageProcessed.width();
227 int originalHeight = imageProcessed.height();
228 int width = optimizeLengthForFft(originalWidth);
229 int height = optimizeLengthForFft(originalHeight);
233 double *image, *sample, *convolution;
234 fftw_complex *imagePrime, *samplePrime;
237 int sampleXCenter, sampleYCenter, sampleXExtent, sampleYExtent;
238 loadImage(imageProcessed,
245 loadSample(samplePointPixels,
254 computeConvolution(imagePrime,
268 dumpToGnuplot(sample,
272 dumpToGnuplot(convolution,
275 "convolution.gnuplot");
280 PointMatchList listCreated;
281 assembleLocalMaxima(convolution,
288 QList<QPoint> pointsCreated;
289 PointMatchList::iterator itr;
290 for (itr = listCreated.begin(); itr != listCreated.end(); itr++) {
293 pointsCreated.push_back (triplet.
point ());
301 releaseImageArray(image);
302 releasePhaseArray(imagePrime);
303 releaseImageArray(sample);
304 releasePhaseArray(samplePrime);
305 releaseImageArray(convolution);
307 return pointsCreated;
310 void PointMatchAlgorithm::loadImage(
const QImage &imageProcessed,
312 const Points &pointsExisting,
316 fftw_complex** imagePrime)
318 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::loadImage";
320 allocateMemory(image,
325 populateImageArray(imageProcessed,
330 removePixelsNearExistingPoints(*image,
337 fftw_plan pImage = fftw_plan_dft_r2c_2d(width,
342 fftw_execute(pImage);
345 void PointMatchAlgorithm::loadSample(
const QList<PointMatchPixel> &samplePointPixels,
349 fftw_complex** samplePrime,
355 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::loadImage";
359 allocateMemory(sample,
364 populateSampleArray(samplePointPixels,
374 fftw_plan pSample = fftw_plan_dft_r2c_2d(width,
379 fftw_execute(pSample);
382 void PointMatchAlgorithm::multiplyMatrices(
int width,
388 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::multiplyMatrices";
390 for (
int x = 0; x < width; x++) {
391 for (
int y = 0; y < height; y++) {
393 int index = FOLD2DINDEX(x, y, height);
395 out [index] [0] = in1 [index] [0] * in2 [index] [0] - in1 [index] [1] * in2 [index] [1];
396 out [index] [1] = in1 [index] [0] * in2 [index] [1] + in1 [index] [1] * in2 [index] [0];
401 int PointMatchAlgorithm::optimizeLengthForFft(
int originalLength)
405 const int INITIAL_CLOSEST_LENGTH = 0;
410 int closestLength = INITIAL_CLOSEST_LENGTH;
411 for (
int power2 = 1; power2 < originalLength; power2 *= 2) {
412 for (
int power3 = 1; power3 < originalLength; power3 *= 3) {
413 for (
int power5 = 1; power5 < originalLength; power5 *= 5) {
414 for (
int power7 = 1; power7 < originalLength; power7 *= 7) {
416 int newLength = power2 * power3 * power5 * power7;
417 if (originalLength <= newLength) {
419 if ((closestLength == INITIAL_CLOSEST_LENGTH) ||
420 (newLength < closestLength)) {
424 closestLength = newLength;
432 if (closestLength == INITIAL_CLOSEST_LENGTH) {
435 closestLength = originalLength;
438 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::optimizeLengthForFft"
439 <<
" originalLength=" << originalLength
440 <<
" newLength=" << closestLength;
442 return closestLength;
445 void PointMatchAlgorithm::populateImageArray(
const QImage &imageProcessed,
450 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::populateImageArray";
454 for (
int x = 0; x < width; x++) {
455 for (
int y = 0; y < height; y++) {
460 (*image) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ?
467 void PointMatchAlgorithm::populateSampleArray(
const QList<PointMatchPixel> &samplePointPixels,
476 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::populateSampleArray";
481 int xMin = width, yMin = height, xMax = 0, yMax = 0;
482 for (i = 0; i < (
unsigned int) samplePointPixels.size(); i++) {
484 int x = (samplePointPixels.at(i)).xOffset();
485 int y = (samplePointPixels.at(i)).yOffset();
486 if (first || (x < xMin))
488 if (first || (x > xMax))
490 if (first || (y < yMin))
492 if (first || (y > yMax))
498 const int border = 1;
507 for (x = 0; x < width; x++) {
508 for (y = 0; y < height; y++) {
509 (*sample) [FOLD2DINDEX(x, y, height)] = PIXEL_OFF;
516 double xSumOn = 0, ySumOn = 0, countOn = 0;
518 for (i = 0; i < (
unsigned int) samplePointPixels.size(); i++) {
521 x = (samplePointPixels.at(i)).xOffset() - xMin;
522 y = (samplePointPixels.at(i)).yOffset() - yMin;
523 ENGAUGE_ASSERT((0 < x) && (x < width));
524 ENGAUGE_ASSERT((0 < y) && (y < height));
526 bool pixelIsOn = samplePointPixels.at(i).pixelIsOn();
528 (*sample) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ? PIXEL_ON : PIXEL_OFF);
538 countOn = qMax (1.0, countOn);
539 *sampleXCenter = (int) (0.5 + xSumOn / countOn);
540 *sampleYCenter = (int) (0.5 + ySumOn / countOn);
543 *sampleXExtent = xMax - xMin + 1;
544 *sampleYExtent = yMax - yMin + 1;
547 void PointMatchAlgorithm::releaseImageArray(
double* array)
549 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::releaseImageArray";
551 ENGAUGE_CHECK_PTR(array);
555 void PointMatchAlgorithm::releasePhaseArray(fftw_complex* arrayPrime)
557 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::releasePhaseArray";
559 ENGAUGE_CHECK_PTR(arrayPrime);
563 void PointMatchAlgorithm::removePixelsNearExistingPoints(
double* image,
566 const Points &pointsExisting,
569 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::removePixelsNearExistingPoints";
571 for (
int i = 0; i < pointsExisting.size(); i++) {
573 int xPoint = pointsExisting.at(i).posScreen().x();
574 int yPoint = pointsExisting.at(i).posScreen().y();
577 int yMin = yPoint - pointSeparation;
580 int yMax = yPoint + pointSeparation;
581 if (imageHeight < yMax)
584 for (
int y = yMin; y < yMax; y++) {
587 int radical = pointSeparation * pointSeparation - (y - yPoint) * (y - yPoint);
590 int xMin = (int) (xPoint - qSqrt((
double) radical));
593 int xMax = xPoint + (xPoint - xMin);
594 if (imageWidth < xMax)
598 for (
int x = xMin; x < xMax; x++) {
600 image [FOLD2DINDEX(x, y, imageHeight)] = PIXEL_OFF;
double maxPointSize() const
Get method for max point size.
Model for DlgSettingsPointMatch and CmdSettingsPointMatch.
Class for filtering image to remove unimportant information.
QPoint point() const
Return (x,y) coordinates as a point.
Representation of one matched point as produced from the point match algorithm.
bool pixelFilteredIsOn(const QImage &image, int x, int y) const
Return true if specified filtered pixel is 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.