← Back to branch summary

~valavanisalex/ubuntu/oneiric/inkscape/inkscape_0.48.1-2ubuntu4

« back to all changes in this revision

Viewing changes to src/2geom/elliptical-arc.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Kees Cook, Ted Gould, Kees Cook
  • Date: 2009-06-24 14:00:43 UTC
  • mfrom: (1.1.8 upstream)
  • Revision ID: james.westby@ubuntu.com-20090624140043-07stp20mry48hqup
Tags: 0.47~pre0-0ubuntu1
* New upstream release

[ Ted Gould ]
* debian/control: Adding libgsl0 and removing version specifics on boost

[ Kees Cook ]
* debian/watch: updated to run uupdate and mangle pre-release versions.
* Dropped patches that have been taken upstream:
  - 01_mips
  - 02-poppler-0.8.3
  - 03-chinese-inkscape
  - 05_fix_latex_patch
  - 06_gcc-4.4
  - 07_cdr2svg
  - 08_skip-bad-utf-on-pdf-import
  - 09_gtk-clist
  - 10_belarussian
  - 11_libpng
  - 12_desktop
  - 13_slider
  - 100_svg_import_improvements
  - 102_sp_pattern_painter_free
  - 103_bitmap_type_print

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/2geom/elliptical-arc.cpp
 
1
/*
 
2
 * SVG Elliptical Arc Class
 
3
 *
 
4
 * Copyright 2008  Marco Cecchetti <mrcekets at gmail.com>
 
5
 *
 
6
 * This library is free software; you can redistribute it and/or
 
7
 * modify it either under the terms of the GNU Lesser General Public
 
8
 * License version 2.1 as published by the Free Software Foundation
 
9
 * (the "LGPL") or, at your option, under the terms of the Mozilla
 
10
 * Public License Version 1.1 (the "MPL"). If you do not alter this
 
11
 * notice, a recipient may use your version of this file under either
 
12
 * the MPL or the LGPL.
 
13
 *
 
14
 * You should have received a copy of the LGPL along with this library
 
15
 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
 
16
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 
17
 * You should have received a copy of the MPL along with this library
 
18
 * in the file COPYING-MPL-1.1
 
19
 *
 
20
 * The contents of this file are subject to the Mozilla Public License
 
21
 * Version 1.1 (the "License"); you may not use this file except in
 
22
 * compliance with the License. You may obtain a copy of the License at
 
23
 * http://www.mozilla.org/MPL/
 
24
 *
 
25
 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
 
26
 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
 
27
 * the specific language governing rights and limitations.
 
28
 */
 
29
 
 
30
 
 
31
#include <2geom/elliptical-arc.h>
 
32
#include <2geom/bezier-curve.h>
 
33
#include <2geom/poly.h>
 
34
#include <2geom/sbasis-math.h>
 
35
 
 
36
#include <cfloat>
 
37
#include <limits>
 
38
 
 
39
 
 
40
 
 
41
 
 
42
namespace Geom
 
