- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# C:/AngioLocalizeNoAuto/MevisLabNoAutoReg/STUMevisApp/Angiolocalize/Sources/ML/STUChamferMatching/amoeha.h

    00001 //----------------------------------------------------------------------------------
    00003 //----------------------------------------------------------------------------------
    00004
    00005 ML_START_NAMESPACE
    00006
    00007 #ifndef __AMOEBA_H
    00008 #define __AMOEBA_H
    00009
    00010
    00011 class Amoeba {
    00012
    00013 public:
    00014         const Doub ftol;
    00015         Int nfunc;
    00016         Int mpts;
    00017         Int ndim;
    00018         Doub fmin;
    00019         VecDoub y;
    00020         MatDoub p;
    00021     Int iterations;
    00022         Int updateFreq;
    00023
    00024         Amoeba(const Doub ftoll, Int iterations, Int updateFreq) : ftol(ftoll), iterations(iterations), updateFreq(updateFreq) {}
    00025
    00026         template <class T>
    00027         VecDoub minimize(VecDoub_I &point, const Doub del, T &func, void (*updateFunc)(void*), bool (*stopFunc)(void*), void *obj)
    00028         {
    00029                 VecDoub dels(point.size(),del);
    00030                 return minimize(point,dels, func, updateFunc, stopFunc, obj);
    00031         }
    00032
    00033         template <class T>
    00034         VecDoub minimize(const VecDoub &point, VecDoub_I &dels, T &func,  void (*updateFunc)(void*), bool (*stopFunc)(void*), void *obj)
    00035         {
    00036                 Int ndim=point.size();
    00037                 MatDoub pp(ndim+1,ndim);
    00038                 for (Int i=0;i<ndim+1;i++) {
    00039                         for (Int j=0;j<ndim;j++)
    00040                                 pp[i][j]=point[j];
    00041                         if (i !=0 ) pp[i][i-1] += dels[i-1];
    00042                 }
    00043                 return minimize(pp,func, updateFunc, stopFunc, obj);
    00044         }
    00045
    00046     template <class T>
    00047         VecDoub minimize(MatDoub_I &pp,  T &func, void (*updateFunc)(void *), bool (*stopFunc)(void*), void *obj)
    00048         {
    00049                 //const Int NMAX=5000;
    00050                 const Int NMAX=iterations;
    00051                 const Doub TINY=1.0e-10;
    00052                 Int ihi,ilo,inhi;
    00053                 mpts=pp.nrows();
    00054                 ndim=pp.ncols();
    00055                 VecDoub psum(ndim),pmin(ndim),x(ndim);
    00056                 p=pp;
    00057                 y.resize(mpts);
    00058                 for (Int i=0;i<mpts;i++) {
    00059                         for (Int j=0;j<ndim;j++)
    00060                                 x[j]=p[i][j];
    00061                         y[i]=func(x, obj);
    00062                 }
    00063                 nfunc=0;
    00064                 get_psum(p,psum);
    00065                 for (;;) {
    00066                         /*
    00067                         if (stopFunc(obj) == true){
    00068                                 updateFunc(obj);
    00069                                 break;
    00070                         }
    00071                         */
    00072                         ilo=0;
    00073                         ihi = y[0]>y[1] ? (inhi=1,0) : (inhi=0,1);
    00074                         for (Int i=0;i<mpts;i++) {
    00075                                 if (y[i] <= y[ilo]) ilo=i;
    00076                                 if (y[i] > y[ihi]) {
    00077                                         inhi=ihi;
    00078                                         ihi=i;
    00079                                 } else if (y[i] > y[inhi] && i != ihi) inhi=i;
    00080                         }
    00081                         Doub rtol=2.0*abs(y[ihi]-y[ilo])/(abs(y[ihi])+abs(y[ilo])+TINY);
    00082                         // ToDO Lina: The current Stop trigger does not work yet
    00083                         // if ((rtol < ftol) || (stopFunc(obj) == true)){
    00084                         if (((rtol < ftol)  == true)|| (nfunc > NMAX)){
    00085                                 SWAP(y[0],y[ilo]);
    00086                                 for (Int i=0;i<ndim;i++) {
    00087                                         SWAP(p[0][i],p[ilo][i]);
    00088                                         pmin[i]=p[0][i];
    00089                                 }
    00090                                 fmin=y[0];
    00091                                 return pmin;
    00092                         }
    00093                         /* ToDO Lina: Does not work yet
    00094                         if (nfunc == updateFreq) {
    00095                                 updateFunc(obj);
    00096
    00097                         }
    00098                         */
    00099
    00100                         nfunc += 2;
    00101                         Doub ytry=amotry(p,y,psum,ihi,-1.0,func, obj);
    00102                         if (ytry <= y[ilo])
    00103                                 ytry=amotry(p,y,psum,ihi,2.0,func, obj);
    00104                         else if (ytry >= y[inhi]) {
    00105                                 Doub ysave=y[ihi];
    00106                                 ytry=amotry(p,y,psum,ihi,0.5,func, obj);
    00107                                 if (ytry >= ysave) {
    00108                                         for (Int i=0;i<mpts;i++) {
    00109                                                 if (i != ilo) {
    00110                                                         for (Int j=0;j<ndim;j++)
    00111                                                                 p[i][j]=psum[j]=0.5*(p[i][j]+p[ilo][j]);
    00112                                                         y[i]=func(psum, obj);
    00113                                                 }
    00114                                         }
    00115                                         nfunc += ndim;
    00116                                         get_psum(p,psum);
    00117                                 }
    00118                         } else --nfunc;
    00119                 }
    00120         }
    00121         inline void get_psum(MatDoub_I &p, VecDoub_O &psum)
    00122         {
    00123                 for (Int j=0;j<ndim;j++) {
    00124                         Doub sum=0.0;
    00125                         for (Int i=0;i<mpts;i++)
    00126                                 sum += p[i][j];
    00127                         psum[j]=sum;
    00128                 }
    00129         }
    00130         template <class T>
    00131         Doub amotry(MatDoub_IO &p, VecDoub_O &y, VecDoub_IO &psum,
    00132                 const Int ihi, const Doub fac, T&func, void *obj)
    00133         {
    00134                 VecDoub ptry(ndim);
    00135                 Doub fac1=(1.0-fac)/ndim;
    00136                 Doub fac2=fac1-fac;
    00137                 for (Int j=0;j<ndim;j++)
    00138                         ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
    00139                 Doub ytry=func(ptry, obj);
    00140                 if (ytry < y[ihi]) {
    00141                         y[ihi]=ytry;
    00142                         for (Int j=0;j<ndim;j++) {
    00143                                 psum[j] += ptry[j]-p[ihi][j];
    00144                                 p[ihi][j]=ptry[j];
    00145                         }
    00146                 }
    00147                 return ytry;
    00148         }
    00149 };
    00150
    00151 #endif
    00152
    00153 ML_END_NAMESPACE

All Classes Namespaces Files Functions Variables Enumerator Friends