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

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);
  }
}

void
Extractor::Go()
{
  double val;

  for (int i=0; i<newcube.dimx; i++) {
    for (int j=0; j<newcube.dimy; j++) {
      for (int k=0; k<newcube.dimz; k++) {
	val=regioncode(i,j,k);
	if (xvflag==0 || val==xval)
	  newcube.SetValue(i,j,k,val);
      }
    }
  }
  newcube.WriteFile();

  return;
}

//.........这里部分代码省略.........
      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
            res=calc_chisquared(bm,dvbm,f_yates);
          if (f_twotailed)
            res.p*=2.0;
          // FIXME??? if we need the p value or z score, get it
          // if doing fdr, stash the p value
          if (f_fdr || pmapname.size()) {
            pmap.SetValue(i,j,k,res.p);
            pvals.push_back(res.p);
          }
          // if doing confidence intervals, stash that
          if (f_zscore)
            statmap.SetValue(i,j,k,res.z);
          else
            statmap.SetValue(i,j,k,res.x2);
          statlookup[bm]=res;
        }
        else {
          // this is a previously encountered pattern
          if (f_fdr && !f_nodup) {
            fdrmask.SetValue(i,j,k,1); 
            pvals.push_back(iter->second.p);
          }
          statmap.SetValue(i,j,k,iter->second.x2);
          if (f_fdr || pmapname.size())
            pmap.SetValue(i,j,k,iter->second.p);
        }
      }
    }
  }
  if (f_non1) {
    cout << "[W] vbcmap: non-0/1 values found in lesion map" << endl;
  }
  if (f_fdr) {
    vector<fdrstat> ffs=calc_multi_fdr_thresh(statmap,pmap,fdrmask,q);
    if (ffs.size()) {
      cout << (format("[I] vbcmap: FDR calculation included %d voxels with p values from %.4f to %.4f\n")
               %ffs[0].nvoxels%ffs[0].low%ffs[0].high).str();
      statmap.AddHeader("# the following thresholds must be exceeded for FDR control");
      vbforeach(fdrstat ff,ffs) {
        if (ff.maxind>=0)
          cout << (format("[I] vbcmap: FDR threhsold for q=%.2f is %.4f\n")%ff.q%ff.statval).str();
        else
          cout << (format("[I] vbcmap: no FDR threhsold could be identified for q=%.2f\n")%ff.q).str();
        string tmps=(format("fdrthresh: %g %g")%ff.q%ff.statval).str().c_str();
        statmap.AddHeader(tmps);
      }
    }
  }