1 #include "ColorFilter.h"
2 #include "DocumentModelPointMatch.h"
3 #include "EngaugeAssert.h"
6 #include "PointMatchAlgorithm.h"
10 #include <QTextStream>
14 #define FOLD2DINDEX(i,j,jmax) ((i)*(jmax)+j)
16 const int PIXEL_OFF = -1;
18 const int PIXEL_ON = 1;
21 m_isGnuplot (isGnuplot)
23 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::PointMatchAlgorithm";
26 void PointMatchAlgorithm::allocateMemory(
double** array,
27 fftw_complex** arrayPrime,
31 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::allocateMemory";
33 *array =
new double [width * height];
34 ENGAUGE_CHECK_PTR(*array);
36 *arrayPrime =
new fftw_complex [width * height];
37 ENGAUGE_CHECK_PTR(*arrayPrime);
40 void PointMatchAlgorithm::assembleLocalMaxima(
double* convolution,
41 PointMatchList& listCreated,
47 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::assembleLocalMaxima";
50 const double SINGLE_PIXEL_CORRELATION = 1.0;
52 for (
int i = 0; i < width; i++) {
53 for (
int j = 0; j < height; j++) {
55 double convIJ = convolution[FOLD2DINDEX(i, j, height)];
58 bool isLocalMax =
true;
59 for (
int iDelta = -1; (iDelta <= 1) && isLocalMax; iDelta++) {
61 int iNeighbor = i + iDelta;
62 if ((0 <= iNeighbor) && (iNeighbor < width)) {
64 for (
int jDelta = -1; (jDelta <= 1) && isLocalMax; jDelta++) {
66 int jNeighbor = j + jDelta;
67 if ((0 <= jNeighbor) && (jNeighbor < height)) {
69 double convNeighbor = convolution[FOLD2DINDEX(iNeighbor, jNeighbor, height)];
70 if (convIJ < convNeighbor) {
74 }
else if (convIJ == convNeighbor) {
77 if ((jDelta < 0) || (jDelta == 0 && iDelta < 0)) {
87 if (isLocalMax && (convIJ > SINGLE_PIXEL_CORRELATION)) {
92 convolution [FOLD2DINDEX(i, j, height)]);
94 listCreated.append(t);
100 void PointMatchAlgorithm::computeConvolution(fftw_complex* imagePrime,
101 fftw_complex* samplePrime,
102 int width,
int height,
103 double** convolution)
105 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::computeConvolution";
107 fftw_complex* convolutionPrime;
109 allocateMemory(convolution,
115 conjugateMatrix(width,
120 multiplyMatrices(width,
127 fftw_plan pConvolution = fftw_plan_dft_c2r_2d(width,
133 fftw_execute (pConvolution);
135 releasePhaseArray(convolutionPrime);
138 void PointMatchAlgorithm::conjugateMatrix(
int width,
140 fftw_complex* matrix)
142 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::conjugateMatrix";
144 ENGAUGE_CHECK_PTR(matrix);
146 for (
int x = 0; x < width; x++) {
147 for (
int y = 0; y < height; y++) {
149 int index = FOLD2DINDEX(x, y, height);
150 matrix [index] [1] = -1.0 * matrix [index] [1];
155 void PointMatchAlgorithm::dumpToGnuplot (
double* convolution,
158 const QString &filename)
const
160 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::dumpToGnuplot";
162 cout <<
"Writing gnuplot file: " << filename.toLatin1().data() <<
"\n";
164 QFile file (filename);
165 if (file.open (QIODevice::WriteOnly | QIODevice::Text)) {
167 QTextStream str (&file);
169 str <<
"# Suggested gnuplot commands:" << endl;
170 str <<
"# set hidden3d" << endl;
171 str <<
"# splot \"" << filename <<
"\" u 1:2:3 with pm3d" << endl;
174 str <<
"# I J Convolution" << endl;
175 for (
int i = 0; i < width; i++) {
176 for (
int j = 0; j < height; j++) {
178 double convIJ = convolution[FOLD2DINDEX(i, j, height)];
179 str << i <<
" " << j <<
" " << convIJ << endl;
189 const QImage &imageProcessed,
191 const Points &pointsExisting)
193 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::findPoints"
194 <<
" samplePointPixels=" << samplePointPixels.count();
197 int originalWidth = imageProcessed.width();
198 int originalHeight = imageProcessed.height();
199 int width = optimizeLengthForFft(originalWidth);
200 int height = optimizeLengthForFft(originalHeight);
204 double *image, *sample, *convolution;
205 fftw_complex *imagePrime, *samplePrime;
208 int sampleXCenter, sampleYCenter, sampleXExtent, sampleYExtent;
209 loadImage(imageProcessed,
216 loadSample(samplePointPixels,
225 computeConvolution(imagePrime,
237 dumpToGnuplot(sample,
241 dumpToGnuplot(convolution,
244 "convolution.gnuplot");
249 PointMatchList listCreated;
250 assembleLocalMaxima(convolution,
259 QList<QPoint> pointsCreated;
260 PointMatchList::iterator itr;
261 for (itr = listCreated.begin(); itr != listCreated.end(); itr++) {
264 pointsCreated.push_back (triplet.
point ());
272 releaseImageArray(image);
273 releasePhaseArray(imagePrime);
274 releaseImageArray(sample);
275 releasePhaseArray(samplePrime);
276 releaseImageArray(convolution);
278 return pointsCreated;
281 void PointMatchAlgorithm::loadImage(
const QImage &imageProcessed,
283 const Points &pointsExisting,
287 fftw_complex** imagePrime)
289 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::loadImage";
291 allocateMemory(image,
296 populateImageArray(imageProcessed,
301 removePixelsNearExistingPoints(*image,
308 fftw_plan pImage = fftw_plan_dft_r2c_2d(width,
313 fftw_execute(pImage);
316 void PointMatchAlgorithm::loadSample(
const QList<PointMatchPixel> &samplePointPixels,
320 fftw_complex** samplePrime,
326 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::loadImage";
330 allocateMemory(sample,
335 populateSampleArray(samplePointPixels,
345 fftw_plan pSample = fftw_plan_dft_r2c_2d(width,
350 fftw_execute(pSample);
353 void PointMatchAlgorithm::multiplyMatrices(
int width,
359 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::multiplyMatrices";
361 for (
int x = 0; x < width; x++) {
362 for (
int y = 0; y < height; y++) {
364 int index = FOLD2DINDEX(x, y, height);
366 out [index] [0] = in1 [index] [0] * in2 [index] [0] - in1 [index] [1] * in2 [index] [1];
367 out [index] [1] = in1 [index] [0] * in2 [index] [1] + in1 [index] [1] * in2 [index] [0];
372 int PointMatchAlgorithm::optimizeLengthForFft(
int originalLength)
376 const int INITIAL_CLOSEST_LENGTH = 0;
381 int closestLength = INITIAL_CLOSEST_LENGTH;
382 for (
int power2 = 1; power2 < originalLength; power2 *= 2) {
383 for (
int power3 = 1; power3 < originalLength; power3 *= 3) {
384 for (
int power5 = 1; power5 < originalLength; power5 *= 5) {
385 for (
int power7 = 1; power7 < originalLength; power7 *= 7) {
387 int newLength = power2 * power3 * power5 * power7;
388 if (originalLength <= newLength) {
390 if ((closestLength == INITIAL_CLOSEST_LENGTH) ||
391 (newLength < closestLength)) {
395 closestLength = newLength;
403 if (closestLength == INITIAL_CLOSEST_LENGTH) {
406 closestLength = originalLength;
409 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::optimizeLengthForFft"
410 <<
" originalLength=" << originalLength
411 <<
" newLength=" << closestLength;
413 return closestLength;
416 void PointMatchAlgorithm::populateImageArray(
const QImage &imageProcessed,
421 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::populateImageArray";
425 for (
int x = 0; x < width; x++) {
426 for (
int y = 0; y < height; y++) {
431 (*image) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ?
438 void PointMatchAlgorithm::populateSampleArray(
const QList<PointMatchPixel> &samplePointPixels,
447 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::populateSampleArray";
452 int xMin = width, yMin = height, xMax = 0, yMax = 0;
453 for (i = 0; i < (
unsigned int) samplePointPixels.size(); i++) {
455 int x = (samplePointPixels.at(i)).xOffset();
456 int y = (samplePointPixels.at(i)).yOffset();
457 if (first || (x < xMin))
459 if (first || (x > xMax))
461 if (first || (y < yMin))
463 if (first || (y > yMax))
469 const int border = 1;
478 for (x = 0; x < width; x++) {
479 for (y = 0; y < height; y++) {
480 (*sample) [FOLD2DINDEX(x, y, height)] = PIXEL_OFF;
487 double xSumOn = 0, ySumOn = 0, countOn = 0;
489 for (i = 0; i < (
unsigned int) samplePointPixels.size(); i++) {
492 x = (samplePointPixels.at(i)).xOffset() - xMin;
493 y = (samplePointPixels.at(i)).yOffset() - yMin;
494 ENGAUGE_ASSERT((0 < x) && (x < width));
495 ENGAUGE_ASSERT((0 < y) && (y < height));
497 bool pixelIsOn = samplePointPixels.at(i).pixelIsOn();
499 (*sample) [FOLD2DINDEX(x, y, height)] = (pixelIsOn ? PIXEL_ON : PIXEL_OFF);
509 countOn = qMax (1.0, countOn);
510 *sampleXCenter = (int) (0.5 + xSumOn / countOn);
511 *sampleYCenter = (int) (0.5 + ySumOn / countOn);
514 *sampleXExtent = xMax - xMin + 1;
515 *sampleYExtent = yMax - yMin + 1;
518 void PointMatchAlgorithm::releaseImageArray(
double* array)
520 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::releaseImageArray";
522 ENGAUGE_CHECK_PTR(array);
526 void PointMatchAlgorithm::releasePhaseArray(fftw_complex* arrayPrime)
528 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::releasePhaseArray";
530 ENGAUGE_CHECK_PTR(arrayPrime);
534 void PointMatchAlgorithm::removePixelsNearExistingPoints(
double* image,
537 const Points &pointsExisting,
540 LOG4CPP_INFO_S ((*mainCat)) <<
"PointMatchAlgorithm::removePixelsNearExistingPoints";
542 for (
int i = 0; i < pointsExisting.size(); i++) {
544 int xPoint = pointsExisting.at(i).posScreen().x();
545 int yPoint = pointsExisting.at(i).posScreen().y();
548 int yMin = yPoint - pointSeparation;
551 int yMax = yPoint + pointSeparation;
552 if (imageHeight < yMax)
555 for (
int y = yMin; y < yMax; y++) {
558 int radical = pointSeparation * pointSeparation - (y - yPoint) * (y - yPoint);
561 int xMin = (int) (xPoint - qSqrt((
double) radical));
564 int xMax = xPoint + (xPoint - xMin);
565 if (imageWidth < xMax)
569 for (
int x = xMin; x < xMax; x++) {
571 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.