C++ volume类代码示例

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

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

//calculates the mask from the mesh by spreading an initial point outside the mesh, and stopping it when the mesh is reached.
volume<short> make_mask_from_mesh(const volume<float> & image, const Mesh& m)
{
  //  cout<<"make_mask_from_mesh begins"<<endl;

  double xdim = (double) image.xdim();
  double ydim = (double) image.ydim();
  double zdim = (double) image.zdim();

  volume<short> mask;
  copyconvert(image,mask);
  
  int xsize = mask.xsize();
  int ysize = mask.ysize();
  int zsize = mask.zsize();
  
  mask = 1;
  
  mask = draw_mesh(mask, m);

  vector<Pt> current;
  current.clear();
  Pt c(0., 0., 0.);
  for (vector<Mpoint *>::const_iterator it=m._points.begin(); it!=m._points.end(); it++)
    c+=(*it)->get_coord();

  c*=(1./m._points.size());
  c.X/=xdim; c.Y/=ydim; c.Z/=zdim;

  current.push_back(c);

  while (!current.empty())
    {
      Pt pc = current.back();
      int x, y, z;
      x=(int) pc.X; y=(int) pc.Y; z=(int) pc.Z;
      //current.remove(pc);
      current.pop_back();
      mask.value(x, y, z) = 0;
      if (0<=x-1 && mask.value(x-1, y, z)==1) current.push_back(Pt(x-1, y, z));
      if (0<=y-1 && mask.value(x, y-1, z)==1) current.push_back(Pt(x, y-1, z));
      if (0<=z-1 && mask.value(x, y, z-1)==1) current.push_back(Pt(x, y, z-1));
      if (xsize>x+1 && mask.value(x+1, y, z)==1) current.push_back(Pt(x+1, y, z));
      if (ysize>y+1 && mask.value(x, y+1, z)==1) current.push_back(Pt(x, y+1, z));
      if (zsize>z+1 && mask.value(x, y, z+1)==1) current.push_back(Pt(x, y, z+1)); 
    }

  //  cout<<"make_mask_from_mesh ends"<<endl;
  return mask;
}

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

//extracts profiles in the area around the eyes
vector<double> t1only_special_extract(const volume<float> & t1, const Pt & point, const Vec & n) {
  vector<double> resul;
  resul.clear();
  bool output = true; 
  const double INNER_DEPTH = 3;
  const double OUTER_DEPTH = 100;
  
  Profile pt1;
  for (double d = -INNER_DEPTH; d < OUTER_DEPTH; d+=.5)
    {
      Pt c = point + d * n;
      double tmp1 = t1.interpolate(c.X, c.Y, c.Z);
      pt1.add (d, tmp1);
    }
  pt1.init_roi();

  //outer skin
  double outskin = pt1.last_point_over(pt1.end(), .2);
  double check = pt1.last_point_over(outskin - 1.5, .2);
  if (outskin - check > 2) output = false;
  
  pt1.set_rroi(outskin);
  
  double inskull = pt1.next_point_under(-INNER_DEPTH, .25);
  if (inskull > 5) inskull = 0;
  
  double outskull = pt1.next_point_over(inskull, .35);
  
  resul.push_back(inskull);
  resul.push_back(outskull);
  resul.push_back(outskin);
  
  return resul;
}

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

void draw_segment(volume<short>& image, const Pt& p1, const Pt& p2)
{
  double xdim = (double) image.xdim();
  double ydim = (double) image.ydim();
  double zdim = (double) image.zdim();
  double mininc = min(xdim,min(ydim,zdim)) * .5;
 
  Vec n = p1 - p2;
  double d = n.norm();
  n.normalize();

  for (double i=0; i<=d; i+=mininc)
    {
      Pt p = p2 + i* n;
      image((int) floor((p.X)/xdim +.5),(int) floor((p.Y)/ydim +.5),(int) floor((p.Z)/zdim +.5)) = 1;
    }
}

