Mercurial > hg > minc-tools

changeset 2577:55b122487957

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fixing saving dti images
author Vladimir S. FONOV <vfonov@ubuntu-64bit.(none)>
date Sun, 22 Jan 2012 06:47:16 +0800
parents f96c8a022cf6
children d4bcf5a71362
files CMakeLists.txt minc4itk/examples/itk_convert.cpp minc4itk/itkMincImageIO.cxx
diffstat 3 files changed, 147 insertions(+), 87 deletions(-) [+]
line wrap: on
line diff
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -55,7 +55,7 @@
   ExternalProject_Add(NETCDF 
     SOURCE_DIR NETCDF
     URL "ftp://ftp.unidata.ucar.edu/pub/netcdf/netcdf-4.0.1.tar.gz"
-    URL_MD5 a251453c5477599f050fa4e593295186
+    URL_MD5 "a251453c5477599f050fa4e593295186"
     BUILD_IN_SOURCE 1
     INSTALL_DIR "${CMAKE_INSTALL_PREFIX}"
     BUILD_COMMAND   make 
@@ -76,7 +76,7 @@
     ExternalProject_Add(HDF5
       SOURCE_DIR HDF5
       URL "http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.7/src/hdf5-1.8.7.tar.gz"
-      URL_MD5 37711d4bcb72997e93d495f97c76c33a
+      URL_MD5 "37711d4bcb72997e93d495f97c76c33a"
       BUILD_IN_SOURCE 1
       INSTALL_DIR     "${CMAKE_INSTALL_PREFIX}"
       BUILD_COMMAND   make 
--- a/minc4itk/examples/itk_convert.cpp
+++ b/minc4itk/examples/itk_convert.cpp
@@ -6,6 +6,8 @@
 #include <itkImageIOFactory.h>
 #include <itkImageIOBase.h>
 #include <itkFlipImageFilter.h>
+#include <itkMetaDataObject.h>
+#include <itkMetaDataDictionary.h>
 
 #include "itkMincImageIOFactory.h"
 #include "itkMincImageIO.h"
@@ -22,6 +24,7 @@
       << "--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
+      << "--show-meta show meta information (if present)"<<std::endl
 ;
 }
 typedef itk::ImageIOBase IOBase;
@@ -33,22 +36,39 @@
                     bool inv_y=false,
                     bool inv_z=false,
                     bool center=false,
-                    bool verbose=false)
+                    bool verbose=false,
+                    bool show_meta=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();
-     
+  
+  if(show_meta)
+  {
+    itk::MetaDataDictionary &thisDic=img->GetMetaDataDictionary();
+    
+    //let's write some meta information if there is any 
+    for(itk::MetaDataDictionary::ConstIterator it=thisDic.Begin();it!=thisDic.End();++it)
+    {
+			itk::MetaDataObjectBase *bs=(*it).second;
+			itk::MetaDataObject<std::string> * str=dynamic_cast<itk::MetaDataObject<std::string> *>(bs);
+			if(str)
+				std::cout<<(*it).first.c_str()<<" = "<< str->GetMetaDataObjectValue().c_str()<<std::endl;
+			else
+				std::cout<<(*it).first.c_str()<<" type: "<< typeid(*bs).name()<<std::endl;
+		}
+  }
+    
+  //TODO: convert metadata to and from minc format for DTI 
   /* WRITING */
   if(verbose)
     std::cout<<"Writing "<<fname<<"..."<<std::endl;
   
