7 #include "ColorFilter.h" 8 #include "DocumentModelPointMatch.h" 9 #include "EngaugeAssert.h" 13 #include "PointMatchAlgorithm.h" 17 #include <QTextStream> 21 #define FOLD2DINDEX(i,j,jmax) ((i)*(jmax)+j) 23 const int PIXEL_OFF = -1;
25 const int PIXEL_ON = 1;
28 m_isGnuplot (isGnuplot)
30 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::PointMatchAlgorithm";
33 void PointMatchAlgorithm::allocateMemory(
double** array,
34 fftw_complex** arrayPrime,
38 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::allocateMemory";
40 *array =
new double [width * height];
41 ENGAUGE_CHECK_PTR(*array);
43 *arrayPrime =
new fftw_complex [width * height];
44 ENGAUGE_CHECK_PTR(*arrayPrime);
47 void PointMatchAlgorithm::assembleLocalMaxima(
double* convolution,
48 PointMatchList& listCreated,
52 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::assembleLocalMaxima";
55 const double SINGLE_PIXEL_CORRELATION = 1.0;
57 for (
int i = 0; i < width; i++) {
58 for (
int j = 0; j < height; j++) {
61 double convIJ = log10 (convolution[FOLD2DINDEX(i, j, height)]);
64 bool isLocalMax =
true;
65 for (
int iDelta = -1; (iDelta <= 1) && isLocalMax; iDelta++) {
67 int iNeighbor = i + iDelta;
68 if ((0 <= iNeighbor) && (iNeighbor < width)) {
70 for (
int jDelta = -1; (jDelta <= 1) && isLocalMax; jDelta++) {
72 int jNeighbor = j + jDelta;
73 if ((0 <= jNeighbor) && (jNeighbor < height)) {
76 double convNeighbor = log10 (convolution[FOLD2DINDEX(iNeighbor, jNeighbor, height)]);
77 if (convIJ < convNeighbor) {
81 }
else if (convIJ == convNeighbor) {
84 if ((jDelta < 0) || (jDelta == 0 && iDelta < 0)) {
95 (convIJ > SINGLE_PIXEL_CORRELATION) ) {
100 convolution [FOLD2DINDEX(i, j, height)]);
102 listCreated.append(t);
108 void PointMatchAlgorithm::computeConvolution(fftw_complex* imagePrime,
109 fftw_complex* samplePrime,
110 int width,
int height,
111 double** convolution,
115 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::computeConvolution";
117 fftw_complex* convolutionPrime;
119 allocateMemory(convolution,
125 conjugateMatrix(width,
130 multiplyMatrices(width,
137 fftw_plan pConvolution = fftw_plan_dft_c2r_2d(width,
143 fftw_execute (pConvolution);
145 releasePhaseArray(convolutionPrime);
149 double *temp =
new double [width * height];
150 ENGAUGE_CHECK_PTR(temp);
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)];
157 for (
int iFrom = 0; iFrom < width; iFrom++) {
158 for (
int jFrom = 0; jFrom < height; jFrom++) {
160 int iTo = (iFrom + sampleXCenter) % width;
161 int jTo = (jFrom + sampleYCenter) % height;
162 (*convolution) [FOLD2DINDEX(iTo, jTo, height)] = temp [FOLD2DINDEX(iFrom, jFrom, height)];
168 void PointMatchAlgorithm::conjugateMatrix(
int width,
170 fftw_complex* matrix)
172 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::conjugateMatrix";
174 ENGAUGE_CHECK_PTR(matrix);
176 for (
int x = 0; x < width; x++) {
177 for (
int y = 0; y < height; y++) {
179 int index = FOLD2DINDEX(x, y, height);
180 matrix [index] [1] = -1.0 * matrix [index] [1];
185 void PointMatchAlgorithm::dumpToGnuplot (
double* convolution,
188 const QString &filename)
const 190 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::dumpToGnuplot";
192 cout << GNUPLOT_FILE_MESSAGE.toLatin1().data() << filename.toLatin1().data() <<
"\n";
194 QFile file (filename);
195 if (file.open (QIODevice::WriteOnly | QIODevice::Text)) {
197 QTextStream str (&file);
199 str <<
"# Suggested gnuplot commands:" << endl;
200 str <<
"# set hidden3d" << endl;
201 str <<
"# splot \"" << filename <<
"\" u 1:2:3 with pm3d" << endl;
204 str <<
"# I J Convolution" << endl;
205 for (
int i = 0; i < width; i++) {
206 for (
int j = 0; j < height; j++) {
208 double convIJ = convolution[FOLD2DINDEX(i, j, height)];
209 str << i <<
" " << j <<
" " << convIJ << endl;
219 const QImage &imageProcessed,
221 const Points &pointsExisting)
223 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::findPoints" 224 <<
" samplePointPixels=" << samplePointPixels.count();
227 int originalWidth = imageProcessed.width();
228 int originalHeight = imageProcessed.height();
229 int width = optimizeLengthForFft(originalWidth);
230 int height = optimizeLengthForFft(originalHeight);
234 double *image, *sample, *convolution;
235 fftw_complex *imagePrime, *samplePrime;
238 int sampleXCenter, sampleYCenter, sampleXExtent, sampleYExtent;
239 loadImage(imageProcessed,
246 loadSample(samplePointPixels,
255 computeConvolution(imagePrime,
269 dumpToGnuplot(sample,
273 dumpToGnuplot(convolution,
276 "convolution.gnuplot");
281 PointMatchList listCreated;
282 assembleLocalMaxima(convolution,
289 QList<QPoint> pointsCreated;
290 PointMatchList::iterator itr;
291 for (itr = listCreated.begin(); itr != listCreated.end(); itr++) {
294 pointsCreated.push_back (triplet.
point ());
302 releaseImageArray(image);
303 releasePhaseArray(imagePrime);
304 releaseImageArray(sample);
305 releasePhaseArray(samplePrime);
306 releaseImageArray(convolution);
308 return pointsCreated;
311 void PointMatchAlgorithm::loadImage(
const QImage &imageProcessed,
313 const Points &pointsExisting,
317 fftw_complex** imagePrime)
319 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::loadImage";
321 allocateMemory(image,
326 populateImageArray(imageProcessed,
331 removePixelsNearExistingPoints(*image,
338 fftw_plan pImage = fftw_plan_dft_r2c_2d(width,
343 fftw_execute(pImage);
346 void PointMatchAlgorithm::loadSample(
const QList<PointMatchPixel> &samplePointPixels,
350 fftw_complex** samplePrime,
356 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::loadImage";
360 allocateMemory(sample,
365 populateSampleArray(samplePointPixels,
375 fftw_plan pSample = fftw_plan_dft_r2c_2d(width,
380 fftw_execute(pSample);
383 void PointMatchAlgorithm::multiplyMatrices(
int width,
389 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::multiplyMatrices";
391 for (
int x = 0; x < width; x++) {
392 for (
int y = 0; y < height; y++) {
394 int index = FOLD2DINDEX(x, y, height);
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];
402 int PointMatchAlgorithm::optimizeLengthForFft(
int originalLength)
406 const int INITIAL_CLOSEST_LENGTH = 0;
411 int closestLength = INITIAL_CLOSEST_LENGTH;
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) {
417 int newLength = power2 * power3 * power5 * power7;
418 if (originalLength <= newLength) {
420 if ((closestLength == INITIAL_CLOSEST_LENGTH) ||
421 (newLength < closestLength)) {
425 closestLength = newLength;
433 if (closestLength == INITIAL_CLOSEST_LENGTH) {
436 closestLength = originalLength;
439 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::optimizeLengthForFft" 440 <<
" originalLength=" << originalLength
441 <<
" newLength=" << closestLength;
443 return closestLength;
446 void PointMatchAlgorithm::populateImageArray(
const QImage &imageProcessed,
451 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::populateImageArray";
455 for (
int x = 0; x < width; x++) {
456 for (
int y = 0; y < height; y++) {
461 (*image) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ?
468 void PointMatchAlgorithm::populateSampleArray(
const QList<PointMatchPixel> &samplePointPixels,
477 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::populateSampleArray";
482 int xMin = width, yMin = height, xMax = 0, yMax = 0;
483 for (i = 0; i < (
unsigned int) samplePointPixels.size(); i++) {
485 int x = (samplePointPixels.at(i)).xOffset();
486 int y = (samplePointPixels.at(i)).yOffset();
487 if (first || (x < xMin))
489 if (first || (x > xMax))
491 if (first || (y < yMin))
493 if (first || (y > yMax))
499 const int border = 1;
508 for (x = 0; x < width; x++) {
509 for (y = 0; y < height; y++) {
510 (*sample) [FOLD2DINDEX(x, y, height)] = PIXEL_OFF;
517 double xSumOn = 0, ySumOn = 0, countOn = 0;
519 for (i = 0; i < (
unsigned int) samplePointPixels.size(); i++) {
522 x = (samplePointPixels.at(i)).xOffset() - xMin;
523 y = (samplePointPixels.at(i)).yOffset() - yMin;
524 ENGAUGE_ASSERT((0 < x) && (x < width));
525 ENGAUGE_ASSERT((0 < y) && (y < height));
527 bool pixelIsOn = samplePointPixels.at(i).pixelIsOn();
529 (*sample) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ? PIXEL_ON : PIXEL_OFF);
539 countOn = qMax (1.0, countOn);
540 *sampleXCenter = (int) (0.5 + xSumOn / countOn);
541 *sampleYCenter = (int) (0.5 + ySumOn / countOn);
544 *sampleXExtent = xMax - xMin + 1;
545 *sampleYExtent = yMax - yMin + 1;
548 void PointMatchAlgorithm::releaseImageArray(
double* array)
550 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::releaseImageArray";
552 ENGAUGE_CHECK_PTR(array);
556 void PointMatchAlgorithm::releasePhaseArray(fftw_complex* arrayPrime)
558 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::releasePhaseArray";
560 ENGAUGE_CHECK_PTR(arrayPrime);
564 void PointMatchAlgorithm::removePixelsNearExistingPoints(
double* image,
567 const Points &pointsExisting,
570 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::removePixelsNearExistingPoints";
572 for (
int i = 0; i < pointsExisting.size(); i++) {
574 int xPoint = pointsExisting.at(i).posScreen().x();
575 int yPoint = pointsExisting.at(i).posScreen().y();
578 int yMin = yPoint - pointSeparation;
581 int yMax = yPoint + pointSeparation;
582 if (imageHeight < yMax)
585 for (
int y = yMin; y < yMax; y++) {
588 int radical = pointSeparation * pointSeparation - (y - yPoint) * (y - yPoint);
591 int xMin = (int) (xPoint - qSqrt((
double) radical));
594 int xMax = xPoint + (xPoint - xMin);
595 if (imageWidth < xMax)
599 for (
int x = xMin; x < xMax; x++) {
601 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.