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

vigra/project2ellipse.hxx VIGRA

00001 
00002 /************************************************************************/
00003 /*                                                                      */
00004 /*                 Copyright 2012 by Frank Lenzen &                     */
00005 /*                                           Ullrich Koethe             */
00006 /*                                                                      */
00007 /*    This file is part of the VIGRA computer vision library.           */
00008 /*    The VIGRA Website is                                              */
00009 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00010 /*    Please direct questions, bug reports, and contributions to        */
00011 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00012 /*        vigra@informatik.uni-hamburg.de                               */
00013 /*                                                                      */
00014 /*    Permission is hereby granted, free of charge, to any person       */
00015 /*    obtaining a copy of this software and associated documentation    */
00016 /*    files (the "Software"), to deal in the Software without           */
00017 /*    restriction, including without limitation the rights to use,      */
00018 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00019 /*    sell copies of the Software, and to permit persons to whom the    */
00020 /*    Software is furnished to do so, subject to the following          */
00021 /*    conditions:                                                       */
00022 /*                                                                      */
00023 /*    The above copyright notice and this permission notice shall be    */
00024 /*    included in all copies or substantial portions of the             */
00025 /*    Software.                                                         */
00026 /*                                                                      */
00027 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00028 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00029 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00030 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00031 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00032 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00033 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00034 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */                
00035 /*                                                                      */
00036 /************************************************************************/
00037 
00038 #ifndef VIGRA_PROJECT2ELLIPSE_HXX
00039 #define VIGRA_PROJECT2ELLIPSE_HXX
00040 #include <iostream>
00041 
00042 namespace vigra{
00043 
00044   namespace detail{
00045 
00046 //---------------------------------------------------------------------------
00047 
00048 void projectEllipse2D_FirstQuad (double &vx, double &vy, double a, double b, const double eps, const int iter_max){                
00049   
00050   double t=0,tmin,tmax,d1,d2,f,x1,y1,l;                      
00051   int i;
00052   tmax=std::max(2*a*vx,2*b*vy);
00053   tmin=-b*b;
00054   d1=a*vx/(tmax+a*a);
00055   d2=b*vy/(tmax+b*b);
00056   f=d1*d1+d2*d2-1;
00057     
00058   for (i=0;i<iter_max;i++){
00059   
00060     t=.5*(tmin+tmax);
00061     d1=a*vx/(t+a*a);
00062     d2=b*vy/(t+b*b);
00063     f=d1*d1+d2*d2-1;
00064     x1=a*vx/(t+a*a);
00065     y1=b*vy/(t+b*b);
00066     l=x1*x1+y1*y1-1;
00067    
00068     if (fabs(l)<eps)
00069       break;
00070     if(f>0)
00071       tmin=t;
00072     else if(f<0)
00073       tmax=t;
00074     else
00075       break;
00076   }
00077   d1=vx;
00078   d2=vy;
00079   vx=a*a*vx/(t+a*a);
00080   vy=b*b*vy/(t+b*b);
00081   d1 = vx - d1;
00082   d2 = vy - d2;
00083   return;
00084 }
00085 
00086 
00087 void projectEllipse2D(double &vx, double &vy, const double _a,const double _b,const double eps,const int max_iter){
00088   
00089   //double err;
00090   double a=_a,b=_b;
00091   
00092   //check if inside ellipse
00093   if (((vx/a)*(vx/a)+(vy/b)*(vy/b))<=1){
00094     return;
00095   }
00096   
00097   // special case of a circle
00098   if (fabs(a-b) < eps){
00099     double l = sqrt(vx*vx+vy*vy);
00100     if (l>(a+b)/2.){
00101       vx=(a+b)/(2*l)*vx;
00102       vy=(a+b)/(2*l)*vy;
00103     }
00104     return;
00105   }
00106   
00107   // reflect vx -> -vx, if necessary
00108   bool x_reflect;
00109   if (vx > eps){
00110     x_reflect = false;
00111   }
00112   else if (vx < -eps){
00113     x_reflect = true;
00114     vx = -vx;
00115   }
00116   else{
00117     x_reflect = false;
00118     vx = 0.0;
00119   }
00120   // reflect vy -> vy = -V if necessary
00121   bool y_reflect;
00122   if (vy > eps){
00123     y_reflect = false;
00124   }
00125   else if (vy < -eps){
00126     y_reflect = true;
00127     vy = -vy;
00128   }
00129   else{
00130     y_reflect = false;
00131     vy = 0.0;
00132   }
00133   
00134   // swap axes if necessary
00135   bool swapped;
00136   if (a >= b){
00137     swapped = false;
00138   }
00139   else{
00140     swapped = true;
00141     std::swap(a,b);
00142     std::swap(vx,vy);
00143   }
00144   if (vx != 0.0){
00145     if (vy != 0.0){
00146       projectEllipse2D_FirstQuad(vx,vy,a,b,eps,max_iter);
00147     }
00148     else{
00149       if (vx < a - b*b/a){
00150         double vx_temp = vx;
00151         vx = a*a*vx/(a*a-b*b);
00152         vy = vy*sqrt(fabs(1.0-vx_temp/a*vx_temp/a));
00153       }
00154       else{
00155         vx = a;
00156         vy = 0.0;
00157       }
00158     }
00159   }
00160   else{
00161     vx = 0.0;
00162     vy = b;
00163   }
00164   if (swapped){
00165     std::swap(vx,vy);
00166   }
00167   if (y_reflect){
00168     vy = -vy;
00169   }
00170   if (x_reflect){
00171     vx = -vx;
00172   }
00173   return;
00174 }
00175   } //closing namespace detail
00176 } //closing namespace vigra
00177 #endif // VIGRA_PROJECT2ELLIPSE_HXX

© 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.9.0 (Tue Nov 6 2012)