C++ volmagick::Volume类代码示例

int main(int argc,char *argv[])
{
  VolMagick::Volume v;
  VolMagick::readVolumeFile(v,argv[1]);

  UniformGrid ugrid(1,1,1);
  ugrid.importVolume(v);
  Point min,max;
  min.x=v.XMin();
  min.y=v.YMin();
  min.z=v.ZMin();
  max.x=v.XMax();
  max.y=v.YMax();
  max.z=v.ZMax();

//generateGridFile(min,max,v.XDim(),v.YDim(),v.ZDim());

  Mesh *mesh=new Mesh();

  float iso_val=atof(argv[2]);

  Mesher mesher;

  mesher.setGrid((Grid &)(ugrid));
  mesher.setMesh(mesh);

  mesher.generateMesh(iso_val);

 // mesh->correctNormals();
  mesher.saveMesh(string(argv[3]),&v);
  delete mesh;
  
  return 0;
}

void getMaterialIds(VolMagick::Volume &vol,vector<int> &matIds)
{
  int xdim=vol.XDim(),
      ydim=vol.YDim(),
      zdim=vol.ZDim();

  set<int> materials;

  for(int i=0;i<xdim;i++)
    for(int j=0;j<ydim;j++)
      for(int k=0;k<zdim;k++)
        materials.insert(vol(i,j,k));

  for(set<int>::iterator sIter=materials.begin();sIter!=materials.end();sIter++)
    if(*sIter!=0)
        matIds.push_back(*sIter);
  
}

void generatePartitionSurfaceByMC(VolMagick::Volume& vol, int &matId, MMHLS::Mesh *mesh)
{
	VolMagick::Volume tempv;
	tempv = vol;

	for(int i=0; i<vol.XDim(); i++)
	 for(int j=0; j<vol.YDim(); j++)
	   for(int k=0; k<vol.ZDim(); k++)
	    {
			if(vol(i,j,k) == matId)
				tempv(i,j,k, 1.0);
			else 
				tempv(i,j,k, 0.0);
		}
	FastContouring::ContourExtractor *contexa;
	contexa = new FastContouring::ContourExtractor;

	contexa->setVolume(tempv);
	FastContouring::TriSurf result = contexa->extractContour(0.5);
	for(int i =0; i<result.verts.size(); i+=3)
	{
		vector<float> vert;
		vert.push_back(result.verts[i]);
		vert.push_back(result.verts[i+1]);
		vert.push_back(result.verts[i+2]);
		mesh->vertexList.push_back(vert);
	}

	for(int i=0; i<result.tris.size(); i+=3)
	{
		vector <unsigned int > t;
		t.push_back(result.tris[i]);
		t.push_back(result.tris[i+1]);
		t.push_back(result.tris[i+2]);
		mesh->faceList.push_back(t);
	}
	delete contexa;

}

int main(int argc, char **argv)
{
  if(argc != 3)
    {
      cerr << "Usage: " << argv[0] << " <input volume file> <output volume file>" << endl;
      return 1;
    }

  try
    {
      VolMagickOpStatus status;
      VolMagick::setDefaultMessenger(&status);

      VolMagick::VolumeFileInfo volinfo;

      volinfo.read(argv[1]);

      VolMagick::createVolumeFile(argv[2],
				  volinfo.boundingBox(),
				  volinfo.dimension(),
				  volinfo.voxelTypes(),
				  volinfo.numVariables(),
				  volinfo.numTimesteps(),
				  volinfo.TMin(),volinfo.TMax());

      for(unsigned int var=0; var<volinfo.numVariables(); var++)
	for(unsigned int time=0; time<volinfo.numTimesteps(); time++)
	  {
	    VolMagick::Volume vol;
	    readVolumeFile(vol,argv[1],var,time);
	    for(VolMagick::uint64 i = 0; i < vol.dimension().size(); i++)
	      vol(i, vol.min() + ((vol.max() - vol(i))/(vol.max() - vol.min()))*(vol.max() - vol.min()));
	    vol.desc(volinfo.name(var));
	    writeVolumeFile(vol,argv[2],var,time);
	  }
    }
  catch(VolMagick::Exception &e)
    {
      cerr << e.what() << endl;
    }
  catch(std::exception &e)
    {
      cerr << e.what() << endl;
    }

  return 0;
}

