kig

centerofcurvature_type.cc

00001 // Copyright (C)  2004  Maurizio Paolini <paolini@dmf.unicatt.it>
00002 
00003 // This program is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU General Public License
00005 // as published by the Free Software Foundation; either version 2
00006 // of the License, or (at your option) any later version.
00007 
00008 // This program is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00011 // GNU General Public License for more details.
00012 
00013 // You should have received a copy of the GNU General Public License
00014 // along with this program; if not, write to the Free Software
00015 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
00016 // 02110-1301, USA.
00017 
00018 #include "centerofcurvature_type.h"
00019 
00020 #include "bogus_imp.h"
00021 #include "conic_imp.h"
00022 #include "cubic_imp.h"
00023 //#include "other_imp.h"
00024 #include "point_imp.h"
00025 //#include "line_imp.h"
00026 
00027 #include "../misc/common.h"
00028 #include "../misc/conic-common.h"
00029 //#include "../misc/calcpaths.h"
00030 #include "../kig/kig_part.h"
00031 #include "../kig/kig_view.h"
00032 
00033 static const char constructcenterofcurvaturepoint[] = "SHOULDNOTBESEEN";
00034 //  I18N_NOOP( "Construct the center of curvature corresponding to this point" );
00035 static const char selectcoc1[] = I18N_NOOP( "Select the curve..." );
00036 static const char selectcoc2[] = I18N_NOOP( "Select a point on the curve..." );
00037 
00038 static const ArgsParser::spec argsspecCocConic[] =
00039 {
00040   { ConicImp::stype(), "SHOULDNOTBESEEN", selectcoc1, false },
00041   { PointImp::stype(), constructcenterofcurvaturepoint, selectcoc2, false }
00042 };
00043 
00044 KIG_INSTANTIATE_OBJECT_TYPE_INSTANCE( CocConicType )
00045 
00046 CocConicType::CocConicType()
00047   : ArgsParserObjectType( "CocConic", argsspecCocConic, 2 )
00048 {
00049 }
00050 
00051 CocConicType::~CocConicType()
00052 {
00053 }
00054 
00055 const CocConicType* CocConicType::instance()
00056 {
00057   static const CocConicType t;
00058   return &t;
00059 }
00060 
00061 ObjectImp* CocConicType::calc( const Args& args, const KigDocument& doc ) const
00062 {
00063   if ( !margsparser.checkArgs( args ) )
00064     return new InvalidImp;
00065 
00066   const ConicImp* conic = static_cast<const ConicImp*>( args[0] );
00067   const Coordinate& p = static_cast<const PointImp*>( args[1] )->coordinate();
00068 
00069   if ( !conic->containsPoint( p, doc ) )
00070     return new InvalidImp;
00071 
00072   double x = p.x;
00073   double y = p.y;
00074   ConicCartesianData data = conic->cartesianData();
00075 //  double aconst = data.coeffs[5];
00076   double ax = data.coeffs[3];
00077   double ay = data.coeffs[4];
00078   double axx = data.coeffs[0];
00079   double axy = data.coeffs[2];
00080   double ayy = data.coeffs[1];
00081 
00082 /*
00083  * mp: we need to compute the normal vector and the curvature
00084  * of the curve.  The curve (conic) is given in implicit form
00085  * f(x,y) = 0;  the gradient of f gives the direction of the
00086  * normal; for the curvature we can use the following formula:
00087  * k = div(grad f/|grad f|)
00088  *
00089  * the hessian matrix has elements [hfxx, hfxy]
00090  *                                 [hfxy, hfyy]
00091  *
00092  * kgf is the curvature multiplied by the norm of grad f
00093  */
00094 
00095   double gradfx = 2*axx*x + axy*y + ax;
00096   double gradfy = axy*x + 2*ayy*y + ay;
00097   Coordinate gradf = Coordinate ( gradfx, gradfy );
00098 
00099   double hfxx = 2*axx;
00100   double hfyy = 2*ayy;
00101   double hfxy = axy;
00102 
00103   double kgf = hfxx + hfyy 
00104      - (hfxx*gradfx*gradfx + hfyy*gradfy*gradfy + 2*hfxy*gradfx*gradfy)
00105        /(gradfx*gradfx + gradfy*gradfy);
00106 
00107   bool ok = true;
00108 
00109   const Coordinate coc = p - 1/kgf*gradf;
00110   
00111   if ( !ok )
00112     return new InvalidImp;
00113   
00114   return new PointImp( coc );
00115 }
00116 
00117 const ObjectImpType* CocConicType::resultId() const
00118 {
00119   return PointImp::stype();
00120 }
00121 
00122 /**** Cubic starts here ****/
00123 
00124 static const ArgsParser::spec argsspecCocCubic[] =
00125 {
00126   { CubicImp::stype(), "SHOULDNOTBESEEN", selectcoc1, false },
00127   { PointImp::stype(), constructcenterofcurvaturepoint, selectcoc2, false }
00128 };
00129 
00130 KIG_INSTANTIATE_OBJECT_TYPE_INSTANCE( CocCubicType )
00131 
00132 CocCubicType::CocCubicType()
00133   : ArgsParserObjectType( "CocCubic", argsspecCocCubic, 2 )
00134 {
00135 }
00136 
00137 CocCubicType::~CocCubicType()
00138 {
00139 }
00140 
00141 const CocCubicType* CocCubicType::instance()
00142 {
00143   static const CocCubicType t;
00144   return &t;
00145 }
00146 
00147 ObjectImp* CocCubicType::calc( const Args& args, const KigDocument& doc ) const
00148 {
00149   if ( !margsparser.checkArgs( args ) )
00150     return new InvalidImp;
00151 
00152   const CubicImp* cubic = static_cast<const CubicImp*>( args[0] );
00153   const Coordinate& p = static_cast<const PointImp*>( args[1] )->coordinate();
00154 
00155   if ( !cubic->containsPoint( p, doc ) )
00156     return new InvalidImp;
00157 
00158   double x = p.x;
00159   double y = p.y;
00160   CubicCartesianData data = cubic->data();
00161 //  double aconst = data.coeffs[0];
00162   double ax = data.coeffs[1];
00163   double ay = data.coeffs[2];
00164   double axx = data.coeffs[3];
00165   double axy = data.coeffs[4];
00166   double ayy = data.coeffs[5];
00167   double axxx = data.coeffs[6];
00168   double axxy = data.coeffs[7];
00169   double axyy = data.coeffs[8];
00170   double ayyy = data.coeffs[9];
00171 
00172   /*
00173    * we use here the same mechanism as for the
00174    * conics, see above
00175    */
00176 
00177   double gradfx = 3*axxx*x*x + 2*axxy*x*y + axyy*y*y + 2*axx*x + axy*y + ax;
00178   double gradfy = axxy*x*x + 2*axyy*x*y + 3*ayyy*y*y + axy*x + 2*ayy*y + ay;
00179   Coordinate gradf = Coordinate ( gradfx, gradfy );
00180 
00181   double hfxx = 6*axxx*x + 2*axxy*y + 2*axx;
00182   double hfyy = 6*ayyy*y + 2*axyy*x + 2*ayy;
00183   double hfxy = 2*axxy*x + 2*axyy*y + axy;
00184 
00185   double kgf = hfxx + hfyy 
00186      - (hfxx*gradfx*gradfx + hfyy*gradfy*gradfy + 2*hfxy*gradfx*gradfy)
00187        /(gradfx*gradfx + gradfy*gradfy);
00188 
00189   bool ok = true;
00190 
00191   const Coordinate coc = p - 1/kgf*gradf;
00192 
00193   if ( !ok )
00194     return new InvalidImp;
00195 
00196   return new PointImp( coc );  
00197 }
00198 
00199 const ObjectImpType* CocCubicType::resultId() const
00200 {
00201   return PointImp::stype();
00202 }
00203 
00204 /**** Curve starts here ****/
00205 
00206 static const ArgsParser::spec argsspecCocCurve[] =
00207 {
00208   { CurveImp::stype(), "SHOULDNOTBESEEN", selectcoc1, false },
00209   { PointImp::stype(), constructcenterofcurvaturepoint, selectcoc2, false }
00210 };
00211 
00212 KIG_INSTANTIATE_OBJECT_TYPE_INSTANCE( CocCurveType )
00213 
00214 CocCurveType::CocCurveType()
00215   : ArgsParserObjectType( "CocCurve", argsspecCocCurve, 2 )
00216 {
00217 }
00218 
00219 CocCurveType::~CocCurveType()
00220 {
00221 }
00222 
00223 const CocCurveType* CocCurveType::instance()
00224 {
00225   static const CocCurveType t;
00226   return &t;
00227 }
00228 
00229 ObjectImp* CocCurveType::calc( const Args& args, const KigDocument& doc ) const
00230 {
00231   if ( !margsparser.checkArgs( args ) )
00232     return new InvalidImp;
00233 
00234   const CurveImp* curve = static_cast<const CurveImp*>( args[0] );
00235   const Coordinate& p = static_cast<const PointImp*>( args[1] )->coordinate();
00236 
00237   if ( !curve->containsPoint( p, doc ) )
00238     return new InvalidImp;
00239 
00240 
00241   const double t = curve->getParam( p, doc );
00242   const double tau0 = 5e-4;
00243   const double sigmasq = 1e-12;
00244   const int maxiter = 20;
00245 
00246   double tau = tau0;
00247   Coordinate gminus, g, gplus, tang, acc, curv, err;
00248   double velsq, curvsq;
00249   double tplus = t + tau;
00250   double tminus = t - tau;
00251   double t0 = t;
00252   if ( tplus > 1 ) {tplus = 1; t0 = 1 - tau; tminus = 1 - 2*tau;}
00253   if ( tminus < 0 ) {tminus = 0; t0 = tau; tplus = 2*tau;}
00254   gminus = curve->getPoint( tminus, doc );
00255   g = curve->getPoint( t0, doc );
00256   gplus = curve->getPoint( tplus, doc );
00257   tang = (gplus - gminus)/(2*tau);
00258   acc = (gminus + gplus - 2*g)/(tau*tau);
00259   velsq = tang.x*tang.x + tang.y*tang.y;
00260   tang = tang/velsq;
00261   Coordinate curvold = acc/velsq - (acc.x*tang.x + acc.y*tang.y)*tang;
00262   curvsq = curvold.x*curvold.x + curvold.y*curvold.y;
00263   curvold = curvold/curvsq;
00264 
00265   for (int i = 0; i < maxiter; i++)
00266   {
00267     tau = tau/2;
00268     tplus = t + tau;
00269     tminus = t - tau;
00270     t0 = t;
00271     if ( tplus > 1 ) {tplus = 1; t0 = 1 - tau; tminus = 1 - 2*tau;}
00272     if ( tminus < 0 ) {tminus = 0; t0 = tau; tplus = 2*tau;}
00273 
00274     gminus = curve->getPoint( tminus, doc );
00275     g = curve->getPoint( t0, doc );
00276     gplus = curve->getPoint( tplus, doc );
00277     tang = (gplus - gminus)/(2*tau);
00278     acc = (gminus + gplus - 2*g)/(tau*tau);
00279     velsq = tang.x*tang.x + tang.y*tang.y;
00280     tang = tang/velsq;
00281     curv = acc/velsq - (acc.x*tang.x + acc.y*tang.y)*tang;
00282     curvsq = curv.x*curv.x + curv.y*curv.y;
00283     curv = curv/curvsq;
00284 
00285     err = (curvold - curv)/3;
00286     /*
00287      * curvsq is the inverse squared of the norm of curvsq
00288      * so this is actually a relative test
00289      * in the end we return an extrapolated value
00290      */
00291     if (err.squareLength() < sigmasq/curvsq)
00292     {
00293       curv = (4*curv - curvold)/3;
00294       return new PointImp( p + curv );
00295     }
00296     curvold = curv;
00297   }
00298   return new InvalidImp;
00299 }
00300 
00301 const ObjectImpType* CocCurveType::resultId() const
00302 {
00303   return PointImp::stype();
00304 }
kig

centerofcurvature_type.cc

kig

kdeedu