void draw_mesh(volume<short>& image, const Mesh &m)
{
  double xdim = (double) image.xdim();
  double ydim = (double) image.ydim();
  double zdim = (double) image.zdim();
  double mininc = min(xdim,min(ydim,zdim)) * .5;

  for (list<Triangle*>::const_iterator i = m._triangles.begin(); i!=m._triangles.end(); i++)
    {
      Vec n = (*(*i)->get_vertice(0) - *(*i)->get_vertice(1));
      double d = n.norm();
      n.normalize();

      for (double j=0; j<=d; j+=mininc)
	{
	  Pt p = (*i)->get_vertice(1)->get_coord()  + j* n;
	  draw_segment(image, p, (*i)->get_vertice(2)->get_coord());
	} 
    }

}

void check_integral(fields &f,
		    linear_integrand_data &d, const volume &v,
		    component cgrid)
{
  double x1 = v.in_direction_min(d.dx);
  double x2 = v.in_direction_max(d.dx);
  double y1 = v.in_direction_min(d.dy);
  double y2 = v.in_direction_max(d.dy);
  double z1 = v.in_direction_min(d.dz);
  double z2 = v.in_direction_max(d.dz);
  
  master_printf("Check %d-dim. %s integral in %s cell with %s integrand...",
		(x2 - x1 > 0) + (y2 - y1 > 0) + (z2 - z1 > 0), 
		component_name(cgrid),
		v.dim == D3 ? "3d" : (v.dim == D2 ? "2d" : 
				       (v.dim == Dcyl ? "cylindrical" 
					: "1d")),
		(d.c == 1.0 && !d.axy && !d.ax && !d.ay && !d.az
		 && !d.axy && !d.ayz && !d.axz) ? "unit" : "linear");
  if (0)
    master_printf("\n... grid_volume (%g,%g,%g) at (%g,%g,%g) with integral (%g, %g,%g,%g, %g,%g,%g, %g)...\n",
		  x2 - x1, y2 - y1, z2 - z1,
		  (x1+x2)/2, (y1+y2)/2, (z1+z2)/2,
		  d.c, d.ax,d.ay,d.az, d.axy,d.ayz,d.axz, d.axyz);

  double sum = real(f.integrate(0, 0, linear_integrand, (void *) &d, v));
  if (fabs(sum - correct_integral(v, d)) > 1e-9 * fabs(sum))
    abort("FAILED: %0.16g instead of %0.16g\n", 
	  (double) sum, correct_integral(v, d));
  master_printf("...PASSED.\n");
}

volume<float> draw_mesh_bis(const volume<float>& image, const Mesh &m)
{
  double xdim = (double) image.xdim();
  double ydim = (double) image.ydim();
  double zdim = (double) image.zdim();
  double mininc = min(xdim,min(ydim,zdim)) * .5;
  volume<float> res = image;
  double max = image.max();
  for (list<Triangle*>::const_iterator i = m._triangles.begin(); i!=m._triangles.end(); i++)
    {
      Vec n = (*(*i)->get_vertice(0) - *(*i)->get_vertice(1));
      double d = n.norm();
      n.normalize();

      for (double j=0; j<=d; j+=mininc)
	{
	  Pt p = (*i)->get_vertice(1)->get_coord()  + j* n;
	  draw_segment_bis(res, p, (*i)->get_vertice(2)->get_coord(), max);
	} 
    }
  return res;
}

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

void vega::data::hexagonal_prismatic_lattice::fill_volume_cell( volume& v, const prismatic_hexagon_node & node )
{
    static const int hex_pixel_offset[8][2] = { {1, 0}, {2, 0}, {0, 1}, {1, 1}, {2, 1}, {3, 1}, {1, 2}, {2, 2}};

    math::vector3d vertex = node.get_vertex();

    for(size_t h=0;h<8;++h)
    {
        math::vector3d pixel;
        int xx = (int)vertex.x + hex_pixel_offset[h][0];
        int yy = (int)vertex.y + hex_pixel_offset[h][1];
        int zz = (int)vertex.z;
        if( xx >= 0 && yy >= 0 && zz >= 0 && xx < v.get_width() && yy < v.get_height() && zz < v.get_depth() )
        {
#ifndef USE_COLOR_MAP
            v.set_voxel(xx, yy, zz, (voxel)(node.density * 255.f));
#else
            v.set_voxel(xx, yy, zz, (voxel)(node.color.A));
#endif
        }
    }
}

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

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