Mercurial > hg > minc-tools
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, ¢er, 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; +};