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

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

// calculates the mean for each roi
void RoiMeanCalculation::calculateMeans(volume<float> &actualvolume)
{
	if (means.size())
	{
		for (int i = 0; i < means.size(); i++) means[i] = 0;

		for (int z = 0; z < reference.zsize(); z++)
			for (int y = 0; y < reference.ysize(); y++)
				for (int x = 0; x < reference.xsize(); x++)
				{
					float value = reference.value(x, y, z);
					if (value != 0)
					{
						int idx = mapping[value];
						if (idx >= -1)
							means[idx] += actualvolume.value(x, y, z);
					}
				}

		for (int i = 0; i < means.size(); i++)
		{
			if (counts[i]) means[i] = means[i] / (float)counts[i];
		}
	}
}

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

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

// transforms an array in a volume
void array2Volume(const char *maskFile, float minValue, vector <double> &weightVector, volume<float> &weightVolume)
{
	volume<float> mask;
    if (maskFile != NULL)
    {
      string Maskfile = maskFile;
      read_volume(mask, Maskfile);
    }

	weightVolume.reinitialize(mask.xsize(), mask.ysize(), mask.zsize(), 0, true);
    weightVolume.copyproperties(mask);
	int i = 0;
	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)
			 {
				 weightVolume.value(x,y,z) = (float) weightVector[i];
				 i++;
			 }
			 else weightVolume.value(x,y,z) = (float) 0.0;
}

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

double step_of_computation(const volume<float> & image, Mesh & m, const double bet_main_parameter, const int pass, const double increase_smoothing, const int iteration_number, double & l, const double t2, const double tm, const double t, const double E,const double F, const double zcog, const double radius, const double local_th=0., const int d1=7, const int d2=3){
  double xdim = image.xdim();
  double ydim = image.ydim();
  double zdim = image.zdim();
  
  if (iteration_number==50 || iteration_number%100 == 0 )
    {
      l = 0;
      int counter = 0;
      for (vector<Mpoint*>::iterator i = m._points.begin(); i!=m._points.end(); i++ )
	{
	  counter++;
	  l += (*i)->medium_distance_of_neighbours();
	}
      l/=counter;
    }
  
  for (vector<Mpoint*>::iterator i = m._points.begin(); i!=m._points.end(); i++)
    {
      Vec sn, st, u1, u2, u3, u;
      double f2, f3=0;
      
      Vec n = (*i)->local_normal();
      Vec dv = (*i)->difference_vector();
      
      double tmp = dv|n;
      sn = n * tmp;
      st = dv - sn;
      
      u1 = st*.5;
      
      double rinv = (2 * fabs(sn|n))/(l*l);
      
      f2 = (1+tanh(F*(rinv - E)))*0.5;
      if (pass > 0)
	if (tmp > 0) {
	  f2*=increase_smoothing;
	  f2 = Min(f2, 1.);
	}
      
      u2 = f2 * sn;
      
      //main term of bet
      {
	double local_t = bet_main_parameter;
	if (local_th != 0.0)
	  {
	    local_t = Min(1., Max(0., bet_main_parameter + local_th*((*i)->get_coord().Z - zcog)/radius));
	  }
	
	double Imin = tm;
	double Imax = t;
	
	Pt p = (*i)->get_coord() + (-1)*n;
	double iv = p.X/xdim + .5, jv = p.Y/ydim +.5, kv = p.Z/zdim +.5; 
	if (image.in_bounds((int)iv,(int) jv,(int) kv))
	  {	
	    double im=image.value((int)iv,(int)jv,(int)kv);
	    Imin = Min(Imin, im);
	    Imax = Max(Imax,im);
	    
	    double nxv=n.X/xdim, nyv=n.Y/ydim, nzv=n.Z/zdim;
	    int i2=(int)(iv-(d1-1)*nxv), j2 =(int) (jv-(d1-1)*nyv), k2 =(int)(kv-(d1-1)*nzv); 
	    if (image.in_bounds(i2, j2, k2))
	      {	
		im=image.value(i2,j2,k2);
		Imin = Min(Imin, im);
		
		for (int gi=2; gi<d1; gi++)
		  {
		    // the following is a quick calc of Pt p = (*i)->get_coord() + (-gi)*n;
		    iv-=nxv; jv-=nyv; kv-=nzv;
		    im = image.value((int) (iv), (int) (jv), (int) (kv));
		    Imin = Min(Imin, im);
		
		    if (gi<d2)
		      Imax = Max(Imax,im);
		  }
		
		Imin = Max (t2, Imin);
		Imax = Min (tm, Imax);	
		
		const double tl = (Imax - t2) * local_t + t2;
		
		if (Imax - t2 > 0)
		  f3=2*(Imin - tl)/(Imax - t2);
		else f3=(Imin - tl)*2;
	      }
	  }
	
      }
      
      f3 *= (normal_max_update_fraction * lambda_fit * l);
      
      u3 = f3 * n;
      
      u = u1 + u2 + u3;
            
      //cout<<"l "<<l<<"u1 "<<((u1*n).norm())<<"u2 "<<(u2|n)<<"u3 "<<(u3|n)<<endl;
      
//.........这里部分代码省略.........

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));
//.........这里部分代码省略.........