From 82adeaee3b57730e036b5f56307156c36f1ccbeb Mon Sep 17 00:00:00 2001
From: Mark Jenkinson <mark@fmrib.ox.ac.uk>
Date: Mon, 10 Apr 2006 16:48:26 +0000
Subject: [PATCH] Lots of changes - backed out VS6 mods and Saad changes and MJ
pre-branch changes but added MJ/DF branch changes
---
cspline.cc | 21 +++++----
cspline.h | 1 -
f2z.cc | 6 +--
kernel.cc | 6 +--
minimize.cc | 16 +++----
miscmaths.cc | 125 ++++++++++++++++++++++++++-------------------------
6 files changed, 87 insertions(+), 88 deletions(-)
diff --git a/cspline.cc b/cspline.cc
index ea2cb6b..80ed10f 100644
--- a/cspline.cc
+++ b/cspline.cc
@@ -11,7 +11,6 @@
#include <string>
#include <iostream>
#include <fstream>
-//#include <unistd.h>
#include "newmatap.h"
#include "newmatio.h"
@@ -72,8 +71,8 @@ namespace MISCMATHS{
}
ColumnVector b(n);
b=0;
- for(int ii=2;ii<=b.Nrows()-1;ii++){
- b(ii)=3*(dx(ii)*dydx(ii-1)+dx(ii-1)*dydx(ii));
+ for(int i=2;i<=b.Nrows()-1;i++){
+ b(i)=3*(dx(i)*dydx(i-1)+dx(i-1)*dydx(i));
}
float x31=nodes(3)-nodes(1),xn=nodes(n)-nodes(n-2);
@@ -94,17 +93,17 @@ namespace MISCMATHS{
ColumnVector d(n-1),c(n-1);
- for(int j2=1;j2<n;j2++){
- d(j2)=(s(j2)+s(j2+1)-2*dydx(j2))/dx(j2);
- c(j2)=(dydx(j2)-s(j2))/dx(j2)-d(j2);
+ for(int j=1;j<n;j++){
+ d(j)=(s(j)+s(j+1)-2*dydx(j))/dx(j);
+ c(j)=(dydx(j)-s(j))/dx(j)-d(j);
}
coefs.ReSize(n-1,4);
- for(int j3=1;j3<n;j3++){
- coefs(j3,1)=vals(j3);
- coefs(j3,2)=s(j3);
- coefs(j3,3)=c(j3);
- coefs(j3,4)=d(j3)/dx(j3);
+ for(int j=1;j<n;j++){
+ coefs(j,1)=vals(j);
+ coefs(j,2)=s(j);
+ coefs(j,3)=c(j);
+ coefs(j,4)=d(j)/dx(j);
}
}
diff --git a/cspline.h b/cspline.h
index 840059c..59cbe8c 100644
--- a/cspline.h
+++ b/cspline.h
@@ -13,7 +13,6 @@
#include <string>
#include <iostream>
#include <fstream>
-//#include <unistd.h>
#include "newmatap.h"
#include "newmatio.h"
diff --git a/f2z.cc b/f2z.cc
index 4c91f5f..b21085a 100644
--- a/f2z.cc
+++ b/f2z.cc
@@ -63,14 +63,14 @@ namespace MISCMATHS {
for(int i = 1; i <= N; i++)
{
// cerr << "i=" << i;
- iter = iter + top*MISCMATHS::pow(f,(-(n+i-1)))/(MISCMATHS::pow(alpha,i)*bot);
+ iter = iter + top* ( std::pow( f,float(-(n+i-1)) ) / ( std::pow(alpha,double(i))*bot ) );
top = top*(n-1+i)*(-1);
bot = bot*(n+m-1+i);
// cerr << "iter=" << iter;
}
- if(iter <= 0) throw RBD_COMMON::Exception("iter negative");
+ if(iter <= 0) throw Exception("iter negative");
float logp = loggam-(m-1)*(::log(1+alpha*f))+::log(iter);
@@ -89,7 +89,7 @@ namespace MISCMATHS {
{
logp=largef2logp(f,d1,d2);
}
- catch(RBD_COMMON::Exception& p_excp)
+ catch(Exception& p_excp)
{
cerr << "Negative iter in F2z::largef2logp" << endl;
return false;
diff --git a/kernel.cc b/kernel.cc
index 48335c5..16c2675 100644
--- a/kernel.cc
+++ b/kernel.cc
@@ -6,8 +6,8 @@
/* CCOPYRIGHT */
-#include "miscmaths.h"
#include "kernel.h"
+#include "miscmaths.h"
namespace MISCMATHS {
@@ -77,8 +77,8 @@ namespace MISCMATHS {
int hw = (w-1)/2; // convert to half-width
// set x between +/- width
float halfnk = (nstore-1.0)/2.0;
- for (int nx=1; nx<=nstore; nx++) {
- float x=(nx-halfnk-1)/halfnk*hw;
+ for (int n=1; n<=nstore; n++) {
+ float x=(n-halfnk-1)/halfnk*hw;
if ( (sincwindowtype=="hanning") || (sincwindowtype=="h") ) {
ker(n) = sincfn(x)*hanning(x,hw);
} else if ( (sincwindowtype=="blackman") || (sincwindowtype=="b") ) {
diff --git a/minimize.cc b/minimize.cc
index 7ef1d0c..e1f3e0f 100644
--- a/minimize.cc
+++ b/minimize.cc
@@ -9,7 +9,7 @@
#include <string>
#include <iostream>
#include <fstream>
-//#include <unistd.h>
+#include <unistd.h>
#include <vector>
#include <algorithm>
#include "newmatap.h"
@@ -223,9 +223,9 @@ void minsearch(ColumnVector& x, const EvalFunction& func, ColumnVector& paramsto
ColumnVector onesn(ntot);
onesn=1;
ColumnVector one2n(ntot),two2np1(ntot);
- for(int ii=1;ii<=ntot;ii++){
- one2n(ii)=ii;
- two2np1(ii)=ii+1;
+ for(int i=1;i<=ntot;i++){
+ one2n(i)=i;
+ two2np1(i)=i+1;
}
// We want to store the best n+1 parameter estimates
@@ -244,14 +244,14 @@ void minsearch(ColumnVector& x, const EvalFunction& func, ColumnVector& paramsto
//perturb each parameter by a bit, and store the cost.
ColumnVector y=x;
- for(int iii=1;iii<=ntot;iii++){
+ for(int i=1;i<=ntot;i++){
// The values of nonvarying parameters should be the same in
// all of the optional param vectors and therefore in all
// combinations of them in the remainder of the code.
- if(paramstovary(iii)==1){
- if(y(iii)!=0){y(iii)=(1+usual_delta)*y(iii);}
- else{y(iii)=(1+zero_term_delta);}
+ if(paramstovary(i)==1){
+ if(y(i)!=0){y(i)=(1+usual_delta)*y(i);}
+ else{y(i)=(1+zero_term_delta);}
en=func.evaluate(y);
func_evals++;
diff --git a/miscmaths.cc b/miscmaths.cc
index c38a84f..92d5804 100644
--- a/miscmaths.cc
+++ b/miscmaths.cc
@@ -8,15 +8,12 @@
// Miscellaneous maths functions
-#define NOMINMAX
-
-#include <cstdlib>
-#include <cmath>
#include "miscmaths.h"
#include "miscprob.h"
#include "stdlib.h"
#include "newmatio.h"
+using namespace std;
namespace MISCMATHS {
@@ -24,14 +21,13 @@ namespace MISCMATHS {
// (version egcs-2.91.57)
// A temporary fix of including the std:: in front of all abs() etc
// has been done for now
-/*
using std::abs;
using std::sqrt;
using std::exp;
using std::log;
// using std::pow;
using std::atan2;
-*/
+
string size(const Matrix& mat)
{
@@ -939,7 +935,7 @@ namespace MISCMATHS {
for (int i=1; i<=3; i++) { params(i+3) = transl(i); }
ColumnVector rotparams(3);
(*rotmat2params)(rotparams,rotmat);
- for (int ii=1; ii<=3; ii++) { params(ii) = rotparams(ii); }
+ for (int i=1; i<=3; i++) { params(i) = rotparams(i); }
return 0;
}
@@ -1031,6 +1027,29 @@ float rms_deviation(const Matrix& affmat1, const Matrix& affmat2,
// helper function - calls nifti, but with FSL default case
+
+Matrix Mat44ToNewmat(mat44 m)
+{
+ Matrix r(4,4);
+
+ for(unsigned short i = 0; i < 4; ++i)
+ for(unsigned short j = 0; j < 4; ++j)
+ r(i+1, j+1) = m.m[i][j];
+
+ return r;
+}
+
+mat44 NewmatToMat44(const Matrix& m)
+{
+ mat44 r;
+
+ for(unsigned short i = 0; i < 4; ++i)
+ for(unsigned short j = 0; j < 4; ++j)
+ r.m[i][j] = m(i+1, j+1);
+
+ return r;
+}
+
void get_axis_orientations(const Matrix& sform_mat, int sform_code,
const Matrix& qform_mat, int qform_code,
int& icode, int& jcode, int& kcode)
@@ -1094,7 +1113,7 @@ float median(const ColumnVector& x)
void cart2sph(const ColumnVector& dir, float& th, float& ph)
{
- float mag=std::sqrt(dir(1)*dir(1)+dir(2)*dir(2)+dir(3)*dir(3));
+ float mag=sqrt(dir(1)*dir(1)+dir(2)*dir(2)+dir(3)*dir(3));
if(mag==0){
ph=M_PI/2;
th=M_PI/2;
@@ -1103,13 +1122,13 @@ void cart2sph(const ColumnVector& dir, float& th, float& ph)
if(dir(1)==0 && dir(2)>=0) ph=M_PI/2;
else if(dir(1)==0 && dir(2)<0) ph=-M_PI/2;
- else if(dir(1)>0) ph=std::atan(dir(2)/dir(1));
- else if(dir(2)>0) ph=std::atan(dir(2)/dir(1))+M_PI;
- else ph=std::atan(dir(2)/dir(1))-M_PI;
+ else if(dir(1)>0) ph=atan(dir(2)/dir(1));
+ else if(dir(2)>0) ph=atan(dir(2)/dir(1))+M_PI;
+ else ph=atan(dir(2)/dir(1))-M_PI;
if(dir(3)==0) th=M_PI/2;
- else if(dir(3)>0) th=std::atan(std::sqrt(dir(1)*dir(1)+dir(2)*dir(2))/dir(3));
- else th=std::atan(std::sqrt(dir(1)*dir(1)+dir(2)*dir(2))/dir(3))+M_PI;
+ else if(dir(3)>0) th=atan(sqrt(dir(1)*dir(1)+dir(2)*dir(2))/dir(3));
+ else th=atan(sqrt(dir(1)*dir(1)+dir(2)*dir(2))/dir(3))+M_PI;
}
}
@@ -1126,7 +1145,7 @@ void cart2sph(const Matrix& dir,ColumnVector& th,ColumnVector& ph)
}
for (int i=1;i<=dir.Ncols();i++) {
- float mag=std::sqrt(dir(1,i)*dir(1,i)+dir(2,i)*dir(2,i)+dir(3,i)*dir(3,i));
+ float mag=sqrt(dir(1,i)*dir(1,i)+dir(2,i)*dir(2,i)+dir(3,i)*dir(3,i));
if(mag==0){
ph(i)=M_PI/2;
th(i)=M_PI/2;
@@ -1134,13 +1153,13 @@ void cart2sph(const Matrix& dir,ColumnVector& th,ColumnVector& ph)
else{
if(dir(1,i)==0 && dir(2,i)>=0) ph(i)=M_PI/2;
else if(dir(1,i)==0 && dir(2,i)<0) ph(i)=-M_PI/2;
- else if(dir(1,i)>0) ph(i)=std::atan(dir(2,i)/dir(1,i));
- else if(dir(2,i)>0) ph(i)=std::atan(dir(2,i)/dir(1,i))+M_PI;
- else ph(i)=std::atan(dir(2,i)/dir(1,i))-M_PI;
+ else if(dir(1,i)>0) ph(i)=atan(dir(2,i)/dir(1,i));
+ else if(dir(2,i)>0) ph(i)=atan(dir(2,i)/dir(1,i))+M_PI;
+ else ph(i)=atan(dir(2,i)/dir(1,i))-M_PI;
if(dir(3,i)==0) th(i)=M_PI/2;
- else if(dir(3,i)>0) th(i)=std::atan(std::sqrt(dir(1,i)*dir(1,i)+dir(2,i)*dir(2,i))/dir(3,i));
- else th(i)=std::atan(std::sqrt(dir(1,i)*dir(1,i)+dir(2,i)*dir(2,i))/dir(3,i))+M_PI;
+ else if(dir(3,i)>0) th(i)=atan(sqrt(dir(1,i)*dir(1,i)+dir(2,i)*dir(2,i))/dir(3,i));
+ else th(i)=atan(sqrt(dir(1,i)*dir(1,i)+dir(2,i)*dir(2,i))/dir(3,i))+M_PI;
}
}
@@ -1172,7 +1191,7 @@ ReturnMatrix repmat(const Matrix &mat, const int rows, const int cols)
Matrix res = mat;
for(int ctr = 1; ctr < cols; ctr++){res |= mat;}
Matrix tmpres = res;
- for(int ctr1 = 1; ctr1 < rows; ctr1++){res &= tmpres;}
+ for(int ctr = 1; ctr < rows; ctr++){res &= tmpres;}
res.Release();
return res;
}
@@ -1422,8 +1441,8 @@ ReturnMatrix stdev(const Matrix& mat, const int dim)
ReturnMatrix gt(const Matrix& mat1,const Matrix& mat2)
{
- int ctrcol = Xmin(mat1.Ncols(),mat2.Ncols());
- int ctrrow = Xmin(mat1.Nrows(),mat2.Nrows());
+ int ctrcol = std::min(mat1.Ncols(),mat2.Ncols());
+ int ctrrow = std::min(mat1.Nrows(),mat2.Nrows());
Matrix res(ctrrow,ctrcol);
res=0.0;
@@ -1442,8 +1461,8 @@ ReturnMatrix gt(const Matrix& mat1,const Matrix& mat2)
ReturnMatrix lt(const Matrix& mat1,const Matrix& mat2)
{
- int ctrcol = Xmin(mat1.Ncols(),mat2.Ncols());
- int ctrrow = Xmin(mat1.Nrows(),mat2.Nrows());
+ int ctrcol = std::min(mat1.Ncols(),mat2.Ncols());
+ int ctrrow = std::min(mat1.Nrows(),mat2.Nrows());
Matrix res(ctrrow,ctrcol);
res=0.0;
@@ -1462,8 +1481,8 @@ ReturnMatrix lt(const Matrix& mat1,const Matrix& mat2)
ReturnMatrix geqt(const Matrix& mat1,const Matrix& mat2)
{
- int ctrcol = Xmin(mat1.Ncols(),mat2.Ncols());
- int ctrrow = Xmin(mat1.Nrows(),mat2.Nrows());
+ int ctrcol = std::min(mat1.Ncols(),mat2.Ncols());
+ int ctrrow = std::min(mat1.Nrows(),mat2.Nrows());
Matrix res(ctrrow,ctrcol);
res=0.0;
@@ -1478,30 +1497,12 @@ ReturnMatrix geqt(const Matrix& mat1,const Matrix& mat2)
res.Release();
return res;
}
-ReturnMatrix geqt(const Matrix& mat,const float a)
-{
- int ncols = mat.Ncols();
- int nrows = mat.Nrows();
- Matrix res(nrows,ncols);
- res=0.0;
-
- for (int ctr1 = 1; ctr1 <= nrows; ctr1++) {
- for (int ctr2 =1; ctr2 <= ncols; ctr2++) {
- if( mat(ctr1,ctr2) >= a){
- res(ctr1,ctr2) = 1.0;
- }
- }
- }
-
- res.Release();
- return res;
-}
ReturnMatrix leqt(const Matrix& mat1,const Matrix& mat2)
{
- int ctrcol = Xmin(mat1.Ncols(),mat2.Ncols());
- int ctrrow = Xmin(mat1.Nrows(),mat2.Nrows());
+ int ctrcol = std::min(mat1.Ncols(),mat2.Ncols());
+ int ctrrow = std::min(mat1.Nrows(),mat2.Nrows());
Matrix res(ctrrow,ctrcol);
res=0.0;
@@ -1520,8 +1521,8 @@ ReturnMatrix leqt(const Matrix& mat1,const Matrix& mat2)
ReturnMatrix eq(const Matrix& mat1,const Matrix& mat2)
{
- int ctrcol = Xmin(mat1.Ncols(),mat2.Ncols());
- int ctrrow = Xmin(mat1.Nrows(),mat2.Nrows());
+ int ctrcol = std::min(mat1.Ncols(),mat2.Ncols());
+ int ctrrow = std::min(mat1.Nrows(),mat2.Nrows());
Matrix res(ctrrow,ctrcol);
res=0.0;
@@ -1540,8 +1541,8 @@ ReturnMatrix eq(const Matrix& mat1,const Matrix& mat2)
ReturnMatrix neq(const Matrix& mat1,const Matrix& mat2)
{
- int ctrcol = Xmin(mat1.Ncols(),mat2.Ncols());
- int ctrrow = Xmin(mat1.Nrows(),mat2.Nrows());
+ int ctrcol = std::min(mat1.Ncols(),mat2.Ncols());
+ int ctrrow = std::min(mat1.Nrows(),mat2.Nrows());
Matrix res(ctrrow,ctrcol);
res=0.0;
@@ -1720,7 +1721,7 @@ void element_mod_n(Matrix& Mat,double n)
int nextpow2(int n)
{
- return (int)pow(2,ceil(std::log(float(n))/std::log(float(2))));
+ return (int)pow(2,ceil(log(float(n))/log(float(2))));
}
void xcorr(const Matrix& p_ts, Matrix& ret, int lag, int p_zeropad)
@@ -1763,10 +1764,10 @@ void xcorr(const Matrix& p_ts, Matrix& ret, int lag, int p_zeropad)
float varx = var(x.Rows(1,sizeTS)).AsScalar();
ret.Column(i) = realifft.Rows(1,lag);
- for(int jj = 1; jj <= lag-1; jj++)
+ for(int j = 1; j <= lag-1; j++)
{
// Correction to make autocorr unbiased and normalised
- ret(jj,i) = ret(jj,i)/((sizeTS-jj)*varx);
+ ret(j,i) = ret(j,i)/((sizeTS-j)*varx);
}
}
}
@@ -1794,15 +1795,15 @@ void detrend(Matrix& p_ts, int p_level)
for(int t = 1; t <= sizeTS; t++)
{
for(int l = 0; l <= p_level; l++)
- a(t,l+1) = std::pow((float)t/sizeTS,l);
+ a(t,l+1) = pow((float)t/sizeTS,l);
}
// Form residual forming matrix R:
Matrix R = Identity(sizeTS)-a*pinv(a);
- for(int tt = 1; tt <= sizeTS; tt++)
+ for(int t = 1; t <= sizeTS; t++)
{
- p_ts.Column(tt) = R*p_ts.Column(tt);
+ p_ts.Column(t) = R*p_ts.Column(t);
}
}
@@ -2038,7 +2039,7 @@ float csevl(const float x, const ColumnVector& cs, const int n)
else
{
const float aux = csevl(8/(Sqr(y))-1, apsics, ntapsi);
- psi = std::log(x) - 0.5/x + aux;
+ psi = log(x) - 0.5/x + aux;
}
return psi;
@@ -2118,14 +2119,14 @@ float csevl(const float x, const ColumnVector& cs, const int n)
cout<<i<<",";
////////////////////
// update phim
- for(int li=1; li <= nevs; li++)
+ for(int l=1; l <= nevs; l++)
{
float b_m0 = 1e10;
float c_m0 = 2;
float c_m = 1.0/2.0 + c_m0;
- float b_m = 1.0/(0.5*(Sqr(B(li))+lambdaB(li))+1.0/b_m0);
- gam_m(li) = b_m*c_m;
+ float b_m = 1.0/(0.5*(Sqr(B(l))+lambdaB(l))+1.0/b_m0);
+ gam_m(l) = b_m*c_m;
}
// OUT(gam_m(1));
@@ -2137,8 +2138,8 @@ float csevl(const float x, const ColumnVector& cs, const int n)
SymmetricMatrix lambda_B(nevs);
lambda_B = 0;
- for(int lj=1; lj <= nevs; lj++)
- lambda_B(lj,lj)=gam_m(lj);
+ for(int l=1; l <= nevs; l++)
+ lambda_B(l,l)=gam_m(l);
SymmetricMatrix tmp = lambda_B + gam_y*ZZ;
lambda_B << tmp;
--
GitLab