43
{
 
44
 
 
45
 
 
46
OptRect EllipticalArc::boundsExact() const
 
47
{
 
48
        std::vector<double> extremes(4);
 
49
        double cosrot = std::cos(rotation_angle());
 
50
        double sinrot = std::sin(rotation_angle());
 
51
        extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
 
52
        extremes[1] = extremes[0] + M_PI;
 
53
        if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;
 
54
        extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
 
55
        extremes[3] = extremes[2] + M_PI;
 
56
        if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;
 
57
 
 
58
 
 
59
        std::vector<double>arc_extremes(4);
 
60
        arc_extremes[0] = initialPoint()[X];
 
61
        arc_extremes[1] = finalPoint()[X];
 
62
        if ( arc_extremes[0] < arc_extremes[1] )
 
63
                std::swap(arc_extremes[0], arc_extremes[1]);
 
64
        arc_extremes[2] = initialPoint()[Y];
 
65
        arc_extremes[3] = finalPoint()[Y];
 
66
        if ( arc_extremes[2] < arc_extremes[3] )
 
67
                std::swap(arc_extremes[2], arc_extremes[3]);
 
68
 
 
69
 
 
70
        if ( start_angle() < end_angle() )
 
71
        {
 
72
                if ( sweep_flag() )
 
73
                {
 
74
                        for ( unsigned int i = 0; i < extremes.size(); ++i )
 
75
                        {
 
76
                                if ( start_angle() < extremes[i] && extremes[i] < end_angle() )
 
77
                                {
 
78
                                        arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
 
79
                                }
 
80
                        }
 
81
                }
 
82
                else
 
83
                {
 
84
                        for ( unsigned int i = 0; i < extremes.size(); ++i )
 
85
                        {
 
86
                                if ( start_angle() > extremes[i] || extremes[i] > end_angle() )
 
87
                                {
 
88
                                        arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
 
89
                                }
 
90
                        }
 
91
                }
 
92
        }
 
93
        else
 
94
        {
 
95
                if ( sweep_flag() )
 
96
                {
 
97
                        for ( unsigned int i = 0; i < extremes.size(); ++i )
 
98
                        {
 
99
                                if ( start_angle() < extremes[i] || extremes[i] < end_angle() )
 
100
                                {
 
101
                                        arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
 
102
                                }
 
103
                        }
 
104
                }
 
105
                else
 
106
                {
 
107
                        for ( unsigned int i = 0; i < extremes.size(); ++i )
 
108
                        {
 
109
                                if ( start_angle() > extremes[i] && extremes[i] > end_angle() )
 
110
                                {
 
111
                                        arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
 
112
                                }
 
113
                        }
 
114
                }
 
115
        }
 
116
 
 
117
        return Rect( Point(arc_extremes[1], arc_extremes[3]) ,
 
118
                             Point(arc_extremes[0], arc_extremes[2]) );
 
119
 
 
120
}
 
121
 
 
122
 
 
123
std::vector<double>
 
124
EllipticalArc::roots(double v, Dim2 d) const
 
125
{
 
126
        if ( d > Y )
 
127
        {
 
128
                THROW_RANGEERROR("dimention out of range");
 
129
        }
 
130
 
 
131
        std::vector<double> sol;
 
132
        if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
 
133
        {
 
134
                if ( center(d) == v )
 
135
                        sol.push_back(0);
 
136
                return sol;
 
137
        }
 
138
 
 
139
        const char* msg[2][2] =
 
140
        {
 
141
                { "d == X; ray(X) == 0; "
 
142
                  "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "
 
143
                  "s should be contained in [-1,1]",
 
144
                  "d == X; ray(Y) == 0; "
 
145
                  "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "
 
146
                  "s should be contained in [-1,1]"
 
147
                },
 
148
                { "d == Y; ray(X) == 0; "
 
149
                  "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "
 
150
                  "s should be contained in [-1,1]",
 
151
                  "d == Y; ray(Y) == 0; "
 
152
                  "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "
 
153
                  "s should be contained in [-1,1]"
 
154
                },
 
155
        };
 
156
 
 
157
        for ( unsigned int dim = 0; dim < 2; ++dim )
 
158
        {
 
159
                if ( are_near(ray(dim), 0) )
 
160
                {
 
161
 
 
162
                        if ( initialPoint()[d] == v && finalPoint()[d] == v )
 
163
                        {
 
164
                                THROW_INFINITESOLUTIONS(0);
 
165
                        }
 
166
                        if ( (initialPoint()[d] < finalPoint()[d])
 
167
                                 && (initialPoint()[d] > v || finalPoint()[d] < v) )
 
168
                        {
 
169
                                return sol;
 
170
                        }
 
171
                        if ( (initialPoint()[d] > finalPoint()[d])
 
172
                                 && (finalPoint()[d] > v || initialPoint()[d] < v) )
 
173
                        {
 
174
                                return sol;
 
175
                        }
 
176
                        double ray_prj;
 
177
                        switch(d)
 
178
                        {
 
179
                                case X:
 
180
                                        switch(dim)
 
181
                                        {
 
182
                                                case X: ray_prj = -ray(Y) * std::sin(rotation_angle());
 
183
                                                                break;
 
184
                                                case Y: ray_prj = ray(X) * std::cos(rotation_angle());
 
185
                                                                break;
 
186
                                        }
 
187
                                        break;
 
188
                                case Y:
 
189
                                        switch(dim)
 
190
                                        {
 
191
                                                case X: ray_prj = ray(Y) * std::cos(rotation_angle());
 
192
                                                                break;
 
193
                                                case Y: ray_prj = ray(X) * std::sin(rotation_angle());
 
194
                                                                break;
 
195
                                        }
 
196
                                        break;
 
197
                        }
 
198
 
 
199
                        double s = (v - center(d)) / ray_prj;
 
200
                        if ( s < -1 || s > 1 )
 
201
                        {
 
202
                                THROW_LOGICALERROR(msg[d][dim]);
 
203
                        }
 
204
                        switch(dim)
 
205
                        {
 
206
                                case X:
 
207
                                        s = std::asin(s); // return a value in [-PI/2,PI/2]
 
208
                                        if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) )  )
 
209
                                        {
 
210
                                                if ( s < 0 ) s += 2*M_PI;
 
211
                                        }
 
212
                                        else
 
213
                                        {
 
214
                                                s = M_PI - s;
 
215
                                                if (!(s < 2*M_PI) ) s -= 2*M_PI;
 
216
                                        }
 
217
                                        break;
 
218
                                case Y:
 
219
                                        s = std::acos(s); // return a value in [0,PI]
 
220
                                        if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )
 
221
                                        {
 
222
                                                s = 2*M_PI - s;
 
223
                                                if ( !(s < 2*M_PI) ) s -= 2*M_PI;
 
224
                                        }
 
225
                                        break;
 
226
                        }
 