-
-  
   if(inv_x||inv_y||inv_z)
   {
     typename itk::FlipImageFilter<ImageType >::FlipAxesArrayType arr;
@@ -95,6 +115,7 @@
 {
   int verbose=0;
   int clobber=0;
+  int show_meta=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
 
@@ -106,6 +127,7 @@
     {"inv-y", no_argument,  &inv_y, 1},
     {"inv-z", no_argument,  &inv_z, 1},
     {"center", no_argument, &center, 1},
+    {"show-meta", no_argument, &show_meta, 1},
     {0, 0, 0, 0}
   };
     
@@ -172,69 +194,69 @@
     if(verbose)
     {
       std::cout<<"dimensions:"<<nd<<std::endl
-          <<"components:"<<nc<<std::endl
-          <<"type:"<<ct_s.c_str()<<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");
-        }
+      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,show_meta);
+	  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,show_meta);
+	  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,show_meta);
+	  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,show_meta);
+	  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,show_meta);
+	  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,show_meta);
+	  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,show_meta);
+	  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,show_meta);
+	  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);
+            load_and_save_image<itk::Image<unsigned char, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::Image<char, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::Image<unsigned short, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::Image<short, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::Image<int, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::Image<unsigned int, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::Image<float, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::Image<double, 4> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             break; 
           default:
             itk::ExceptionObject("Unsupported component type");
@@ -245,28 +267,28 @@
         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);
+            load_and_save_image<itk::VectorImage<unsigned char, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::VectorImage<char, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::VectorImage<unsigned short, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::VectorImage<short, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::VectorImage<int, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::VectorImage<unsigned int, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::VectorImage<float, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             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);
+            load_and_save_image<itk::VectorImage<double, 3> >(io,output.c_str(),inv_x,inv_y,inv_z,center,verbose,show_meta);
             break; 
           default:
             itk::ExceptionObject("Unsupported component type");
--- a/minc4itk/itkMincImageIO.cxx
+++ b/minc4itk/itkMincImageIO.cxx
@@ -421,6 +421,7 @@
       bool is_signed=false;
       bool have_vectors=false;
       bool have_time=false;
+      int dim_no=GetNumberOfDimensions();
       
       switch(GetComponentType())
       {
@@ -474,18 +475,24 @@
       std::vector<int> dimmap(5,-1);
       minc_info info;
       
-      if(GetNumberOfComponents()>1 && GetNumberOfComponents()<=3 ) 
+      if(GetNumberOfComponents()>1)
       {
-        have_vectors=true;
-        have_time=false;
-      } else if(GetNumberOfComponents()>3||GetNumberOfDimensions()>3) {
-        have_vectors=false;
-        have_time=true;
+	if(GetNumberOfComponents()<=3 ) 
+	{
+	  have_vectors=true;
+	  have_time=false;
+	} else if(GetNumberOfComponents()>3) {
+	  have_vectors=false;
+	  have_time=true;
+	}
+	dim_no++;
+      } else if(GetNumberOfDimensions()>3) {
+	have_vectors=false;
+	have_time=true;
       }
       
       if(itk::ExposeMetaData(thisDic,"dimorder",dimorder))
       {
-        //dimmap.resize(dimorder.size());
         for(int i=0,j=0;i<dimorder.size();i++)
         {
           if(dimorder[i]==MIvector_dimension && have_vectors)
@@ -507,29 +514,28 @@
         if(have_vectors)
         {
           dimmap[0]=0;
-        } else if(have_time) {
-          have_vectors=false;
-          have_time=true;
-          dimmap[GetNumberOfDimensions()+1]=GetNumberOfDimensions();
-        }
-        dimmap[1]=(have_vectors?1:0)+0;
-        dimmap[2]=(have_vectors?1:0)+1;
-        dimmap[3]=(have_vectors?1:0)+2;
+        } 
+        for(int i=0;i<dim_no-(have_vectors?1:0);i++)
+	  dimmap[1+i]=(have_vectors?1:0)+i;
       }
       
-      int dim_no=GetNumberOfDimensions()+(have_vectors||have_time?1:0);
       info.resize(dim_no);
       
-/*      std::cout<<"dimmap:";
+      #ifdef _DEBUG
+      std::cout<<"Number of components:"<<GetNumberOfComponents()<<std::endl;
+      std::cout<<"Number of dimensions:"<<GetNumberOfDimensions()<<std::endl;
+      std::cout<<"info.size()="<<info.size()<<std::endl;
+      std::cout<<"dimmap:";
       for(int i=0;i<5;i++)
       {
         std::cout<<dimmap[i]<<",";
       }
-      std::cout<<std::endl;*/
+      std::cout<<std::endl;
+      #endif  //_DEBUG
       
       for(int i=0;i<GetNumberOfDimensions();i++)
       {
-        int _i=dimmap[i+1];//GetNumberOfDimensions()-i-1+(have_vectors?1:0);
+        int _i=dimmap[i+1];
         if(_i<0) 
         {
           throw ExceptionObject(__FILE__, __LINE__,"Internal error");
@@ -555,6 +561,7 @@
             break;
         }
       }
+      
       if(GetNumberOfDimensions()==3) //we are only rotating 3D volumes
       {
         vnl_vector< double> start(3);
@@ -572,7 +579,24 @@
           int _i=dimmap[i+1];
           info[_i].start=start[i];
         }
+      } else {
+	std::cout<<"Warning: maging rotation matrix for "<<GetNumberOfDimensions()<<"D volume"<<std::endl;
+	std::cout<<"Origin:";
+	for(int i=0;i<GetNumberOfDimensions();i++)
+	  std::cout<<GetOrigin(i)<<",";
+	std::cout<<std::endl;
+	
+	std::cout<<"Rotation matrix:"<<std::endl;
+	for(int i=0;i<GetNumberOfDimensions();i++)
+	{
+	  std::cout<<"\t"<<i<<":";
+	  for(int j=0;j<GetNumberOfDimensions();j++)
+	    std::cout<<GetDirection(i)[j]<<",";
+	  std::cout<<std::endl;
+	}
+	std::cout<<std::endl;
       }
