Engauge Digitizer  2
 All Classes Files Functions Variables Enumerations Enumerator Friends Pages
GridClassifier.cpp
1 #include "ColorFilter.h"
2 #include "Correlation.h"
3 #include "DocumentModelCoords.h"
4 #include "EngaugeAssert.h"
5 #include "GridClassifier.h"
6 #include <iostream>
7 #include "Logger.h"
8 #include <QDebug>
9 #include <QFile>
10 #include <QImage>
11 #include "QtToString.h"
12 #include "Transformation.h"
13 
14 int GridClassifier::NUM_PIXELS_PER_HISTOGRAM_BINS = 1;
15 double GridClassifier::PEAK_HALF_WIDTH = 4;
16 int GridClassifier::MIN_STEP_PIXELS = 4 * GridClassifier::PEAK_HALF_WIDTH; // Step includes down ramp, flat part, up ramp
17 const QString GNUPLOT_DELIMITER ("\t");
18 
19 // We set up the picket fence with binStart arbitrarily set close to zero. Peak is
20 // not exactly at zero since we want to include the left side of the first peak.
21 int GridClassifier::BIN_START_UNSHIFTED = GridClassifier::PEAK_HALF_WIDTH;
22 
23 using namespace std;
24 
26 {
27 }
28 
29 int GridClassifier::binFromCoordinate (double coord,
30  double coordMin,
31  double coordMax) const
32 {
33  ENGAUGE_ASSERT (coordMin < coordMax);
34  ENGAUGE_ASSERT (coordMin <= coord);
35  ENGAUGE_ASSERT (coord <= coordMax);
36 
37  int bin = 0.5 + (m_numHistogramBins - 1.0) * (coord - coordMin) / (coordMax - coordMin);
38 
39  return bin;
40 }
41 
42 void GridClassifier::classify (bool isGnuplot,
43  const QPixmap &originalPixmap,
44  const Transformation &transformation,
45  int &countX,
46  double &startX,
47  double &stepX,
48  int &countY,
49  double &startY,
50  double &stepY)
51 {
52  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::classify";
53 
54  QImage image = originalPixmap.toImage ();
55 
56  m_numHistogramBins = image.width() / NUM_PIXELS_PER_HISTOGRAM_BINS;
57 
58  double xMin, xMax, yMin, yMax;
59  double binStartX, binStepX, binStartY, binStepY;
60 
61  m_binsX = new double [m_numHistogramBins];
62  m_binsY = new double [m_numHistogramBins];
63 
64  computeGraphCoordinateLimits (image,
65  transformation,
66  xMin,
67  xMax,
68  yMin,
69  yMax);
70  initializeHistogramBins ();
71  populateHistogramBins (image,
72  transformation,
73  xMin,
74  xMax,
75  yMin,
76  yMax);
77  searchStartStepSpace (isGnuplot,
78  m_binsX,
79  "x",
80  xMin,
81  xMax,
82  startX,
83  stepX,
84  binStartX,
85  binStepX);
86  searchStartStepSpace (isGnuplot,
87  m_binsY,
88  "y",
89  yMin,
90  yMax,
91  startY,
92  stepY,
93  binStartY,
94  binStepY);
95  searchCountSpace (m_binsX,
96  binStartX,
97  binStepX,
98  countX);
99  searchCountSpace (m_binsY,
100  binStartY,
101  binStepY,
102  countY);
103 
104  delete m_binsX;
105  delete m_binsY;
106 }
107 
108 void GridClassifier::computeGraphCoordinateLimits (const QImage &image,
109  const Transformation &transformation,
110  double &xMin,
111  double &xMax,
112  double &yMin,
113  double &yMax)
114 {
115  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::computeGraphCoordinateLimits";
116 
117  // Project every pixel onto the x axis, and onto the y axis. One set of histogram bins will be
118  // set up along each of the axes. The range of bins will encompass every pixel in the image, and no more.
119 
120  QPointF posGraphTL, posGraphTR, posGraphBL, posGraphBR;
121  transformation.transformScreenToRawGraph (QPointF (0, 0) , posGraphTL);
122  transformation.transformScreenToRawGraph (QPointF (image.width(), 0) , posGraphTR);
123  transformation.transformScreenToRawGraph (QPointF (0, image.height()) , posGraphBL);
124  transformation.transformScreenToRawGraph (QPointF (image.width(), image.height()), posGraphBR);
125 
126  // Compute x and y ranges for setting up the histogram bins
127  if (transformation.modelCoords().coordsType() == COORDS_TYPE_CARTESIAN) {
128 
129  // For affine cartesian coordinates, we only need to look at the screen corners
130  xMin = qMin (qMin (qMin (posGraphTL.x(), posGraphTR.x()), posGraphBL.x()), posGraphBR.x());
131  xMax = qMax (qMax (qMax (posGraphTL.x(), posGraphTR.x()), posGraphBL.x()), posGraphBR.x());
132  yMin = qMin (qMin (qMin (posGraphTL.y(), posGraphTR.y()), posGraphBL.y()), posGraphBR.y());
133  yMax = qMax (qMax (qMax (posGraphTL.y(), posGraphTR.y()), posGraphBL.y()), posGraphBR.y());
134 
135  } else {
136 
137  // For affine polar coordinates, easiest approach is to assume the full circle
138  xMin = 0.0;
139  xMax = transformation.modelCoords().thetaPeriod();
140  yMin = transformation.modelCoords().originRadius();
141  yMax = qMax (qMax (qMax (posGraphTL.y(), posGraphTR.y()), posGraphBL.y()), posGraphBR.y());
142 
143  }
144 
145  ENGAUGE_ASSERT (xMin < xMax);
146  ENGAUGE_ASSERT (yMin < yMax);
147 }
148 
149 double GridClassifier::coordinateFromBin (int bin,
150  double coordMin,
151  double coordMax) const
152 {
153  ENGAUGE_ASSERT (bin < m_numHistogramBins);
154  ENGAUGE_ASSERT (coordMin < coordMax);
155 
156  return coordMin + (coordMax - coordMin) * (double) bin / ((double) m_numHistogramBins - 1.0);
157 }
158 
159 void GridClassifier::copyVectorToVector (const double from [],
160  double to []) const
161 {
162  for (int bin = 0; bin < m_numHistogramBins; bin++) {
163  to [bin] = from [bin];
164  }
165 }
166 
167 void GridClassifier::dumpGnuplotCoordinate (const QString &coordinateLabel,
168  double corr,
169  const double *bins,
170  double coordinateMin,
171  double coordinateMax,
172  int binStart,
173  int binStep) const
174 {
175  QString filename = QString ("gridclassifier_%1_corr%2_startMax%3_stepMax%4.gnuplot")
176  .arg (coordinateLabel)
177  .arg (corr, 8, 'f', 3, '0')
178  .arg (binStart)
179  .arg (binStep);
180 
181  cout << "Writing gnuplot file: " << filename.toLatin1().data() << "\n";
182 
183  QFile fileDump (filename);
184  fileDump.open (QIODevice::WriteOnly | QIODevice::Text);
185  QTextStream strDump (&fileDump);
186 
187  int bin;
188 
189  // For consistent scaling, get the max bin count
190  int binCountMax = 0;
191  for (bin = 0; bin < m_numHistogramBins; bin++) {
192  if (bins [bin] > binCountMax) {
193  binCountMax = qMax ((double) binCountMax,
194  bins [bin]);
195  }
196  }
197 
198  // Get picket fence
199  double *picketFence = new double [m_numHistogramBins];
200  loadPicketFence (picketFence,
201  binStart,
202  binStep,
203  0,
204  false);
205 
206  // Header
207  strDump << "bin"
208  << GNUPLOT_DELIMITER << "coordinate"
209  << GNUPLOT_DELIMITER << "binCount"
210  << GNUPLOT_DELIMITER << "startStep"
211  << GNUPLOT_DELIMITER << "picketFence" << "\n";
212 
213  // Data, with one row per coordinate value
214  for (bin = 0; bin < m_numHistogramBins; bin++) {
215 
216  double coordinate = coordinateFromBin (bin,
217  coordinateMin,
218  coordinateMax);
219  double startStepValue (((bin - binStart) % binStep == 0) ? 1 : 0);
220  strDump << bin
221  << GNUPLOT_DELIMITER << coordinate
222  << GNUPLOT_DELIMITER << bins [bin]
223  << GNUPLOT_DELIMITER << binCountMax * startStepValue
224  << GNUPLOT_DELIMITER << binCountMax * picketFence [bin] << "\n";
225  }
226 
227  delete [] picketFence;
228 }
229 
230 void GridClassifier::dumpGnuplotCorrelations (const QString &coordinateLabel,
231  double valueMin,
232  double valueMax,
233  const double signalA [],
234  const double signalB [],
235  const double correlations [])
236 {
237  QString filename = QString ("gridclassifier_%1_correlations.gnuplot")
238  .arg (coordinateLabel);
239 
240  cout << "Writing gnuplot file: " << filename.toLatin1().data() << "\n";
241 
242  QFile fileDump (filename);
243  fileDump.open (QIODevice::WriteOnly | QIODevice::Text);
244  QTextStream strDump (&fileDump);
245 
246  int bin;
247 
248  // Compute max values so curves can be normalized to the same heights
249  double signalAMax = 1, signalBMax = 1, correlationsMax = 1;
250  for (bin = 0; bin < m_numHistogramBins; bin++) {
251  if (bin == 0 || signalA [bin] > signalAMax) {
252  signalAMax = signalA [bin];
253  }
254  if (bin == 0 || signalB [bin] > signalBMax) {
255  signalBMax = signalB [bin];
256  }
257  if (bin == 0 || correlations [bin] > correlationsMax) {
258  correlationsMax = correlations [bin];
259  }
260  }
261 
262  // Output normalized curves
263  for (int bin = 0; bin < m_numHistogramBins; bin++) {
264 
265  strDump << coordinateFromBin (bin,
266  valueMin,
267  valueMax)
268  << GNUPLOT_DELIMITER << signalA [bin] / signalAMax
269  << GNUPLOT_DELIMITER << signalB [bin] / signalBMax
270  << GNUPLOT_DELIMITER << correlations [bin] / correlationsMax << "\n";
271  }
272 }
273 
274 void GridClassifier::initializeHistogramBins ()
275 {
276  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::initializeHistogramBins";
277 
278  for (int bin = 0; bin < m_numHistogramBins; bin++) {
279  m_binsX [bin] = 0;
280  m_binsY [bin] = 0;
281  }
282 }
283 
284 void GridClassifier::loadPicketFence (double picketFence [],
285  int binStart,
286  int binStep,
287  int count,
288  bool isCount) const
289 {
290  const double PEAK_HEIGHT = 1.0; // Arbitrary height, since normalization will counteract any change to this value
291 
292  // Count argument is optional. Note that binStart already excludes PEAK_HALF_WIDTH bins at start,
293  // and we also exclude PEAK_HALF_WIDTH bins at the end
294  ENGAUGE_ASSERT (binStart >= PEAK_HALF_WIDTH);
295  if (!isCount) {
296  count = 1 + (m_numHistogramBins - binStart - PEAK_HALF_WIDTH) / binStep;
297  }
298 
299  // Bins that we need to worry about
300  int binStartMinusHalfWidth = binStart - PEAK_HALF_WIDTH;
301  int binStopPlusHalfWidth = (binStart + (count - 1) * binStep) + PEAK_HALF_WIDTH;
302 
303  // To normalize, we compute the area under the picket fence. Constraint is
304  // areaUnnormalized + NUM_HISTOGRAM_BINS * normalizationOffset = 0
305  double areaUnnormalized = count * PEAK_HEIGHT * PEAK_HALF_WIDTH;
306  double normalizationOffset = -1.0 * areaUnnormalized / m_numHistogramBins;
307 
308  for (int bin = 0; bin < m_numHistogramBins; bin++) {
309 
310  // Outside of the peaks, bin is small negative so correlation with unit function is zero. In other
311  // words, the function is normalized
312  picketFence [bin] = normalizationOffset;
313 
314  if ((binStartMinusHalfWidth <= bin) &&
315  (bin <= binStopPlusHalfWidth)) {
316 
317  // Closest peak
318  int ordinalClosestPeak = (int) ((bin - binStart + binStep / 2) / binStep);
319  int binClosestPeak = binStart + ordinalClosestPeak * binStep;
320 
321  // Distance from closest peak is used to define an isosceles triangle
322  int distanceToClosestPeak = qAbs (bin - binClosestPeak);
323 
324  if (distanceToClosestPeak < PEAK_HALF_WIDTH) {
325 
326  // Map 0 to PEAK_HALF_WIDTH to 1 to 0
327  picketFence [bin] = 1.0 - (double) distanceToClosestPeak / PEAK_HALF_WIDTH + normalizationOffset;
328 
329  }
330  }
331  }
332 }
333 
334 void GridClassifier::populateHistogramBins (const QImage &image,
335  const Transformation &transformation,
336  double xMin,
337  double xMax,
338  double yMin,
339  double yMax)
340 {
341  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::populateHistogramBins";
342 
343  ColorFilter filter;
344  QRgb rgbBackground = filter.marginColor (&image);
345 
346  for (int x = 0; x < image.width(); x++) {
347  for (int y = 0; y < image.height(); y++) {
348 
349  QColor pixel = image.pixel (x, y);
350 
351  // Skip pixels with background color
352  if (!filter.colorCompare (rgbBackground,
353  pixel.rgb ())) {
354 
355  // Add this pixel to histograms
356  QPointF posGraph;
357  transformation.transformScreenToRawGraph (QPointF (x, y), posGraph);
358 
359  if (transformation.modelCoords().coordsType() == COORDS_TYPE_POLAR) {
360 
361  // If out of the 0 to period range, the theta value must shifted by the period to get into that range
362  while (posGraph.x() < xMin) {
363  posGraph.setX (posGraph.x() + transformation.modelCoords().thetaPeriod());
364  }
365  while (posGraph.x() > xMax) {
366  posGraph.setX (posGraph.x() - transformation.modelCoords().thetaPeriod());
367  }
368  }
369 
370  int binX = binFromCoordinate (posGraph.x(), xMin, xMax);
371  int binY = binFromCoordinate (posGraph.y(), yMin, yMax);
372 
373  ENGAUGE_ASSERT (0 <= binX);
374  ENGAUGE_ASSERT (0 <= binY);
375  ENGAUGE_ASSERT (binX < m_numHistogramBins);
376  ENGAUGE_ASSERT (binY < m_numHistogramBins);
377 
378  // Roundoff error in log scaling may let bin go just outside legal range
379  binX = qMin (binX, m_numHistogramBins - 1);
380  binY = qMin (binY, m_numHistogramBins - 1);
381 
382  ++m_binsX [binX];
383  ++m_binsY [binY];
384  }
385  }
386  }
387 }
388 
389 void GridClassifier::searchCountSpace (double bins [],
390  double binStart,
391  double binStep,
392  int &countMax)
393 {
394  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::searchCountSpace"
395  << " start=" << binStart
396  << " step=" << binStep;
397 
398  // Loop though the space of possible counts
399  Correlation correlation (m_numHistogramBins);
400  double *picketFence = new double [m_numHistogramBins];
401  double corr, corrMax;
402  bool isFirst = true;
403  int countStop = 1 + (m_numHistogramBins - binStart) / binStep;
404  for (int count = 2; count <= countStop; count++) {
405 
406  loadPicketFence (picketFence,
407  binStart,
408  binStep,
409  count,
410  true);
411 
412  correlation.correlateWithoutShift (m_numHistogramBins,
413  bins,
414  picketFence,
415  corr);
416  if (isFirst || (corr > corrMax)) {
417  countMax = count;
418  corrMax = corr;
419  }
420 
421  isFirst = false;
422  }
423 
424  free (picketFence);
425 }
426 
427 void GridClassifier::searchStartStepSpace (bool isGnuplot,
428  double bins [],
429  const QString &coordinateLabel,
430  double valueMin,
431  double valueMax,
432  double &start,
433  double &step,
434  double &binStartMax,
435  double &binStepMax)
436 {
437  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::searchStartStepSpace";
438 
439  // Correlations are tracked for logging
440  double *signalA = new double [m_numHistogramBins];
441  double *signalB = new double [m_numHistogramBins];
442  double *correlations = new double [m_numHistogramBins];
443  double *correlationsMax = new double [m_numHistogramBins];
444 
445  // Loop though the space of possible gridlines using the independent variables (start,step).
446  Correlation correlation (m_numHistogramBins);
447  double *picketFence = new double [m_numHistogramBins];
448  int binStart;
449  double corr, corrMax;
450  bool isFirst = true;
451 
452  // We do not explicitly search(=loop) through binStart here, since Correlation::correlateWithShift will take
453  // care of that for us
454 
455  // Step search starts out small, and stops at value that gives count substantially greater than 2
456  for (int binStep = MIN_STEP_PIXELS; binStep < m_numHistogramBins / 4; binStep++) {
457 
458  loadPicketFence (picketFence,
459  BIN_START_UNSHIFTED,
460  binStep,
461  PEAK_HALF_WIDTH,
462  false);
463 
464  correlation.correlateWithShift (m_numHistogramBins,
465  bins,
466  picketFence,
467  binStart,
468  corr,
469  correlations);
470  if (isFirst || (corr > corrMax)) {
471 
472  binStartMax = binStart + BIN_START_UNSHIFTED + 1; // Compensate for the shift performed inside loadPicketFence
473  binStepMax = binStep;
474  corrMax = corr;
475  copyVectorToVector (bins, signalA);
476  copyVectorToVector (picketFence, signalB);
477  copyVectorToVector (correlations, correlationsMax);
478 
479  // Output a gnuplot file. We should see the correlation values consistently increasing
480  if (isGnuplot) {
481 
482  dumpGnuplotCoordinate(coordinateLabel,
483  corr,
484  bins,
485  valueMin,
486  valueMax,
487  binStart,
488  binStep);
489  }
490  }
491 
492  isFirst = false;
493  }
494 
495  // Convert from bins back to graph coordinates
496  start = coordinateFromBin (binStartMax,
497  valueMin,
498  valueMax);
499  double next = coordinateFromBin (binStartMax + binStepMax,
500  valueMin,
501  valueMax);
502  step = next - start;
503 
504  if (isGnuplot) {
505  dumpGnuplotCorrelations (coordinateLabel,
506  valueMin,
507  valueMax,
508  signalA,
509  signalB,
510  correlationsMax);
511  }
512 
513  free (signalA);
514  free (signalB);
515  free (correlations);
516  free (correlationsMax);
517  free (picketFence);
518 }
void transformScreenToRawGraph(const QPointF &coordScreen, QPointF &coordGraph) const
Transform from cartesian pixel screen coordinates to cartesian/polar graph coordinates.
double originRadius() const
Get method for origin radius in polar mode.
Fast cross correlation between two functions.
Definition: Correlation.h:8
Class for filtering image to remove unimportant information.
Definition: ColorFilter.h:12
DocumentModelCoords modelCoords() const
Get method for DocumentModelCoords.
double thetaPeriod() const
Return the period of the theta value for polar coordinates, consistent with CoordThetaUnits.
Affine transformation between screen and graph coordinates, based on digitized axis points...
QRgb marginColor(const QImage *image) const
Identify the margin color of the image, which is defined as the most common color in the four margins...
Definition: ColorFilter.cpp:52
GridClassifier()
Single constructor.
CoordsType coordsType() const
Get method for coordinates type.
bool colorCompare(QRgb rgb1, QRgb rgb2) const
See if the two color values are close enough to be considered to be the same.
Definition: ColorFilter.cpp:13
Classify the grid pattern in an original image.
void classify(bool isGnuplot, const QPixmap &originalPixmap, const Transformation &transformation, int &countX, double &startX, double &stepX, int &countY, double &startY, double &stepY)
Classify the specified image, and return the most probably x and y grid settings. ...