int main(int argc, char **argv)
{
  if(argc < 4)
    {
      std:: cerr << "make templete volume with initial value:" << std::endl;
      std:: cerr << "Usage: " << argv[0] << "<ref. vol> <set val(float)> <out vol>" << std::endl;
      return 1;
    }

  try
    {
      VolMagick::Volume outputVol;

      VolMagick::VolumeFileInfo volinfo;
      volinfo.read(argv[1]);

      float setval = atof( argv[2] );

      outputVol.voxelType(volinfo.voxelType());
      outputVol.dimension(volinfo.dimension());
      outputVol.boundingBox(volinfo.boundingBox());
   
      int dimx = outputVol.XDim();
      int dimy = outputVol.YDim();
      int dimz = outputVol.ZDim();
      fprintf( stderr, "dim: %d %d %d\n", dimx, dimy, dimz );

      for( int kz = 0; kz<dimz; kz++)
        for( int jy = 0; jy<dimy; jy++)
          for( int ix = 0; ix<dimx; ix++)
 	    outputVol(ix,jy,kz, setval);
      
      VolMagick::createVolumeFile(outputVol, argv[argc-1]);

    }

  catch(VolMagick::Exception &e)
    {
      std:: cerr << e.what() << std::endl;
    }
  catch(std::exception &e)
    {
      std::cerr << e.what() << std::endl;
    }

  return 0;
}

int main(int argc, char **argv)
{
  if(argc < 4)
    {
      std:: cerr << 
	"Usage: " << argv[0] << 
	"  first volume - second volume, output volume, threshold " << std::endl;
	      return 1;
    }

  try
    {
      VolMagick::Volume inputVol;

      VolMagick::Volume inputVol2;
      
      VolMagick::Volume outputVol;

      VolMagick::readVolumeFile(inputVol,argv[1]); ///first argument is input volume
      
      VolMagick::readVolumeFile(inputVol2, argv[2]); /// second argument is mask volume
      
	  float threshold = atof(argv[4]);

      VolMagick::VolumeFileInfo volinfo1;
      volinfo1.read(argv[1]);
      std::cout << volinfo1.filename() << ":" <<std::endl;

      VolMagick::VolumeFileInfo volinfo2;
      volinfo2.read(argv[2]);

      std::cout << volinfo2.filename() << ":" <<std::endl;
      
      std::cout<<"minVol1 , maxVol1: "<<volinfo1.min()<<" "<<volinfo1.max()<<std::endl;;

      std::cout<<"minVol2 , maxVol2: "<<volinfo2.min()<<" "<<volinfo2.max()<<std::endl;

      
      
      outputVol.voxelType(inputVol.voxelType());
      outputVol.dimension(inputVol.dimension());
      outputVol.boundingBox(inputVol.boundingBox());
      
      
      std::cout<<"voxeltype "<<inputVol.voxelType()<<std::endl;
      	
	 // int step = 1;
   //   if(argc==5) step= atoi(argv[argc-1]);

	  outputVol= inputVol;

      
      for( int kz = 0; kz<inputVol2.ZDim(); kz++)
	{
	    for( int jy = 0; jy<inputVol2.YDim(); jy++)
	      for( int ix = 0; ix<inputVol2.XDim(); ix++)
		{
	//	  if(inputVol2(ix, jy, kz)<= inputVol(ix, jy, kz)) 
	//	  {
		     outputVol(ix, jy, kz, inputVol(ix,jy,kz)-inputVol2(ix,jy,kz));
					
	//	  }
	//	  else outputVol(ix,jy,kz, inputVol.min());
			  
		}
	}	
      

      VolMagick::createVolumeFile(outputVol, argv[3]);





    }

  catch(VolMagick::Exception &e)
    {
      std:: cerr << e.what() << std::endl;
    }
  catch(std::exception &e)
    {
      std::cerr << e.what() << std::endl;
    }

  return 0;
}

int main(int argc, char **argv)
{
  if(argc < 4)
    {
      std:: cerr << 
	"Usage: " << argv[0] << 

"  inputmesh inputvol outputMesh [range]\n";
      return 1;
    }

  try
    {

	  geometry geom;
	  cvcraw_geometry::read(geom, string(argv[1]));

      VolMagick::Volume inputVol;


      VolMagick::readVolumeFile(inputVol,argv[2]); ///input volume
     
      double range = 0; bool gotRange = false;
      if( argc > 4 ) {
         range = atof( argv[4] );
         gotRange = true;
std::cout << "got range: " << range << std::endl;
      }
      
	  float x,y,z;

	  float value; 
	  
	  color_t color;

	  COLOUR c;

        if( !gotRange ) {
          // compute range, max absolute value
	  for(int i=0; i< geom.points.size(); i++)
	  {
	  		x = geom.points[i][0];
			y = geom.points[i][1];
			z = geom.points[i][2];


	  		value = inputVol.interpolate(x,y,z);
			if( range  < fabs( value ) ) range = fabs( value );
	   }
	}
	 
std::cout<<"range: "<< range << std::endl;
	  for(int i=0; i< geom.points.size(); i++)
	  {
	  		x = geom.points[i][0];
			y = geom.points[i][1];
			z = geom.points[i][2];


	  		value = inputVol.interpolate(x,y,z);
  	
	        c=GetColour(value, range);
			
	        color[0] = c.r;
			color[1] = c.g;
			color[2] = c.b;

			geom.colors.push_back(color);
	  }

	cvcraw_geometry::write(geom, string(argv[3]));

     std::cout<< "done. "<<  std::endl;

    }

  catch(VolMagick::Exception &e)
    {
      std:: cerr << e.what() << std::endl;
    }
  catch(std::exception &e)
    {
      std::cerr << e.what() << std::endl;
    }

  return 0;
}

