C++ volume::xsize方法代码示例

// transforms a 4D Volume in a SVM samples file, based on a mask
void saveSVMFile(volume4D <float> &volSamples, volume<float> &mask, const char *outputFileName, float minValue, vector <int > &indexes, vector <int> &classes)
{
	FILE *f;
	f = fopen(outputFileName, "wt+");
	if (f != NULL)
	{
		int i, t;
		for (int h = 0; h < indexes.size(); h++)
		{
			// picking the right indexes
			t = indexes[h] - 1;
			i = 0;
			// saving the class value. Remember indexes size is different from classes. classes has the same sime of the number of volumes
			fprintf(f, "%d ", classes[t]);
			for (int z = 0; z < mask.zsize(); z++)
				for (int y = 0; y < mask.ysize(); y++)
					for (int x = 0; x < mask.xsize(); x++)
						if (mask.value(x, y, z) > minValue)
						{
							i++;
							// writing each voxel value in svm format
							fprintf(f, "%d:%f ", i, volSamples.value(x, y, z, t));
						}
			fprintf(f, "\n");
		}
		fclose(f);
	}
}

FnirtFileWriter::FnirtFileWriter(const string&              fname,
				 const volume<float>&       fieldx,                         
			         const volume<float>&       fieldy,                         
			         const volume<float>&       fieldz)
{
  volume<float>  tmp(fieldx.xsize(),fieldx.ysize(),fieldx.zsize());
  tmp.copyproperties(fieldx);
  Matrix         aff = IdentityMatrix(4);

  common_field_construction(fname,tmp,fieldx,fieldy,fieldz,aff);
}

void basisfield::Set(const volume<float>& pfield)
{
  if (int(FieldSz_x()) != pfield.xsize() || int(FieldSz_y()) != pfield.ysize() || int(FieldSz_z()) != pfield.zsize()) {
    throw BasisfieldException("basisfield::Set:: Matrix size mismatch beween basisfield class and supplied field");
  }
  if (Vxs_x() != pfield.xdim() || Vxs_y() != pfield.ydim() || Vxs_z() != pfield.zdim()) {
    throw BasisfieldException("basisfield::Set:: Voxel size mismatch beween basisfield class and supplied field");
  }

  volume<float>   volume_of_ones(pfield.xsize(),pfield.ysize(),pfield.zsize());
  volume_of_ones.copyproperties(pfield);
  volume_of_ones = 1.0;

  double lambda = 0.001;
  ColumnVector y = Jte(pfield,0);
  boost::shared_ptr<MISCMATHS::BFMatrix>  XtX = JtJ(volume_of_ones);
  boost::shared_ptr<MISCMATHS::BFMatrix>  BeEn = BendEnergyHess();
  XtX->AddToMe(*BeEn,lambda);
  ColumnVector coef_roof = XtX->SolveForx(y,SYM_POSDEF,1e-6,500);
  SetCoef(coef_roof);
}

volume<float> inside_mesh(const volume<float> & image, const Mesh& m)
{
  volume<float> res = image;
  int xsize = image.xsize();
  int ysize = image.ysize();
  int zsize = image.zsize();
  volume<short> inside = make_mask_from_mesh(image, m);
  for (int k=0; k<zsize; k++)
    for (int j=0; j<ysize; j++)
      for (int i=0; i<xsize; i++)
	res.value(i, j, k) = (1-inside.value(i, j, k)) * image.value(i, j, k);
  return res;
}

void _volume2Sample(svm_model *model, volume<float> &volSample, volume<float> &mask, int sampleSize, float minValue, svm_node * &sample)
{	
	sample=(struct svm_node *) malloc((sampleSize+1)*sizeof(struct svm_node));
    
	int i=0;
    for(int z=0; z < volSample.zsize();z++)
      for(int y=0; y < volSample.ysize();y++)
	     for(int x=0; x < volSample.xsize();x++)
         if (mask.value(x,y,z) > minValue)
		   {
             sample[i].index = (i+1);
		       sample[i].value = volSample.value(x,y,z);
				 i++;
			}
    sample[i].index = -1;
}

