C++ Cube::SetVolume方法代码示例

int
tes1_read_vol(Tes &ts,Cube &cb,int t)
{
  gzFile fp;
  string keyword;
  tokenlist args;
  int cnt;
  
  // Tes::ReadHeader should already have read the header, which
  // includes the mask
  if (!ts.header_valid)
    return (100);
  if (t<0 || t>ts.dimt-1)
    return 101;

  fp=gzopen(ts.GetFileName().c_str(),"r");
  if (!fp)
    return (100);

  // skip the header and mask and advance to our image position
  gzseek(fp,ts.offset+(t*ts.datasize),SEEK_SET);
  
  cb.SetVolume(ts.dimx,ts.dimy,ts.dimz,ts.datatype);
  if (!cb.data)
    return 102;
  int index=0;
  for (int k=0; k<ts.dimz; k++) {
    for (int j=0; j<ts.dimy; j++) {
      for (int i=0; i<ts.dimx; i++) {
	if (ts.mask[index]) {
	  cnt=gzread(fp,cb.data+(ts.datasize*index),ts.datasize);
	  if (cnt!=ts.datasize) {
	    gzclose(fp);
	    return 103;
	  }
	  gzseek(fp,ts.datasize*(ts.dimt-1),SEEK_CUR);
	}
	index++;
      }
    }
  }
  gzclose(fp);
  if (my_endian() != ts.filebyteorder)
    cb.byteswap();
  if (ts.f_scaled) {
    if (ts.datatype==vb_byte || ts.datatype==vb_short || ts.datatype==vb_long)
      cb.convert_type(vb_float);
    cb*=ts.scl_slope;
    cb+=ts.scl_inter;
  }
  return(0);  // no error!
}

void
Extractor::init(tokenlist &args)
{
  // defaults
  xvflag=0;
  newdim[0]=41;
  newdim[1]=51;
  newdim[2]=27;

  newcube.voxsize[0]=3.75;
  newcube.voxsize[1]=3.75;
  newcube.voxsize[2]=5.00;

  newcube.origin[0]=21;
  newcube.origin[1]=31;
  newcube.origin[2]=10;

  original.ReadFile(args[0]);

  if (args.size()==1)
    newcube.SetFileName((string)"m"+original.GetFileName());

  if (args.size()==2 || args.size()==3) {
    newcube.SetFileName(args[1]);
    if (args.size()==3) {
      xval=strtol(args[2]);
      xvflag=1;
    }
  }

  if (args.size()==11 || args.size()==12) {
    newcube.SetFileName(args[1]);
    newdim[0]=strtol(args[2]);
    newdim[1]=strtol(args[3]);
    newdim[2]=strtol(args[4]);
    
    newcube.voxsize[0]=strtod(args[5]);
    newcube.voxsize[1]=strtod(args[6]);
    newcube.voxsize[2]=strtod(args[7]);
    
    newcube.origin[0]=strtol(args[8]);
    newcube.origin[1]=strtol(args[9]);
    newcube.origin[2]=strtol(args[10]);
    if (args.size()==12) {
      xval=strtol(args[11]);
      xvflag=1;
    }
  }
  newcube.SetVolume(newdim[0],newdim[1],newdim[2],original.datatype);
}

void
do_sumrfx(const string &fname,const string &mname,double thresh,vector<VBRegion> &rlist)
{
  int i,j,k;
  Tes rfx;
  Cube mask;

  rfx.ReadFile(fname);
  if (!rfx.data_valid) {
    printf("sumrfx: invalid 4D file %s.\n",fname.c_str());
    return;
  }
  // specified mask, intersect with tes mask
  if (mname.size()) {
    mask.ReadFile(mname);
    for (i=0; i<mask.dimx; i++)
      for (j=0; j<mask.dimy; j++)
	for (k=0; k<mask.dimz; k++)
	  if (rfx.GetMaskValue(i,j,k)<0.5)
	    mask.SetValue(i,j,k,0);
  }
  // default mask
  else {
    mask.SetVolume(rfx.dimx,rfx.dimy,rfx.dimz,vb_byte);
    for (i=0; i<mask.dimx; i++)
      for (j=0; j<mask.dimy; j++)
	for (k=0; k<mask.dimz; k++)
	  mask.SetValue(i,j,k,rfx.GetMaskValue(i,j,k));
  }
  if (!mask.data_valid) {
    printf("sumrfx: invalid mask file %s.\n",mname.c_str());
    return;
  }

  if (mask.dimx != rfx.dimx || mask.dimy != rfx.dimy || mask.dimz != rfx.dimz) {
    printf("sumrfx: mask and rfx volume dimensions don't match\n");
    return;
  }

  printf("sumrfx: thresh %.2f mask file: %s\n",thresh,mname.c_str());

  for (i=0; i<rfx.dimt; i++) {
    Cube mycube;
    rfx.getCube(i,mycube);
    if (!mycube.data_valid) continue;
    do_sumcube(rfx.GetFileName(),i,mycube,mask,thresh,rlist);
  }
}

