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