Engauge Digitizer  2
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 
18 using namespace std;
19 
21 {
22 }
23 
24 int GridClassifier::binFromCoordinate (double coord,
25  double coordMin,
26  double coordMax) const
27 {
28  ENGAUGE_ASSERT (coordMin < coordMax);
29  ENGAUGE_ASSERT (coordMin <= coord);
30  ENGAUGE_ASSERT (coord <= coordMax);
31 
32  int bin = 0.5 + (m_numHistogramBins - 1.0) * (coord - coordMin) / (coordMax - coordMin);
33 
34  return bin;
35 }
36 
37 void GridClassifier::classify (bool isGnuplot,
38  const QPixmap &originalPixmap,
39  const Transformation &transformation,
40  int &countX,
41  double &startX,
42  double &stepX,
43  int &countY,
44  double &startY,
45  double &stepY)
46 {
47  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::classify";
48 
49  QImage image = originalPixmap.toImage ();
50 
51  m_numHistogramBins = image.width() / NUM_PIXELS_PER_HISTOGRAM_BINS;
52 
53  double xMin, xMax, yMin, yMax;
54  double binStartX, binStepX, binStartY, binStepY;
55 
56  m_binsX = new double [m_numHistogramBins];
57  m_binsY = new double [m_numHistogramBins];
58 
59  computeGraphCoordinateLimits (image,
60  transformation,
61  xMin,
62  xMax,
63  yMin,
64  yMax);
65  initializeHistogramBins ();
66  populateHistogramBins (image,
67  transformation,
68  xMin,
69  xMax,
70  yMin,
71  yMax);
72  searchStartStepSpace (isGnuplot,
73  xMin,
74  xMax,
75  yMin,
76  yMax,
77  startX,
78  stepX,
79  startY,
80  stepY,
81  binStartX,
82  binStepX,
83  binStartY,
84  binStepY);
85  searchCountSpace (m_binsX,
86  binStartX,
87  binStepX,
88  countX);
89  searchCountSpace (m_binsY,
90  binStartY,
91  binStepY,
92  countY);
93 
94  delete m_binsX;
95  delete m_binsY;
96 }
97 
98 void GridClassifier::computeGraphCoordinateLimits (const QImage &image,
99  const Transformation &transformation,
100  double &xMin,
101  double &xMax,
102  double &yMin,
103  double &yMax)
104 {
105  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::computeGraphCoordinateLimits";
106 
107  // Project every pixel onto the x axis, and onto the y axis. One set of histogram bins will be
108  // set up along each of the axes. The range of bins will encompass every pixel in the image, and no more.
109 
110  QPointF posGraphTL, posGraphTR, posGraphBL, posGraphBR;
111  transformation.transformScreenToRawGraph (QPointF (0, 0) , posGraphTL);
112  transformation.transformScreenToRawGraph (QPointF (image.width(), 0) , posGraphTR);
113  transformation.transformScreenToRawGraph (QPointF (0, image.height()) , posGraphBL);
114  transformation.transformScreenToRawGraph (QPointF (image.width(), image.height()), posGraphBR);
115 
116  // Compute x and y ranges for setting up the histogram bins
117  if (transformation.modelCoords().coordsType() == COORDS_TYPE_CARTESIAN) {
118 
119  // For affine cartesian coordinates, we only need to look at the screen corners
120  xMin = qMin (qMin (qMin (posGraphTL.x(), posGraphTR.x()), posGraphBL.x()), posGraphBR.x());
121  xMax = qMax (qMax (qMax (posGraphTL.x(), posGraphTR.x()), posGraphBL.x()), posGraphBR.x());
122  yMin = qMin (qMin (qMin (posGraphTL.y(), posGraphTR.y()), posGraphBL.y()), posGraphBR.y());
123  yMax = qMax (qMax (qMax (posGraphTL.y(), posGraphTR.y()), posGraphBL.y()), posGraphBR.y());
124 
125  } else {
126 
127  // For affine polar coordinates, easiest approach is to assume the full circle
128  xMin = 0.0;
129  xMax = transformation.modelCoords().thetaPeriod();
130  yMin = transformation.modelCoords().originRadius();
131  yMax = qMax (qMax (qMax (posGraphTL.y(), posGraphTR.y()), posGraphBL.y()), posGraphBR.y());
132 
133  }
134 
135  ENGAUGE_ASSERT (xMin < xMax);
136  ENGAUGE_ASSERT (yMin < yMax);
137 }
138 
139 double GridClassifier::coordinateFromBin (int bin,
140  double coordMin,
141  double coordMax) const
142 {
143  ENGAUGE_ASSERT (bin < m_numHistogramBins);
144  ENGAUGE_ASSERT (coordMin < coordMax);
145 
146  return coordMin + (coordMax - coordMin) * (double) bin / ((double) m_numHistogramBins - 1.0);
147 }
148 
149 void GridClassifier::dumpGnuplotCoordinate (const QString &filename,
150  const double *bins,
151  double coordinateMin,
152  double coordinateMax,
153  int binStart,
154  int binStep) const
155 {
156  cout << "Writing gnuplot file: " << filename.toLatin1().data() << "\n";
157 
158  QFile fileDump (filename);
159  fileDump.open (QIODevice::WriteOnly | QIODevice::Text);
160  QTextStream strDump (&fileDump);
161 
162  int bin;
163  const QString DELIMITER ("\t");
164 
165  // For consistent scaling, get the max bin count
166  int binCountMax = 0;
167  for (bin = 0; bin < m_numHistogramBins; bin++) {
168  if (bins [bin] > binCountMax) {
169  binCountMax = qMax ((double) binCountMax,
170  bins [bin]);
171  }
172  }
173 
174  // Get picket fence
175  double *picketFence = new double [m_numHistogramBins];
176  loadPicketFence (picketFence,
177  binStart,
178  binStep,
179  0,
180  false);
181 
182  // Header
183  strDump << "bin"
184  << DELIMITER << "coordinate"
185  << DELIMITER << "binCount"
186  << DELIMITER << "startStep"
187  << DELIMITER << "picketFence" << "\n";
188 
189  // Data, with one row per coordinate value
190  for (bin = 0; bin < m_numHistogramBins; bin++) {
191 
192  double coordinate = coordinateFromBin (bin,
193  coordinateMin,
194  coordinateMax);
195  double startStepValue (((bin - binStart) % binStep == 0) ? 1 : 0);
196  strDump << bin
197  << DELIMITER << coordinate
198  << DELIMITER << bins [bin]
199  << DELIMITER << binCountMax * startStepValue
200  << DELIMITER << binCountMax * picketFence [bin] << "\n";
201  }
202 
203  delete picketFence;
204 }
205 
206 void GridClassifier::initializeHistogramBins ()
207 {
208  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::initializeHistogramBins";
209 
210  for (int bin = 0; bin < m_numHistogramBins; bin++) {
211  m_binsX [bin] = 0;
212  m_binsY [bin] = 0;
213  }
214 }
215 
216 void GridClassifier::loadPicketFence (double picketFence [],
217  int binStart,
218  int binStep,
219  int count,
220  bool isCount) const
221 {
222  const double PEAK_HEIGHT = 1.0; // Arbitrary height, since normalization will counteract any change to this value
223 
224  // Count argument is optional. Note that binStart already excludes PEAK_HALF_WIDTH bins at start,
225  // and we also exclude PEAK_HALF_WIDTH bins at the end
226  ENGAUGE_ASSERT (binStart >= PEAK_HALF_WIDTH);
227  if (!isCount) {
228  count = 1 + (m_numHistogramBins - binStart - PEAK_HALF_WIDTH) / binStep;
229  }
230 
231  // Bins that we need to worry about
232  int binStartMinusHalfWidth = binStart - PEAK_HALF_WIDTH;
233  int binStopPlusHalfWidth = (binStart + (count - 1) * binStep) + PEAK_HALF_WIDTH;
234 
235  // To normalize, we compute the area under the picket fence. Constraint is
236  // areaUnnormalized + NUM_HISTOGRAM_BINS * normalizationOffset = 0
237  double areaUnnormalized = count * PEAK_HEIGHT * PEAK_HALF_WIDTH;
238  double normalizationOffset = -1.0 * areaUnnormalized / m_numHistogramBins;
239 
240  for (int bin = 0; bin < m_numHistogramBins; bin++) {
241 
242  // Outside of the peaks, bin is small negative so correlation with unit function is zero. In other
243  // words, the function is normalized
244  picketFence [bin] = normalizationOffset;
245 
246  if ((binStartMinusHalfWidth <= bin) &&
247  (bin <= binStopPlusHalfWidth)) {
248 
249  // Closest peak
250  int ordinalClosestPeak = (int) ((bin - binStart + binStep / 2) / binStep);
251  int binClosestPeak = binStart + ordinalClosestPeak * binStep;
252 
253  // Distance from closest peak is used to define an isosceles triangle
254  int distanceToClosestPeak = qAbs (bin - binClosestPeak);
255 
256  if (distanceToClosestPeak < PEAK_HALF_WIDTH) {
257 
258  // Map 0 to PEAK_HALF_WIDTH to 1 to 0
259  picketFence [bin] = 1.0 - (double) distanceToClosestPeak / PEAK_HALF_WIDTH + normalizationOffset;
260 
261  }
262  }
263  }
264 }
265 
266 void GridClassifier::populateHistogramBins (const QImage &image,
267  const Transformation &transformation,
268  double xMin,
269  double xMax,
270  double yMin,
271  double yMax)
272 {
273  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::populateHistogramBins";
274 
275  ColorFilter filter;
276  QRgb rgbBackground = filter.marginColor (&image);
277 
278  for (int x = 0; x < image.width(); x++) {
279  for (int y = 0; y < image.height(); y++) {
280 
281  QColor pixel = image.pixel (x, y);
282 
283  // Skip pixels with background color
284  if (!filter.colorCompare (rgbBackground,
285  pixel.rgb ())) {
286 
287  // Add this pixel to histograms
288  QPointF posGraph;
289  transformation.transformScreenToRawGraph (QPointF (x, y), posGraph);
290 
291  if (transformation.modelCoords().coordsType() == COORDS_TYPE_POLAR) {
292 
293  // If out of the 0 to period range, the theta value must shifted by the period to get into that range
294  while (posGraph.x() < xMin) {
295  posGraph.setX (posGraph.x() + transformation.modelCoords().thetaPeriod());
296  }
297  while (posGraph.x() > xMax) {
298  posGraph.setX (posGraph.x() - transformation.modelCoords().thetaPeriod());
299  }
300  }
301 
302  int binX = binFromCoordinate (posGraph.x(), xMin, xMax);
303  int binY = binFromCoordinate (posGraph.y(), yMin, yMax);
304 
305  ENGAUGE_ASSERT (0 <= binX);
306  ENGAUGE_ASSERT (0 <= binY);
307  ENGAUGE_ASSERT (binX < m_numHistogramBins);
308  ENGAUGE_ASSERT (binY < m_numHistogramBins);
309 
310  // Roundoff error in log scaling may let bin go just outside legal range
311  binX = qMin (binX, m_numHistogramBins - 1);
312  binY = qMin (binY, m_numHistogramBins - 1);
313 
314  ++m_binsX [binX];
315  ++m_binsY [binY];
316  }
317  }
318  }
319 }
320 
321 void GridClassifier::searchCountSpace (double bins [],
322  double binStart,
323  double binStep,
324  int &countMax)
325 {
326  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::searchCountSpace"
327  << " start=" << binStart
328  << " step=" << binStep;
329 
330  // Loop though the space of possible counts
331  Correlation correlation (m_numHistogramBins);
332  double picketFence [m_numHistogramBins];
333  double corr, corrMax;
334  bool isFirst = true;
335  int countStop = 1 + (m_numHistogramBins - binStart) / binStep;
336  for (int count = 2; count <= countStop; count++) {
337 
338  loadPicketFence (picketFence,
339  binStart,
340  binStep,
341  count,
342  true);
343 
344  correlation.correlateWithoutShift (m_numHistogramBins,
345  bins,
346  picketFence,
347  corr);
348  if (isFirst || (corr > corrMax)) {
349  countMax = count;
350  corrMax = corr;
351  }
352 
353  isFirst = false;
354  }
355 }
356 
357 void GridClassifier::searchStartStepSpace (bool isGnuplot,
358  double xMin,
359  double xMax,
360  double yMin,
361  double yMax,
362  double &startX,
363  double &stepX,
364  double &startY,
365  double &stepY,
366  double &binStartXMax,
367  double &binStepXMax,
368  double &binStartYMax,
369  double &binStepYMax)
370 {
371  LOG4CPP_INFO_S ((*mainCat)) << "GridClassifier::searchStartStepSpace";
372 
373  // Loop though the space of possible gridlines using the independent variables (start,step).
374  Correlation correlation (m_numHistogramBins);
375  double picketFence [m_numHistogramBins];
376  int binStartX, binStartY;
377  double corrX, corrY, corrXMax, corrYMax;
378  bool isFirst = true;
379 
380  // We do not explicitly search(=loop) through binStart here, since Correlation::correlateWithShift will take
381  // care of that for us. So, we set up the picket fence with binStart arbitrarily set to zero
382  const int BIN_START_UNSHIFTED = PEAK_HALF_WIDTH;
383 
384  // Step search starts out small, and stops at value that gives count substantially greater than 2
385  for (int binStep = MIN_STEP_PIXELS; binStep < m_numHistogramBins / 4; binStep++) {
386 
387  loadPicketFence (picketFence,
388  BIN_START_UNSHIFTED,
389  binStep,
390  PEAK_HALF_WIDTH,
391  false);
392 
393  correlation.correlateWithShift (m_numHistogramBins,
394  m_binsX,
395  picketFence,
396  binStartX,
397  corrX);
398  correlation.correlateWithShift (m_numHistogramBins,
399  m_binsY,
400  picketFence,
401  binStartY,
402  corrY);
403  if (isFirst || (corrX > corrXMax)) {
404  binStartXMax = binStartX;
405  binStepXMax = binStep;
406  corrXMax = corrX;
407 
408  // Output a gnuplot file. We should see the correlation values consistently increasing
409  if (isGnuplot) {
410  QString filenameGnuplotX = QString ("gridclassifier_x_corr%1_startMax%2_stepMax%3.gnuplot")
411  .arg (corrX, 8, 'f', 3, '0')
412  .arg (binStartXMax)
413  .arg (binStepXMax);
414  dumpGnuplotCoordinate(filenameGnuplotX,
415  m_binsX,
416  xMin,
417  xMax,
418  binStartX,
419  binStep);
420  }
421  }
422 
423  if (isFirst || (corrY > corrYMax)) {
424  binStartYMax = binStartY;
425  binStepYMax = binStep;
426  corrYMax = corrY;
427 
428  // Output a gnuplot file. We should see the correlation values consistently increasing
429  if (isGnuplot) {
430  QString filenameGnuplotY = QString ("gridclassifier_y_corr%1_startMax%2_stepMax%3.gnuplot")
431  .arg (corrY, 8, 'f', 3, '0')
432  .arg (binStartYMax)
433  .arg (binStepYMax);
434  dumpGnuplotCoordinate(filenameGnuplotY,
435  m_binsY,
436  yMin,
437  yMax,
438  binStartY,
439  binStep);
440  }
441  }
442 
443  isFirst = false;
444  }
445 
446  // Convert from bins back to graph coordinates
447  startX = coordinateFromBin (binStartXMax,
448  xMin,
449  xMax);
450  startY = coordinateFromBin (binStartYMax,
451  yMin,
452  yMax);
453  double nextX = coordinateFromBin (binStartXMax + binStepXMax,
454  xMin,
455  xMax);
456  double nextY = coordinateFromBin (binStartYMax + binStepYMax,
457  yMin,
458  yMax);
459 
460  stepX = nextX - startX;
461  stepY = nextY - startY;
462 }
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:7
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. ...