changeset 2558:274a2a6b7cb7

Added minc4itk examples
author Vladimir S. FONOV <vladimir.fonov@gmail.com>
date Thu, 08 Dec 2011 19:16:52 -0500
parents 0d08a9897649
children 1ad066b39b7d
files minc4itk/examples/itk_convert.cpp minc4itk/examples/itk_distance.cpp minc4itk/examples/itk_dti.cpp minc4itk/examples/itk_resample.cpp minc4itk/examples/volume_2_csv.cpp
diffstat 5 files changed, 989 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
new file mode 100644
--- /dev/null
+++ b/minc4itk/examples/itk_convert.cpp
@@ -0,0 +1,287 @@
+#include <iostream>
+
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <itkImageIOFactory.h>
+#include <itkImageIOBase.h>
+#include <itkFlipImageFilter.h>
+
+#include "itkMincImageIOFactory.h"
+#include "itkMincImageIO.h"
+
+#include <getopt.h>
+
+void show_usage (const char *name)
+{
+  std::cerr 
+      << "Usage: "<<name<<" <input> <output> " << std::endl
+      << "--clobber clobber output files"<<std::endl
+      << "--verbose be verbose"<<std::endl
+      << "--inv-x invert X axis"<<std::endl
+      << "--inv-y invert Y axis"<<std::endl
+      << "--inv-z invert Z axis"<<std::endl
+      << "--center set origin to the center of the image"<<std::endl
+;
+}
+typedef itk::ImageIOBase IOBase;
+typedef itk::SmartPointer<IOBase> IOBasePointer;
+
+template<class ImageType> void load_and_save_image(IOBase* base,
+                    const char *fname,
+                    bool inv_x=false,
+                    bool inv_y=false,
+                    bool inv_z=false,
+                    bool center=false,
+                    bool verbose=false)
+{
+  typename itk::ImageFileReader<ImageType >::Pointer reader = itk::ImageFileReader<ImageType >::New();
+  typename itk::FlipImageFilter<ImageType >::Pointer flip=itk::FlipImageFilter<ImageType >::New();
+  reader->SetImageIO(base);
+  reader->SetFileName(base->GetFileName());
+  reader->Update();
+  
+  typename ImageType::Pointer img=reader->GetOutput();
+     
+  /* WRITING */
+  if(verbose)
+    std::cout<<"Writing "<<fname<<"..."<<std::endl;
+  
+
+  
+  if(inv_x||inv_y||inv_z)
+  {
+    typename itk::FlipImageFilter<ImageType >::FlipAxesArrayType arr;
+    arr[0]=inv_x;
+    arr[1]=inv_y;
+    arr[2]=inv_z;
+    flip->SetFlipAxes(arr);
+    flip->SetInput(img);
+    flip->Update();
+    img=flip->GetOutput();
+  }
+  
+  if(center)//move origin to the center of the image
+  {
+    typename ImageType::RegionType r=img->GetLargestPossibleRegion();
+    std::vector<double> corner[3];
+    
+    typename ImageType::IndexType idx;
+    typename ImageType::PointType c;
+    
+    idx[0]=r.GetIndex()[0]+r.GetSize()[0]/2.0;
+    idx[1]=r.GetIndex()[1]+r.GetSize()[1]/2.0;
+    idx[2]=r.GetIndex()[2]+r.GetSize()[2]/2.0;
+    
+    img->TransformIndexToPhysicalPoint(idx,c);
+    
+    typename ImageType::PointType org=img->GetOrigin();
+    
+    org[0]-=c[0];
+    org[1]-=c[1];
+    org[2]-=c[2];
+    
+    img->SetOrigin(org);
+  }  
+  
+  typename itk::ImageFileWriter< ImageType >::Pointer writer = itk::ImageFileWriter<ImageType >::New();
+  writer->SetFileName(fname);
+  writer->SetInput( img );
+  writer->Update();
+  
+}
+
+int main(int argc,char **argv)
+{
+  int verbose=0;
+  int clobber=0;
+  int inv_x=0,inv_y=0,inv_z=0,center=0;
+  //char *history = time_stamp(argc, argv); //maybe we should free it afterwards
+
+  static struct option long_options[] = { 
+    {"verbose", no_argument, &verbose, 1},
+    {"quiet", no_argument, &verbose, 0},
+    {"clobber", no_argument, &clobber, 1},
+    {"inv-x", no_argument,  &inv_x, 1},
+    {"inv-y", no_argument,  &inv_y, 1},
+    {"inv-z", no_argument,  &inv_z, 1},
+    {"center", no_argument, &center, 1},
+    {0, 0, 0, 0}
+  };
+    
+  int c;
+  for (;;)
+  {
+    /* getopt_long stores the option index here. */
+    int option_index = 0;
+
+    c = getopt_long (argc, argv, "", long_options, &option_index);
+
+    /* Detect the end of the options. */
+    if (c == -1)
+      break;
+
+    switch (c)
+    {
+      case 0:
+        break;
+      case '?':
+        /* getopt_long already printed an error message. */
+      default:
+        show_usage (argv[0]);
+        return 1;
+    }
+  }
+  
+  if((argc - optind)<2)
+  {
+    show_usage(argv[0]);
+    return 1;
+  }
+  std::string input=argv[optind];
+  std::string output=argv[optind+1];
+  
+  if (!clobber && !access (output.c_str(), F_OK))
+  {
+    std::cerr << output.c_str () << " Exists!" << std::endl;
+    return 1;
+  }
+ 
+  try
+  {
+    //registering the MINC_IO factory
+    itk::ObjectFactoryBase::RegisterFactory(itk::MincImageIOFactory::New());
+    /* READING */
+    if(verbose) 
+      std::cout<<"Reading "<<input.c_str()<<"..."<<std::endl;
+    
+    //try to figure out what we have got
+    IOBasePointer io = itk::ImageIOFactory::CreateImageIO(input.c_str(), itk::ImageIOFactory::ReadMode );
+    
+    if(!io)
+      throw itk::ExceptionObject("Unsupported image file type");
+    
+    io->SetFileName(input.c_str());
+    io->ReadImageInformation();
+
+    size_t nd = io->GetNumberOfDimensions();
+    size_t nc = io->GetNumberOfComponents();
+    itk::ImageIOBase::IOComponentType  ct = io->GetComponentType();
+    std::string ct_s = io->GetComponentTypeAsString(ct);
+    
+    if(verbose)
+    {
+      std::cout<<"dimensions:"<<nd<<std::endl
+          <<"components:"<<nc<<std::endl
+          <<"type:"<<ct_s.c_str()<<std::endl;
+    }
+    
+    if(nd==3 && nc==1)
+    {
+      if(verbose) std::cout<<"Writing 3D image..."<<std::endl;
+        switch(ct)
+        {
+          case itk::ImageIOBase::UCHAR :
+            load_and_save_image<itk::Image<unsigned char, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::CHAR :
+            load_and_save_image<itk::Image<char, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::USHORT :
+            load_and_save_image<itk::Image<unsigned short, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::SHORT :
+            load_and_save_image<itk::Image<short, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::INT :
+            load_and_save_image<itk::Image<int, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          case itk::ImageIOBase::UINT:
+            load_and_save_image<itk::Image<unsigned int, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          case itk::ImageIOBase::FLOAT :
+            load_and_save_image<itk::Image<float, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          case itk::ImageIOBase::DOUBLE:
+            load_and_save_image<itk::Image<double, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          default:
+            itk::ExceptionObject("Unsupported component type");
+        }
+    } else if((nd==4 && nc==1)) {
+      if(verbose) std::cout<<"Writing 4D image..."<<std::endl;
+        switch(ct)
+        {
+          case itk::ImageIOBase::UCHAR:
+            load_and_save_image<itk::Image<unsigned char, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::CHAR:
+            load_and_save_image<itk::Image<char, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::USHORT:
+            load_and_save_image<itk::Image<unsigned short, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::SHORT:
+            load_and_save_image<itk::Image<short, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::INT:
+            load_and_save_image<itk::Image<int, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          case itk::ImageIOBase::UINT:
+            load_and_save_image<itk::Image<unsigned int, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          case itk::ImageIOBase::FLOAT:
+            load_and_save_image<itk::Image<float, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          case itk::ImageIOBase::DOUBLE:
+            load_and_save_image<itk::Image<double, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          default:
+            itk::ExceptionObject("Unsupported component type");
+        }
+    } else if((nd==3 && nc>1)) {
+        if(verbose) std::cout<<"Writing multicomponent 3D image..."<<std::endl;
+
+        switch(ct)
+        {
+          case itk::ImageIOBase::UCHAR:
+            load_and_save_image<itk::VectorImage<unsigned char, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::CHAR:
+            load_and_save_image<itk::VectorImage<char, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::USHORT:
+            load_and_save_image<itk::VectorImage<unsigned short, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::SHORT:
+            load_and_save_image<itk::VectorImage<short, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break;
+          case itk::ImageIOBase::INT:
+            load_and_save_image<itk::VectorImage<int, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          case itk::ImageIOBase::UINT:
+            load_and_save_image<itk::VectorImage<unsigned int, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          case itk::ImageIOBase::FLOAT:
+            load_and_save_image<itk::VectorImage<float, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          case itk::ImageIOBase::DOUBLE:
+            load_and_save_image<itk::VectorImage<double, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose);
+            break; 
+          default:
+            itk::ExceptionObject("Unsupported component type");
+        }
+    } else {
+        throw itk::ExceptionObject("Unsupported number of dimensions");
+    }
+    
+  }
+  
+  catch( itk::ExceptionObject & err )
+  {
+    std::cerr << "ExceptionObject caught !" << std::endl;
+    std::cerr << err << std::endl;
+    return 2;
+  }
+  return 0;
+}
new file mode 100644
--- /dev/null
+++ b/minc4itk/examples/itk_distance.cpp
@@ -0,0 +1,139 @@
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       : 
+@DESCRIPTION: 
+@COPYRIGHT  :
+              Copyright 2006 Vladimir Fonov, McConnell Brain Imaging Centre, 
+              Montreal Neurological Institute, McGill University.
+              Permission to use, copy, modify, and distribute this
+              software and its documentation for any purpose and without
+              fee is hereby granted, provided that the above copyright
+              notice appear in all copies.  The author and McGill University
+              make no representations about the suitability of this
+              software for any purpose.  It is provided "as is" without
+              express or implied warranty.
+---------------------------------------------------------------------------- */
+#include <unistd.h>
+#include <getopt.h>
+#include <iostream>
+#include <itkDanielssonDistanceMapImageFilter.h>
+#include <itkSignedDanielssonDistanceMapImageFilter.h>
+#include <time_stamp.h>    // for creating minc style history entry
+#include "itkMincImageIOFactory.h"
+#include "itkMincImageIO.h"
+#include "minc_helpers.h"
+
+//#include <minc_wrappers.h>
+
+using namespace minc;
+using namespace std;
+void show_usage(const char *name)
+{
+  std::cerr 
+    << "Usage: "<<name<<" <input> <output> " << endl
+    << "--verbose be verbose "    << endl
+    << "--clobber clobber output files"<<endl
+    << "--signed produce signed distance map"<<endl;
+}
+
+int main (int argc, char **argv)
+{
+  
+  int verbose=1;
+  double sigma=0.5;
+  double keep=1.0;
+  int order=5;
+  int approx=0;
+  int ss=0;
+  int clobber=0;
+  char *history = time_stamp(argc, argv); 
+
+  //int voxel_neibourhood=5;
+  static struct option long_options[] = { 
+    {"clobber", no_argument, &clobber, 1},
+    {"verbose", no_argument, &verbose, 1},
+    {"quiet", no_argument, &verbose, 0},
+    {"signed",no_argument, &ss, 1},
+    {0, 0, 0, 0}
+  };
+  
+  for (;;) {
+    /* getopt_long stores the option index here. */
+    int option_index = 0;
+
+    int c = getopt_long (argc, argv, "vq", long_options, &option_index);
+
+    /* Detect the end of the options. */
+    if (c == -1) break;
+
+    switch (c)
+    {
+    //  case 'n':
+    //    voxel_neibourhood=atoi(optarg);break;
+      case 0:
+        break;
+      case '?':
+        /* getopt_long already printed an error message. */
+      default:
+        show_usage (argv[0]);
+        return 1;
+    }
+  }
+  if ((argc - optind) < 2) {
+    show_usage (argv[0]);
+    return 1;
+  }
+  std::string input_f=argv[optind],  out_f=argv[optind+1];
+  
+  // check if the file already present
+  if (!clobber && !access (out_f.c_str (), F_OK))
+  {
+    cerr << out_f.c_str () << " Exists!" << endl;
+    return 1;
+  }
+    
+  try
+  {
+    itk::ObjectFactoryBase::RegisterFactory(itk::MincImageIOFactory::New());
+    itk::ImageFileReader<minc::mask3d >::Pointer reader = itk::ImageFileReader<minc::mask3d >::New();
+    
+    //initializing the reader
+    reader->SetFileName(input_f.c_str());
+    reader->Update();
+    
+    minc::mask3d::Pointer input=reader->GetOutput();
+    minc::image3d::Pointer output;
+    
+    typedef itk::DanielssonDistanceMapImageFilter< minc::mask3d, minc::image3d >
+            DistanceMapFilter;
+    typedef itk::SignedDanielssonDistanceMapImageFilter< minc::mask3d, minc::image3d >
+        SignedDistanceMapFilter;
+    
+    if(ss)
+    {
+      SignedDistanceMapFilter::Pointer dist(SignedDistanceMapFilter::New());
+      dist->SetInput(input);
+      dist->Update();
+      output=dist->GetOutput();
+    } else {
+      DistanceMapFilter::Pointer dist(DistanceMapFilter::New());
+      dist->SetInput(input);
+      dist->Update();
+      output=dist->GetOutput();
+    }
+    
+    minc::copy_metadata(output,input);
+    minc::append_history(output,history);
+    free(history);
+    
+    itk::ImageFileWriter< minc::image3d >::Pointer writer = itk::ImageFileWriter<minc::image3d >::New();
+    writer->SetFileName(out_f.c_str());
+    writer->SetInput( output );
+    writer->Update();
+    
+  } catch (const minc::generic_error & err) {
+    cerr << "Got an error at:" << err.file () << ":" << err.line () << endl;
+    cerr << err.msg()<<endl;
+    return 1;
+  }
+  return 0;
+}
new file mode 100644
--- /dev/null
+++ b/minc4itk/examples/itk_dti.cpp
@@ -0,0 +1,148 @@
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       :  itk_dti
+@DESCRIPTION:  an example of processing DTI information
+@COPYRIGHT  :
+              Copyright 2011 Vladimir Fonov, McConnell Brain Imaging Centre, 
+              Montreal Neurological Institute, McGill University.
+              Permission to use, copy, modify, and distribute this
+              software and its documentation for any purpose and without
+              fee is hereby granted, provided that the above copyright
+              notice appear in all copies.  The author and McGill University
+              make no representations about the suitability of this
+              software for any purpose.  It is provided "as is" without
+              express or implied warranty.
+---------------------------------------------------------------------------- */
+#include <iostream>
+
+#include <itkMetaDataObject.h>
+#include <itkVectorImage.h>
+
+#include <itkImage.h>
+
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <itkImageIOFactory.h>
+#include <itkDiffusionTensor3DReconstructionImageFilter.h>
+#include <itkTensorFractionalAnisotropyImageFilter.h>
+
+#include "itkMincImageIOFactory.h"
+#include "itkMincImageIO.h"
+
+typedef itk::VectorImage<float, 3>  DTIImageType;
+typedef itk::DiffusionTensor3DReconstructionImageFilter<float,float,float> DtiReconstructionFilter;
+typedef itk::Image< itk::DiffusionTensor3D< float >, 3 > TensorImage;
+typedef DtiReconstructionFilter::GradientDirectionContainerType Gradients;
+typedef itk::Image<float,3> faImage;
+typedef itk::TensorFractionalAnisotropyImageFilter<TensorImage,faImage > FAFilter;
+
+
+// a helper function for minc reading
+template <class T> typename T::Pointer load_minc(const char *file)
+{
+  typedef itk::MincImageIO ImageIOType;
+  ImageIOType::Pointer minc2ImageIO = ImageIOType::New();
+     
+  typename itk::ImageFileReader<T>::Pointer reader = itk::ImageFileReader<T>::New();
+    
+  reader->SetFileName(file);
+  reader->SetImageIO( minc2ImageIO );
+  reader->Update();
+    
+  return reader->GetOutput();
+}
+
+// a helper function for minc writing
+template <class T> void save_minc(const char *file,typename T::ConstPointer img)
+{
+  typedef itk::MincImageIO ImageIOType;
+  ImageIOType::Pointer minc2ImageIO = ImageIOType::New();
+     
+  typename itk::ImageFileWriter< T >::Pointer writer = itk::ImageFileWriter<T>::New();
+  writer->SetFileName(file);
+  writer->SetImageIO( minc2ImageIO );
+  writer->SetInput( img );
+  writer->Update();
+}
+
+
+int main(int argc,char **argv)
+{
+  if(argc<3)
+  {
+    std::cerr<<"Usage:"<<argv[0]<<" <in.mnc> <out_tensor.mnc> <out_eign.mnc> "<<std::endl;
+    return 1;
+  }
+  
+  try
+  {
+    std::cout<<"Reading "<<argv[1]<<"..."<<std::endl;
+    
+    typedef itk::MincImageIO ImageIOType;
+    ImageIOType::Pointer minc2ImageIO = ImageIOType::New();
+    DTIImageType::Pointer img=load_minc<DTIImageType>(argv[1]);
+    DtiReconstructionFilter::Pointer filter=DtiReconstructionFilter::New();
+    filter->SetNumberOfThreads(1); //as per request in documentation
+    
+    //extracting the parameters of acquisition from the metadata
+    typedef std::vector<double> double_vector;
+    double_vector bvalues,direction_x,direction_y,direction_z;
+
+    //making sure that all vcrtors contain the same number of parameters (just in case)
+    if(!itk::ExposeMetaData<double_vector>( img->GetMetaDataDictionary() , "acquisition:bvalues",bvalues) ||
+        !itk::ExposeMetaData<double_vector>( img->GetMetaDataDictionary() , "acquisition:direction_x",direction_x) ||
+        !itk::ExposeMetaData<double_vector>( img->GetMetaDataDictionary() , "acquisition:direction_y",direction_y) ||
+        !itk::ExposeMetaData<double_vector>( img->GetMetaDataDictionary() , "acquisition:direction_z",direction_z))
+    {
+      std::cerr<<"Image doesn't have information on DTI gradients, can't process!"<<std::endl;
+      return 2;
+    }
+    if(bvalues.size()!=direction_x.size() || 
+       bvalues.size()!=direction_y.size() ||
+       bvalues.size()!=direction_z.size() )
+    {
+      std::cerr<<"Different number of components of gradients"<<std::endl;
+        return 2;
+    }
+    std::cout<<"Found "<<bvalues.size()<<" gradient directions"<<std::endl;
+    
+    //converting metadata representation to the format used by DiffusionTensor3DReconstructionImageFilter
+    Gradients::Pointer gradients=Gradients::New();
+    gradients->resize(direction_x.size());
+    double bval=0;
+    //copying values one-by-one
+    for(int i=0;i<bvalues.size();i++)
+      if(bval<bvalues[i]) bval=bvalues[i];
+    
+    for(int i=0;i<direction_x.size();i++)
+    {
+      (*gradients)[i][0]=direction_x[i]*sqrt(bvalues[i]/bval);
+      (*gradients)[i][1]=direction_y[i]*sqrt(bvalues[i]/bval);
+      (*gradients)[i][2]=direction_z[i]*sqrt(bvalues[i]/bval);
+    }
+    
+    std::cout<<"Calculating tensor..."<<std::endl;
+    filter->SetGradientImage(gradients,img);
+    filter->SetBValue(bval);
+    filter->Update();
+    
+    std::cout<<"Writing "<<argv[2]<<"..."<<std::endl;
+    save_minc<TensorImage>(argv[2],filter->GetOutput() );
+    std::cout<<"Calculating FA ..."<<std::endl;
+    
+    FAFilter::Pointer fa=FAFilter::New();
+    fa->SetInput(filter->GetOutput());
+    fa->SetNumberOfThreads(1); //just in case
+    fa->Update();
+    
+    std::cout<<"Writing "<<argv[3]<<"..."<<std::endl;
+    save_minc<faImage>(argv[3],fa->GetOutput() );
+    
+  }
+  catch( itk::ExceptionObject & err )
+  {
+    std::cerr << "ExceptionObject caught !" << std::endl;
+    std::cerr << err << std::endl;
+    return 2;
+  }
+  return 0;
+}
new file mode 100644
--- /dev/null
+++ b/minc4itk/examples/itk_resample.cpp
@@ -0,0 +1,259 @@
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       :  itk_resample
+@DESCRIPTION:  an example of using spline itnerpolation with MINC xfm transform
+@COPYRIGHT  :
+              Copyright 2011 Vladimir Fonov, McConnell Brain Imaging Centre, 
+              Montreal Neurological Institute, McGill University.
+              Permission to use, copy, modify, and distribute this
+              software and its documentation for any purpose and without
+              fee is hereby granted, provided that the above copyright
+              notice appear in all copies.  The author and McGill University
+              make no representations about the suitability of this
+              software for any purpose.  It is provided "as is" without
+              express or implied warranty.
+---------------------------------------------------------------------------- */
+#include <iostream>
+#include <fstream>
+#include <getopt.h>
+#include <vector>
+#include <valarray>
+#include <math.h>
+#include <limits>
+#include <unistd.h>
+
+#include <itkResampleImageFilter.h>
+#include <itkAffineTransform.h>
+#include <itkNearestNeighborInterpolateImageFunction.h>
+#include <itkBSplineInterpolateImageFunction.h>
+#include <minc_general_transform.h>
+
+#include <unistd.h>
+#include <getopt.h>
+#include <time_stamp.h>    // for creating minc style history entry
+
+
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <itkImageIOFactory.h>
+
+#include "itkMincImageIOFactory.h"
+#include "itkMincImageIO.h"
+#include "minc_helpers.h"
+
+//typedef itk::MincImageIO ImageIOType;
+typedef itk::BSplineInterpolateImageFunction< minc::image3d, double, double >  InterpolatorType;
+typedef itk::ResampleImageFilter<minc::image3d, minc::image3d> FilterType;
+typedef minc::XfmTransform<double,3,3>  TransformType;
+
+using namespace  std;
+
+void show_usage (const char * prog)
+{
+  std::cerr 
+    << "Usage: "<<prog<<" <input> <output.mnc> " << std::endl
+    << "--clobber overwrite files"    << std::endl
+    << "--like <example> (default behaviour analogous to use_input_sampling)"<<std::endl
+    << "--transform <xfm_transform> "<<std::endl
+    << "--order <n> spline order, default 2 "<<std::endl
+    << "--uniformize <step> - will make a volume with uniform step size and no direction cosines" << std::endl
+    << "--invert_transform  - apply inverted transform"<<std::endl;
+}
+
+void generate_uniform_sampling(FilterType* flt, const minc::image3d* img,double step)
+{
+  //obtain physical coordinats of all image corners
+  minc::image3d::RegionType r=img->GetLargestPossibleRegion();
+  std::vector<double> corner[3];
+  for(int i=0;i<8;i++)
+  {
+    minc::image3d::IndexType idx;
+    minc::image3d::PointType c;
+    idx[0]=r.GetIndex()[0]+r.GetSize()[0]*(i%2);
+    idx[1]=r.GetIndex()[1]+r.GetSize()[1]*((i/2)%2);
+    idx[2]=r.GetIndex()[2]+r.GetSize()[2]*((i/4)%2);
+    img->TransformIndexToPhysicalPoint(idx,c);
+    for(int j=0;j<3;j++)
+      corner[j].push_back(c[j]);
+  }
+  minc::image3d::IndexType start;
+  FilterType::SizeType size;
+  FilterType::OriginPointType org;
+  minc::image3d::SpacingType spc;
+  spc.Fill(step);
+  for(int j=0;j<3;j++)
+  {
+    std::sort(corner[j].begin(),corner[j].end());
+    size[j]=ceil((corner[j][7]-corner[j][0])/step);
+    org[j]=corner[j][0];
+  }
+  minc::image3d::DirectionType identity;
+  identity.SetIdentity();
+  flt->SetOutputDirection(identity);
+  start.Fill(0);
+  flt->SetOutputStartIndex(start);
+  flt->SetSize(size);
+  flt->SetOutputOrigin(org);
+  flt->SetOutputSpacing(spc);
+}
+
+int main (int argc, char **argv)
+{
+  int verbose=0, clobber=0,skip_grid=0;
+  int order=2;
+  std::string like_f,xfm_f,output_f,input_f;
+  double uniformize=0.0;
+  int invert=0;
+  char *history = time_stamp(argc, argv); 
+  
+  static struct option long_options[] = {
+		{"verbose", no_argument,       &verbose, 1},
+		{"quiet",   no_argument,       &verbose, 0},
+		{"clobber", no_argument,       &clobber, 1},
+		{"like",    required_argument, 0, 'l'},
+		{"transform",    required_argument, 0, 't'},
+    {"order",    required_argument, 0, 'o'},
+    {"uniformize",    required_argument, 0, 'u'},
+    {"invert_transform", no_argument, &invert, 1},
+		{0, 0, 0, 0}
+		};
+  
+  for (;;) {
+      /* getopt_long stores the option index here. */
+      int option_index = 0;
+
+      int c = getopt_long (argc, argv, "vqcl:t:o:u:", long_options, &option_index);
+
+      /* Detect the end of the options. */
+      if (c == -1) break;
+
+      switch (c)
+			{
+			case 0:
+				break;
+			case 'v':
+				cout << "Version: 0.1" << endl;
+				return 0;
+      case 'l':
+        like_f=optarg; break;
+      case 't':
+        xfm_f=optarg; break;
+      case 'o':
+        order=atoi(optarg);break;
+      case 'u':
+        uniformize=atof(optarg);break;
+			case '?':
+				/* getopt_long already printed an error message. */
+			default:
+				show_usage (argv[0]);
+				return 1;
+			}
+    }
+
+	if ((argc - optind) < 2) {
+		show_usage(argv[0]);
+		return 1;
+	}
+  input_f=argv[optind];
+  output_f=argv[optind+1];
+  
+  if (!clobber && !access (output_f.c_str (), F_OK))
+  {
+    std::cerr << output_f.c_str () << " Exists!" << std::endl;
+    return 1;
+  }
+  
+	try
+  {
+    itk::ObjectFactoryBase::RegisterFactory(itk::MincImageIOFactory::New());
+    itk::ImageFileReader<minc::image3d >::Pointer reader = itk::ImageFileReader<minc::image3d >::New();
+    
+    //initializing the reader
+    reader->SetFileName(input_f.c_str());
+    reader->Update();
+    
+		minc::image3d::Pointer in=reader->GetOutput();
+
+		FilterType::Pointer filter  = FilterType::New();
+    
+    //creating coordinate transformation objects
+		TransformType::Pointer transform = TransformType::New();
+    if(!xfm_f.empty())
+    {
+      //reading a minc style xfm file
+      transform->OpenXfm(xfm_f.c_str());
+      if(!invert) transform->Invert(); //should be inverted by default to walk through target space
+      filter->SetTransform( transform );
+    }
+
+    //creating the interpolator
+		InterpolatorType::Pointer interpolator = InterpolatorType::New();
+		interpolator->SetSplineOrder(order);
+		filter->SetInterpolator( interpolator );
+		filter->SetDefaultPixelValue( 0 );
+    
+    //this is for processing using batch system
+    filter->SetNumberOfThreads(1);
+    
+    minc::image3d::Pointer like=0;
+    if(!like_f.empty())
+    {
+      itk::ImageFileReader<minc::image3d >::Pointer reader = itk::ImageFileReader<minc::image3d >::New();
+      reader->SetFileName(like_f.c_str());
+      reader->Update();
+      if(uniformize!=0.0)
+      {
+        generate_uniform_sampling(filter,reader->GetOutput(),uniformize);
+      } else {
+        filter->SetOutputParametersFromImage(reader->GetOutput());
+        filter->SetOutputDirection(reader->GetOutput()->GetDirection());
+      }
+      like=reader->GetOutput();
+      like->DisconnectPipeline();
+    }
+    else
+    {
+      if(uniformize!=0.0)
+      {
+        generate_uniform_sampling(filter,in,uniformize);
+      } else {
+        //we are using original sampling
+        filter->SetOutputParametersFromImage(in);
+        filter->SetOutputDirection(in->GetDirection());
+      }
+    }
+    
+		filter->SetInput(in);
+    filter->Update();
+    //copy the metadate information, for some reason it is not preserved
+    //filter->GetOutput()->SetMetaDataDictionary(reader->GetOutput()->GetMetaDataDictionary());
+    minc::image3d::Pointer out=filter->GetOutput();
+    minc::copy_metadata(out,in);
+    minc::append_history(out,history);
+    free(history);
+    
+    //correct dimension order
+    if(like.IsNotNull())
+      minc::copy_dimorder(out,like);
+    
+    //generic file writer
+    itk::ImageFileWriter< minc::image3d >::Pointer writer = itk::ImageFileWriter<minc::image3d >::New();
+    writer->SetFileName(output_f.c_str());
+     
+    writer->SetInput( out );
+    //writer->UseInputMetaDataDictionaryOn();
+    
+    writer->Update();
+    
+		return 0;
+	} catch (const minc::generic_error & err) {
+    cerr << "Got an error at:" << err.file () << ":" << err.line () << endl;
+    return 1;
+  }
+  catch( itk::ExceptionObject & err )
+  {
+    std::cerr << "ExceptionObject caught !" << std::endl;
+    std::cerr << err << std::endl;
+    return 2;
+  }
+	return 0;
+};
new file mode 100644
--- /dev/null
+++ b/minc4itk/examples/volume_2_csv.cpp
@@ -0,0 +1,156 @@
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       :  volume_2_csv
+@DESCRIPTION:  an example of converting a minc file volume to a text format
+@COPYRIGHT  :
+              Copyright 2011 Vladimir Fonov, McConnell Brain Imaging Centre, 
+              Montreal Neurological Institute, McGill University.
+              Permission to use, copy, modify, and distribute this
+              software and its documentation for any purpose and without
+              fee is hereby granted, provided that the above copyright
+              notice appear in all copies.  The author and McGill University
+              make no representations about the suitability of this
+              software for any purpose.  It is provided "as is" without
+              express or implied warranty.
+---------------------------------------------------------------------------- */
+#include <iostream>
+#include <fstream>
+#include <getopt.h>
+#include <vector>
+#include <valarray>
+#include <math.h>
+#include <limits>
+#include <unistd.h>
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <itkImageIOFactory.h>
+
+#include "itkMincImageIOFactory.h"
+#include "itkMincImageIO.h"
+
+#include "minc_helpers.h"
+
+
+using namespace  std;
+using namespace  minc;
+
+void show_usage (const char * prog)
+{
+  std::cerr<<"Usage:"<<prog<<" in.mnc out.csv --mask <mask> --clobber --terse"<<endl; 
+}
+
+typedef itk::MincImageIO ImageIOType;
+
+
+int main (int argc, char **argv)
+{
+  std::string mask_f;
+  int verbose=0,clobber=0;
+  int terse=0;
+  static struct option long_options[] = {
+		{"verbose", no_argument,       &verbose, 1},
+		{"quiet",   no_argument,       &verbose, 0},
+		{"clobber", no_argument,       &clobber, 1},
+		{"mask",    required_argument,   0, 'm'},
+    {"terse",   no_argument,       &terse, 1},
+		{0, 0, 0, 0}
+		};
+  
+  for (;;) {
+      /* getopt_long stores the option index here. */
+      int option_index = 0;
+
+      int c = getopt_long (argc, argv, "m:vd:k:i:", long_options, &option_index);
+
+      /* Detect the end of the options. */
+      if (c == -1) break;
+
+      switch (c)
+			{
+			case 0:
+				break;
+			case 'v':
+				cout << "Version: 0.1" << endl;
+				return 0;
+      case 'm':
+        mask_f=optarg; break;
+			case '?':
+				/* getopt_long already printed an error message. */
+			default:
+				show_usage (argv[0]);
+				return 1;
+			}
+    }
+
+	if((argc - optind) < 2) {
+		show_usage (argv[0]);
+		return 1;
+	}
+  if (!clobber && !access(argv[optind+1], F_OK))
+  {
+    cerr << argv[optind+1] << " Exists!" << endl;
+    return 1;
+  }
+  
+	try
+  {
+		minc::mask3d::Pointer  mask(minc::mask3d::New());
+    
+    //creating a minc reader
+    ImageIOType::Pointer minc2ImageIO = ImageIOType::New();
+    
+    itk::ImageFileReader<minc::image3d >::Pointer reader = itk::ImageFileReader<minc::image3d >::New();
+    //initializing the reader
+    reader->SetFileName(argv[optind]);
+    reader->SetImageIO( minc2ImageIO );
+    reader->Update();
+    
+		minc::image3d::Pointer img=reader->GetOutput();
+    
+    
+    if(!mask_f.empty())
+    {
+      itk::ImageFileReader<minc::mask3d >::Pointer reader = itk::ImageFileReader<minc::mask3d >::New();
+      //initializing the reader
+      reader->SetFileName(mask_f.c_str());
+      reader->SetImageIO( minc2ImageIO );
+      reader->Update();
+      mask=reader->GetOutput();
+    }
+    const minc::mask3d::RegionType& reg=mask->GetLargestPossibleRegion();
+  
+    ofstream out(argv[optind+1]);
+    if(out.bad())
+      REPORT_ERROR ("can't open file");
+
+    image3d_iterator it(img,img->GetLargestPossibleRegion());
+    if(!terse) out<<"x,y,z,I"<<endl;
+    if(out.bad())
+      REPORT_ERROR ("can't write to file");
+    
+    mask3d::IndexType idx;
+    
+    for(it.GoToBegin();!it.IsAtEnd();++it) {
+      tag_point p;
+      img->TransformIndexToPhysicalPoint(it.GetIndex(),p);
+      
+      if(!mask_f.empty())
+      {
+        mask->TransformPhysicalPointToIndex(p,idx);
+        if(!reg.IsInside(idx) || !mask->GetPixel(idx))
+            continue;
+      }
+      
+      if(terse)
+        out<<it.Value()<<endl;
+      else
+      out<<p[0]<<","<<p[1]<<","<<p[2]<<","<<it.Value()<<endl;
+      if(out.bad())
+        REPORT_ERROR ("can't write to file");
+    }
+		return 0;
+	} catch (const minc::generic_error & err) {
+    cerr << "Got an error at:" << err.file () << ":" << err.line () << endl;
+    return 1;
+  }
+	return 0;
+};