int main(int argc, char **argv)
{
  if(argc < 6)
    {
      std:: cerr << 
	"Usage: " << argv[0] << 

	
	"  <first volume>  <second volume> < threshold> <output volume left> <output volume masked> \n";
	std::cerr<<"second volume (son) is subset of first volume (mother), threshold is in [0,255], output masked is the volume in first volume with > threshold in the second volume. output volume left is the first volume - output volume masked. " << std::endl;

      return 1;
    }

  try
    {
      VolMagick::Volume inputVol;

      VolMagick::Volume inputVol2;

      
      
      VolMagick::readVolumeFile(inputVol,argv[1]); ///first argument is input volume
      
      VolMagick::readVolumeFile(inputVol2, argv[2]); /// second argument is mask volume
      
	  float threshold = atof(argv[3]);


	  VolMagick::Volume  outputLeft;
	  VolMagick::Volume  outputMasked;

      if(!inputVol2.boundingBox().isWithin(inputVol.boundingBox()) )
	  	throw VolMagick::SubVolumeOutOfBounds("The mask volume bounding box must be within the mother volume's bounding box.");



      outputLeft.voxelType(inputVol.voxelType());
      outputLeft.dimension(inputVol.dimension());
      outputLeft.boundingBox(inputVol.boundingBox());
      
      outputMasked.voxelType(inputVol.voxelType());
      outputMasked.dimension(inputVol.dimension());
      outputMasked.boundingBox(inputVol.boundingBox());
    
   

	  float x, y, z;
	  float fvalue;
	  float minvalue = (float)inputVol2.min();
	  float maxvalue = (float)inputVol2.max();

 	  for( int kz = 0; kz<inputVol.ZDim(); kz++)
	  {
	  	for( int jy = 0; jy<inputVol.YDim(); jy++)
	    {
			for( int ix = 0; ix<inputVol.XDim(); ix++)
		   	{
		     	z = inputVol.ZMin() + kz*inputVol.ZSpan();
			 	y = inputVol.YMin() + jy*inputVol.YSpan();
				x = inputVol.XMin() + ix*inputVol.XSpan();
				if((z >= (float)inputVol2.ZMin()) && (z <= (float)inputVol2.ZMax()) && ( y >=(float)inputVol2.YMin()) && (y <= (float)inputVol2.YMax()) &&(x >= (float)inputVol2.XMin()) && (x <= (float)inputVol2.XMax()) )
				{
				  fvalue = 255.0*(inputVol2.interpolate(x,y,z)-minvalue)/(maxvalue-minvalue);
				  if(fvalue > threshold)
				  {
				  	outputMasked(ix,jy,kz, inputVol(ix,jy,kz));
					outputLeft(ix,jy,kz, inputVol.min());
				  }	  
				  else
				  {
				  	outputMasked(ix,jy,kz, inputVol.min());
					outputLeft(ix,jy,kz, inputVol(ix,jy,kz));
				  }
				}
				else
				{
					outputMasked(ix,jy,kz, inputVol.min());
					outputLeft(ix,jy,kz, inputVol(ix,jy,kz));
				}
		     }		   
		  }	
		}



	VolMagick::createVolumeFile(outputLeft, argv[4]);
	VolMagick::createVolumeFile(outputMasked, argv[5]);

  
    std::cout<<"done!"<<std::endl;


    }

  catch(VolMagick::Exception &e)
    {
      std:: cerr << e.what() << std::endl;
    }
  catch(std::exception &e)
//.........这里部分代码省略.........

int main(int argc, char **argv) {
  // Initialize log4cplus
  std::ifstream testfile("log4cplus.properties");
  if (testfile) {
    testfile.close();
    log4cplus::PropertyConfigurator::doConfigure("log4cplus.properties");
  }
  else {
    log4cplus::BasicConfigurator::doConfigure();
  }

  if(argc < 9) {
    std::cerr << "Usage: " << argv[0] << " <first volume>" <<
	"<xmin> <ymin> <zmin> <xmax> <ymax> <zmax>  <output volume> \n";
    return 1;
  }

  try {
    VolMagick::Volume inputVol;
    VolMagick::readVolumeFile(inputVol,argv[1]); ///first argument is input volume
	
    VolMagick::Volume  outputVol;
	
    VolMagick::VolumeFileInfo volinfo1;
    volinfo1.read(argv[1]);
    std::cout << volinfo1.filename() << ":" <<std::endl;

    std::cout << "minVol1 , maxVol1: " << volinfo1.min() << " " << 
        volinfo1.max() << std::endl;;
 
    float span[3];
    span[0]=inputVol.XSpan(); 
    span[1]=inputVol.YSpan(); 
    span[2]=inputVol.ZSpan();
	
    VolMagick::BoundingBox bbox;
    bbox.minx = atof(argv[2]);
    bbox.miny = atof(argv[3]);
    bbox.minz = atof(argv[4]);
    bbox.maxx = atof(argv[5]);
    bbox.maxy = atof(argv[6]);
    bbox.maxz = atof(argv[7]);
	
    VolMagick::Dimension dim;
    dim.xdim = (int) ((bbox.maxx-bbox.minx)/ span[0])+1;
    dim.ydim = (int) ((bbox.maxy-bbox.miny)/ span[1])+1;
    dim.zdim = (int) ((bbox.maxz-bbox.minz)/ span[2])+1;

    bbox.maxx = bbox.minx + (dim.xdim-1)*span[0];
    bbox.maxy = bbox.miny + (dim.ydim-1)*span[1];
    bbox.maxz = bbox.minz + (dim.zdim-1)*span[2];
 	
    outputVol.voxelType(inputVol.voxelType());
    outputVol.dimension(dim);
    outputVol.boundingBox(bbox);
      
    float temp;
    int i, j, k;
    int xsh = (int)((bbox.minx - inputVol.XMin())/span[0]);
    int ysh = (int)((bbox.miny - inputVol.YMin())/span[1]);
    int zsh = (int)((bbox.minz - inputVol.ZMin())/span[2]);

    for( int kz = 0; kz<outputVol.ZDim(); kz++) {
      for( int jy = 0; jy<outputVol.YDim(); jy++) {
        for( int ix = 0; ix<outputVol.XDim(); ix++) {
          i = ix + xsh;
          j = jy + ysh;
          k = kz + zsh;
          if (i<0 || i >= inputVol.XDim()|| j<0 || j>= inputVol.YDim() || k <0 || k>=inputVol.ZDim()) {
            outputVol(ix, jy, kz, 0.0);
          } else {
            outputVol(ix,jy,kz, inputVol(i,j,k));
          }
        }
      }
    }
    std::cout << "creating volume file..." << std::endl;
    VolMagick::createVolumeFile(outputVol, argv[8]);
    std::cout << "done!" << std::endl;
  } catch(VolMagick::Exception &e) {
    std:: cerr << e.what() << std::endl;
  } catch(std::exception &e) {
    std::cerr << e.what() << std::endl;
  }

  return 0;
}

void GeometryInterface::setViewerState()
{
  using namespace std;
  using namespace boost;

  // static log4cplus::Logger logger = log4cplus::getFuncLogger();
  static log4cplus::Logger logger = FUNCTION_LOGGER;

  //so it responds to volumeInterfaceViewer.* property changes
  _geometryInterfaceViewer->setObjectName("geometryInterfaceViewer");

  CVC_NAMESPACE::PropertyMap properties;
  properties["geometryInterfaceViewer.rendering_mode"] = "colormapped";
  properties["geometryInterfaceViewer.shaded_rendering_enabled"] = "false";
  properties["geometryInterfaceViewer.draw_bounding_box"] = "true";
  properties["geometryInterfaceViewer.draw_subvolume_selector"] = "false";
  // arand, 6-14-2011: changed default below to 0.5
  properties["geometryInterfaceViewer.volume_rendering_quality"] = "0.5"; //[0.0,1.0]
  properties["geometryInterfaceViewer.volume_rendering_near_plane"] = "0.0";
  properties["geometryInterfaceViewer.projection_mode"] = "perspective";
  properties["geometryInterfaceViewer.draw_corner_axis"] = "true";
  properties["geometryInterfaceViewer.draw_geometry"] = "true";
  properties["geometryInterfaceViewer.draw_volumes"] = "true";
  properties["geometryInterfaceViewer.clip_geometry"] = "true";
  properties["geometryInterfaceViewer.volumes"] = "";
  properties["geometryInterfaceViewer.geometry_line_width"] = "1.2";
  properties["geometryInterfaceViewer.background_color"] = "#000000";

  //stereo related properties
  properties["geometryInterfaceViewer.io_distance"] = "0.062";
  properties["geometryInterfaceViewer.physical_distance_to_screen"] = "2.0";
  properties["geometryInterfaceViewer.physical_screen_width"] = "1.8";
  properties["geometryInterfaceViewer.focus_distance"] = "1000.0";

  //viewers_transfer_function should be a boost::shared_array<float>
  //on the data map.
  properties["geometryInterfaceViewer.transfer_function"] = "none";

  cvcapp.addProperties(properties);

  //set up the volume bounding box to encompass the geometry
  vector<cvcraw_geometry::cvcgeom_t> geoms =
    _geometryInterfaceViewer->getGeometriesFromDatamap();
  LOG4CPLUS_TRACE(logger, "geoms.size() == " << geoms.size());
  // cvcapp.log(3,str(format("%s :: geoms.size() == %s\n")
  //                  % BOOST_CURRENT_FUNCTION
  //                  % geoms.size()));

  if(geoms.empty())
    {
      //just use the default bounding box of a volume
      _geometryInterfaceViewer->colorMappedVolume(VolMagick::Volume());
      _geometryInterfaceViewer->rgbaVolumes(vector<VolMagick::Volume>());
    }
  else
    {
      VolMagick::BoundingBox bbox;
      cvcraw_geometry::cvcgeom_t initial_geom;
      
      //find a non empty geom to supply an initial bbox
      BOOST_FOREACH(cvcraw_geometry::cvcgeom_t geom, geoms)
        if(!geom.empty())
          {
            initial_geom = geom;
            break;
          }
      
      cvcraw_geometry::cvcgeom_t::point_t minpt;
      cvcraw_geometry::cvcgeom_t::point_t maxpt;
      if(initial_geom.empty())
        {
          //no non empty geometry so lets just use the default box
          _geometryInterfaceViewer->colorMappedVolume(VolMagick::Volume());
          _geometryInterfaceViewer->rgbaVolumes(vector<VolMagick::Volume>());
        }
      else
        {
          minpt = initial_geom.min_point();
          maxpt = initial_geom.max_point();
          // arand: slightly enlarging
          double eps0 = (maxpt[0]-minpt[0])/20.0;
          double eps1 = (maxpt[1]-minpt[1])/20.0;
          double eps2 = (maxpt[2]-minpt[2])/20.0;
          
          bbox = VolMagick::BoundingBox(minpt[0]-eps0,minpt[1]-eps1,minpt[2]-eps2,
                                        maxpt[0]+eps0,maxpt[1]+eps1,maxpt[2]+eps2);
        }
      
      //build a bounding box that encompasses all bounding boxes of the geometries
      BOOST_FOREACH(cvcraw_geometry::cvcgeom_t geom, geoms)
        {
          if(geom.empty()) continue;
          minpt = geom.min_point();
          maxpt = geom.max_point();
          VolMagick::BoundingBox geobox(minpt[0],minpt[1],minpt[2],
                                        maxpt[0],maxpt[1],maxpt[2]);
          
          cvcapp.log(5,str(boost::format("geobox: (%f,%f,%f) (%f,%f,%f)")
                           % geobox.minx % geobox.miny % geobox.minz
                           % geobox.maxx % geobox.maxy % geobox.maxz));
//.........这里部分代码省略.........

int main(int argc, char **argv)
{
  if(argc != 3)
    {
      cerr << "Usage: " << argv[0] << " <input volume file> <output volume file>" << endl;
      return 1;
    }

  try
    {
#ifndef OUT_OF_CORE
      cerr << "In-core convert" << endl;
      VolMagick::Volume vol;

      //TODO: read/write a slice at a time instead of reading the whole volume in memory then writing it out...
      VolMagick::readVolumeFile(vol,argv[1]/*,var,time*/);
      //VolMagick::writeVolumeFile(vol,argv[2]/*,var,time*/);
      VolMagick::createVolumeFile(vol,argv[2]);
#else
      cerr << "Out-of-core convert" << endl;
      VolMagick::VolumeFileInfo volinfo;
      volinfo.read(argv[1]);

	  //VolMagick::createVolumeFile in Utlity.h 
      VolMagick::createVolumeFile(argv[2],volinfo); 


	 // cout<<"convert volinfo:" << volinfo.boundingBox().minx <<" " << volinfo.boundingBox().maxx<< endl;
      //read in slice by slice
      for(unsigned int k = 0; k < volinfo.ZDim(); k++)
	{
	  for(unsigned int var=0; var<volinfo.numVariables(); var++)
	    for(unsigned int time=0; time<volinfo.numTimesteps(); time++)
	      {
		VolMagick::Volume vol;
		readVolumeFile(vol,argv[1],
			       var,time,
			       0,0,k,
			       VolMagick::Dimension(volinfo.XDim(),volinfo.YDim(),1));

		vol.desc(volinfo.name(var));
		writeVolumeFile(vol,argv[2],
				var,time,
				0,0,k);
	      }
	  fprintf(stderr,"Converting: %5.2f %%\r",(((float)k)/((float)((int)(volinfo.ZDim()-1))))*100.0);
	}
      fprintf(stderr,"\n");
#endif
    }
  catch(VolMagick::Exception &e)
    {
      cerr << e.what() << endl;
    }
  catch(std::exception &e)
    {
      cerr << e.what() << endl;
    }

  return 0;
}

int main(int argc, char **argv)
{
  if(argc < 4)
    {
      std:: cerr << 
	"Usage: " << argv[0] << 

	
	"<input vol 1> <input vol 2> <output vol> \n";

      return 1;
    }

  try
    {
      VolMagick::Volume inputVol;

      VolMagick::Volume inputVol2;
      
      VolMagick::Volume outputVol;

      VolMagick::VolumeFileInfo volinfo;
      volinfo.read(argv[1]);


      VolMagick::readVolumeFile(inputVol,argv[1]); ///first argument is input volume
      
      VolMagick::readVolumeFile(inputVol2, argv[2]); /// second argument is mask volume
      

      VolMagick::VolumeFileInfo volinfo1;
      volinfo1.read(argv[1]);
      std::cout << volinfo1.filename() << ":" <<std::endl;

      VolMagick::VolumeFileInfo volinfo2;
      volinfo2.read(argv[2]);
      std::cout << volinfo2.filename() << ":" <<std::endl;
      
      std::cout<<"minVol1 , maxVol1: "<<volinfo1.min()<<" "<<volinfo1.max()<<std::endl;;

      std::cout<<"minVol2 , maxVol2: "<<volinfo2.min()<<" "<<volinfo2.max()<<std::endl;

           

      outputVol.voxelType(inputVol.voxelType());
      outputVol.dimension(inputVol.dimension());
      outputVol.boundingBox(inputVol.boundingBox());
      
      
      std::cout<<"voxeltype "<<inputVol.voxelType()<<std::endl;
      

      double min=100000000.0;
      double max=-100000000.0;
      
      for( int kz = 0; kz<inputVol.ZDim(); kz++)
	{
	  std::cout<<kz<<"..";
	    for( int jy = 0; jy<inputVol.YDim(); jy++)
	      for( int ix = 0; ix<inputVol.XDim(); ix++)
		{		  
		  outputVol(ix,jy,kz, inputVol(ix,jy,kz)+inputVol2(ix,jy,kz));

		  if (inputVol(ix,jy,kz)+inputVol2(ix,jy,kz) < min) {
		    min = inputVol(ix,jy,kz)+inputVol2(ix,jy,kz);
		  }

		  if (inputVol(ix,jy,kz)+inputVol2(ix,jy,kz) >max) {
		    max = inputVol(ix,jy,kz)+inputVol2(ix,jy,kz);
		  }



		}
	}	
      

      std::cout << std::endl << "New Min/Max " << min << " " << max << std::endl;

      VolMagick::createVolumeFile(argv[3],volinfo);

      VolMagick::writeVolumeFile(outputVol, argv[3]);

    }

  catch(VolMagick::Exception &e)
    {
      std:: cerr << e.what() << std::endl;
    }
  catch(std::exception &e)
    {
      std::cerr << e.what() << std::endl;
    }

  return 0;
}

/*using namespace std;

class VolMagickOpStatus : public VolMagick::VoxelOperationStatusMessenger
{
public:
  void start(const VolMagick::Voxels *vox, Operation op, VolMagick::uint64 numSteps) const
  {
    _numSteps = numSteps;
  }

  void step(const VolMagick::Voxels *vox, Operation op, VolMagick::uint64 curStep) const
  {
    const char *opStrings[] = { "CalculatingMinMax", "CalculatingMin", "CalculatingMax",
				"SubvolumeExtraction", "Fill", "Map", "Resize", "Composite",
				"BilateralFilter", "ContrastEnhancement"};

    fprintf(stderr,"%s: %5.2f %%\r",opStrings[op],(((float)curStep)/((float)((int)(_numSteps-1))))*100.0);
  }

  void end(const VolMagick::Voxels *vox, Operation op) const
  {
    printf("\n");
  }

private:
  mutable VolMagick::uint64 _numSteps;
};

typedef boost::tuple<double, double, double> Color;
*/
int main(int argc, char **argv)
{
  if(argc < 3)
    {
      std:: cerr << 
	"Usage: " << argv[0] << "[vol 1] [vol 2] [out vol] [threshold]";
      std:: cerr <<
	"first volume + second volume = output volume. The voxels are remapped giving two different intensity ranges for the volumes, and threshold will be used for the second volume if it was set \n";

      return 1;
    }

  try
    {
      VolMagick::Volume inputVol;

      VolMagick::Volume inputVol2;
      
      VolMagick::Volume outputVol;

      VolMagick::readVolumeFile(inputVol,argv[1]); ///first argument is input volume
      
      VolMagick::readVolumeFile(inputVol2, argv[2]); /// second argument is mask volume
      

      VolMagick::VolumeFileInfo volinfo1;
      volinfo1.read(argv[1]);
      std::cout << volinfo1.filename() << ":" <<std::endl;

      VolMagick::VolumeFileInfo volinfo2;
      volinfo2.read(argv[2]);
      std::cout << volinfo2.filename() << ":" <<std::endl;
      
      std::cout<<"minVol1 , maxVol1: "<<volinfo1.min()<<" "<<volinfo1.max()<<std::endl;;

      std::cout<<"minVol2 , maxVol2: "<<volinfo2.min()<<" "<<volinfo2.max()<<std::endl;

      
      double thold = 0.0001;
      if( argc > 3) thold = atof( argv[4] );
      if( thold < 0.0001 ) thold = 0.0001;

      outputVol.voxelType(inputVol.voxelType());
      outputVol.dimension(inputVol.dimension());
      outputVol.boundingBox(inputVol.boundingBox());
      
      
      std::cout<<"voxeltype "<<inputVol.voxelType()<<std::endl;
      double min1 = volinfo1.min(); 
      double range1 = volinfo1.max()-volinfo1.min(); 
      double min2 = volinfo2.min(); 
      double range2 = volinfo2.max()-volinfo2.min(); 
      
      for( int kz = 0; kz<inputVol.ZDim(); kz++)
	{
	  std::cout<<kz<<"..";
	    for( int jy = 0; jy<inputVol.YDim(); jy++)
	      for( int ix = 0; ix<inputVol.XDim(); ix++)
		{
		  
		  outputVol(ix,jy,kz, 0);
		  
		  if( (inputVol(ix,jy,kz)- min1 ) > 0.0001 )
		    
		    outputVol(ix,jy,kz, ( (inputVol(ix,jy,kz)-min1)/range1)*128.0 );
		  
		  if( (inputVol2(ix,jy,kz)-min2) > thold)
		    
		    outputVol(ix,jy,kz, ( (inputVol2(ix,jy,kz)-min2)/range2)*128.0 + 128.0);
		  
//.........这里部分代码省略.........

int main(int argc, char **argv)
{
  if(argc < 6)
    {
      std:: cerr << 
	"Usage: " << argv[0] << 

	
	" <first volume>  <xmin> <ymin> <zmin> <output volume> \n";

      return 1;
    }

  try
    {
      VolMagick::Volume inputVol;


      
      
      VolMagick::readVolumeFile(inputVol,argv[1]); ///first argument is input volume
      
	
	  VolMagick::Volume  outputVol;

	
      VolMagick::VolumeFileInfo volinfo1;
      volinfo1.read(argv[1]);
      std::cout << volinfo1.filename() << ":" <<std::endl;

      
      std::cout<<"minVol1 , maxVol1: "<<volinfo1.min()<<" "<<volinfo1.max()<<std::endl;;
 
	  float span[3];
	  span[0]=inputVol.XSpan(); 
	  span[1]=inputVol.YSpan(); 
	  span[2]=inputVol.ZSpan();
	
	  VolMagick::BoundingBox bbox;
	  bbox.minx = atof(argv[2]);
	  bbox.miny = atof(argv[3]);
	  bbox.minz = atof(argv[4]);
	  bbox.maxx = bbox.minx + volinfo1.XMax()-volinfo1.XMin();
	  bbox.maxy = bbox.miny + volinfo1.YMax()-volinfo1.YMin();
	  bbox.maxz = bbox.minz + volinfo1.ZMax()-volinfo1.ZMin();
	
/*	  VolMagick::Dimension dim;
	  dim.xdim = (int) ((bbox.maxx-bbox.minx)/ span[0])+1;
	  dim.ydim = (int) ((bbox.maxy-bbox.miny)/ span[1])+1;
	  dim.zdim = (int) ((bbox.maxz-bbox.minz)/ span[2])+1;
*/

 	
      outputVol.voxelType(inputVol.voxelType());
      outputVol.dimension(volinfo1.dimension());
      outputVol.boundingBox(bbox);
      
      
/*	 float temp;
	 int i, j, k;
     int xsh = (int)((bbox.minx - inputVol.XMin())/span[0]);
     int ysh = (int)((bbox.miny - inputVol.YMin())/span[1]);
     int zsh = (int)((bbox.minz - inputVol.ZMin())/span[2]);
*/	 

	  for( int kz = 0; kz<outputVol.ZDim(); kz++)
	   for( int jy = 0; jy<outputVol.YDim(); jy++)
	     for( int ix = 0; ix<outputVol.XDim(); ix++)
		    {
		//		i = ix + xsh;
		//		j = jy + ysh;
		//		k = kz + zsh;
		//		cout <<"ijk "<< i <<" " << j <<" " <<k <<endl;
		//		if(i<0 || i >= inputVol.XDim()|| j<0 || j>= inputVol.YDim() || k <0 || k>=inputVol.ZDim()) outputVol(ix, jy, kz, 0.0);
		//		else
				outputVol(ix,jy,kz, inputVol(ix,jy,kz)  );
			}
	//		stringstream s;
	//		s<<argv[3]<< num <<".rawiv";
	      VolMagick::createVolumeFile(outputVol, argv[5]);
    

    std::cout<<"done!"<<std::endl;


    }

  catch(VolMagick::Exception &e)
    {
      std:: cerr << e.what() << std::endl;
    }
  catch(std::exception &e)
    {
      std::cerr << e.what() << std::endl;
    }

  return 0;
}

int main(int argc, char **argv)
{
  if(argc < 3)
    {
      std:: cerr << 
	"Usage: " << argv[0] << 

	
	"  <first volume>  <output volume> [ bool 0/1 ].   \n"; 
	cerr<<"bool 1 keeps the original volume, bool 0 only the reflection, default 0."  << endl;

      return 1;
    }

  try
    {
      VolMagick::Volume inputVol;

      
      VolMagick::Volume outputVol;

      VolMagick::readVolumeFile(inputVol,argv[1]); ///first argument is input volume
      
      VolMagick::VolumeFileInfo volinfo1;
      volinfo1.read(argv[1]);
      std::cout << volinfo1.filename() << ":" <<std::endl;

       
      std::cout<<"minVol1 , maxVol1: "<<volinfo1.min()<<" "<<volinfo1.max()<<std::endl;;

  
      VolMagick::BoundingBox bbox;
	  bbox.minx = inputVol.XMin();
	  bbox.maxx = inputVol.XMax();
	  bbox.miny = inputVol.YMin();
	  bbox.maxy = inputVol.YMax();
	  bbox.minz = inputVol.ZMin();
	  bbox.maxz = inputVol.ZMax();

	  VolMagick::Dimension dim;
	  dim.xdim = inputVol.XDim();
	  dim.ydim = inputVol.YDim();
	  dim.zdim = inputVol.ZDim();


//	 cout<<bbox.minz <<" " << bbox.maxz<<" "<< bbox.maxy <<endl;
//	 cout<<dim.zdim <<" " << dim.ydim << endl;
      
      outputVol.voxelType(inputVol.voxelType());
      outputVol.dimension(dim);
      outputVol.boundingBox(bbox);

	
      
   
      //Works for GroEL4.2
      for( int kz = 0; kz<outputVol.ZDim(); kz++)
	   for( int jy = 0; jy<inputVol.YDim(); jy++)
	      for( int ix = 0; ix<inputVol.XDim(); ix++)
		       {
			        float temp=0.0;
					if(inputVol.ZDim() >= kz && inputVol.ZDim() - kz < inputVol.ZDim()
					&& inputVol.YDim() >= jy && inputVol.YDim() - jy < inputVol.YDim() )
					temp = inputVol(ix, inputVol.YDim()-jy, inputVol.ZDim()-kz);
	  				outputVol(ix, jy, kz, temp);
			   }	
	 
	 if(argc==4 && atoi(argv[3])== 1)
	 {
      for( int kz = 0; kz<inputVol.ZDim(); kz++)
	   for( int jy = 0; jy<inputVol.YDim(); jy++)
	      for( int ix = 0; ix<inputVol.XDim(); ix++)
		       {
			   		outputVol(ix, jy, kz, outputVol(ix,jy,kz)+inputVol(ix,jy,kz));
			   }		  
	 }
		  
			  
		 
	      

      VolMagick::createVolumeFile(outputVol, argv[2]);


      cout<<"done!"<<endl;

    }

  catch(VolMagick::Exception &e)
    {
      std:: cerr << e.what() << std::endl;
    }
  catch(std::exception &e)
    {
      std::cerr << e.what() << std::endl;
    }

  return 0;
}