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

void
cmp3d(Cube &c1,Cube &c2,Cube &mask,string spacer)
{
  if (!(c1.dimsequal(c2))) {
    printf("[I] vbcmp: %sdimensions don't match, can't compare.\n",spacer.c_str());
    return;
  }
  if (mask.data && !(c1.dimsequal(mask))) {
    printf("[I] vbcmp: %smask dimensions don't agree with images\n",spacer.c_str());
    return;
  }
  double totaldiff,maxdiff=0.0,v1,v2,diff;
  uint32 diffcount,voxelcount;
  int i,j,k;

  diffcount=0;
  totaldiff=0.0;
  int maxx=0,maxy=0,maxz=0;
  voxelcount=c1.dimx * c1.dimy * c1.dimz;
  // pre-count mask
  if (mask.data) {
    int newvoxelcount=0;
    for (size_t i=0; i<voxelcount; i++)
      if (mask.testValue(i))
        newvoxelcount++;
    voxelcount=newvoxelcount;
  }
  VB_Vector vec1(voxelcount),vec2(voxelcount);
  uint64 ind=0;
  for (i=0; i<c1.dimx; i++) {
    for (j=0; j<c1.dimy; j++) {
      for (k=0; k<c1.dimz; k++) {
        if (mask.data)
          if (!(mask.testValue(i,j,k))) continue;
        v1=c1.GetValue(i,j,k);
        v2=c2.GetValue(i,j,k);
        // to support correlation, smush into vectors:
        vec1[ind]=v1; vec2[ind]=v2;  ind++;        
        if (v1 != v2) {
          if (isnan(v1) && isnan(v2)) continue;
          if (isinf(v1) && isinf(v2)) continue;
          diffcount++;
          diff=fabs(v1-v2);
          totaldiff += diff;
          if (diff>maxdiff) {
            maxdiff=diff;
            maxx=i;
            maxy=j;
            maxz=k;
          }
        }
      }
    }
  }
  if (!diffcount) {
    printf("[I] vbcmp: %sthe data are identical\n",spacer.c_str());
  }
  else {
    cout << format("[I] vbcmp: %sdifferent voxels: %ld of %ld (%.0f%%)\n")
      % spacer % diffcount % voxelcount % ((diffcount*100.0)/voxelcount);
    cout << format("[I] vbcmp: %smean difference: %.8f\n") % spacer % (totaldiff/diffcount);
    cout << format("[I] vbcmp: %smax difference: %.8f (%d,%d,%d)\n")
      % spacer % maxdiff % maxx % maxy % maxz;
    cout << format("[I] vbcmp: %scorrelation between images: %g\n")%spacer%correlation(vec1,vec2);
  }
}

void
do_sumcube(string fname,int index,Cube &mycube,Cube &mask,double thresh,vector<VBRegion> &rlist)
{ 
  Cube tmpmask=mask; 
  int i,j,k;

  printf("%s-%02d\n",fname.c_str(),index);
  // eliminate from mask any sub-threshold voxels
  for (i=0; i<mycube.dimx; i++) {
    for (j=0; j<mycube.dimy; j++) {
      for (k=0; k<mycube.dimz; k++) {
        if (mycube.GetValue(i,j,k) < thresh)
          tmpmask.SetValue(i,j,k,0);
      }
    }
  }

  // now find regions
  vector<VBRegion> regionlist;
  regionlist=findregions(tmpmask,vb_gt,0.0);

  // now that we have a list of regions we can go through them and produce some stats
  double grandsum=0.0;
  double val;
  int grandcnt=0;
  for (i=0; i<(int)regionlist.size(); i++) {
    double sum=0.0,x=0.0,y=0.0,z=0.0;
    int cnt=0;
    for (VI myvox=regionlist[i].begin(); myvox!=regionlist[i].end(); myvox++) {
      if (val==0.0 && nozeroflag)
        continue;
      sum+=val;
      grandsum+=val;
      x+=myvox->second.x;  y+=myvox->second.y;  z+=myvox->second.z;
      cnt++;
      grandcnt++;
    }
    if (detailflag)
      printf("   region %d: %d voxels, mean %.2f (%d,%d,%d)\n",i,cnt,sum/(double)cnt,
             (int)x/cnt,(int)y/cnt,(int)z/cnt);
  }
  printf("  sum: %.3f  count: %d  average(sum/cnt): %.6f\n",grandsum,grandcnt,grandsum/(double)grandcnt);
  
  // now do the region comparisons
  if (rlist.size()) {
    // build a vb_vector for each regionlist
    vector<VB_Vector> myvecs;
    for (i=0; i<(int)rlist.size(); i++) {
      vector<double> vals;
      for (VI myvox=rlist[i].begin(); myvox!=rlist[i].end(); myvox++) {
        val=mycube.GetValue(myvox->second.x,myvox->second.y,myvox->second.z);
        if (val==0.0 && nozeroflag)
          continue;
        if (val<thresh)
          continue;
        vals.push_back(val);
      }
      VB_Vector tmpv(vals.size());
      for (j=0; j<(int)vals.size(); j++)
        tmpv[j]=vals[j];
      myvecs.push_back(tmpv);
    }
    for (i=0; i<(int)rlist.size()-1; i++) {
      for (j=i+1; j<(int)rlist.size(); j++) {
        double tvalue=ttest(myvecs[i],myvecs[j]);
        cout << "  region comparison: " << rlist[i].name << " vs. " << rlist[j].name << endl;
	cout << "    means: " << myvecs[i].getVectorMean() << " vs. " << myvecs[j].getVectorMean() << endl;
	cout << "    voxels: " << myvecs[i].getLength() << " vs. " << myvecs[j].getLength() << endl;
	cout << "    t value= " << tvalue << endl;
      }
    }
  }
}