// returns in `output` the best voxels of `region`
void RegionExtraction::regionBestVoxels(RoiMeanCalculation &reference, volume<float>&values, volume<float>&output, int region, int regionSize, float percentage)
{
   vector<roiPoint> roi;
   roi.resize(regionSize);
   greaterRoiPoint greaterFirst;
   int index=0;
   
   // Filter the voxels of the region `region`
   for(int z=0;z<values.zsize();z++)
      for(int y=0;y<values.ysize();y++)
         for(int x=0;x<values.xsize();x++)
         {
            // get the voxel intensity in reference
            int voxelRegion = (int) reference.voxelValues(x,y,z);
            
            // if is the chosen region, records the voxel values (T values or other voxel value)
            if (region == voxelRegion)
            {
               roi[index].value = values.value(x,y,z);
               roi[index].roiValue = region;
               
               roi[index].x=x;
               roi[index].y=y;
               roi[index].z=z;
               index++;
            }
            else values.value(x,y,z)=(float)0.0;
         }
   
   // sorts the vector of values in descending order
   std::sort(roi.begin(), roi.end(), greaterFirst);
   
   // calculates the cut Index. Remember the voxels descending order
   int cutIndex = (int) (roi.size() * percentage + 0.5);
   
   // recording the result in `output`
   for (int j=0; j<=cutIndex;j++)
   {
      output.value(roi[j].x, roi[j].y, roi[j].z) = roi[j].roiValue;
   }
}

void basisfield::AsVolume(volume<float>& vol, FieldIndex fi)
{
  if (int(FieldSz_x()) != vol.xsize() || int(FieldSz_y()) != vol.ysize() || int(FieldSz_z()) != vol.zsize()) {
    throw BasisfieldException("basisfield::AsVolume:: Matrix size mismatch beween field and volume");
  }
  if (Vxs_x() != vol.xdim() || Vxs_y() != vol.ydim() || Vxs_z() != vol.zdim()) {
    throw BasisfieldException("basisfield::AsVolume:: Voxel size mismatch beween field and volume");
  }
  if (!coef) {throw BasisfieldException("basisfield::AsVolume: Coefficients undefined");}

  if (!UpToDate(fi)) {Update(fi);}

  const boost::shared_ptr<NEWMAT::ColumnVector> tmptr = Get(fi);
  int vindx=0;
  for (unsigned int k=0; k<FieldSz_z(); k++) {
    for (unsigned int j=0; j<FieldSz_y(); j++) {
      for (unsigned int i=0; i<FieldSz_x(); i++) {
        vol(i,j,k) = tmptr->element(vindx++);
      }
    }
  }
}

void FnirtFileWriter::common_field_construction(const string&            fname,
                                                const volume<float>&     ref,
                                                const volume<float>&     fieldx,
                                                const volume<float>&     fieldy,
                                                const volume<float>&     fieldz, 
                                                const Matrix&            aff)
{
  volume4D<float>   fields(ref.xsize(),ref.ysize(),ref.zsize(),3);
  fields.copyproperties(ref); 

  Matrix M;
  bool   add_affine = false;
  if (add_affine = ((aff-IdentityMatrix(4)).MaximumAbsoluteValue() > 1e-6)) { // Add affine part to fields
    M = (aff.i() - IdentityMatrix(4))*ref.sampling_mat();
  }

  if (samesize(ref,fieldx,true)) { // If ref is same size as the original field
    fields[0] = fieldx; fields[1] = fieldy; fields[2] = fieldz;
    fields.copyproperties(ref);    // Put qform/sform and stuff back.
    if (add_affine) {
      ColumnVector xv(4), xo(4);
      int zs = ref.zsize(), ys = ref.ysize(), xs = ref.xsize();
      xv(4) = 1.0;
      for (int z=0; z<zs; z++) {
        xv(3) = double(z);
        for (int y=0; y<ys; y++) {
          xv(2) = double(y);
          for (int x=0; x<xs; x++) {
            xv(1) = double(x);
            xo = M*xv;
            fields(x,y,z,0) += xo(1);
            fields(x,y,z,1) += xo(2);
            fields(x,y,z,2) += xo(3);
	  }
        }
      }
    } 
  }
  else {
    fieldx.setextrapolationmethod(extraslice);
    fieldy.setextrapolationmethod(extraslice);
    fieldz.setextrapolationmethod(extraslice);
    Matrix R2F = fieldx.sampling_mat().i() * ref.sampling_mat();
    ColumnVector xv(4), xo(4), xr(4);
    int zs = ref.zsize(), ys = ref.ysize(), xs = ref.xsize();
    xv(4) = 1.0;
    for (int z=0; z<zs; z++) {
      xv(3) = double(z);
      for (int y=0; y<ys; y++) {
        xv(2) = double(y);
        for (int x=0; x<xs; x++) {
          xv(1) = double(x);
          xr = R2F*xv;
          fields(x,y,z,0) = fieldx.interpolate(xr(1),xr(2),xr(3));
          fields(x,y,z,1) = fieldy.interpolate(xr(1),xr(2),xr(3));
          fields(x,y,z,2) = fieldz.interpolate(xr(1),xr(2),xr(3));
          if (add_affine) {
            xo = M*xv;
            fields(x,y,z,0) += xo(1);
            fields(x,y,z,1) += xo(2);
            fields(x,y,z,2) += xo(3);
	  }
        }
      }
    }
  }

  fields.set_intent(FSL_FNIRT_DISPLACEMENT_FIELD,fields.intent_param(0),fields.intent_param(1),fields.intent_param(2));
  fields.setDisplayMaximum(0.0);
  fields.setDisplayMinimum(0.0);

  // Save resulting field
  save_volume4D(fields,fname);
}

