File indexing completed on 2024-12-22 04:14:22

 /*
  * SPDX-FileCopyrightText: 2010 Geoffry Song <goffrie@gmail.com>
  *
  *  SPDX-License-Identifier: LGPL-2.0-or-later
  */
 
 #include "Ellipse.h"
 #include <cmath>
 
 #include <iostream>
 
 Ellipse::Ellipse() : a(-1), b(-1)
 {
 }
 Ellipse::Ellipse(const QPointF& _p1, const QPointF& _p2, const QPointF& _p3) : p1(_p1), p2(_p2), p3(_p3) {
     changeMajor();
 }
 
 Ellipse::~Ellipse()
 {
 }
 
 bool Ellipse::set(const QPointF& m1, const QPointF& m2, const QPointF& p)
 {
     bool changedMajor = m1 != p1 || m2 != p2,
          changedMinor = !changedMajor && p != p3;
     p1 = m1;
     p2 = m2;
     p3 = p;
     if (changedMajor) {
         return changeMajor();
     } else if (changedMinor) {
         return changeMinor();
     } else {
         return a > 0 && b > 0;
     }
 }
 
 QPointF Ellipse::project(const QPointF& pt) const
 {
     if (a <= 0 || b <= 0) return pt; // not a valid ellipse
     QPointF p = matrix.map(pt);
     /*
      * intersect line from (0,0) to p with the ellipse in canonical position
      * the equation of the line is y = py/px x
      * the equation of the ellipse is x^2/a^2 + y^2/b^2 = 1
      * x=(a*b*px)/sqrt(a^2*py^2+b^2*px^2)
      * y=(a*b*py)/sqrt(a^2*py^2+b^2*px^2)
      */
     const qreal divisor = sqrt(a * a * p.y() * p.y() + b * b * p.x() * p.x());
     if (divisor <= 0) return inverse.map(QPointF(a, 0)); // give up
     const qreal ab = a * b, factor = 1.0 / divisor;
     QPointF ep(ab * p.x() * factor, ab * p.y() * factor);
     return inverse.map(ep);
 /*    return inverse.map(closest(matrix.map(pt)));*/
 }
 
 inline QPointF rotate90(const QPointF& p) {
     return QPointF(p.y(), -p.x());
 }
 
 QRectF Ellipse::boundingRect() const
 {
     const QPointF d = rotate90((p2 - p1) * 0.5 * b / a);
     const QPointF pts[4] = {
         p1 + d,
         p1 - d,
         p2 + d,
         p2 - d
     };
     QRectF ret;
     for (int i = 0; i < 4; ++i) {
         ret = ret.united(QRectF(pts[i], QSizeF(0.0001, 0.0001)));
     }
     return ret;
 }
 
 inline qreal sqrlength(const QPointF& vec)
 {
     return vec.x() * vec.x() + vec.y() * vec.y();
 }
 inline qreal length(const QPointF& vec)
 {
     return sqrt(vec.x() * vec.x() + vec.y() * vec.y());
 }
 
 bool Ellipse::changeMajor()
 {
     a = length(p1 - p2) * 0.5;
     
     /*
      * calculate transformation matrix
      * x' = m11*x + m21*y + dx
      * y' = m22*y + m12*x + dy
      * m11 = m22, m12 = -m21 (rotations and translations only)
      * x' = m11*x - m12*y + dx
      * y' = m11*y + m12*x + dy
      * 
      * then, transforming (x1, y1) to (x1', y1') and (x2, y2) to (x2', y2'):
      * 
      * m11 = (y2*y2' + y1 * (y1'-y2') - y2*y1' + x2*x'2 - x1*x'2 + (x1-x2)*x1')
      *       ------------------------------------------------------------------
      *                 (y2^2 - 2*y1*y2 + y1^2 + x2^2 - 2*x1*x2 + x1^2)
      * m12 = -(x1*(y2'-y1') - x2*y2' + x2*y1' + x2'*y2 - x1'*y2 + (x1'-x2')*y1)
      *       ------------------------------------------------------------------
      *                 (y2^2 - 2*y1*y2 + y1^2 + x2^2 - 2*x1*x2 + x1^2)
      * dx = (x1*(-y2*y2' + y2*y1' - x2*x2') + y1*( x2*y2' - x2*y1' - x2'*y2 - x1'*y2) + x2'*y1^2 + x1^2*x2' + x1'*(y2^2 + x2^2 - x1*x2))
      *      ----------------------------------------------------------------------------------------------------------------------------
      *                 (y2^2 - 2*y1*y2 + y1^2 + x2^2 - 2*x1*x2 + x1^2)
      * dy = (x1*(-x2*y2' - x2*y1' + x2'*y2) + y1*(-y2*y2' - y2*y1' - x2*x2' + x2*x1') + y2'*y1^2 + x1^2*y2' + y1'(y2^2 + x2^2) - x1*x1'*y2)
      *      -------------------------------------------------------------------------------------------------------------------------------
      *                 (y2^2 - 2*y1*y2 + y1^2 + x2^2 - 2*x1*x2 + x1^2)
      * 
      * in our usage, to move the ellipse into canonical position:
      * 
      * (x1, y1) = p1
      * (x2, y2) = p2
      * (x1', y1') = (-a, 0)
      * (x2', y2') = (a, 0)
      */
     
     const qreal
         x1 = p1.x(),
         x2 = p2.x(),
         y1 = p1.y(),
         y2 = p2.y(),
         x1p = -a,
         x2p = a,
         x1sqr = x1 * x1,
         x2sqr = x2 * x2,
         y1sqr = y1 * y1,
         y2sqr = y2 * y2,
         factor = 1.0 / (x1sqr + y1sqr + x2sqr + y2sqr - 2.0 * y1 * y2 - 2.0 * x1 * x2),
         m11 = (x2*x2p - x1*x2p + (x1-x2)*x1p) * factor,
         m12 = -(x2p*y2 - x1p*y2 + (x1p-x2p)*y1) * factor,
         dx = (x1*(-x2*x2p) + y1*(-x2p*y2 - x1p*y2) + x2p*y1sqr + x1sqr*x2p + x1p*(y2sqr + x2sqr - x1*x2)) * factor,
         dy = (x1*(x2p*y2) + y1*(-x2*x2p + x2*x1p) - x1*x1p*y2) * factor;
     
     matrix = QTransform(m11, m12, -m12, m11, dx, dy);
     inverse = matrix.inverted();
 
     return changeMinor();
 }
 
 bool Ellipse::changeMinor()
 {
     QPointF p = matrix.map(p3);
 
     /*
      * ellipse formula:
      * x^2/a^2 + y^2/b^2 = 1
      * b = sqrt(y^2 / (1 - x^2/a^2))
      */
     const qreal
         asqr = a * a,
         xsqr = p.x() * p.x(),
         ysqr = p.y() * p.y(),
         divisor = (1.0 - xsqr / asqr);
     if (divisor <= 0) {
         // division by zero!
         b = -1;
         return false;
     }
     b = sqrt(ysqr / divisor);
     return true;
 }
 
 bool Ellipse::setMajor1(const QPointF& p)
 {
     p1 = p;
     return changeMajor();
 }
 bool Ellipse::setMajor2(const QPointF& p) {
     p2 = p;
     return changeMajor();
 }
 bool Ellipse::setPoint(const QPointF& p)
 {
     p3 = p;
     return changeMinor();
 }