Engauge Digitizer 2
Loading...
Searching...
No Matches
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
13const double PI = 3.1415926535;
14
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
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
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
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
138 m10, m11, m12,
139 m20, m21, m22);
141 -1.0, m11, m12,
142 0.0 , m21, m22);
144 m10, -1.0, m12,
145 m20, 0.0 , m22);
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));
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
169QRgb 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
182QRgb 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
195QRgb pixelRGB8(const QImage &image8Bit, int x, int y)
196{
197 int tableIndex = *(image8Bit.scanLine(y) + x);
198 return image8Bit.color(tableIndex);
199}
200
201QRgb 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
212 double yToProject,
213 double xStart,
214 double yStart,
215 double xStop,
216 double yStop,
217 double *xProjection,
218 double *yProjection,
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);
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);
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
254
255 // Bring projection point to inside line
258
259 } else if (s > 1) {
260
262 (*yProjection - yStop) * (*yProjection - yStop));
264 (yToProject - yStop) * (yToProject - yStop));
265
266 // Bring projection point to inside line
269
270 } else {
271
274
275 // Projected point is aleady inside line
277
278 }
279}
280
281void 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
297void 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
319void 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
330void setPixelRGB32(QImage &image32Bit, int x, int y, QRgb q)
331{
332 int* p = (int *)image32Bit.scanLine(y) + x;
333 *p = signed (q);
334}
const int INNER_RADIUS_MIN
QRgb pixelRGB8(const QImage &image8Bit, int x, int y)
Get pixel method for 8 bit depth.
Definition mmsubs.cpp:195
double angleFromVectorToVector(const QPointF &vFrom, const QPointF &vTo)
Angle between two vectors. Direction is positive when rotation is about +z vector,...
Definition mmsubs.cpp:32
QRgb pixelRGB32(const QImage &image32Bit, int x, int y)
Get pixel method for 32 bit depth.
Definition mmsubs.cpp:201
QRgb pixelRGB(const QImage &image, int x, int y)
Get pixel method for any bit depth.
Definition mmsubs.cpp:169
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
void setPixelRGB(QImage &image, int x, int y, QRgb q)
Set pixel method for any bit depth.
Definition mmsubs.cpp:281
void setPixelRGB1(QImage &image1Bit, int x, int y, QRgb q)
Set pixel method for one bit depth.
Definition mmsubs.cpp:297
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
const double PI
Definition mmsubs.cpp:13
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,...
Definition mmsubs.cpp:52
QRgb pixelRGB1(const QImage &image1Bit, int x, int y)
Get pixel method for one bit depth.
Definition mmsubs.cpp:182
void setPixelRGB32(QImage &image32Bit, int x, int y, QRgb q)
Set pixel method for 32 bit depth.
Definition mmsubs.cpp:330