double
Extractor::regioncode(int x,int y,int z)
{
  vector<value> valuelist;
  double x0,x1,y0,y1,z0,z1;
  int xx0,xx1,yy0,yy1,zz0,zz1;
  double val;
  int added;

  // xyz are the newcube coordinates of the voxel we'd like.  subtract
  // the origin of the new cube and divide by newvox to get that in mm
  // relative to the origin.

  x0=(double)(x-newcube.origin[0])*newcube.voxsize[0];
  x1=(double)((x+1)-newcube.origin[0])*newcube.voxsize[0];
  y0=(double)(y-newcube.origin[1])*newcube.voxsize[1];
  y1=(double)((y+1)-newcube.origin[1])*newcube.voxsize[1];
  z0=(double)(z-newcube.origin[2])*newcube.voxsize[2];
  z1=(double)((z+1)-newcube.origin[2])*newcube.voxsize[2];

  // cout << "x: " << x0 << ":" << x1 << endl;
  // cout << "y: " << y0 << ":" << y1 << endl;
  // cout << "z: " << z0 << ":" << z1 << endl;

  // now then divide by oldvox and subtract from original origin to
  // get original cube voxel.

  xx0=original.origin[0]+(int)(x0/original.voxsize[0]);
  xx1=original.origin[0]+(int)(x1/original.voxsize[0]);
  yy0=original.origin[1]+(int)(y0/original.voxsize[1]);
  yy1=original.origin[1]+(int)(y1/original.voxsize[1]);
  zz0=original.origin[2]+(int)(z0/original.voxsize[2]);
  zz1=original.origin[2]+(int)(z1/original.voxsize[2]);

  // cout << "coords: " << xx0 << ":" << xx1;
  // cout << ", " << yy0 << ":" << yy1;
  // cout << ", " << zz0 << ":" << zz1 << endl;
   
  for (int i=xx0; i<xx1; i++) {
    if (i<0 || i>=original.dimx)
      continue;
    for (int j=yy0; j<yy1; j++) {
      if (j<0 || j>=original.dimy)
	continue;
      for (int k=zz0; k<zz1; k++) {
	if (k<0 || k>=original.dimz)
	  continue;
	val=original.GetValue(i,j,k);
	added=0;
	for (int m=0; m<(int)valuelist.size(); m++) {
	  if (valuelist[m].val==val) {
	    valuelist[m].cnt++;
	    added=1;
	    break;
	  }
	}
	if (added==0) {
	  value tmp;
	  tmp.val=val;
	  tmp.cnt=1;
	  valuelist.push_back(tmp);
	}
	break;
      }
    }
  }
  // cout << "voxels: " << (xx1-xx0)*(yy1-yy0)*(zz1-zz0) << " ";
  // cout << "vals: " << valuelist.size() << " ";
  int commoncnt=0;
  double commonval=0.0;
  for (int i=0; i<(int)valuelist.size(); i++) {
    if (valuelist[i].cnt > commoncnt) {
      commonval=valuelist[i].val;
      commoncnt=valuelist[i].cnt;
    }
  }
  // cout << "common: " << commonval << "(" << commoncnt << ")" << endl;
  return (commonval);
}