227
 
 
228
                        //std::cerr << "s = " << rad_to_deg(s);
 
229
                        s = map_to_01(s);
 
230
                        //std::cerr << " -> t: " << s << std::endl;
 
231
                        if ( !(s < 0 || s > 1) )
 
232
                                sol.push_back(s);
 
233
                        return sol;
 
234
                }
 
235
        }
 
236
 
 
237
        double rotx, roty;
 
238
        switch(d)
 
239
        {
 
240
                case X:
 
241
                        rotx = std::cos(rotation_angle());
 
242
                        roty = -std::sin(rotation_angle());
 
243
                        break;
 
244
                case Y:
 
245
                        rotx = std::sin(rotation_angle());
 
246
                        roty = std::cos(rotation_angle());
 
247
                        break;
 
248
        }
 
249
        double rxrotx = ray(X) * rotx;
 
250
        double c_v = center(d) - v;
 
251
 
 
252
        double a = -rxrotx + c_v;
 
253
        double b = ray(Y) * roty;
 
254
        double c = rxrotx + c_v;
 
255
        //std::cerr << "a = " << a << std::endl;
 
256
        //std::cerr << "b = " << b << std::endl;
 
257
        //std::cerr << "c = " << c << std::endl;
 
258
 
 
259
        if ( are_near(a,0) )
 
260
        {
 
261
                sol.push_back(M_PI);
 
262
                if ( !are_near(b,0) )
 
263
                {
 
264
                        double s = 2 * std::atan(-c/(2*b));
 
265
                        if ( s < 0 ) s += 2*M_PI;
 
266
                        sol.push_back(s);
 
267
                }
 
268
        }
 
269
        else
 
270
        {
 
271
                double delta = b * b - a * c;
 
272
                //std::cerr << "delta = " << delta << std::endl;
 
273
                if ( are_near(delta, 0) )
 
274
                {
 
275
                        double s = 2 * std::atan(-b/a);
 
276
                        if ( s < 0 ) s += 2*M_PI;
 
277
                        sol.push_back(s);
 
278
                }
 
279
                else if ( delta > 0 )
 
280
                {
 
281
                        double sq = std::sqrt(delta);
 
282
                        double s = 2 * std::atan( (-b - sq) / a );
 
283
                        if ( s < 0 ) s += 2*M_PI;
 
284
                        sol.push_back(s);
 
285
                        s = 2 * std::atan( (-b + sq) / a );
 
286
                        if ( s < 0 ) s += 2*M_PI;
 
287
                        sol.push_back(s);
 
288
                }
 
289
        }
 
290
 
 
291
        std::vector<double> arc_sol;
 
292
        for (unsigned int i = 0; i < sol.size(); ++i )
 
293
        {
 
294
                //std::cerr << "s = " << rad_to_deg(sol[i]);
 
295
                sol[i] = map_to_01(sol[i]);
 
296
                //std::cerr << " -> t: " << sol[i] << std::endl;
 
297
                if ( !(sol[i] < 0 || sol[i] > 1) )
 
298
                        arc_sol.push_back(sol[i]);
 
299
        }
 
300
        return arc_sol;
 
301
 
 
302
 
 
303
//      return SBasisCurve(toSBasis()).roots(v, d);
 
304
}
 
305
 
 
306
// D(E(t,C),t) = E(t+PI/2,O)
 
307
Curve* EllipticalArc::derivative() const
 
308
{
 
309
        EllipticalArc* result = new EllipticalArc(*this);
 
310
        result->m_center[X] = result->m_center[Y] = 0;
 
311
        result->m_start_angle += M_PI/2;
 
312
        if( !( result->m_start_angle < 2*M_PI ) )
 
313
        {
 
314
                result->m_start_angle -= 2*M_PI;
 
315
        }
 
316
        result->m_end_angle += M_PI/2;
 
317
        if( !( result->m_end_angle < 2*M_PI ) )
 
318
        {
 
319
                result->m_end_angle -= 2*M_PI;
 
320
        }
 
321
        result->m_initial_point = result->pointAtAngle( result->start_angle() );
 
322
        result->m_final_point = result->pointAtAngle( result->end_angle() );
 
323
        return result;
 
324
 
 
325
}
 
326
 
 
327
std::vector<Point>
 
328
EllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
 
