pt_matrix.cc

/*  Copyright (C) 2004 University of Oxford  */

/*  CCOPYRIGHT  */

#include "pt_matrix.h"

using namespace std;
using namespace NEWIMAGE;
using namespace TRACT;
using namespace Utilities;
using namespace PARTICLE;
using namespace TRACTVOLS;
using namespace mesh;

void matrix2(){
  
  probtrackOptions& opts =probtrackOptions::getInstance();
  Log& logger = LogSingleton::getInstance();

  ////////////////////////////
  //  Log& logger = LogSingleton::getInstance();
  //  logger.makeDir(opts.logdir.value(),"logfile",true,false);
  ////////////////////////////////////
  float xst,yst,zst,x,y,z;
  int nparticles=opts.nparticles.value();
  int nsteps=opts.nsteps.value();
  ///////////////////////////////////
  volume<int> mask;
  read_volume(mask,opts.maskfile.value());  
  volume<int> skipmask;
  if(opts.skipmask.value()!="") read_volume(skipmask,opts.skipmask.value());
  float lcrat=5;
  volume4D<float> loopcheck;
  if(opts.loopcheck.value()){
    loopcheck.reinitialize(int(ceil(mask.xsize()/lcrat)+1),int(ceil(mask.ysize()/lcrat)+1),int(ceil(mask.zsize()/lcrat)+1),3);
    loopcheck=0;
  }
  //////////////////////////////////////////////
  // Segmented Volumes                        //
  //////////////////////////////////////////////
 //   vector<string> masknames;
//    read_masks(masknames,opts.cortexfile.value());
//    vector<volume<int> > cortex_masks;
//    vector<volume<int> > thal_segs;
//    volume<int> tmpcort;
//    cerr<<"Number of masks "<<masknames.size()<<endl;
//    for( unsigned int m = 0; m < masknames.size(); m++ ){
//      cerr<<"Reading "<<masknames[m]<<endl;
//      read_volume(tmpcort,masknames[m]);
//      cortex_masks.push_back(tmpcort);
//      tmpcort=0;
//      thal_segs.push_back(tmpcort);
//  }


  volume<int> Seeds;
  read_volume(Seeds,opts.seedfile.value());

  volume<int> RUBBISH;
  if(opts.rubbishfile.value()!=""){
    read_volume(RUBBISH,opts.rubbishfile.value());
  }

  volume<int> prob;
  volume<int> beenhere;
  prob=Seeds;prob=0;

  int numseeds=0;
  for(int Wz=Seeds.minz();Wz<=Seeds.maxz();Wz++){
    for(int Wy=Seeds.miny();Wy<=Seeds.maxy();Wy++){