+      
       //here we assume that we had a grid file
       if(have_vectors)
       {
@@ -584,36 +608,50 @@
         info[_i].have_dir_cos=false;
       } else if(have_time){
         int _i=dimmap[4];
-        info[_i].length=GetNumberOfComponents();
+        info[_i].length=(GetNumberOfComponents()==1)?GetDimensions(3):GetNumberOfComponents();
         double tstep=1;
         double tstart=0;
-        itk::ExposeMetaData(thisDic,"tstep",tstep);
-        itk::ExposeMetaData(thisDic,"tstart",tstart);
-        
-        info[_i].step=tstep;
-        info[_i].start=tstart;
+	
+	if( itk::ExposeMetaData(thisDic,"tstep",tstep) &&
+	    itk::ExposeMetaData(thisDic,"tstart",tstart) )
+	{
+	  info[_i].step=tstep;
+	  info[_i].start=tstart;
+	} else if( GetNumberOfComponents()==1 ) {
+	  info[_i].step=GetSpacing(3);
+	  info[_i].start=GetOrigin(3);
+	} else {
+	  info[_i].step=1.0; //TODO: show somehow differently that we don't have info?
+	  info[_i].start=0;
+	}
+	
         info[_i].dim=dim_info::DIM_TIME;
-        info[_i].have_dir_cos=false;
+        info[_i].have_dir_cos=false; //TODO: what to do here?
       }
       
-/*      std::cout<<"info:";
+      #ifdef _DEBUG
+      std::cout<<"info:";
       for(int i=0;i<info.size();i++)
       {
         switch(info[i].dim)
         {
           case dim_info::DIM_VEC:
-            std::cout<<"vector_dimension,";break;
+            std::cout<<"vector_dimension";break;
           case dim_info::DIM_Z:
-            std::cout<<"zspace,";break;
+            std::cout<<"zspace";break;
           case dim_info::DIM_Y:
-            std::cout<<"yspace,";break;
+            std::cout<<"yspace";break;
           case dim_info::DIM_X:
-            std::cout<<"xspace,";break;
-          default: std::cout<<"Unknown! ";break;
-            
+            std::cout<<"xspace";break;
+	  case dim_info::DIM_TIME:
+	    std::cout<<"time";break;
+          default: 
+	    std::cout<<"unknown";break;
         }
+        std::cout<<"("<<info[i].length<<"),";
       }
-      std::cout<<std::endl;      */
+      std::cout<<std::endl;
+      #endif //_DEBUG
       
       //TODO:  shuffle info based on the dimnames
       _wrt=new minc_1_writer;