329
{
 
330
    unsigned int nn = n+1; // nn represents the size of the result vector.
 
331
        std::vector<Point> result;
 
332
        result.reserve(nn);
 
333
        double angle = map_unit_interval_on_circular_arc(t, start_angle(),
 
334
                                                                 end_angle(), sweep_flag());
 
335
        EllipticalArc ea(*this);
 
336
        ea.m_center = Point(0,0);
 
337
        unsigned int m = std::min(nn, 4u);
 
338
        for ( unsigned int i = 0; i < m; ++i )
 
339
        {
 
340
                result.push_back( ea.pointAtAngle(angle) );
 
341
                angle += M_PI/2;
 
342
                if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;
 
343
        }
 
344
        m = nn / 4;
 
345
        for ( unsigned int i = 1; i < m; ++i )
 
346
        {
 
347
                for ( unsigned int j = 0; j < 4; ++j )
 
348
                        result.push_back( result[j] );
 
349
        }
 
350
        m = nn - 4 * m;
 
351
        for ( unsigned int i = 0; i < m; ++i )
 
352
        {
 
353
                result.push_back( result[i] );
 
354
        }
 
355
        if ( !result.empty() ) // nn != 0
 
356
                result[0] = pointAtAngle(angle);
 
357
        return result;
 
358
}
 
359
 
 
360
D2<SBasis> EllipticalArc::toSBasis() const
 
361
{
 
362
    // the interval of parametrization has to be [0,1]
 
363
    Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );
 
364
    Linear param(start_angle(), et);
 
365
    Coord cos_rot_angle = std::cos(rotation_angle());
 
366
    Coord sin_rot_angle = std::sin(rotation_angle());
 
367
    // order = 4 seems to be enough to get a perfect looking elliptical arc
 
368
    // should it be choosen in function of the arc length anyway ?
 
369
    // or maybe a user settable parameter: toSBasis(unsigned int order) ?
 
370
    SBasis arc_x = ray(X) * cos(param,4);
 
371
    SBasis arc_y = ray(Y) * sin(param,4);
 
372
    D2<SBasis> arc;
 
373
    arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));
 
374
    arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));
 
375
    return arc;
 
376
}
 
377
 
 
378
 
 
379
bool EllipticalArc::containsAngle(Coord angle) const
 
380
{
 
381
        if ( sweep_flag() )
 
382
                if ( start_angle() < end_angle() )
 
383
                        return ( !( angle < start_angle() || angle > end_angle() ) );
 
384
                else
 
385
                        return ( !( angle < start_angle() && angle > end_angle() ) );
 
386
        else
 
387
                if ( start_angle() > end_angle() )
 
388
                        return ( !( angle > start_angle() || angle < end_angle() ) );
 
389
                else
 
390
                        return ( !( angle > start_angle() && angle < end_angle() ) );
 
391
}
 
392
 
 
393
 
 
394
double EllipticalArc::valueAtAngle(Coord t, Dim2 d) const
 
395
{
 
396
    double sin_rot_angle = std::sin(rotation_angle());
 
397
    double cos_rot_angle = std::cos(rotation_angle());
 
398
    if ( d == X )
 
399
    {
 
400
        return    ray(X) * cos_rot_angle * std::cos(t)
 
401
                - ray(Y) * sin_rot_angle * std::sin(t)
 
402
                + center(X);
 
403
    }
 
404
    else if ( d == Y )
 
405
    {
 
406
        return    ray(X) * sin_rot_angle * std::cos(t)
 
407
                + ray(Y) * cos_rot_angle * std::sin(t)
 
408
                + center(Y);
 
409
    }
 
410
    THROW_RANGEERROR("dimension parameter out of range");
 
411
}
 
412
 
 
413
 
 
414
Curve* EllipticalArc::portion(double f, double t) const
 
415
{
 
416
        if (f < 0) f = 0;
 
417
        if (f > 1) f = 1;
 
418
        if (t < 0) t = 0;
 
419
        if (t > 1) t = 1;
 
420
        if ( are_near(f, t) )
 
421
        {
 
422
                EllipticalArc* arc = new EllipticalArc();
 
423
                arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);
 
424
                arc->m_start_angle = arc->m_end_angle = m_start_angle;
 
425
                arc->m_rot_angle = m_rot_angle;
 
426
                arc->m_sweep = m_sweep;
 
427
                arc->m_large_arc = m_large_arc;
 
428
        }
 
429
    EllipticalArc* arc = new EllipticalArc( *this );
 
430
    arc->m_initial_point = pointAt(f);
 
431
    arc->m_final_point = pointAt(t);
 
432
    double sa = sweep_flag() ? sweep_angle() : -sweep_angle();
 
433
    arc->m_start_angle = m_start_angle + sa * f;
 
434
    if ( !(arc->m_start_angle < 2*M_PI) )
 
435
        arc->m_start_angle -= 2*M_PI;
 
436
    if ( arc->m_start_angle < 0 )
 
437
        arc->m_start_angle += 2*M_PI;
 
438
    arc->m_end_angle = m_start_angle + sa * t;
 
