diff --git a/histogram.cc b/histogram.cc
index 5fe931f0c361179a67d94f6236f1c386e044d629..149b6e50908043eec4f74963945e7fa628135fa5 100644
--- a/histogram.cc
+++ b/histogram.cc
@@ -14,13 +14,13 @@ using namespace NEWMAT;
namespace MISCMATHS {
- float Histogram::getPercentile(float perc)
+ float Histogram::getPercentile(float perc)
{
if(histogram.Nrows()==0) generate();
-
+
generateCDF();
float percentile=getValue(1);
-
+
for (int i=2;i<=CDF.Nrows();i++){
if(CDF(i)>perc){
double diff=(CDF(i)-perc)/(CDF(i)-CDF(i-1));
@@ -35,9 +35,9 @@ namespace MISCMATHS {
void Histogram::generate()
{
Tracer ts("Histogram::generate");
-
+
int size = sourceData.Nrows();
-
+
if(calcRange)
{
// calculate range automatically
@@ -50,21 +50,21 @@ namespace MISCMATHS {
histMin=sourceData(i);
}
}
-
+
// zero histogram
histogram.ReSize(bins);
histogram=0;
-
+
// create histogram; the MIN is so that the maximum value falls in the
// last valid bin, not the (last+1) bin
for(int i=1; i<=size; i++)
- {
+ {
histogram(getBin(sourceData(i)))++;
}
datapoints=size;
}
-
+
void Histogram::generate(ColumnVector exclude)
{
Tracer ts("Histogram::generate");
@@ -81,13 +81,13 @@ namespace MISCMATHS {
histMin=sourceData(i);
}
}
-
-
+
+
histogram.ReSize(bins);
histogram=0;
-
+
datapoints=size;
-
+
for(int i=1; i<=size; i++)
{
if(exclude(i)>0){
@@ -95,7 +95,7 @@ namespace MISCMATHS {
}else datapoints--;
}
-
+
}
@@ -108,19 +108,19 @@ namespace MISCMATHS {
// smooth in i direction
newhist=0;
- ColumnVector kernel(3);
+ ColumnVector kernel(3);
// corresponds to Gaussian with sigma=0.8 voxels
// kernel(1)=0.5;
- // kernel(2)=0.2283;
+ // kernel(2)=0.2283;
// kernel(3)=0.0219;
// corresponds to Gaussian with sigma=0.6 voxels
// kernel(1)=0.6638;
- // kernel(2)=0.1655;
+ // kernel(2)=0.1655;
// kernel(3)=0.0026;
//gauss(0.5,5,1)
kernel(1)=0.7866;
- kernel(2)=0.1065;
+ kernel(2)=0.1065;
kernel(3)=0.0003;
for(int i=1; i<=bins; i++)
@@ -136,7 +136,7 @@ namespace MISCMATHS {
if(i>2)
{
val+=kernel(3)*(histogram(i-2));
- norm+=kernel(3);
+ norm+=kernel(3);
}
if(i<bins)
{
@@ -146,7 +146,7 @@ namespace MISCMATHS {
if(i<bins-1)
{
val+=kernel(3)*(histogram(i+2));
- norm+=kernel(3);
+ norm+=kernel(3);
}
val/=norm;
@@ -188,7 +188,7 @@ namespace MISCMATHS {
void Histogram::match(Histogram &target){
-
+
int bin, newbin;
double cdfval, val,dist;
ColumnVector CDF_ref;
@@ -200,52 +200,52 @@ namespace MISCMATHS {
olddata=sourceData;
histnew=0; dist=0;
if(exclusion.Nrows()==0){cout << " Histogram::match no excl " << endl; exclusion.ReSize(sourceData.Nrows()); exclusion=1;}
-
+
float binwidthtarget=(target.getHistMax() - target.getHistMin())/target.getNumBins();
-
-
+
+
for (int i=1;i<=sourceData.Nrows();i++){
-
+
if(exclusion(i)>0){
bin=getBin(sourceData(i));
-
+
cdfval=CDF(bin);
newbin=1;
-
+
val=sourceData(i)-histMin;
-
+
if(bin==target.getNumBins()){
newbin=bin;
}else if( (cdfval- CDF_ref(bin)) > 1e-20){
-
-
+
+
for (int j=bin;j<target.getNumBins();j++){
-
+
if((cdfval - CDF_ref(j)) > 1e-20 && (cdfval - CDF_ref(j+1)) <= 1e-20 ){
newbin=j+1;
dist=(cdfval-CDF_ref(j))/(CDF_ref(j+1)-CDF_ref(j)); // find distance between current bins as proportion of bin width
break;
}
}
-
+
}
else{
-
+
//search backwards
-
-
-
+
+
+
int j=bin-1;
while (j>0){
-
+
// if(sourceData(i)<1.2 && sourceData(i)>1)
if((cdfval - CDF_ref(1)) < -1e-20 && fabs(CDF_ref(bin) - CDF_ref(1)) < 1e-20){newbin=j+1; break; }
else if ((cdfval - CDF_ref(1)) > -1e-20 && fabs(CDF_ref(bin) - CDF_ref(1)) < 1e-20){newbin=j+1; break; }
else if((cdfval - CDF_ref(j)) > 1e-20 && (cdfval - CDF_ref(j+1)) <= 1e-20 ){
newbin=j+1;
-
+
dist=(cdfval-CDF_ref(j))/(CDF_ref(j+1)-CDF_ref(j)); // find distance between current bins as proportion of bin width
break;
}
@@ -253,51 +253,17 @@ namespace MISCMATHS {
}
//search forwards
}
-
-
+
+
sourceData(i)=target.getHistMin() + (newbin-1)*binwidthtarget+dist*binwidthtarget;
histnew(newbin)++;
-
+
if(sourceData(i) < target.getHistMin()) sourceData(i) = target.getHistMin();
if(sourceData(i) > target.getHistMax()) sourceData(i) = target.getHistMax();
}
}
}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+}
diff --git a/kernel.cc b/kernel.cc
index 5f438e1bc9dcc71e65466961b44ef0ccc3bc9d1d..60381b1714ef6163629693a8e93084d01015f1dd 100644
--- a/kernel.cc
+++ b/kernel.cc
@@ -13,7 +13,7 @@ using namespace NEWMAT;
namespace MISCMATHS {
- set<kernelstorage*, kernelstorage::comparer> kernel::existingkernels;
+ set<kernelstorage*, kernelstorage::comparer> kernel::existingkernels;
//////// Support function /////////
@@ -28,42 +28,42 @@ namespace MISCMATHS {
dn -= n;
if (n>(kernel.Nrows()-1)) return 0.0;
if (n<1) return 0.0;
-
+
return kernel(n)*(1.0-dn) + kernel(n+1)*dn;
}
-
+
inline bool in_bounds(const ColumnVector& data, int index)
{ return ( (index>=1) && (index<=data.Nrows())); }
-
+
inline bool in_bounds(const ColumnVector& data, float index)
{ return ( ((int)floor(index)>=1) && ((int)ceil(index)<=data.Nrows())); }
-
+
float sincfn(float x)
{
if (fabs(x)<1e-7) { return 1.0-fabs(x); }
float y=M_PI*x;
return sin(y)/y;
}
-
+
float hanning(float x, int w)
{ // w is half-width
- if (fabs(x)>w)
+ if (fabs(x)>w)
return 0.0;
else
return (0.5 + 0.5 *cos(M_PI*x/w));
}
-
+
float blackman(float x, int w)
{ // w is half-width
- if (fabs(x)>w)
+ if (fabs(x)>w)
return 0.0;
else
return (0.42 + 0.5 *cos(M_PI*x/w) + 0.08*cos(2.0*M_PI*x/w));
}
-
+
float rectangular(float x, int w)
{ // w is half-width
- if (fabs(x)>w)
+ if (fabs(x)>w)
return 0.0;
else
return 1.0;
@@ -94,8 +94,8 @@ namespace MISCMATHS {
return ker;
}
-
- kernel sinckernel(const string& sincwindowtype, int w, int nstore)
+
+ kernel sinckernel(const string& sincwindowtype, int w, int nstore)
{
kernel sinck;
sinck = sinckernel(sincwindowtype,w,w,w,nstore);
@@ -119,7 +119,7 @@ namespace MISCMATHS {
kx = sinckernel1D(sincwindowtype,wx,nstore);
ky = sinckernel1D(sincwindowtype,wy,nstore);
kz = sinckernel1D(sincwindowtype,wz,nstore);
-
+
sinckern.setkernel(kx,ky,kz,hwx,hwy,hwz);
return sinckern;
}
@@ -140,13 +140,13 @@ namespace MISCMATHS {
return extrapval;
}
-
+
// basic trilinear call
float interpolate_1d(const ColumnVector& data, const float index)
{
float interpval;
int low_bound = (int)floor(index);
- int high_bound = (int)ceil(index);
+ int high_bound = (int)ceil(index);
if (in_bounds(data, index))
interpval = data(low_bound) + (index - low_bound)*(data(high_bound) - data(low_bound));
@@ -156,7 +156,7 @@ namespace MISCMATHS {
return interpval;
}
-
+
//////// Spline Support /////////
float hermiteinterpolation_1d(const ColumnVector& data, int p1, int p4, float t)
@@ -165,7 +165,7 @@ namespace MISCMATHS {
// inputs: points P_1, P_4; tangents R_1, R_4; interpolation index t;
float retval, r1 = 0.0, r4 = 0.0;
-
+
if (!in_bounds(data,p1) || !in_bounds(data,p4)) {
cerr << "Hermite Interpolation - ERROR: One or more indicies lie outside the data range. Returning ZERO" << endl;
retval = 0.0;
@@ -176,17 +176,17 @@ namespace MISCMATHS {
retval = data(p1);
} else if (t == 1.0) {
retval = data(p4);
- */
+ */
} else {
r1 = 0.5 * (extrapolate_1d(data, p1) - extrapolate_1d(data, p1 - 1)) + 0.5 * (extrapolate_1d(data, p1 + 1) - extrapolate_1d(data, p1));// tangent @ P_1
r4 = 0.5 * (extrapolate_1d(data, p4) - extrapolate_1d(data, p4 - 1)) + 0.5 * (extrapolate_1d(data, p4 + 1) - extrapolate_1d(data, p4));// tangent @ P_4
-
+
float t2 = t*t; float t3 = t2*t;
retval = (2*t3 - 3*t2 + 1)*data(p1) + (-2*t3 + 3*t2)*data(p4) + (t3 - 2*t2 + t)*r1 + (t3 - t2)*r4;
}
// cerr << "p1, p4, t, r1, r4 = " << p1 << ", " << p4 << ", " << t << ", " << r1 << ", " << r4 << endl;
-
+
return retval;
}
@@ -199,14 +199,14 @@ namespace MISCMATHS {
// kernel half-width (i.e. range is +/- w)
int ix0;
ix0 = (int) floor(index);
-
+
int wx(widthx);
vector<float> storex(2*wx+1);
- for (int d=-wx; d<=wx; d++)
+ for (int d=-wx; d<=wx; d++)
storex[d+wx] = kernelval((index-ix0+d),wx,userkernel);
float convsum=0.0, interpval=0.0, kersum=0.0;
-
+
int xj;
for (int x1=ix0-wx; x1<=ix0+wx; x1++) {
if (in_bounds(data, x1)) {
@@ -224,7 +224,7 @@ namespace MISCMATHS {
}
return interpval;
}
-
+
////// Kernel wrappers //////
float kernelinterpolation_1d(const ColumnVector& data, float index)
@@ -239,4 +239,3 @@ namespace MISCMATHS {
return kernelinterpolation_1d(data.t(), index, userkernel, 7);
}
}
-
diff --git a/optimise.cc b/optimise.cc
index 22c308ed75fb6b1d6ea878153f280a11bf34865f..80005ce750f13f624c28068ae5487926bf028920 100644
--- a/optimise.cc
+++ b/optimise.cc
@@ -28,7 +28,7 @@ namespace MISCMATHS {
using std::exp;
using std::log;
- bool estquadmin(float &xnew, float x1, float xmid, float x2,
+ bool estquadmin(float &xnew, float x1, float xmid, float x2,
float y1, float ymid, float y2)
{
// Finds the estimated quadratic minimum's position
@@ -64,7 +64,7 @@ namespace MISCMATHS {
}
return xnew;
}
-
+
float nextpt(float x1, float xmid, float x2, float y1, float ymid, float y2)
@@ -82,10 +82,10 @@ namespace MISCMATHS {
return xnew;
}
-
- void findinitialbound(float &x1, float &xmid, float &x2,
- float &y1, float &ymid, float &y2,
+
+ void findinitialbound(float &x1, float &xmid, float &x2,
+ float &y1, float &ymid, float &y2,
float (*func)(const ColumnVector &),
const ColumnVector &unitdir, const ColumnVector &pt)
{
@@ -109,15 +109,15 @@ namespace MISCMATHS {
y2 = (*func)(x2*unitdir + pt);
while (ymid > y2) { // note: must maintain y1 >= ymid
-
- // cout << " <" << Min(x1,x2) << "," << xmid
+
+ // cout << " <" << Min(x1,x2) << "," << xmid
// << "," << Max(x1,x2) << ">" << endl;
maxx2 = xmid + maxextrap*(x2 - xmid);
quadok = estquadmin(newx2,x1,xmid,x2,y1,ymid,y2);
if ((!quadok) || ((newx2 - x1)*dir<0) || ((newx2 - maxx2)*dir>0)) {
newx2 = xmid + extrapolationfactor*(x2-x1);
}
-
+
newy2 = (*func)(newx2*unitdir + pt);
if ((newx2 - xmid)*(newx2 - x1)<0) { // newx2 is between x1 and xmid
@@ -141,19 +141,19 @@ namespace MISCMATHS {
x2 = newx2; y2 = newy2;
}
}
-
+
}
if ( (y2<ymid) || (y1<ymid) ) {
cerr << "findinitialbound failed to bracket: current triplet is" << endl;
}
}
-
- float optimise1d(ColumnVector &pt, const ColumnVector dir,
- const ColumnVector tol, int &iterations_done,
+
+ float optimise1d(ColumnVector &pt, const ColumnVector dir,
+ const ColumnVector tol, int &iterations_done,
float (*func)(const ColumnVector &), int max_iter,
- float &init_value, float boundguess)
+ float &init_value, float boundguess)
{
// Golden Search Routine
// Must pass in the direction vector in N-space (dir), the initial
@@ -247,10 +247,10 @@ namespace MISCMATHS {
-
- float optimise(ColumnVector &pt, int numopt, const ColumnVector &tol,
- float (*func)(const ColumnVector &), int &iterations_done,
- int max_iter, const ColumnVector& boundguess,
+
+ float optimise(ColumnVector &pt, int numopt, const ColumnVector &tol,
+ float (*func)(const ColumnVector &), int &iterations_done,
+ int max_iter, const ColumnVector& boundguess,
const string type)
{
// Note that numopt can be less than pt.Nrows() - e.g. 6 dof optimisation
@@ -291,12 +291,12 @@ namespace MISCMATHS {
if (avtol < 1.0) break;
// if continuing then change the directions if using Powell's method
- if (type=="powell")
+ if (type=="powell")
{
// find direction of maximal change
int bestm=1;
- for (int m=1; m<=numopt; m++) {
- if (deltaf(m)<deltaf(bestm)) bestm=m;
+ for (int m=1; m<=numopt; m++) {
+ if (deltaf(m)<deltaf(bestm)) bestm=m;
}
fend=fval;
fextrap=(*func)(initpt + 2*(pt-initpt));
@@ -317,7 +317,7 @@ namespace MISCMATHS {
}
}
}
-
+
}
// cout << endl << "Major iterations = " << it << endl;
iterations_done = littot;
@@ -326,11 +326,3 @@ namespace MISCMATHS {
}
-
-
-
-
-
-
-
-
diff --git a/rungekutta.cc b/rungekutta.cc
index 9e7042d4b8f240e4e28175528711442d31832eef..50b92b224e3977042147cb5ec76f17566611ad4a 100644
--- a/rungekutta.cc
+++ b/rungekutta.cc
@@ -14,34 +14,34 @@ using namespace NEWMAT;
namespace MISCMATHS {
void rk(ColumnVector& ret, const ColumnVector& y, const ColumnVector& dy, float x, float h, const Derivative& deriv,const ColumnVector& paramvalues)
-{
- Tracer tr("rk");
+{
+ Tracer tr("rk");
float hh=h*0.5;
float xh=x+hh;
-
+
//first step
ColumnVector yt=y+hh*dy;
-
+
//second step
ColumnVector dyt = deriv.evaluate(xh,yt,paramvalues);
yt=y+hh*dyt;
-
+
//third step
ColumnVector dym = deriv.evaluate(xh,yt,paramvalues);
yt=y+h*dym;
dym=dym+dyt;
-
+
//fourth step
dyt = deriv.evaluate(x+h,yt,paramvalues);
-
+
//addup
ret = y+h*(dy+dyt+2*dym)/6;
}
void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry, float eps, const Derivative& deriv,const ColumnVector& paramvalues)
{
- Tracer tr("rkqc");
+ Tracer tr("rkqc");
float xsav = x;
ColumnVector dysav = dy;
@@ -53,12 +53,12 @@ void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry,
while(true)
{
// take 2 1/2 step sizes
-
+
// first 1/2 step
float hh=h*0.5;
rk(ytemp,ysav,dysav,xsav,hh,deriv,paramvalues);
-
+
// second 1/2 step
x=xsav+hh;
dy = deriv.evaluate(x,ytemp,paramvalues);
@@ -68,32 +68,32 @@ void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry,
// take large step size
rk(ytemp,ysav,dysav,xsav,h,deriv,paramvalues);
-
+
// eval accuracy
float errmax = 0.0;
for(int i=1; i<=y.Nrows(); i++)
{
//errmax=max(abs((y-ytemp)./y));
-
+
float tmp = fabs((y(i)-ytemp(i))/y(i));
if(tmp > errmax) errmax = tmp;
}
errmax=errmax/eps;
-
- if(errmax <=1.0)
+
+ if(errmax <=1.0)
{
// step OK, compute step size for next step
hdid=h;
-
+
if(errmax>6e-4)
hnext=h*std::exp(-0.2*std::log(errmax));
else
hnext=4*h;
-
+
break;
}
- else
+ else
{
// step too large,
h=h*std::exp(-0.25*std::log(errmax));
@@ -105,7 +105,7 @@ void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry,
void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVector& ystart, float x1, float x2, float eps, float hmin, const Derivative& deriv,const ColumnVector& paramvalues)
{
- Tracer tr("runge_kutta");
+ Tracer tr("runge_kutta");
int MAXSTEP=1000;
@@ -124,20 +124,20 @@ void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVec
yp = 0;
int kout=1;
-
+
ColumnVector dy;
for(int k=1; k <= MAXSTEP; k++)
- {
+ {
dy = deriv.evaluate(x,y,paramvalues);
// store results:
xp(kout)=x;
yp.Row(kout)=y;
hp(kout)=h;
-
+
kout=kout+1;
-
+
// stop overshoot of step past x2:
if((x+h-x2)*(x+h-x1)>0) h=x2-x;
@@ -150,7 +150,7 @@ void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVec
yp.Row(kout)=y;
hp(kout)=h;
//kout=kout+1;
-
+
xp = xp.Rows(1,kout);
yp = yp.Rows(1,kout);
@@ -160,8 +160,8 @@ void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVec
{
if(hnext<=hmin) cerr << "step size too small" << endl;
h=hnext;
- }
-
+ }
+
}
cerr << "too many steps" << endl;
}