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

vigra/project2ellipse.hxx |

    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)