439
    if ( !(arc->m_end_angle < 2*M_PI) )
 
440
        arc->m_end_angle -= 2*M_PI;
 
441
    if ( arc->m_end_angle < 0 )
 
442
        arc->m_end_angle += 2*M_PI;
 
443
    if ( f > t ) arc->m_sweep = !sweep_flag();
 
444
    if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )
 
445
        arc->m_large_arc = false;
 
446
    return arc;
 
447
}
 
448
 
 
449
// NOTE: doesn't work with 360 deg arcs
 
450
void EllipticalArc::calculate_center_and_extreme_angles()
 
451
{
 
452
    // as stated in the svg standard the rotation angle parameter
 
453
    // value must be modulo 2*PI
 
454
    m_rot_angle = std::fmod(m_rot_angle, 2*M_PI);
 
455
    if (m_rot_angle < 0) m_rot_angle += 2*M_PI;
 
456
 
 
457
    if ( are_near(initialPoint(), finalPoint()) )
 
458
    {
 
459
        if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
 
460
        {
 
461
                m_start_angle = m_end_angle = 0;
 
462
                m_center = initialPoint();
 
463
                return;
 
464
        }
 
465
        else
 
466
        {
 
467
                THROW_RANGEERROR("initial and final point are the same");
 
468
        }
 
469
    }
 
470
        if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
 
471
        { // but initialPoint != finalPoint
 
472
                THROW_RANGEERROR(
 
473
                        "there is no ellipse that satisfies the given constraints: "
 
474
                        "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"
 
475
                );
 
476
        }
 
477
        if ( are_near(ray(Y), 0) )
 
478
        {
 
479
                Point v = initialPoint() - finalPoint();
 
480
                if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )
 
481
                {
 
482
                        double angle = std::atan2(v[Y], v[X]);
 
483
                        if (angle < 0) angle += 2*M_PI;
 
484
                        if ( are_near( angle, rotation_angle() ) )
 
485
                        {
 
486
                                m_start_angle = 0;
 
487
                                m_end_angle = M_PI;
 
488
                                m_center = v/2 + finalPoint();
 
489
                                return;
 
490
                        }
 
491
                        angle -= M_PI;
 
492
                        if ( angle < 0 ) angle += 2*M_PI;
 
493
                        if ( are_near( angle, rotation_angle() ) )
 
494
                        {
 
495
                                m_start_angle = M_PI;
 
496
                                m_end_angle = 0;
 
497
                                m_center = v/2 + finalPoint();
 
498
                                return;
 
499
                        }
 
500
                        THROW_RANGEERROR(
 
501
                                "there is no ellipse that satisfies the given constraints: "
 
502
                                "ray(Y) == 0 "
 
503
                                "and slope(initialPoint - finalPoint) != rotation_angle "
 
504
                                "and != rotation_angle + PI"
 
505
                        );
 
506
                }
 
507
                if ( L2sq(v) > 4*ray(X)*ray(X) )
 
508
                {
 
509
                        THROW_RANGEERROR(
 
510
                                "there is no ellipse that satisfies the given constraints: "
 
511
                                "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"
 
512
                        );
 
513
                }
 
514
                else
 
515
                {
 
516
                        THROW_RANGEERROR(
 
517
                                "there is infinite ellipses that satisfy the given constraints: "
 
518
                                "ray(Y) == 0  and distance(initialPoint, finalPoint) < 2*ray(X)"
 
519
                        );
 
520
                }
 
521
 
 
522
        }
 
523
 
 
524
        if ( are_near(ray(X), 0) )
 
525
        {
 
526
                Point v = initialPoint() - finalPoint();
 
527
                if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )
 
528
                {
 
529
                        double angle = std::atan2(v[Y], v[X]);
 
530
                        if (angle < 0) angle += 2*M_PI;
 
531
                        double rot_angle = rotation_angle() + M_PI/2;
 
532
                        if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;
 
533
                        if ( are_near( angle, rot_angle ) )
 
534
                        {
 
535
                                m_start_angle = M_PI/2;
 
536
                                m_end_angle = 3*M_PI/2;
 
537
                                m_center = v/2 + finalPoint();
 
538
                                return;
 
539
                        }
 
540
                        angle -= M_PI;
 
541
                        if ( angle < 0 ) angle += 2*M_PI;
 
542
                        if ( are_near( angle, rot_angle ) )
 
543
                        {
 
544
                                m_start_angle = 3*M_PI/2;
 
545
                                m_end_angle = M_PI/2;
 
546
                                m_center = v/2 + finalPoint();
 
547
                                return;
 
548
                        }
 
549
                        THROW_RANGEERROR(
 
550
                                "there is no ellipse that satisfies the given constraints: "
 
551
                                "ray(X) == 0 "
 
552
                                "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "
 
553
                                "and != rotation_angle + (3/2)*PI"
 
554
                        );
 
555
                }
 
