diff --git a/diff_pvm.cc b/diff_pvm.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f8e1366898bf55fe10e5894a0ed36c4ddb5fc671
--- /dev/null
+++ b/diff_pvm.cc
@@ -0,0 +1,394 @@
+/* Diffusion Partial Volume Model
+
+ Tim Behrens - FMRIB Image Analysis Group
+
+ Copyright (C) 2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <strstream>
+#define WANT_STREAM
+#define WANT_MATH
+// #include "newmatap.h"
+// #include "newmatio.h"
+#include <string>
+#include <math.h>
+#include "utils/log.h"
+#include "diff_pvmoptions.h"
+#include "utils/tracer_plus.h"
+#include "miscmaths/miscprob.h"
+#include "stdlib.h"
+#include "bint/model.h"
+#include "bint/lsmcmcmanager.h"
+#include "bint/lslaplacemanager.h"
+
+using namespace Bint;
+using namespace Utilities;
+using namespace NEWMAT;
+using namespace MISCMATHS;
+
+
+
+const float maxfloat=1e10;
+const float minfloat=1e-10;
+const float maxlogfloat=23;
+const float minlogfloat=-23;
+const int maxint=1000000000;
+
+
+inline float min(float a,float b){
+ return a<b ? a:b;}
+inline float max(float a,float b){
+ return a>b ? a:b;}
+inline Matrix Anis()
+{
+ Matrix A(3,3);
+ A << 1 << 0 << 0
+ << 0 << 0 << 0
+ << 0 << 0 << 0;
+ return A;
+}
+
+inline Matrix Is()
+{
+ Matrix I(3,3);
+ I << 1 << 0 << 0
+ << 0 << 1 << 0
+ << 0 << 0 << 1;
+ return I;
+}
+
+inline ColumnVector Cross(const ColumnVector& A,const ColumnVector& B)
+{
+ ColumnVector res(3);
+ res << A(2)*B(3)-A(3)*B(2)
+ << A(3)*B(1)-A(1)*B(3)
+ << A(1)*B(2)-B(1)*A(2);
+ return res;
+}
+
+inline Matrix Cross(const Matrix& A,const Matrix& B)
+{
+ Matrix res(3,1);
+ res << A(2,1)*B(3,1)-A(3,1)*B(2,1)
+ << A(3,1)*B(1,1)-A(1,1)*B(3,1)
+ << A(1,1)*B(2,1)-B(1,1)*A(2,1);
+ return res;
+}
+
+float mod(float a, float b){
+ while(a>b){a=a-b;}
+ while(a<0){a=a+b;}
+ return a;
+}
+
+
+Matrix form_Amat(const Matrix& r,const Matrix& b)
+{
+ Matrix A(r.Ncols(),7);
+ Matrix tmpvec(3,1), tmpmat;
+
+ for( int i = 1; i <= r.Ncols(); i++){
+ tmpvec << r(1,i) << r(2,i) << r(3,i);
+ tmpmat = tmpvec*tmpvec.t()*b(1,i);
+ A(i,1) = tmpmat(1,1);
+ A(i,2) = 2*tmpmat(1,2);
+ A(i,3) = 2*tmpmat(1,3);
+ A(i,4) = tmpmat(2,2);
+ A(i,5) = 2*tmpmat(2,3);
+ A(i,6) = tmpmat(3,3);
+ A(i,7) = 1;
+ }
+ return A;
+}
+inline SymmetricMatrix vec2tens(ColumnVector& Vec){
+ SymmetricMatrix tens(3);
+ tens(1,1)=Vec(1);
+ tens(2,1)=Vec(2);
+ tens(3,1)=Vec(3);
+ tens(2,2)=Vec(4);
+ tens(3,2)=Vec(5);
+ tens(3,3)=Vec(6);
+ return tens;
+}
+
+inline void cart2sph(const ColumnVector& dir, float& th, float& ph)
+{
+ 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;
+ }
+ else{
+
+ 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=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=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;
+ }
+}
+
+
+
+void cart2sph(const Matrix& dir,ColumnVector& th,ColumnVector& ph)
+{
+ for (int i=1;i<=dir.Ncols();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;
+ }
+ 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)=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)=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;
+
+ }
+ }
+}
+
+
+
+
+
+class Diff_pvmModel : public ForwardModel
+ {
+ public:
+
+ Diff_pvmModel(const Matrix& pbvecs,const Matrix& pbvals,int pdebuglevel)
+ : ForwardModel(pdebuglevel), r(pbvecs) , b(pbvals), alpha(pbvals.Ncols()), beta(pbvals.Ncols()), debuglevel(pdebuglevel)
+
+ {
+ Amat=form_Amat(r,b);
+ cart2sph(r,alpha,beta);
+ }
+
+ ~Diff_pvmModel(){}
+
+ virtual void setparams(float pstd);
+ ReturnMatrix nonlinearfunc(const ColumnVector& paramvalues) const;
+ void initialise(const ColumnVector& S);
+
+
+ protected:
+
+ const Matrix& r;
+ const Matrix& b;
+ ColumnVector alpha;
+ ColumnVector beta;
+ Matrix Amat;
+ int debuglevel;
+};
+
+void Diff_pvmModel::setparams(float pstd)
+ {
+ Tracer_Plus tr("Diff_pvmModel::setdata");
+ if(debuglevel>2){
+ cout << "Diff_pvmModel::setparams"<<endl;
+ }
+ clear_params();
+
+ SinPrior thtmp(1);
+ add_param("th",0.2,thtmp,true);
+ UnifPrior phtmp(0.2,2000*M_PI);
+ add_param("ph",0,phtmp,true);
+ UnifPrior ftmp(0,1);
+ add_param("f",0.5,ftmp,true);
+ GammaPrior dtmp(4,1.0/0.0003); //test this out,
+ add_param("d",0.005,dtmp,true);
+ UnifPrior S0tmp(0,100000);
+ add_param("S0",10000,S0tmp,true);//false);
+
+ }
+
+ReturnMatrix Diff_pvmModel::nonlinearfunc(const ColumnVector& paramvalues) const
+ {
+ Tracer_Plus trace("Diff_pvmModel::nonlinearfunc");
+ if(debuglevel>2){
+ cout << "Diff_pvmModel::nonlinearfunc"<<endl;
+ cout<<paramvalues<<endl;
+ }
+
+ float th=paramvalues(1);
+ float ph=paramvalues(2);
+ float f=paramvalues(3);
+ float D=paramvalues(4);
+ float S0=paramvalues(5);
+ // cout <<" nlf "<<S0<<endl;
+
+ ColumnVector ret(b.Ncols());
+ float angtmp;
+ for (int i = 1; i <= ret.Nrows(); i++){
+ angtmp=cos(ph-beta(i))*sin(alpha(i))*sin(th) + cos(alpha(i))*cos(th);
+ angtmp=angtmp*angtmp;
+ // cout <<angtmp<<endl;
+ // cout <<i<<endl;
+ ret(i)=S0*(f*exp(-D*b(1,i)*angtmp)+(1-f)*exp(-D*b(1,i)));
+ }
+
+
+ if(debuglevel>2){
+ cout<<ret<<endl;
+ cout <<"done"<<endl;
+ }
+ ret.Release();
+ return ret;
+ }
+
+
+
+void Diff_pvmModel::initialise(const ColumnVector& S){
+
+ Tracer_Plus trace("Diff_pvmModel::initialise");
+ if(debuglevel>2){
+ cout << "Diff_pvmModel::initialise"<<endl;
+ }
+
+
+ ColumnVector logS(S.Nrows()),tmp(S.Nrows()),Dvec(7),dir(3);
+ SymmetricMatrix tens; //Basser's Diffusion Tensor;
+ DiagonalMatrix Dd; //eigenvalues
+ Matrix Vd; //eigenvectors
+ float mDd,fsquared;
+ float th,ph,f,D,S0;
+
+ for ( int i = 1; i <= S.Nrows(); i++)
+ {
+ if(S(i)>0){
+ logS(i)=log(S(i));
+ }
+ else{
+ logS(i)=0;
+ }
+ }
+ Dvec = -pinv(Amat)*logS;
+ if( Dvec(7) > -maxlogfloat ){
+ S0=exp(-Dvec(7));
+ }
+ else{
+ S0=S.MaximumAbsoluteValue();
+ }
+
+ for ( int i = 1; i <= S.Nrows(); i++)
+ {
+ if(S0<S.Sum()/S.Nrows()){ S0=S.MaximumAbsoluteValue(); }
+ logS(i)=(S(i)/S0)>0.01 ? log(S(i)):log(0.01*S0);
+ }
+ Dvec = -pinv(Amat)*logS;
+ S0=exp(-Dvec(7));
+ if(S0<S.Sum()/S.Nrows()){ S0=S.Sum()/S.Nrows(); }
+ tens = vec2tens(Dvec);
+ EigenValues(tens,Dd,Vd);
+ mDd = Dd.Sum()/Dd.Nrows();
+ int maxind = Dd(1) > Dd(2) ? 1:2; //finding maximum eigenvalue
+ maxind = Dd(maxind) > Dd(3) ? maxind:3;
+ dir << Vd(1,maxind) << Vd(2,maxind) << Vd(3,maxind);
+ cart2sph(dir,th,ph);
+ th= mod(th,M_PI);
+ ph= mod(ph,2*M_PI);
+ D = Dd(maxind);
+
+ float numer=(1.5*(Dd(1)-mDd)*(Dd(1)-mDd)+(Dd(2)-mDd)*(Dd(2)-mDd)+(Dd(3)-mDd)*(Dd(3)-mDd));
+ float denom=(Dd(1)*Dd(1)+Dd(2)*Dd(2)+Dd(3)*Dd(3));
+ if(denom>0) fsquared=numer/denom;
+ else fsquared=0;
+ if(fsquared>0){f=sqrt(fsquared);}
+ else{f=0;}
+ if(f>=0.95) f=0.95;
+ if(f<=0.001) f=0.001;
+ cout<<"S0 "<<S0<<endl;
+ cout<<"S1 "<<S(1)<<endl;
+ getparam(0).setinitvalue(th);
+ getparam(1).setinitvalue(ph);
+ getparam(2).setinitvalue(f);
+ getparam(3).setinitvalue(D);
+ getparam(4).setinitvalue(S0);
+}
+
+
+int main(int argc, char *argv[])
+{
+ try{
+
+ // Setup logging:
+ Log& logger = LogSingleton::getInstance();
+
+ // parse command line - will output arguments to logfile
+ Diff_pvmOptions& opts = Diff_pvmOptions::getInstance();
+ opts.parse_command_line(argc, argv, logger);
+
+ srand(Diff_pvmOptions::getInstance().seed.value());
+
+ if(opts.debuglevel.value()==1)
+ Tracer_Plus::setrunningstackon();
+
+ if(opts.timingon.value())
+ Tracer_Plus::settimingon();
+
+ // read data
+ VolumeSeries data;
+ data.read(opts.datafile.value());
+ int ntpts = data.tsize();
+ Matrix bvecs = read_ascii_matrix(opts.bvecsfile.value());
+ Matrix bvals = read_ascii_matrix(opts.bvalsfile.value());
+ // mask:
+ Volume mask;
+ mask.read(opts.maskfile.value());
+ mask.threshold(1e-16);
+
+ // threshold using mask:
+ data.setPreThresholdPositions(mask.getPreThresholdPositions());
+ data.thresholdSeries();
+
+ cout << "ntpts=" << ntpts << endl;
+ cout << "nvoxels=" << mask.getVolumeSize() << endl;
+
+ Diff_pvmModel model(bvecs,bvals,Diff_pvmOptions::getInstance().debuglevel.value());
+
+ LSMCMCManager lsmcmc(Diff_pvmOptions::getInstance(),model,data,mask);
+ LSLaplaceManager lslaplace(Diff_pvmOptions::getInstance(),model,data,mask,(Diff_pvmOptions::getInstance().inference.value()=="full"));
+
+
+ if(Diff_pvmOptions::getInstance().inference.value()=="mcmc")
+ {
+ lsmcmc.setup();
+ lsmcmc.run();
+ lsmcmc.save();
+ }
+ else
+ {
+ lslaplace.setup();
+ lslaplace.run();
+ lslaplace.save();
+ }
+
+ if(opts.timingon.value())
+ Tracer_Plus::dump_times(logger.getDir());
+
+ cout << endl << "Log directory was: " << logger.getDir() << endl;
+ }
+ catch(Exception& e)
+ {
+ cerr << endl << e.what() << endl;
+ }
+ catch(X_OptionError& e)
+ {
+ cerr << endl << e.what() << endl;
+ }
+
+ return 0;
+}
diff --git a/diff_pvmoptions.cc b/diff_pvmoptions.cc
new file mode 100644
index 0000000000000000000000000000000000000000..fcef3cd310a5d7ec959f0564ebe15f10fafe8238
--- /dev/null
+++ b/diff_pvmoptions.cc
@@ -0,0 +1,38 @@
+/* meanoptions.cc
+
+ Mark Woolrich - FMRIB Image Analysis Group
+
+ Copyright (C) 2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#define WANT_STREAM
+#define WANT_MATH
+
+#include <iostream.h>
+#include <fstream.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "diff_pvmoptions.h"
+#include "utils/log.h"
+#include "utils/tracer_plus.h"
+
+using namespace Utilities;
+
+namespace Bint {
+
+ Diff_pvmOptions* Diff_pvmOptions::gopt = NULL;
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/diff_pvmoptions.h b/diff_pvmoptions.h
new file mode 100644
index 0000000000000000000000000000000000000000..b4440da681ce2cc2328c08d2b6f66eb187b8debd
--- /dev/null
+++ b/diff_pvmoptions.h
@@ -0,0 +1,80 @@
+/* meanoptions.h
+
+ Mark Woolrich - FMRIB Image Analysis Group
+
+ Copyright (C) 2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#if !defined(Diff_pvmOptions_h)
+#define Diff_pvmOptions_h
+
+#include <string>
+#include <iostream.h>
+#include <fstream.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "utils/options.h"
+#include "utils/log.h"
+#include "bint/bintoptions.h"
+
+using namespace Utilities;
+
+namespace Bint {
+
+class Diff_pvmOptions : public BintOptions {
+ public:
+ static Diff_pvmOptions& getInstance();
+ ~Diff_pvmOptions() { delete gopt; }
+
+ Option<string> bvecsfile;
+ Option<string> bvalsfile;
+
+
+ private:
+ Diff_pvmOptions();
+ const Diff_pvmOptions& operator=(Diff_pvmOptions&);
+ Diff_pvmOptions(Diff_pvmOptions&);
+
+ static Diff_pvmOptions* gopt;
+
+};
+
+ inline Diff_pvmOptions& Diff_pvmOptions::getInstance(){
+ if(gopt == NULL)
+ gopt = new Diff_pvmOptions();
+
+ return *gopt;
+ }
+
+ inline Diff_pvmOptions::Diff_pvmOptions() :
+ BintOptions("diff_pvm", "diff_pvm --verbose\n"),
+ bvecsfile(string("-r,--bvecs"),"bvecs",
+ string("gradient directions"),
+ true, requires_argument),
+ bvalsfile(string("-b,--bvals"),"bvals",
+ string("b values"),
+ true, requires_argument)
+
+ {
+ try {
+ options.add(bvecsfile);
+ options.add(bvalsfile);
+ }
+ catch(X_OptionError& e) {
+ options.usage();
+ cerr << endl << e.what() << endl;
+ }
+ catch(std::exception &e) {
+ cerr << e.what() << endl;
+ }
+
+ }
+}
+
+#endif
+
+
+
+
+
diff --git a/fdt.tcl b/fdt.tcl
index af746a80c2f1abf876e8f4288bd9014fd3bc99bf..4dc464b3a2a4af19e982565db0e1289ac2eedd88 100644
--- a/fdt.tcl
+++ b/fdt.tcl
@@ -816,7 +816,7 @@ proc fdt:dialog { w tclstartupfile } {
collapsible frame $w.advanced -title "Advanced Options"
- set probtrack(nsteps) 1000
+ set probtrack(nsteps) 2000
tixControl $w.advanced.nsteps -label "Maximum number of steps" \
-variable probtrack(nsteps) -step 10 -min 2 \
-selectmode immediate -options { entry.width 6 }
diff --git a/probtrack.cc b/probtrack.cc
index 5bc485cbd44c44496b0313f553e4a4cc75b498d6..7e41d4777408a004cc7b51a98d2b3a9a24ded7d3 100644
--- a/probtrack.cc
+++ b/probtrack.cc
@@ -1765,7 +1765,7 @@ void seeds_to_targets()
path(it+(direc-1)*nsteps/2,1)=round(part.x());
path(it+(direc-1)*nsteps/2,2)=round(part.y());
- path(it+(direc-1)*nsteps/2,3)=round(part.z()); //stopring path in DTI space here
+ path(it+(direc-1)*nsteps/2,3)=round(part.z()); //stopping path in DTI space here
for(unsigned int m=0;m<masknames.size();m++){
if(target_masks[m](x_s,y_s,z_s)>0 && flags[m]==0){
@@ -1892,7 +1892,8 @@ void meshtrack(){
ColumnVector th_ph_f;
Mesh mseeds;
- int ftype=mseeds.load(opts.seedfile.value());
+ int ftype=mseeds.load(opts.meshfile.value());
+ mseeds.load_fs_label(opts.seedfile.value());
ColumnVector mni_origin(3),fs_coord_mm(3),xyz_dti,xyz_seeds,dim_seeds(3);
dim_seeds<<prob.xdim()<<prob.ydim()<<prob.zdim(); //In seedref space if exists. Else in dti space
mni_origin << 92 << 128 << 37;
@@ -2028,7 +2029,7 @@ void meshtrack(){
// save_volume(prob,thisout);
}
} //Close Seed number Loop
- mseeds.save(opts.outfile.value(),ftype);
+ mseeds.save_fs_label(logger.AppendDir(opts.outfile.value()));
}
@@ -2057,6 +2058,8 @@ int main ( int argc, char **argv ){
matrix2();
else if(opts.mode.value()=="maskmatrix")
maskmatrix();
+ else if(opts.mode.value()=="meshtrack")
+ meshtrack();
else{
cout <<"Invalid setting for option mode -- try setting mode=help"<<endl;
}
@@ -2077,3 +2080,4 @@ int main ( int argc, char **argv ){
+
diff --git a/probtrackOptions.h b/probtrackOptions.h
index feb13905c6703ed625806cd2f19dea084a9c225d..658bbe0afd2a00633688db9eaff72550ebab57ca 100644
--- a/probtrackOptions.h
+++ b/probtrackOptions.h
@@ -40,6 +40,7 @@ class probtrackOptions {
Option<string> skipmask;
Option<string> seedref;
Option<string> mask2;
+ Option<string> meshfile;
Option<string> lrmask;
Option<string> logdir;
Option<bool> forcedir;
@@ -113,6 +114,9 @@ class probtrackOptions {
mask2(string("--mask2"), string(""),
string("second mask in both 2 mask modes."),
false, requires_argument),
+ meshfile(string("--mesh"), string(""),
+ string(""),
+ false, requires_argument),
lrmask(string("--lrmask"), string(""),
string("low resolution binary brain mask for stroring connectivity distribution in matrix2 mode."),
false, requires_argument),
@@ -125,7 +129,7 @@ class probtrackOptions {
nparticles(string("-P,--nsamples"), 10000,
string("Number of samples"),
false, requires_argument),
- nsteps(string("-S,--nsteps"), 1000,
+ nsteps(string("-S,--nsteps"), 2000,
string("Number of steps per sample"),
false, requires_argument),
c_thr(string("-c,--cthr"), 0.2,
@@ -160,6 +164,7 @@ class probtrackOptions {
options.add(targetfile);
options.add(skipmask);
options.add(mask2);
+ options.add(meshfile);
options.add(lrmask);
options.add(seedref);
options.add(logdir);
diff --git a/replacevols.cc b/replacevols.cc
index dba7e33f867920e982a57f85b10d9a4efd36b0db..37fed99c70111215fcf550290e4bba7960e8b4ac 100644
--- a/replacevols.cc
+++ b/replacevols.cc
@@ -94,16 +94,16 @@ int main ( int argc, char **argv ){
}
- for(int j=0;j<avgs[0].size();j++){//loop over volume numbers
+ for(unsigned int j=0;j<avgs[0].size();j++){//loop over volume numbers
//Next loop is within volume number over averages just
// Working out which ones to replace and which to keep.
vector<int> repthis,keepthis;
- for(int i=0;i<avgs.size();i++){ //loop over averages
+ for(unsigned int i=0;i<avgs.size();i++){ //loop over averages
bool replaced=false;
- for(int r=0;r<repvols.size();r++){// loop over things to be replaced
+ for(unsigned int r=0;r<repvols.size();r++){// loop over things to be replaced
if(avgs[i][j]==repvols[r]){
replaced=true;
repthis.push_back(avgs[i][j]);
@@ -119,12 +119,12 @@ int main ( int argc, char **argv ){
if(repthis.size()>0){
cerr<<"Replacing volumes: ";
- for(int r=0;r<repthis.size();r++){
+ for(unsigned int r=0;r<repthis.size();r++){
cerr<<repthis[r]<<" ";
}
cerr <<endl;
cerr<<"with the average of volumes: ";
- for(int r=0;r<keepthis.size();r++){
+ for(unsigned int r=0;r<keepthis.size();r++){
cerr<<keepthis[r]<<" ";
}
cerr<<endl;
@@ -135,7 +135,7 @@ int main ( int argc, char **argv ){
volume<float> tmp;
tmp=data4D[keepthis[0] ];
- for(int n=1;n<keepthis.size();n++){
+ for(unsigned int n=1;n<keepthis.size();n++){
tmp=tmp+data4D[keepthis[n] ];
}
tmp=tmp/keepthis.size(); //Average of all non-replaced ones.
@@ -143,7 +143,7 @@ int main ( int argc, char **argv ){
//Next loop replaces all the ones to be replaced with this average
- for(int n=0;n<repthis.size();n++){
+ for(unsigned int n=0;n<repthis.size();n++){
data4D[repthis[n] ]=tmp; //replacing.
}
diff --git a/sausages.cc b/sausages.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e5bc456238c1d7c93530ab23f49c0b72b4bf990c
--- /dev/null
+++ b/sausages.cc
@@ -0,0 +1,73 @@
+#include <iostream>
+#include <fstream>
+#include "newimage/newimageall.h"
+#include <vector>
+using namespace std;
+using namespace NEWIMAGE;
+using namespace NEWMAT;
+
+
+
+
+
+
+
+int main ( int argc, char **argv ){
+ if(argc<6){
+ cerr<<"usage: sausages <refvol> <coordfile> <output> <start1> <end1> ... <startn> <endn>"<<endl;
+ exit(0);
+ }
+
+
+
+ volume<int> out;
+ read_volume(out,argv[1]);
+ out*=0;
+
+ volume<int> coord;
+ read_volume(coord,argv[2]);
+ vector<pair<int,int> > chunks;
+ if(argc/2 != argc/2.0f ){
+ cerr<<"your chunks are not in pairs"<<endl;
+ return -1;
+ }
+ pair<int,int> tmppair;
+ for(int i=4;i<argc;i+=2){
+ int tmp=atoi(argv[i]);
+ if(tmp>coord.xsize()){
+ cerr<<"There aren't "<<argv[i]<<" elements in "<<argv[2]<<endl;
+ return -1;
+ }
+ else{
+ tmppair.first=tmp;
+ }
+ tmp=atoi(argv[i+1]);
+ if(tmp>coord.xsize()){
+ cerr<<"There aren't "<<argv[i+1]<<" elements in "<<argv[2]<<endl;
+ return -1;
+ }
+ else{
+ tmppair.second=tmp;
+ }
+ chunks.push_back(tmppair);
+ }
+
+ for(int chunkno=0;chunkno < chunks.size();chunkno++){
+ for(int i=chunks[chunkno].first;i<=chunks[chunkno].second;i++){
+ out(coord(i,0,0),coord(i,1,0),coord(i,2,0))=chunkno+1;
+ }
+ }
+
+ save_volume_dtype(out,argv[3],dtype(string(argv[1])));
+ return 0;
+}
+
+
+
+
+
+
+
+
+
+