bet_parameters adjust_initial_mesh(const volume<float> & image, Mesh& m, const double & rad = 0., const double xpara=0.,  const double ypara=0.,  const double zpara=0.)
{
  bet_parameters bp;
  double xdim = image.xdim();
  double ydim = image.ydim();
  double zdim = image.zdim();
  double t2, t98, t;

  //computing t2 && t98
  //  cout<<"computing robust min && max begins"<<endl;

  bp.min = image.min();
  bp.max = image.max();

  t2 = image.robustmin();
  t98 = image.robustmax();
  //t2=32.;
  //t98=16121.;
  
  //  cout<<"computing robust min && max ends"<<endl;
  
  t = t2 + .1*(t98 - t2);
  bp.t98 = t98;
  bp.t2 = t2;
  bp.t = t;
  //  cout<<"t2 "<<t2<<" t98 "<<t98<<" t "<<t<<endl;
  
  //  cout<<"computing center && radius begins"<<endl;
  
  //finds the COG
  Pt center(0, 0, 0);
  double counter = 0;
  if (xpara == 0. & ypara==0. & zpara==0.)
    {
      double tmp = t - t2;
      for (int k=0; k<image.zsize(); k++)
	for (int j=0; j<image.ysize(); j++)
	  for (int i=0; i<image.xsize(); i++)
	    {
	      double c = image(i, j, k ) - t2;
	      if (c > tmp)
		{
		  c = min(c, t98 - t2);   
		  counter+=c;
		  center +=  Pt(c*i*xdim, c*j*ydim, c*k*zdim);
		}
	    }
      center=Pt(center.X/counter, center.Y/counter, center.Z/counter);
      //cout<<counter<<endl;
      //  cout<<"cog "<<center.X<<" "<<center.Y<<" "<<center.Z<<endl;
    }
  else center=Pt(xpara, ypara, zpara);
  
  bp.cog = center;

  if (rad == 0.)
    {
      double radius=0;
      counter=0;
      double scale=xdim*ydim*zdim;
      for (int k=0; k<image.zsize(); k++)
	for (int j=0; j<image.ysize(); j++)
	  for (int i=0; i<image.xsize(); i++)
	    {
	      double c = image(i, j, k);
	      if (c > t)
		{
		  counter+=1;
		}
	    }
      radius = pow (.75 * counter*scale/M_PI, 1.0/3.0);
      //      cout<<radius<<endl;
      bp.radius = radius;
    } 
  else (bp.radius = rad);

  m.translation(center.X, center.Y, center.Z);
  m.rescale (bp.radius/2, center);

  //  cout<<"computing center && radius ends"<<endl;

  //computing tm
  //  cout<<"computing tm begins"<<endl;
  vector<double> vm;
  for (int k=0; k<image.zsize(); k++)
    for (int j=0; j<image.ysize(); j++)
      for (int i=0; i<image.xsize(); i++)
	{
	  double d = image.value(i, j, k);
	  Pt p(i*xdim, j*ydim, k*zdim);
	  if (d > t2 && d < t98 && ((p - center)|(p - center)) < bp.radius * bp.radius)
	    vm.push_back(d);
	}

  int med = (int) floor(vm.size()/2.);
  //  cout<<"before sort"<<endl;
  nth_element(vm.begin(), vm.begin() + med - 1, vm.end());
  //partial_sort(vm.begin(), vm.begin() + med + 1, vm.end());
  //double tm = vm[med];
  double tm=(*max_element(vm.begin(), vm.begin() + med));
//.........这里部分代码省略.........