556
                if ( L2sq(v) > 4*ray(Y)*ray(Y) )
 
557
                {
 
558
                        THROW_RANGEERROR(
 
559
                                "there is no ellipse that satisfies the given constraints: "
 
560
                                "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"
 
561
                        );
 
562
                }
 
563
                else
 
564
                {
 
565
                        THROW_RANGEERROR(
 
566
                                "there is infinite ellipses that satisfy the given constraints: "
 
567
                                "ray(X) == 0  and distance(initialPoint, finalPoint) < 2*ray(Y)"
 
568
                        );
 
569
                }
 
570
 
 
571
        }
 
572
 
 
573
    double sin_rot_angle = std::sin(rotation_angle());
 
574
    double cos_rot_angle = std::cos(rotation_angle());
 
575
 
 
576
    Point sp = sweep_flag() ? initialPoint() : finalPoint();
 
577
    Point ep = sweep_flag() ? finalPoint() : initialPoint();
 
578
 
 
579
    Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,
 
580
             -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,
 
581
              0,                      0 );
 
582
    Matrix im = m.inverse();
 
583
    Point sol = (ep - sp) * im;
 
584
    double half_sum_angle = std::atan2(-sol[X], sol[Y]);
 
585
    double half_diff_angle;
 
586
    if ( are_near(std::fabs(half_sum_angle), M_PI/2) )
 
587
    {
 
588
        double anti_sgn_hsa = (half_sum_angle > 0) ? -1 : 1;
 
589
        double arg = anti_sgn_hsa * sol[X] / 2;
 
590
        // if |arg| is a little bit > 1 acos returns nan
 
591
        if ( are_near(arg, 1) )
 
592
            half_diff_angle = 0;
 
593
        else if ( are_near(arg, -1) )
 
594
            half_diff_angle = M_PI;
 
595
        else
 
596
        {
 
597
                if ( !(-1 < arg && arg < 1) )
 
598
                        THROW_RANGEERROR(
 
599
                                "there is no ellipse that satisfies the given constraints"
 
600
                        );
 
601
            // assert( -1 < arg && arg < 1 );
 
602
            // if it fails
 
603
            // => there is no ellipse that satisfies the given constraints
 
604
            half_diff_angle = std::acos( arg );
 
605
        }
 
606
 
 
607
        half_diff_angle = M_PI/2 - half_diff_angle;
 
608
    }
 
609
    else
 
610
    {
 
611
        double  arg = sol[Y] / ( 2 * std::cos(half_sum_angle) );
 
612
        // if |arg| is a little bit > 1 asin returns nan
 
613
        if ( are_near(arg, 1) )
 
614
            half_diff_angle = M_PI/2;
 
615
        else if ( are_near(arg, -1) )
 
616
            half_diff_angle = -M_PI/2;
 
617
        else
 
618
        {
 
619
                if ( !(-1 < arg && arg < 1) )
 
620
                        THROW_RANGEERROR(
 
621
                                "there is no ellipse that satisfies the given constraints"
 
622
                        );
 
623
            // assert( -1 < arg && arg < 1 );
 
624
            // if it fails
 
625
            // => there is no ellipse that satisfies the given constraints
 
626
            half_diff_angle = std::asin( arg );
 
627
        }
 
628
    }
 
629
 
 
630
    if (   ( m_large_arc && half_diff_angle > 0 )
 
631
        || (!m_large_arc && half_diff_angle < 0 ) )
 
632
    {
 
633
        half_diff_angle = -half_diff_angle;
 
634
    }
 
635
    if ( half_sum_angle < 0 ) half_sum_angle += 2*M_PI;
 
636
    if ( half_diff_angle < 0 ) half_diff_angle += M_PI;
 
637
 
 
638
    m_start_angle = half_sum_angle - half_diff_angle;
 
639
    m_end_angle =  half_sum_angle + half_diff_angle;
 
640
    // 0 <= m_start_angle, m_end_angle < 2PI
 
641
    if ( m_start_angle < 0 ) m_start_angle += 2*M_PI;
 
642
    if( !(m_end_angle < 2*M_PI) ) m_end_angle -= 2*M_PI;
 
643
    sol[0] = std::cos(m_start_angle);
 
644
    sol[1] = std::sin(m_start_angle);
 
645
    m_center = sp - sol * m;
 
646
    if ( !sweep_flag() )
 
647
    {
 
648
        double angle = m_start_angle;
 
649
        m_start_angle = m_end_angle;
 
650
        m_end_angle = angle;
 
651
    }
 
652
}
 
653
 
 
654
Coord EllipticalArc::map_to_02PI(Coord t) const
 