//.........这里部分代码省略.........
  if (perm_index>-1) {
    VBMatrix vm(perm_mat,0,0,perm_index,perm_index);
    perm_order=vm.GetColumn(0);
    VB_Vector tmp(depvar.size());
    for (uint32 i=0; i<depvar.size(); i++)
      tmp[i]=depvar[(int)perm_order[i]];
    depvar=tmp;
  }
  // convert dv to bitmask
  bitmask dvbm;
  dvbm.resize(depvar.size());
  dvbm.clear();
  for (size_t i=0; i<depvar.size(); i++) {
    if (!f_flip && fabs(depvar[i])>FLT_MIN)
      dvbm.set(i);
    if (f_flip && !(fabs(depvar[i])>FLT_MIN))
      dvbm.set(i);
  }
  
  string partstring;
  if (nparts>1)
    partstring="_part_"+strnum(part);

  bitmask bm;
  bm.resize(ts.dimt);
  map<bitmask,x2val> statlookup;
  map<bitmask,x2val>::iterator iter;
  Cube statmap(ts.dimx,ts.dimy,ts.dimz,vb_float);
  Cube pmap;
  Tes cimap;
  x2val res;
  Cube fdrmask;
  if (f_fdr || pmapname.size()) {
    pmap.SetVolume(ts.dimx,ts.dimy,ts.dimz,vb_float);
    fdrmask=mask;
    fdrmask.zero();
  }
  vector<double> pvals;   // all p vals used in fdr calculation
  bool f_non1=0;
  int16 val;
  for (int i=0; i<ts.dimx; i++) {
    for (int j=0; j<ts.dimy; j++) {
      for (int k=0; k<ts.dimz; k++) {
        if (!mask.testValue(i,j,k))
          continue;
        for (int m=0; m<ts.dimt; m++) {
          val=ts.getValue<int16>(i,j,k,m);
          if (val) {
            bm.set(m);
            if (val!=1) f_non1=1;
          }
          else
            bm.unset(m);
        }
        if (bm.count()<minlesions)
          continue;
        if (ts.dimt-bm.count()<2)
          continue;
        iter=statlookup.find(bm);
        if (iter==statlookup.end()) {
          // this is a new pattern
          if (f_fdr)
            fdrmask.SetValue(i,j,k,1);
          if (f_fisher)
            res=calc_fisher(bm,dvbm);
          else

int
VBSim::Go(int argc,char *argv[])
{
  tokenlist args;
  rng=NULL;
  dimx=0;
  dimy=0;
  dimz=0;
  dimt=0;
  n_mean=10.0;
  n_variance=5.0;
  n_fwhm=0.0;
  float vx=1.0;
  float vy=1.0;
  float vz=1.0;
  float vt=2000;
  uint32 rngseed=VBRandom();

  args.Transfer(argc,argv);

  if (args.size() == 0) {
    vbsim_help();
    exit(0);
  }

  for (size_t i=0; i<args.size(); i++) {
    // -d x y z for dims
    if (args[i]=="-d" && i<args.size()-4) {
      dimx=strtol(args[i+1]);
      dimy=strtol(args[i+2]);
      dimz=strtol(args[i+3]);
      dimt=strtol(args[i+4]);
      i+=4;
    }
    // -c for loading anatomy
    else if (args[i]=="-c" && i<args.size()-1) {
      anatname=args[i+1];
      i++;
    }
    // -z for voxel sizes
    else if (args[i]=="-z" && i<args.size()-3) {
      vx=strtod(args[++i]);
      vy=strtod(args[++i]);
      vz=strtod(args[++i]);
      vt=strtod(args[++i]);
    }
    // -n for per-volume noise
    else if (args[i]=="-n" && i<args.size()-3) {
      n_mean=strtod(args[i+1]);
      n_variance=strtod(args[i+2]);
      n_fwhm=strtod(args[i+3]);
      i+=3;
    }
    else if (args[i]=="-s" && i<args.size()-1)
      rngseed=strtol(args[++i]);
    else if (args[i]=="-o" && i<args.size()-1) {
      outfile=args[i+1];
      i++;
    }
    else {
      printf("[E] vbsim: unrecognized argument %s\n",args(i));
      return 140;
    }
  }
  if (dimx<1 || dimy<0 || dimz<0 || dimt<0) {
    printf("[E] vbsim: bad dimensions\n");
    return 110;
  }

  // initialize RNG
  rng=gsl_rng_alloc(gsl_rng_mt19937);
  assert(rng);
  gsl_rng_set(rng,rngseed);

  // FIXME tell the user here what we're doing

  // CREATE ANATOMY (constant image)
  if (anatname.size()) {
    if (anat.ReadFile(anatname)) {
      printf("[E] vbsim: couldn't read %s\n",anatname.c_str());
      return 101;
    }
    if (anat.dimx!=dimx || anat.dimy!=dimy || anat.dimz != dimz) {
      printf("[E] vbsim: %s doesn't match your dimensions\n",anatname.c_str());
      return 102;
    }
  }
  else {
    // create volume according to dimensions, add random noise if requested
    anat.SetVolume(dimx,dimy,dimz,vb_float);
  }

  // SPECIAL CASE-- vecs if dimy=dimz=dimt==0
  if (dimt==0 && dimy==0 && dimz==0) {
    printf("[I] vbsim: creating a vector with %d elements with N(%g,%g)\n",dimx,
           n_mean,sqrt(n_variance));
    if (outfile=="")
      outfile="data.ref";
    VB_Vector vv(dimx);
    
//.........这里部分代码省略.........