21#define FOLD2DINDEX(i,j,jmax) ((i)*(jmax)+j)
28 m_isGnuplot (isGnuplot)
33void PointMatchAlgorithm::allocateMemory(
double** array,
34 fftw_complex** arrayPrime,
40 *array =
new double [unsigned (width * height)];
43 *arrayPrime =
new fftw_complex [unsigned (width * height)];
47void PointMatchAlgorithm::assembleLocalMaxima(
double* convolution,
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) ) {
98 PointMatchTriplet t (i,
102 listCreated.append(t);
108void PointMatchAlgorithm::computeConvolution(fftw_complex* imagePrime,
109 fftw_complex* samplePrime,
110 int width,
int height,
111 double** convolution,
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 [unsigned (width * height)];
152 for (
int i = 0; i < width; i++) {
153 for (
int j = 0; j < height; j++) {
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;
168void PointMatchAlgorithm::conjugateMatrix(
int width,
170 fftw_complex* matrix)
176 for (
int x = 0; x < width; x++) {
177 for (
int y = 0; y < height; y++) {
180 matrix [index] [1] = -1.0 * matrix [index] [1];
185void PointMatchAlgorithm::dumpToGnuplot (
double* convolution,
188 const QString &filename)
const
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)
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");
282 assembleLocalMaxima(convolution,
286 std::sort (listCreated.begin(),
290 QList<QPoint> pointsCreated;
291 PointMatchList::iterator itr;
292 for (itr = listCreated.begin(); itr != listCreated.end(); itr++) {
295 pointsCreated.push_back (triplet.
point ());
303 releaseImageArray(image);
304 releasePhaseArray(imagePrime);
305 releaseImageArray(sample);
306 releasePhaseArray(samplePrime);
307 releaseImageArray(convolution);
309 return pointsCreated;
312void PointMatchAlgorithm::loadImage(
const QImage &imageProcessed,
314 const Points &pointsExisting,
318 fftw_complex** imagePrime)
322 allocateMemory(image,
327 populateImageArray(imageProcessed,
332 removePixelsNearExistingPoints(*image,
339 fftw_plan pImage = fftw_plan_dft_r2c_2d(width,
344 fftw_execute(pImage);
347void PointMatchAlgorithm::loadSample(
const QList<PointMatchPixel> &samplePointPixels,
351 fftw_complex** samplePrime,
361 allocateMemory(sample,
366 populateSampleArray(samplePointPixels,
376 fftw_plan pSample = fftw_plan_dft_r2c_2d(width,
381 fftw_execute(pSample);
384void PointMatchAlgorithm::multiplyMatrices(
int width,
392 for (
int x = 0; x < width; x++) {
393 for (
int y = 0; y < height; y++) {
397 out [index] [0] = in1 [index] [0] * in2 [index] [0] - in1 [index] [1] * in2 [index] [1];
398 out [index] [1] = in1 [index] [0] * in2 [index] [1] + in1 [index] [1] * in2 [index] [0];
403int PointMatchAlgorithm::optimizeLengthForFft(
int originalLength)
407 const int INITIAL_CLOSEST_LENGTH = 0;
412 int closestLength = INITIAL_CLOSEST_LENGTH;
413 for (
int power2 = 1; power2 < originalLength; power2 *= 2) {
414 for (
int power3 = 1; power3 < originalLength; power3 *= 3) {
415 for (
int power5 = 1; power5 < originalLength; power5 *= 5) {
416 for (
int power7 = 1; power7 < originalLength; power7 *= 7) {
418 int newLength = power2 * power3 * power5 * power7;
419 if (originalLength <= newLength) {
421 if ((closestLength == INITIAL_CLOSEST_LENGTH) ||
422 (newLength < closestLength)) {
426 closestLength = newLength;
434 if (closestLength == INITIAL_CLOSEST_LENGTH) {
437 closestLength = originalLength;
441 <<
" originalLength=" << originalLength
442 <<
" newLength=" << closestLength;
444 return closestLength;
447void PointMatchAlgorithm::populateImageArray(
const QImage &imageProcessed,
455 ColorFilter colorFilter;
456 for (
int x = 0; x < width; x++) {
457 for (
int y = 0; y < height; y++) {
462 (*image) [
FOLD2DINDEX(x, y, height)] = (pixelIsOn ?
469void PointMatchAlgorithm::populateSampleArray(
const QList<PointMatchPixel> &samplePointPixels,
483 int xMin = width, yMin = height, xMax = 0, yMax = 0;
484 for (i = 0; i < static_cast<unsigned int> (samplePointPixels.size()); i++) {
486 int x = samplePointPixels.at(
signed (i)).xOffset();
487 int y = samplePointPixels.at(
signed (i)).yOffset();
488 if (first || (x < xMin))
490 if (first || (x > xMax))
492 if (first || (y < yMin))
494 if (first || (y > yMax))
500 const int border = 1;
509 for (x = 0; x < width; x++) {
510 for (y = 0; y < height; y++) {
518 double xSumOn = 0, ySumOn = 0, countOn = 0;
520 for (i = 0; i < static_cast<unsigned int> (samplePointPixels.size()); i++) {
523 x = (samplePointPixels.at(
signed (i))).xOffset() - xMin;
524 y = (samplePointPixels.at(
signed (i))).yOffset() - yMin;
528 bool pixelIsOn = samplePointPixels.at(
signed (i)).pixelIsOn();
540 countOn = qMax (1.0, countOn);
541 *sampleXCenter = qFloor (0.5 + xSumOn / countOn);
542 *sampleYCenter = qFloor (0.5 + ySumOn / countOn);
545 *sampleXExtent = xMax - xMin + 1;
546 *sampleYExtent = yMax - yMin + 1;
549void PointMatchAlgorithm::releaseImageArray(
double* array)
557void PointMatchAlgorithm::releasePhaseArray(fftw_complex* arrayPrime)
565void PointMatchAlgorithm::removePixelsNearExistingPoints(
double* image,
568 const Points &pointsExisting,
573 for (
int i = 0; i < pointsExisting.size(); i++) {
575 int xPoint = qFloor (pointsExisting.at(signed (i)).posScreen().x());
576 int yPoint = qFloor (pointsExisting.at(signed (i)).posScreen().y());
579 int yMin = yPoint - pointSeparation;
582 int yMax = yPoint + pointSeparation;
583 if (imageHeight < yMax)
586 for (
int y = yMin; y < yMax; y++) {
589 int radical = pointSeparation * pointSeparation - (y - yPoint) * (y - yPoint);
592 int xMin = qFloor (xPoint - qSqrt(
double (radical)));
595 int xMax = xPoint + (xPoint - xMin);
596 if (imageWidth < xMax)
600 for (
int x = xMin; x < xMax; x++) {
#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
#define FOLD2DINDEX(i, j, jmax)
QList< PointMatchTriplet > PointMatchList
bool pixelFilteredIsOn(const QImage &image, int x, int y) const
Return true if specified filtered pixel is on.
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)
const QString GNUPLOT_FILE_MESSAGE