[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

project2ellipse.hxx
1
2/************************************************************************/
3/* */
4/* Copyright 2012 by Frank Lenzen & */
5/* Ullrich Koethe */
6/* */
7/* This file is part of the VIGRA computer vision library. */
8/* The VIGRA Website is */
9/* http://hci.iwr.uni-heidelberg.de/vigra/ */
10/* Please direct questions, bug reports, and contributions to */
11/* ullrich.koethe@iwr.uni-heidelberg.de or */
12/* vigra@informatik.uni-hamburg.de */
13/* */
14/* Permission is hereby granted, free of charge, to any person */
15/* obtaining a copy of this software and associated documentation */
16/* files (the "Software"), to deal in the Software without */
17/* restriction, including without limitation the rights to use, */
18/* copy, modify, merge, publish, distribute, sublicense, and/or */
19/* sell copies of the Software, and to permit persons to whom the */
20/* Software is furnished to do so, subject to the following */
21/* conditions: */
22/* */
23/* The above copyright notice and this permission notice shall be */
24/* included in all copies or substantial portions of the */
25/* Software. */
26/* */
27/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
28/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
29/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
30/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
31/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
32/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
33/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
34/* OTHER DEALINGS IN THE SOFTWARE. */
35/* */
36/************************************************************************/
37
38#ifndef VIGRA_PROJECT2ELLIPSE_HXX
39#define VIGRA_PROJECT2ELLIPSE_HXX
40#include <iostream>
41
42namespace vigra{
43
44 namespace detail{
45
46//---------------------------------------------------------------------------
47
48inline void projectEllipse2D_FirstQuad (double &vx, double &vy, double a, double b, const double eps, const int iter_max){
49
50 double t=0,tmin,tmax,d1,d2,f,x1,y1,l;
51 int i;
52 tmax=std::max(2*a*vx,2*b*vy);
53 tmin=-b*b;
54 d1=a*vx/(tmax+a*a);
55 d2=b*vy/(tmax+b*b);
56 f=d1*d1+d2*d2-1;
57
58 for (i=0;i<iter_max;i++){
59
60 t=.5*(tmin+tmax);
61 d1=a*vx/(t+a*a);
62 d2=b*vy/(t+b*b);
63 f=d1*d1+d2*d2-1;
64 x1=a*vx/(t+a*a);
65 y1=b*vy/(t+b*b);
66 l=x1*x1+y1*y1-1;
67
68 if (fabs(l)<eps)
69 break;
70 if(f>0)
71 tmin=t;
72 else if(f<0)
73 tmax=t;
74 else
75 break;
76 }
77 d1=vx;
78 d2=vy;
79 vx=a*a*vx/(t+a*a);
80 vy=b*b*vy/(t+b*b);
81 d1 = vx - d1;
82 d2 = vy - d2;
83 return;
84}
85
86
87inline void projectEllipse2D(double &vx, double &vy, const double _a,const double _b,const double eps,const int max_iter){
88
89 //double err;
90 double a=_a,b=_b;
91
92 //check if inside ellipse
93 if (((vx/a)*(vx/a)+(vy/b)*(vy/b))<=1){
94 return;
95 }
96
97 // special case of a circle
98 if (fabs(a-b) < eps){
99 double l = sqrt(vx*vx+vy*vy);
100 if (l>(a+b)/2.){
101 vx=(a+b)/(2*l)*vx;
102 vy=(a+b)/(2*l)*vy;
103 }
104 return;
105 }
106
107 // reflect vx -> -vx, if necessary
108 bool x_reflect;
109 if (vx > eps){
110 x_reflect = false;
111 }
112 else if (vx < -eps){
113 x_reflect = true;
114 vx = -vx;
115 }
116 else{
117 x_reflect = false;
118 vx = 0.0;
119 }
120 // reflect vy -> vy = -V if necessary
121 bool y_reflect;
122 if (vy > eps){
123 y_reflect = false;
124 }
125 else if (vy < -eps){
126 y_reflect = true;
127 vy = -vy;
128 }
129 else{
130 y_reflect = false;
131 vy = 0.0;
132 }
133
134 // swap axes if necessary
135 bool swapped;
136 if (a >= b){
137 swapped = false;
138 }
139 else{
140 swapped = true;
141 std::swap(a,b);
142 std::swap(vx,vy);
143 }
144 if (vx != 0.0){
145 if (vy != 0.0){
146 projectEllipse2D_FirstQuad(vx,vy,a,b,eps,max_iter);
147 }
148 else{
149 if (vx < a - b*b/a){
150 double vx_temp = vx;
151 vx = a*a*vx/(a*a-b*b);
152 vy = vy*sqrt(fabs(1.0-vx_temp/a*vx_temp/a));
153 }
154 else{
155 vx = a;
156 vy = 0.0;
157 }
158 }
159 }
160 else{
161 vx = 0.0;
162 vy = b;
163 }
164 if (swapped){
165 std::swap(vx,vy);
166 }
167 if (y_reflect){
168 vy = -vy;
169 }
170 if (x_reflect){
171 vx = -vx;
172 }
173 return;
174}
175 } //closing namespace detail
176} //closing namespace vigra
177#endif // VIGRA_PROJECT2ELLIPSE_HXX
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1