655
{
 
656
    if ( sweep_flag() )
 
657
    {
 
658
        Coord angle = start_angle() + sweep_angle() * t;
 
659
        if ( !(angle < 2*M_PI) )
 
660
            angle -= 2*M_PI;
 
661
        return angle;
 
662
    }
 
663
    else
 
664
    {
 
665
        Coord angle = start_angle() - sweep_angle() * t;
 
666
        if ( angle < 0 ) angle += 2*M_PI;
 
667
        return angle;
 
668
    }
 
669
}
 
670
 
 
671
Coord EllipticalArc::map_to_01(Coord angle) const
 
672
{
 
673
        return map_circular_arc_on_unit_interval(angle, start_angle(),
 
674
                                                         end_angle(), sweep_flag());
 
675
}
 
676
 
 
677
 
 
678
std::vector<double> EllipticalArc::
 
679
allNearestPoints( Point const& p, double from, double to ) const
 
680
{
 
681
        if ( from > to ) std::swap(from, to);
 
682
        if ( from < 0 || to > 1 )
 
683
        {
 
684
                THROW_RANGEERROR("[from,to] interval out of range");
 
685
        }
 
686
        std::vector<double> result;
 
687
        if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) )  || are_near(from, to) )
 
688
        {
 
689
                result.push_back(from);
 
690
                return result;
 
691
        }
 
692
        else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
 
693
        {
 
694
                LineSegment seg(pointAt(from), pointAt(to));
 
695
                Point np = seg.pointAt( seg.nearestPoint(p) );
 
696
                if ( are_near(ray(Y), 0) )
 
697
                {
 
698
                        if ( are_near(rotation_angle(), M_PI/2)
 
699
                                 || are_near(rotation_angle(), 3*M_PI/2) )
 
700
                        {
 
701
                                result = roots(np[Y], Y);
 
702
                        }
 
703
                        else
 
704
                        {
 
705
                                result = roots(np[X], X);
 
706
                        }
 
707
                }
 
708
                else
 
709
                {
 
710
                        if ( are_near(rotation_angle(), M_PI/2)
 
711
                                 || are_near(rotation_angle(), 3*M_PI/2) )
 
712
                        {
 
713
                                result = roots(np[X], X);
 
714
                        }
 
715
                        else
 
716
                        {
 
717
                                result = roots(np[Y], Y);
 
718
                        }
 
719
                }
 
720
                return result;
 
721
        }
 
722
        else if ( are_near(ray(X), ray(Y)) )
 
723
        {
 
724
                Point r = p - center();
 
725
                if ( are_near(r, Point(0,0)) )
 
726
                {
 
727
                        THROW_INFINITESOLUTIONS(0);
 
728
                }
 
729
                // TODO: implement case r != 0
 
730
//              Point np = ray(X) * unit_vector(r);
 
731
//              std::vector<double> solX = roots(np[X],X);
 
732
//              std::vector<double> solY = roots(np[Y],Y);
 
733
//              double t;
 
734
//              if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))
 
735
//              {
 
736
//                      t = solX[0];
 
737
//              }
 
738
//              else
 
739
//              {
 
740
//                      t = solX[1];
 
741
//              }
 
742
//              if ( !(t < from || t > to) )
 
743
//              {
 
744
//                      result.push_back(t);
 
745
//              }
 
746
//              else
 
747
//              {
 
748
//
 
749
//              }
 
750
        }
 
751
 
 
752
        // solve the equation <D(E(t),t)|E(t)-p> == 0
 
753
        // that provides min and max distance points
 
754
        // on the ellipse E wrt the point p
 
755
        // after the substitutions:
 
756
        // cos(t) = (1 - s^2) / (1 + s^2)
 
757
        // sin(t) = 2t / (1 + s^2)
 
758
        // where s = tan(t/2)
 
759
        // we get a 4th degree equation in s
 
760
        /*
 
761
         *      ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +
 
762
         *      ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +
 
763
         *      2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +
 
764
         *      2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])
 
765
         */
 
766
 
 
767
        Point p_c = p - center();
 
768
        double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
 
769
        double cosrot = std::cos( rotation_angle() );
 
770
        double sinrot = std::sin( rotation_angle() );
 
771
        double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
 
772
        Poly coeff;
 
773
        coeff.resize(5);
 
774
        coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );
 
775
        coeff[3] = 2 * ( rx2_ry2 + expr1 );
 
776
        coeff[2] = 0;
 
777
        coeff[1] = 2 * ( -rx2_ry2 + expr1 );
 
778
        coeff[0] = -coeff[4];
 
779
 
 
780
//      for ( unsigned int i = 0; i < 5; ++i )
 
781
//              std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;
 
782
 
 
783
        std::vector<double> real_sol;
 
784
        // gsl_poly_complex_solve raises an error
 
785
        // if the leading coefficient is zero
 
786
        if ( are_near(coeff[4], 0) )
 
