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;
Model for DlgSettingsPointMatch and CmdSettingsPointMatch.
bool pixelFilteredIsOn(const QImage &image, int x, int y) const
Return true if specified filtered pixel is on.
Class for filtering image to remove unimportant information.
Representation of one matched point as produced from the point match algorithm.
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.
QPoint point() const
Return (x,y) coordinates as a point.
double maxPointSize() const
Get method for max point size.