00001
00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00003 #pragma implementation
00004 #endif
00005
00006
00007
00008
00009
00010
00011
00012 #include "vnl_cost_function.h"
00013 #include <vcl_cassert.h>
00014
00015 static bool f_calling_compute;
00016
00017 void vnl_cost_function::compute(vnl_vector<double> const& x, double *f, vnl_vector<double>* g)
00018 {
00019 if (f) *f = this->f(x);
00020 if (g) this->gradf(x, *g);
00021 }
00022
00023
00024 double vnl_cost_function::f(vnl_vector<double> const& x)
00025 {
00026
00027 if (f_calling_compute)
00028 assert(!"vnl_cost_function: RECURSION");
00029 double f;
00030 f_calling_compute = true;
00031 this->compute(x, &f, 0);
00032 f_calling_compute = false;
00033 return f;
00034 }
00035
00036
00037 void vnl_cost_function::gradf(vnl_vector<double> const& x, vnl_vector<double>& g)
00038 {
00039 if (f_calling_compute)
00040 assert(!"vnl_cost_function: RECURSION");
00041 f_calling_compute = true;
00042 this->compute(x, 0, &g);
00043 f_calling_compute = false;
00044 }
00045
00046
00047 void vnl_cost_function::fdgradf(vnl_vector<double> const& x,
00048 vnl_vector<double> & gradient,
00049 double stepsize )
00050 {
00051 vnl_vector<double> tx = x;
00052 double h = stepsize;
00053 for (int i = 0; i < dim; ++i) {
00054
00055 double tplus = x[i] + h;
00056 tx[i] = tplus;
00057 double fplus = this->f(tx);
00058
00059 double tminus = x[i] - h;
00060 tx[i] = tminus;
00061 double fminus = this->f(tx);
00062
00063 gradient[i] = (fplus - fminus) / (tplus - tminus);
00064 tx[i] = x[i];
00065 }
00066 }
00067
00068 vnl_vector<double> vnl_cost_function::gradf(vnl_vector<double> const& x)
00069 {
00070 vnl_vector<double> g(dim);
00071 this->gradf(x, g);
00072 return g;
00073 }
00074
00075 vnl_vector<double> vnl_cost_function::fdgradf(vnl_vector<double> const& x)
00076 {
00077 vnl_vector<double> g(dim);
00078 this->fdgradf(x, g);
00079 return g;
00080 }