22#define FOLD2DINDEX(i,j,jmax) ((i)*(jmax)+j)
29 m_isGnuplot (isGnuplot)
34void PointMatchAlgorithm::allocateMemory(
double** array,
35 fftw_complex** arrayPrime,
41 *array =
new double [unsigned (width * height)];
44 *arrayPrime =
new fftw_complex [unsigned (width * height)];
48void PointMatchAlgorithm::assembleLocalMaxima(
double* convolution,
56 const double SINGLE_PIXEL_CORRELATION = 1.0;
58 for (
int i = 0; i < width; i++) {
59 for (
int j = 0; j < height; j++) {
62 double convIJ = log10 (convolution[
FOLD2DINDEX(i, j, height)]);
65 bool isLocalMax =
true;
66 for (
int iDelta = -1; (iDelta <= 1) && isLocalMax; iDelta++) {
68 int iNeighbor = i + iDelta;
69 if ((0 <= iNeighbor) && (iNeighbor < width)) {
71 for (
int jDelta = -1; (jDelta <= 1) && isLocalMax; jDelta++) {
73 int jNeighbor = j + jDelta;
74 if ((0 <= jNeighbor) && (jNeighbor < height)) {
77 double convNeighbor = log10 (convolution[
FOLD2DINDEX(iNeighbor, jNeighbor, height)]);
78 if (convIJ < convNeighbor) {
82 }
else if (convIJ == convNeighbor) {
85 if ((jDelta < 0) || (jDelta == 0 && iDelta < 0)) {
96 (convIJ > SINGLE_PIXEL_CORRELATION) ) {
99 PointMatchTriplet t (i,
103 listCreated.append(t);
109void PointMatchAlgorithm::computeConvolution(fftw_complex* imagePrime,
110 fftw_complex* samplePrime,
111 int width,
int height,
112 double** convolution,
118 fftw_complex* convolutionPrime;
120 allocateMemory(convolution,
126 conjugateMatrix(width,
131 multiplyMatrices(width,
138 fftw_plan pConvolution = fftw_plan_dft_c2r_2d(width,
144 fftw_execute (pConvolution);
146 releasePhaseArray(convolutionPrime);
150 double *temp =
new double [unsigned (width * height)];
153 for (
int i = 0; i < width; i++) {
154 for (
int j = 0; j < height; j++) {
158 for (
int iFrom = 0; iFrom < width; iFrom++) {
159 for (
int jFrom = 0; jFrom < height; jFrom++) {
161 int iTo = (iFrom + sampleXCenter) % width;
162 int jTo = (jFrom + sampleYCenter) % height;
169void PointMatchAlgorithm::conjugateMatrix(
int width,
171 fftw_complex* matrix)
177 for (
int x = 0; x < width; x++) {
178 for (
int y = 0; y < height; y++) {
181 matrix [index] [1] = -1.0 * matrix [index] [1];
186void PointMatchAlgorithm::dumpToGnuplot (
double* convolution,
189 const QString &filename)
const
195 QFile file (filename);
196 if (file.open (QIODevice::WriteOnly | QIODevice::Text)) {
198 QTextStream str (&file);
206 for (
int i = 0; i < width; i++) {
207 for (
int j = 0; j < height; j++) {
209 double convIJ = convolution[
FOLD2DINDEX(i, j, height)];
220 const QImage &imageProcessed,
222 const Points &pointsExisting)
225 <<
" samplePointPixels=" << samplePointPixels.count();
228 int originalWidth = imageProcessed.width();
229 int originalHeight = imageProcessed.height();
230 int width = optimizeLengthForFft(originalWidth);
231 int height = optimizeLengthForFft(originalHeight);
235 double *image, *sample, *convolution;
236 fftw_complex *imagePrime, *samplePrime;
239 int sampleXCenter, sampleYCenter, sampleXExtent, sampleYExtent;
240 loadImage(imageProcessed,
247 loadSample(samplePointPixels,
256 computeConvolution(imagePrime,
270 dumpToGnuplot(sample,
274 dumpToGnuplot(convolution,
277 "convolution.gnuplot");
283 assembleLocalMaxima(convolution,
287 std::sort (listCreated.begin(),
291 QList<QPoint> pointsCreated;
292 PointMatchList::iterator itr;
293 for (itr = listCreated.begin(); itr != listCreated.end(); itr++) {
296 pointsCreated.push_back (triplet.
point ());
304 releaseImageArray(image);
305 releasePhaseArray(imagePrime);
306 releaseImageArray(sample);
307 releasePhaseArray(samplePrime);
308 releaseImageArray(convolution);
310 return pointsCreated;
313void PointMatchAlgorithm::loadImage(
const QImage &imageProcessed,
315 const Points &pointsExisting,
319 fftw_complex** imagePrime)
323 allocateMemory(image,
328 populateImageArray(imageProcessed,
333 removePixelsNearExistingPoints(*image,
340 fftw_plan pImage = fftw_plan_dft_r2c_2d(width,
345 fftw_execute(pImage);
348void PointMatchAlgorithm::loadSample(
const QList<PointMatchPixel> &samplePointPixels,
352 fftw_complex** samplePrime,
362 allocateMemory(sample,
367 populateSampleArray(samplePointPixels,
377 fftw_plan pSample = fftw_plan_dft_r2c_2d(width,
382 fftw_execute(pSample);
385void PointMatchAlgorithm::multiplyMatrices(
int width,
393 for (
int x = 0; x < width; x++) {
394 for (
int y = 0; y < height; y++) {
398 out [index] [0] = in1 [index] [0] * in2 [index] [0] - in1 [index] [1] * in2 [index] [1];
399 out [index] [1] = in1 [index] [0] * in2 [index] [1] + in1 [index] [1] * in2 [index] [0];
404int PointMatchAlgorithm::optimizeLengthForFft(
int originalLength)
408 const int INITIAL_CLOSEST_LENGTH = 0;
413 int closestLength = INITIAL_CLOSEST_LENGTH;
414 for (
int power2 = 1; power2 < originalLength; power2 *= 2) {
415 for (
int power3 = 1; power3 < originalLength; power3 *= 3) {
416 for (
int power5 = 1; power5 < originalLength; power5 *= 5) {
417 for (
int power7 = 1; power7 < originalLength; power7 *= 7) {
419 int newLength = power2 * power3 * power5 * power7;
420 if (originalLength <= newLength) {
422 if ((closestLength == INITIAL_CLOSEST_LENGTH) ||
423 (newLength < closestLength)) {
427 closestLength = newLength;
435 if (closestLength == INITIAL_CLOSEST_LENGTH) {
438 closestLength = originalLength;
442 <<
" originalLength=" << originalLength
443 <<
" newLength=" << closestLength;
445 return closestLength;
448void PointMatchAlgorithm::populateImageArray(
const QImage &imageProcessed,
456 ColorFilter colorFilter;
457 for (
int x = 0; x < width; x++) {
458 for (
int y = 0; y < height; y++) {
463 (*image) [
FOLD2DINDEX(x, y, height)] = (pixelIsOn ?
470void PointMatchAlgorithm::populateSampleArray(
const QList<PointMatchPixel> &samplePointPixels,
484 int xMin = width, yMin = height, xMax = 0, yMax = 0;
485 for (i = 0; i < static_cast<unsigned int> (samplePointPixels.size()); i++) {
487 int x = samplePointPixels.at(
signed (i)).xOffset();
488 int y = samplePointPixels.at(
signed (i)).yOffset();
489 if (first || (x < xMin))
491 if (first || (x > xMax))
493 if (first || (y < yMin))
495 if (first || (y > yMax))
501 const int border = 1;
510 for (x = 0; x < width; x++) {
511 for (y = 0; y < height; y++) {
519 double xSumOn = 0, ySumOn = 0, countOn = 0;
521 for (i = 0; i < static_cast<unsigned int> (samplePointPixels.size()); i++) {
524 x = (samplePointPixels.at(
signed (i))).xOffset() - xMin;
525 y = (samplePointPixels.at(
signed (i))).yOffset() - yMin;
529 bool pixelIsOn = samplePointPixels.at(
signed (i)).pixelIsOn();
541 countOn = qMax (1.0, countOn);
542 *sampleXCenter = qFloor (0.5 + xSumOn / countOn);
543 *sampleYCenter = qFloor (0.5 + ySumOn / countOn);
546 *sampleXExtent = xMax - xMin + 1;
547 *sampleYExtent = yMax - yMin + 1;
550void PointMatchAlgorithm::releaseImageArray(
double* array)
558void PointMatchAlgorithm::releasePhaseArray(fftw_complex* arrayPrime)
566void PointMatchAlgorithm::removePixelsNearExistingPoints(
double* image,
569 const Points &pointsExisting,
574 for (
int i = 0; i < pointsExisting.size(); i++) {
576 int xPoint = qFloor (pointsExisting.at(signed (i)).posScreen().x());
577 int yPoint = qFloor (pointsExisting.at(signed (i)).posScreen().y());
580 int yMin = yPoint - pointSeparation;
583 int yMax = yPoint + pointSeparation;
584 if (imageHeight < yMax)
587 for (
int y = yMin; y < yMax; y++) {
590 int radical = pointSeparation * pointSeparation - (y - yPoint) * (y - yPoint);
593 int xMin = qFloor (xPoint - qSqrt(
double (radical)));
596 int xMax = xPoint + (xPoint - xMin);
597 if (imageWidth < xMax)
601 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.
static QTextStream & endl(QTextStream &stream)
End of line.
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