Engauge Digitizer  2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mmsubs.cpp
Go to the documentation of this file.
1 /******************************************************************************************************
2  * (C) 2014 markummitchell@github.com. This file is part of Engauge Digitizer, which is released *
3  * under GNU General Public License version 2 (GPLv2) or (at your option) any later version. See file *
4  * LICENSE or go to gnu.org/licenses for details. Distribution requires prior written permission. *
5  ******************************************************************************************************/
6 
7 #include "mmsubs.h"
8 #include <QImage>
9 #include <QPointF>
10 #include <qmath.h>
11 #include <QTransform>
12 
13 const double PI = 3.1415926535;
14 
15 double angleBetweenVectors (const QPointF &v1,
16  const QPointF &v2)
17 {
18  double v1Mag = qSqrt (v1.x() * v1.x() + v1.y() * v1.y());
19  double v2Mag = qSqrt (v2.x() * v2.x() + v2.y() * v2.y());
20 
21  double angle = 0;
22  if ((v1Mag > 0) || (v2Mag > 0)) {
23 
24  double cosArg = (v1.x() * v2.x() + v1.y() * v2.y()) / (v1Mag * v2Mag);
25  cosArg = qMin (qMax (cosArg, -1.0), 1.0);
26  angle = qAcos (cosArg);
27  }
28 
29  return angle;
30 }
31 
32 double angleFromVectorToVector (const QPointF &vFrom,
33  const QPointF &vTo)
34 {
35  double angleFrom = qAtan2 (vFrom.y(), vFrom.x());
36  double angleTo = qAtan2 (vTo.y() , vTo.x());
37 
38  // Rotate both angles to put from vector along x axis. Note that angleFrom-angleFrom is zero,
39  // and angleTo-angleFrom is -pi to +pi radians
40  double angleSeparation = angleTo - angleFrom;
41 
42  while (angleSeparation < -1.0 * PI) {
43  angleSeparation += 2.0 * PI;
44  }
45  while (angleSeparation > PI) {
46  angleSeparation -= 2.0 * PI;
47  }
48 
49  return angleSeparation;
50 }
51 
52 void ellipseFromParallelogram (double xTL,
53  double yTL,
54  double xTR,
55  double yTR,
56  double xBR,
57  double yBR,
58  double &angleRadians,
59  double &aAligned,
60  double &bAligned)
61 {
62  // Given input describing a parallelogram centered (for simplicity) on the origin,
63  // with three successive corners at (xTL,yTL) and (xTR,yTR) and two other implicit corners
64  // given by symmetry at (-xTL,-yTL) and (-xTR,-yTR), this computes the inscribed ellipse
65  // and returns it as an angle rotation of an axis-aligned ellipse with (x/a)^2+(y/b)^2=1,
66  // also centered on the origin. Great reference is arxiv:0808.0297v1 by A Horwitz.
67  //
68  // Translations to/from the origin must be handled externally since they would
69  // dramatically increase the complexity of the code in this function, and they are
70  // so easily handled before/after calling this method.
71  //
72  // Ellipse will go through midpoints of the 4 parallelogram sides. Interestingly,
73  // the resulting ellipse will NOT in general be aligned with the axes
74 
75  double xT = (xTL + xTR) / 2.0;
76  double yT = (yTL + yTR) / 2.0;
77  double xR = (xTR + xBR) / 2.0;
78  double yR = (yTR + yBR) / 2.0;
79  //
80  // Math:
81  // Ellipse equation: A x^2 + B y^2 + 2 C x y + D x + E y + F = 0
82  //
83  // 6 equations and 6 unknowns (A,B,C,D,E,F)
84  // 1) Passes through (xT,yT)
85  // 2) Passes through (xR,yR)
86  // 3) Slope at top midpoint is the constant mT which comes from y = mT x + bT
87  // 4) Slope at right midpoint is the constant mR which comes from y = mR x + bR
88  // 5+6) Symmetry through the origin means replacing (x,y) by (-x,-y) should leave the
89  // ellipse equation unchanged.
90  // A x^2 + B y^2 + 2 C x y + D x + E y + F = A (-x)^2 + B (-y)^2 + 2 C (-x) (-y) + D (-x) + E (-y) + F
91  // D x + E y = D (-x) + E (-y)
92  // which implies both D and E are zero. This one equation counts as two equations (one for x and one for y)
93 
94  // Taking differentials of ellipse equation
95  // dx (A x + C y) + dy (B y + C x) = 0
96  // dy/dx = -(A x + C y) / (B y + C x) = m
97  // This gives the following set of 4 equations for 4 unknowns
98  // A xT^2 + B yT^2 + 2 C xT yT + F = 0
99  // A xR^2 + B yR^2 + 2 C xR yR + F = 0
100  // A xT + B mT yT + C (mT xT + yT) = 0
101  // A xR + B mR yR + C (mR xR + yR) = 0
102  // but we need to move terms without x and y to the right or else A=B=C=F=0 will be the solution.
103  // At this point we realize that scaling is arbitrary, so divide out F
104  // A xT^2 + B yT^2 + 2 C xT yT = -1
105  // A xR^2 + B yR^2 + 2 C xR yR = -1
106  // A xT + B mT yT + C (mT xT + yT) = 0
107  // Then we apply Kramers Rule to solve for A, B and C
108 
109  double m00 = xT * xT;
110  double m01 = yT * yT;
111  double m02 = 2.0 * xT * yT;
112  double m10 = xR * xR;
113  double m11 = yR * yR;
114  double m12 = 2.0 * xR * yR;
115  double m20 = 0;
116  double m21 = 0;
117  double m22 = 0;
118  // We pick either the top or right side, whichever has the smaller slope to prevent divide by zero error
119  // |mT| = |yTR - yTL| / |xTR - xTL| versus |mR| = |yTR - yBR| / |xTR - xBR|
120  // |yTR - yTL| * |xTR - xBR| versus |yTR - yBR| * |xTR - xTL|
121  if (qAbs (yTR - yTL) * qAbs (xTR - xBR) < qAbs (yTR - yBR) * qAbs (xTR - xTL)) {
122  // Top slope is less so we use it
123  double mT = (yTR - yTL) / (xTR - xTL);
124  //double bT = yTL - mT * xTL;
125  m20 = xT;
126  m21 = mT * yT;
127  m22 = mT * xT + yT;
128  } else {
129  // Right slope is less so we use it
130  double mR = (yTR - yBR) / (xTR - xBR);
131  //double bR = yTR - mR * xTR;
132  m20 = xR;
133  m21 = mR * yR;
134  m22 = mR * xR + yR;
135  }
136 
137  QTransform denominator (m00, m01, m02,
138  m10, m11, m12,
139  m20, m21, m22);
140  QTransform numeratorA (-1.0, m01, m02,
141  -1.0, m11, m12,
142  0.0 , m21, m22);
143  QTransform numeratorB (m00, -1.0, m02,
144  m10, -1.0, m12,
145  m20, 0.0 , m22);
146  QTransform numeratorC (m00, m01, -1.0,
147  m10, m11, -1.0,
148  m20, m21, 0.0 );
149  double A = numeratorA.determinant () / denominator.determinant ();
150  double B = numeratorB.determinant () / denominator.determinant ();
151  double C = numeratorC.determinant () / denominator.determinant ();
152  double F = 1.0;
153 
154  // Semimajor and semiminor axes are from equations 1.1 and 1.2 in the arXiv reference, with
155  // D and E terms set to zero
156  double numerator = 4.0 * F * C * C - 4.0 * A * B * F;
157  double denominatorMinus = 2.0 * (A * B - C * C) * (A + B + qSqrt ((B - A) * (B - A) + 4 * C * C));
158  double denominatorPlus = 2.0 * (A * B - C * C) * (A + B - qSqrt ((B - A) * (B - A) + 4 * C * C));
159  aAligned = qSqrt (numerator / denominatorMinus);
160  bAligned = qSqrt (numerator / denominatorPlus);
161  // Angle is from equation 1.3 in the arXiv reference
162  if (qAbs (2.0 * C) > 10000.0 * qAbs (A - B)) {
163  angleRadians = 90.0;
164  } else {
165  angleRadians = 0.5 * qAtan (2 * C / (A - B));
166  }
167 }
168 
169 QRgb pixelRGB(const QImage &image, int x, int y)
170 {
171  switch (image.depth())
172  {
173  case 1:
174  return pixelRGB1(image, x, y);
175  case 8:
176  return pixelRGB8(image, x, y);
177  default:
178  return pixelRGB32(image, x, y);
179  }
180 }
181 
182 QRgb pixelRGB1(const QImage &image1Bit, int x, int y)
183 {
184  unsigned int bit;
185  if (image1Bit.format () == QImage::Format_MonoLSB) {
186  bit = *(image1Bit.scanLine (y) + (x >> 3)) & (1 << (x & 7));
187  } else {
188  bit = *(image1Bit.scanLine (y) + (x >> 3)) & (1 << (7 - (x & 7)));
189  }
190 
191  int tableIndex = ((bit == 0) ? 0 : 1);
192  return image1Bit.color(tableIndex);
193 }
194 
195 QRgb pixelRGB8(const QImage &image8Bit, int x, int y)
196 {
197  int tableIndex = *(image8Bit.scanLine(y) + x);
198  return image8Bit.color(tableIndex);
199 }
200 
201 QRgb pixelRGB32(const QImage &image32Bit, int x, int y)
202 {
203  // QImage::scanLine documentation says:
204  // 1) Cast return value to QRgb* (which is 32 bit) which hides platform-specific byte orders
205  // 2) Scanline data is aligned on 32 bit boundary
206  // unsigned int* p = (unsigned int *) image32Bit.scanLine(y) + x;
207  const QRgb *p = reinterpret_cast<const QRgb *> (image32Bit.scanLine(y)) + x;
208  return *p;
209 }
210 
211 void projectPointOntoLine(double xToProject,
212  double yToProject,
213  double xStart,
214  double yStart,
215  double xStop,
216  double yStop,
217  double *xProjection,
218  double *yProjection,
219  double *projectedDistanceOutsideLine,
220  double *distanceToLine)
221 {
222  double s;
223  if (qAbs (yStart - yStop) > qAbs (xStart - xStop)) {
224 
225  // More vertical than horizontal. Compute slope and intercept of y=slope*x+yintercept line through (xToProject, yToProject)
226  double slope = (xStop - xStart) / (yStart - yStop);
227  double yintercept = yToProject - slope * xToProject;
228 
229  // Intersect projected point line (slope-intercept form) with start-stop line (parametric form x=(1-s)*x1+s*x2, y=(1-s)*y1+s*y2)
230  s = (slope * xStart + yintercept - yStart) /
231  (yStop - yStart + slope * (xStart - xStop));
232 
233  } else {
234 
235  // More horizontal than vertical. Compute slope and intercept of x=slope*y+xintercept line through (xToProject, yToProject)
236  double slope = (yStop - yStart) / (xStart - xStop);
237  double xintercept = xToProject - slope * yToProject;
238 
239  // Intersect projected point line (slope-intercept form) with start-stop line (parametric form x=(1-s)*x1+s*x2, y=(1-s)*y1+s*y2)
240  s = (slope * yStart + xintercept - xStart) /
241  (xStop - xStart + slope * (yStart - yStop));
242 
243  }
244 
245  *xProjection = (1.0 - s) * xStart + s * xStop;
246  *yProjection = (1.0 - s) * yStart + s * yStop;
247 
248  if (s < 0) {
249 
250  *projectedDistanceOutsideLine = qSqrt ((*xProjection - xStart) * (*xProjection - xStart) +
251  (*yProjection - yStart) * (*yProjection - yStart));
252  *distanceToLine = qSqrt ((xToProject - xStart) * (xToProject - xStart) +
253  (yToProject - yStart) * (yToProject - yStart));
254 
255  // Bring projection point to inside line
256  *xProjection = xStart;
257  *yProjection = yStart;
258 
259  } else if (s > 1) {
260 
261  *projectedDistanceOutsideLine = qSqrt ((*xProjection - xStop) * (*xProjection - xStop) +
262  (*yProjection - yStop) * (*yProjection - yStop));
263  *distanceToLine = qSqrt ((xToProject - xStop) * (xToProject - xStop) +
264  (yToProject - yStop) * (yToProject - yStop));
265 
266  // Bring projection point to inside line
267  *xProjection = xStop;
268  *yProjection = yStop;
269 
270  } else {
271 
272  *distanceToLine = qSqrt ((xToProject - *xProjection) * (xToProject - *xProjection) +
273  (yToProject - *yProjection) * (yToProject - *yProjection));
274 
275  // Projected point is aleady inside line
276  *projectedDistanceOutsideLine = 0.0;
277 
278  }
279 }
280 
281 void setPixelRGB(QImage &image, int x, int y, QRgb q)
282 {
283  switch (image.depth())
284  {
285  case 1:
286  setPixelRGB1(image, x, y, q);
287  return;
288  case 8:
289  setPixelRGB8(image, x, y, q);
290  return;
291  case 32:
292  setPixelRGB32(image, x, y, q);
293  return;
294  }
295 }
296 
297 void setPixelRGB1(QImage &image1Bit, int x, int y, QRgb q)
298 {
299  for (int index = 0; index < image1Bit.colorCount(); index++) {
300  if (q == image1Bit.color(index))
301  {
302  if (image1Bit.format () == QImage::Format_MonoLSB)
303  {
304  *(image1Bit.scanLine (y) + (x >> 3)) &= ~(1 << (x & 7));
305  if (index > 0)
306  *(image1Bit.scanLine (y) + (x >> 3)) |= index << (x & 7);
307  }
308  else
309  {
310  *(image1Bit.scanLine (y) + (x >> 3)) &= ~(1 << (7 - (x & 7)));
311  if (index > 0)
312  *(image1Bit.scanLine (y) + (x >> 3)) |= index << (7 - (x & 7));
313  }
314  return;
315  }
316  }
317 }
318 
319 void setPixelRGB8(QImage &image8Bit, int x, int y, QRgb q)
320 {
321  for (int index = 0; index < image8Bit.colorCount(); index++) {
322  if (q == image8Bit.color(index))
323  {
324  *(image8Bit.scanLine(y) + x) = static_cast<uchar> (index);
325  return;
326  }
327  }
328 }
329 
330 void setPixelRGB32(QImage &image32Bit, int x, int y, QRgb q)
331 {
332  int* p = (int *)image32Bit.scanLine(y) + x;
333  *p = signed (q);
334 }
QRgb pixelRGB8(const QImage &image8Bit, int x, int y)
Get pixel method for 8 bit depth.
Definition: mmsubs.cpp:195
void setPixelRGB32(QImage &image32Bit, int x, int y, QRgb q)
Set pixel method for 32 bit depth.
Definition: mmsubs.cpp:330
void projectPointOntoLine(double xToProject, double yToProject, double xStart, double yStart, double xStop, double yStop, double *xProjection, double *yProjection, double *projectedDistanceOutsideLine, double *distanceToLine)
Find the projection of a point onto a line segment such that the line through the point and its proje...
Definition: mmsubs.cpp:211
void setPixelRGB(QImage &image, int x, int y, QRgb q)
Set pixel method for any bit depth.
Definition: mmsubs.cpp:281
QRgb pixelRGB(const QImage &image, int x, int y)
Get pixel method for any bit depth.
Definition: mmsubs.cpp:169
void ellipseFromParallelogram(double xTL, double yTL, double xTR, double yTR, double xBR, double yBR, double &angleRadians, double &aAligned, double &bAligned)
Calculate ellipse parameters that is incribed in a parallelogram centered at the origin, given three successive corners of that parallelogram.
Definition: mmsubs.cpp:52
QRgb pixelRGB32(const QImage &image32Bit, int x, int y)
Get pixel method for 32 bit depth.
Definition: mmsubs.cpp:201
const double PI
void setPixelRGB8(QImage &image8Bit, int x, int y, QRgb q)
Set pixel method for 8 bit depth.
Definition: mmsubs.cpp:319
double angleBetweenVectors(const QPointF &v1, const QPointF &v2)
Angle between two vectors. Direction is unimportant, so result is between 0 to pi radians...
Definition: mmsubs.cpp:15
QRgb pixelRGB1(const QImage &image1Bit, int x, int y)
Get pixel method for one bit depth.
Definition: mmsubs.cpp:182
double angleFromVectorToVector(const QPointF &vFrom, const QPointF &vTo)
Angle between two vectors. Direction is positive when rotation is about +z vector, so result is betwen -pi to pi radians.
Definition: mmsubs.cpp:32
void setPixelRGB1(QImage &image1Bit, int x, int y, QRgb q)
Set pixel method for one bit depth.
Definition: mmsubs.cpp:297