787
        {
 
788
                real_sol.push_back(0);
 
789
                if ( !are_near(coeff[3], 0) )
 
790
                {
 
791
                        double sq = -coeff[1] / coeff[3];
 
792
                        if ( sq > 0 )
 
793
                        {
 
794
                                double s = std::sqrt(sq);
 
795
                                real_sol.push_back(s);
 
796
                                real_sol.push_back(-s);
 
797
                        }
 
798
                }
 
799
        }
 
800
        else
 
801
        {
 
802
                real_sol = solve_reals(coeff);
 
803
        }
 
804
//      else
 
805
//      {
 
806
//              double sol[8];
 
807
//              gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc(5);
 
808
//              gsl_poly_complex_solve(coeff, 5, w, sol );
 
809
//              gsl_poly_complex_workspace_free(w);
 
810
//
 
811
//              for ( unsigned int i = 0; i < 4; ++i )
 
812
//              {
 
813
//                      if ( sol[2*i+1] == 0 ) real_sol.push_back(sol[2*i]);
 
814
//              }
 
815
//      }
 
816
 
 
817
        for ( unsigned int i = 0; i < real_sol.size(); ++i )
 
818
        {
 
819
                real_sol[i] = 2 * std::atan(real_sol[i]);
 
820
                if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;
 
821
        }
 
822
        // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0
 
823
        // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity
 
824
        if ( (real_sol.size() % 2) != 0 )
 
825
        {
 
826
                real_sol.push_back(M_PI);
 
827
        }
 
828
 
 
829
        double mindistsq1 = std::numeric_limits<double>::max();
 
830
        double mindistsq2 = std::numeric_limits<double>::max();
 
831
        double dsq;
 
832
        unsigned int mi1, mi2;
 
833
        for ( unsigned int i = 0; i < real_sol.size(); ++i )
 
834
        {
 
835
                dsq = distanceSq(p, pointAtAngle(real_sol[i]));
 
836
                if ( mindistsq1 > dsq )
 
837
                {
 
838
                        mindistsq2 = mindistsq1;
 
839
                        mi2 = mi1;
 
840
                        mindistsq1 = dsq;
 
841
                        mi1 = i;
 
842
                }
 
843
                else if ( mindistsq2 > dsq )
 
844
                {
 
845
                        mindistsq2 = dsq;
 
846
                        mi2 = i;
 
847
                }
 
848
        }
 
849
 
 
850
        double t = map_to_01( real_sol[mi1] );
 
851
        if ( !(t < from || t > to) )
 
852
        {
 
853
                result.push_back(t);
 
854
        }
 
855
 
 
856
        bool second_sol = false;
 
857
        t = map_to_01( real_sol[mi2] );
 
858
        if ( real_sol.size() == 4 && !(t < from || t > to) )
 
859
        {
 
860
        if ( result.empty() || are_near(mindistsq1, mindistsq2) )
 
861
        {
 
862
                result.push_back(t);
 
863
                second_sol = true;
 
864
        }
 
865
        }
 
866
 
 
867
        // we need to test extreme points too
 
868
        double dsq1 = distanceSq(p, pointAt(from));
 
869
        double dsq2 = distanceSq(p, pointAt(to));
 
870
        if ( second_sol )
 
871
        {
 
872
                if ( mindistsq2 > dsq1 )
 
873
                {
 
874
                        result.clear();
 
875
                        result.push_back(from);
 
876
                        mindistsq2 = dsq1;
 
877
                }
 
878
                else if ( are_near(mindistsq2, dsq) )
 
879
                {
 
880
                        result.push_back(from);
 
881
                }
 
882
                if ( mindistsq2 > dsq2 )
 
883
                {
 
884
                        result.clear();
 
885
                        result.push_back(to);
 
886
                }
 
887
                else if ( are_near(mindistsq2, dsq2) )
 
888
                {
 
889
                        result.push_back(to);
 
890
                }
 
891
 
 
892
        }
 
893
        else
 
894
        {
 
895
                if ( result.empty() )
 
896
                {
 
897
                        if ( are_near(dsq1, dsq2) )
 
898
                        {
 
899
                                result.push_back(from);
 
900
                                result.push_back(to);
 
901
                        }
 
902
                        else if ( dsq2 > dsq1 )
 
903
                        {
 
904
                                result.push_back(from);
 
905
                        }
 
906
                        else
 
907
                        {
 
908
                                result.push_back(to);
 
909
                        }
 
910
                }
 
911
        }
 
912
 
 
913
        return result;
 
914
}
 
915
 
 
916
 
 
917
} // end namespace Geom
 
918
 
 
919
 
 
920
/*
 
921
  Local Variables:
 
922
  mode:c++
 
923
  c-file-style:"stroustrup"
 
924
  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 
925
  indent-tabs-mode:nil
 
926
  fill-column:99
 
927
  End:
 
928
*/
 
929
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
 
930
 
 
931