src/2geom/ellipse.cpp

   1 /*
   2  * Ellipse Curve
   3  *
   4  * Authors:
   5  *      Marco Cecchetti <mrcekets at gmail.com>
   6  *
   7  * Copyright 2008  authors
   8  *
   9  * This library is free software; you can redistribute it and/or
  10  * modify it either under the terms of the GNU Lesser General Public
  11  * License version 2.1 as published by the Free Software Foundation
  12  * (the "LGPL") or, at your option, under the terms of the Mozilla
  13  * Public License Version 1.1 (the "MPL"). If you do not alter this
  14  * notice, a recipient may use your version of this file under either
  15  * the MPL or the LGPL.
  16  *
  17  * You should have received a copy of the LGPL along with this library
  18  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  20  * You should have received a copy of the MPL along with this library
  21  * in the file COPYING-MPL-1.1
  22  *
  23  * The contents of this file are subject to the Mozilla Public License
  24  * Version 1.1 (the "License"); you may not use this file except in
  25  * compliance with the License. You may obtain a copy of the License at
  26  * http://www.mozilla.org/MPL/
  27  *
  28  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  29  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  30  * the specific language governing rights and limitations.
  31  */
  32
  33
  34 #include <2geom/ellipse.h>
  35 #include <2geom/svg-elliptical-arc.h>
  36 #include <2geom/numeric/fitting-tool.h>
  37 #include <2geom/numeric/fitting-model.h>
  38
  39
  40 namespace Geom
  41 {
  42
  43 void Ellipse::set(double A, double B, double C, double D, double E, double F)
  44 {
  45     double den = 4*A*C - B*B;
  46     if ( den == 0 )
  47     {
  48         THROW_LOGICALERROR("den == 0, while computing ellipse centre");
  49     }
  50     m_centre[X] = (B*E - 2*C*D) / den;
  51     m_centre[Y] = (B*D - 2*A*E) / den;
  52
  53     // evaluate the a coefficient of the ellipse equation in normal form
  54     // E(x,y) = a*(x-cx)^2 + b*(x-cx)*(y-cy) + c*(y-cy)^2 = 1
  55     // where b = a*B , c = a*C, (cx,cy) == centre
  56     double num =   A * sqr(m_centre[X])
  57                  + B * m_centre[X] * m_centre[Y]
  58                  + C * sqr(m_centre[Y])
  59                  - F;
  60
  61
  62     //evaluate ellipse rotation angle
  63     double rot = std::atan2( -B, -(A - C) )/2;
  64 //      std::cerr << "rot = " << rot << std::endl;
  65     bool swap_axes = false;
  66     if ( are_near(rot, 0) ) rot = 0;
  67     if ( are_near(rot, M_PI/2)  || rot < 0 )
  68     {
  69         swap_axes = true;
  70     }
  71
  72     // evaluate the length of the ellipse rays
  73     double cosrot = std::cos(rot);
  74     double sinrot = std::sin(rot);
  75     double cos2 = cosrot * cosrot;
  76     double sin2 = sinrot * sinrot;
  77     double cossin = cosrot * sinrot;
  78
  79     den = A * cos2 + B * cossin + C * sin2;
  80     if ( den == 0 )
  81     {
  82         THROW_LOGICALERROR("den == 0, while computing 'rx' coefficient");
  83     }
  84     double rx2 =  num/den;
  85     if ( rx2 < 0 )
  86     {
  87         THROW_LOGICALERROR("rx2 < 0, while computing 'rx' coefficient");
  88     }
  89     double rx = std::sqrt(rx2);
  90
  91     den = C * cos2 - B * cossin + A * sin2;
  92     if ( den == 0 )
  93     {
  94         THROW_LOGICALERROR("den == 0, while computing 'ry' coefficient");
  95     }
  96     double ry2 =  num/den;
  97     if ( ry2 < 0 )
  98     {
  99         THROW_LOGICALERROR("ry2 < 0, while computing 'rx' coefficient");
 100     }
 101     double ry = std::sqrt(ry2);
 102
 103     // the solution is not unique so we choose always the ellipse
 104     // with a rotation angle between 0 and PI/2
 105     if ( swap_axes ) std::swap(rx, ry);
 106     if (    are_near(rot,  M_PI/2)
 107          || are_near(rot, -M_PI/2)
 108          || are_near(rx, ry)       )
 109     {
 110         rot = 0;
 111     }
 112     else if ( rot < 0 )
 113     {
 114         rot += M_PI/2;
 115     }
 116
 117     m_ray[X] = rx;
 118     m_ray[Y] = ry;
 119     m_angle = rot;
 120 }
 121
 122
 123 std::vector<double> Ellipse::implicit_form_coefficients() const
 124 {
 125     if (ray(X) == 0 || ray(Y) == 0)
 126     {
 127         THROW_LOGICALERROR("a degenerate ellipse doesn't own an implicit form");
 128     }
 129
 130     std::vector<double> coeff(6);
 131     double cosrot = std::cos(rot_angle());
 132     double sinrot = std::sin(rot_angle());
 133     double cos2 = cosrot * cosrot;
 134     double sin2 = sinrot * sinrot;
 135     double cossin = cosrot * sinrot;
 136     double invrx2 = 1 / (ray(X) * ray(X));
 137     double invry2 = 1 / (ray(Y) * ray(Y));
 138
 139     coeff[0] = invrx2 * cos2 + invry2 * sin2;
 140     coeff[1] = 2 * (invrx2 - invry2) * cossin;
 141     coeff[2] = invrx2 * sin2 + invry2 * cos2;
 142     coeff[3] = -(2 * coeff[0] * center(X) + coeff[1] * center(Y));
 143     coeff[4] = -(2 * coeff[2] * center(Y) + coeff[1] * center(X));
 144     coeff[5] = coeff[0] * center(X) * center(X)
 145              + coeff[1] * center(X) * center(Y)
 146              + coeff[2] * center(Y) * center(Y)
 147              - 1;
 148     return coeff;
 149 }
 150
 151
 152 void Ellipse::set(std::vector<Point> const& points)
 153 {
 154     size_t sz = points.size();
 155     if (sz < 5)
 156     {
 157         THROW_RANGEERROR("fitting error: too few points passed");
 158     }
 159     NL::LFMEllipse model;
 160     NL::least_squeares_fitter<NL::LFMEllipse> fitter(model, sz);
 161
 162     for (size_t i = 0; i < sz; ++i)
 163     {
 164         fitter.append(points[i]);
 165     }
 166     fitter.update();
 167
 168     NL::Vector z(sz, 0.0);
 169     model.instance(*this, fitter.result(z));
 170 }
 171
 172
 173 SVGEllipticalArc
 174 Ellipse::arc(Point const& initial, Point const& inner, Point const& final,
 175              bool _svg_compliant)
 176 {
 177     Point sp_cp = initial - center();
 178     Point ep_cp = final   - center();
 179     Point ip_cp = inner   - center();
 180
 181     double angle1 = angle_between(sp_cp, ep_cp);
 182     double angle2 = angle_between(sp_cp, ip_cp);
 183     double angle3 = angle_between(ip_cp, ep_cp);
 184
 185     bool large_arc_flag = true;
 186     bool sweep_flag = true;
 187
 188     if ( angle1 > 0 )
 189     {
 190         if ( angle2 > 0 && angle3 > 0 )
 191         {
 192             large_arc_flag = false;
 193             sweep_flag = true;
 194         }
 195         else
 196         {
 197             large_arc_flag = true;
 198             sweep_flag = false;
 199         }
 200     }
 201     else
 202     {
 203         if ( angle2 < 0 && angle3 < 0 )
 204         {
 205             large_arc_flag = false;
 206             sweep_flag = false;
 207         }
 208         else
 209         {
 210             large_arc_flag = true;
 211             sweep_flag = true;
 212         }
 213     }
 214
 215     SVGEllipticalArc ea( initial, ray(X), ray(Y), rot_angle(),
 216                       large_arc_flag, sweep_flag, final, _svg_compliant);
 217     return ea;
 218 }
 219
 220 Ellipse Ellipse::transformed(Matrix const& m) const
 221 {
 222     double cosrot = std::cos(rot_angle());
 223     double sinrot = std::sin(rot_angle());
 224     Matrix A(  ray(X) * cosrot, ray(X) * sinrot,
 225               -ray(Y) * sinrot, ray(Y) * cosrot,
 226                0,               0                );
 227     Point new_center = center() * m;
 228     Matrix M = m.without_translation();
 229     Matrix AM = A * M;
 230     if ( are_near(AM.det(), 0) )
 231     {
 232         double angle;
 233         if (AM[0] != 0)
 234         {
 235             angle = std::atan2(AM[2], AM[0]);
 236         }
 237         else if (AM[1] != 0)
 238         {
 239             angle = std::atan2(AM[3], AM[1]);
 240         }
 241         else
 242         {
 243             angle = M_PI/2;
 244         }
 245         Point V(std::cos(angle), std::sin(angle));
 246         V *= AM;
 247         double rx = L2(V);
 248         angle = atan2(V);
 249         return Ellipse(new_center[X], new_center[Y], rx, 0, angle);
 250     }
 251
 252     std::vector<double> coeff = implicit_form_coefficients();
 253     Matrix Q( coeff[0],   coeff[1]/2,
 254               coeff[1]/2, coeff[2],
 255               0,          0   );
 256
 257     Matrix invm = M.inverse();
 258     Q = invm * Q ;
 259     std::swap( invm[1], invm[2] );
 260     Q *= invm;
 261     Ellipse e(Q[0], 2*Q[1], Q[3], 0, 0, -1);
 262     e.m_centre = new_center;
 263
 264     return e;
 265 }
 266
 267 Ellipse::Ellipse(Geom::Circle const &c)
 268 {
 269     m_centre = c.center();
 270     m_ray[X] = m_ray[Y] = c.ray();
 271 }
 272
 273 }  // end namespace Geom
 274
 275
 276 /*
 277   Local Variables:
 278   mode:c++
 279   c-file-style:"stroustrup"
 280   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 281   indent-tabs-mode:nil
 282   fill-column:99
 283   End:
 284 */
 285 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
 286
 287