changeset 2556:7cdc8b84f038

added ezminc examples and test
author Vladimir S. FONOV <vladimir.fonov@gmail.com>
date Thu, 08 Dec 2011 19:06:27 -0500
parents 425fa26827ab
children 0d08a9897649
files ezminc/examples/fuzzy_volume_similarity.cpp ezminc/examples/trilinear_resample.cpp ezminc/examples/volume_avg.cpp ezminc/examples/volume_msq_dist.cpp ezminc/examples/volume_similarity.cpp ezminc/tests/minc_com.cpp ezminc/tests/minc_com_4d.cpp ezminc/tests/minc_rw_test.cpp ezminc/tests/minc_rw_test2.cpp ezminc/tests/minc_rw_test_4d.cpp ezminc/tests/minc_rw_test_simple.cpp ezminc/tests/transformpoint.cpp
diffstat 12 files changed, 1255 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
new file mode 100644
--- /dev/null
+++ b/ezminc/examples/fuzzy_volume_similarity.cpp
@@ -0,0 +1,209 @@
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       :  fuzzy_volume_similarity
+@DESCRIPTION:  an example of implimentation of fuzzy volume similarity 
+@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 "minc_1_rw.h"
+#include <iostream>
+#include "minc_1_simple.h"
+#include <getopt.h>
+#include <algorithm>
+
+using namespace minc;
+
+void show_usage (const char * prog)
+{
+  std::cout<<"Program calculates fuzzy volume similarity metrics"<<std::endl
+          <<"based on :  William R. Crum, Oscar Camara, and Derek L. G. Hill"
+          <<"\"Generalized Overlap Measures for Evaluation and Validation in Medical Image Analysis \""
+          <<" IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 25, NO. 11, NOVEMBER 2006"<<std::endl
+          <<"http://dx.doi.org/10.1109/TMI.2006.880587"<<std::endl<<std::endl
+          <<"Usage: "<<prog<<" <input1.mnc> <input2.mnc> [--verbose --mask <mask.mnc>]"<<std::endl;
+  
+}
+
+template<class T>class find_min_max
+{
+  public:
+  T _min;
+  T _max;
+  bool _initialized;
+  int _count;
+  
+  find_min_max():_initialized(false),_min(0),_max(0),_count(0)
+  {
+  }
+  
+  void operator()(const T& v)
+  {
+    if(_initialized)
+    {
+      if(v<_min) _min=v;
+      if(v>_max) _max=v;
+    } else {
+      _initialized=true;
+      _min=v;
+      _max=v;
+    }
+    _count++;
+  }
+  
+  const T& min(void) const {
+    return _min;
+  }
+  
+  const T& max(void) const{
+    return _max;
+  }
+  
+  int count(void) const {
+    return _count;
+  }
+};
+
+int main(int argc,char **argv)
+{
+  int verbose=0;
+  std::string mask_f;
+  static struct option long_options[] = {
+    {"verbose", no_argument,             &verbose, 1},
+    {"quiet",   no_argument,             &verbose, 0},
+    {"mask",     required_argument,              0,'m'},
+    
+    {0, 0, 0, 0}
+  };
+  
+  for (;;) {
+    /* getopt_long stores the option index here. */
+    int option_index = 0;
+
+    int c = getopt_long (argc, argv, "vqm:", long_options, &option_index);
+
+    /* Detect the end of the options. */
+    if (c == -1) break;
+
+    switch (c)
+    {
+      case 0:
+        break;
+      case 'v':
+        std::cout << "Version: 0.1" << std::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;
+  }
+  
+  try
+  {
+    minc_1_reader rdr1;
+    rdr1.open(argv[optind]);
+    
+    minc_1_reader rdr2;
+    rdr2.open(argv[optind+1]);
+    
+    if(rdr1.dim_no()!=rdr2.dim_no() )
+    {
+      std::cerr<<"Different number of dimensions!"<<std::endl;
+      return 1;
+    }
+    unsigned long size=1;
+    
+    for(int i=0;i<5;i++)
+    {
+      if(rdr1.ndim(i)!=rdr2.ndim(i))
+        std::cerr<<"Different dimensions length! "<<std::endl;
+      
+      if(rdr1.ndim(i)>0) size*=rdr1.ndim(i);
+    }
+    
+    for(int i=0;i<5;i++)
+    {
+      if(rdr1.nspacing(i)!=rdr2.nspacing(i) )
+        std::cerr<<"Different step size! "<<std::endl;
+    }
+    
+    std::vector<unsigned char> mask(size,1);
+    if(!mask_f.empty())
+    {
+      minc_1_reader rdr_m;
+      rdr_m.open(mask_f.c_str());
+      
+      for(int i=0;i<5;i++)
+      {
+        if(rdr1.ndim(i)!=rdr_m.ndim(i))
+          std::cerr<<"Different mask dimensions length! "<<std::endl;
+        
+        if(rdr1.nspacing(i)!=rdr_m.nspacing(i) )
+          std::cerr<<"Different mask step size! "<<std::endl;
+      }
+      rdr_m.setup_read_byte();
+      load_standard_volume<unsigned char>(rdr_m,&mask[0]);
+    }
+    
+    rdr1.setup_read_double();
+    rdr2.setup_read_double();
+    
+    std::vector<double> buffer1(size),buffer2(size);
+    
+    load_standard_volume<double>(rdr1,&buffer1[0]);
+    load_standard_volume<double>(rdr2,&buffer2[0]);
+    
+    find_min_max<double> f1,f2;
+    
+    //get min and max
+    f1=std::for_each(buffer1.begin(), buffer1.end(), f1);
+    f2=std::for_each(buffer2.begin(), buffer2.end(), f2);
+    
+    if( verbose )
+    {
+      std::cout<<"Volume 1 min="<<(int)f1.min()<<" max="<<(int)f1.max()<<" count="<<f1.count()<<std::endl;
+      std::cout<<"Volume 2 min="<<(int)f2.min()<<" max="<<(int)f2.max()<<" count="<<f2.count()<<std::endl;
+    }
+
+    unsigned char low= std::min<double>(f1.min(), f2.min());
+    unsigned char hi = std::max<double>(f1.max(), f2.max());
+
+    int v1=0,v2=0;
+    double a=0.0,b=0.0,c=0.0,d=0.0;
+    for(int i=0; i<size ; i++ )
+    {
+      
+      if(mask[i]>0)
+      {
+        a+=std::min<double>(buffer1[i],buffer2[i]);
+        b+=std::max<double>(buffer1[i],buffer2[i]);
+      }
+    }
+
+    std::cout.precision(10);
+    std::cout<<a/b<<std::endl;
+
+    
+  } catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  return 0;
+}
new file mode 100644
--- /dev/null
+++ b/ezminc/examples/trilinear_resample.cpp
@@ -0,0 +1,129 @@
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       :  trilinear_resample
+@DESCRIPTION:  an example of using trilinear resampling algorithm
+@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 "minc_1_rw.h"
+#include <iostream>
+#include "minc_1_simple.h"
+#include "minc_1_simple_rw.h"
+#include <getopt.h>
+#include <stdlib.h>
+#include <math.h>
+//based on http://www3.interscience.wiley.com/journal/121496529/abstract
+using namespace minc;
+
+void show_usage (const char * prog)
+{
+  std::cerr<<"Usage: "<<prog<<" <input.mnc> <output.mnc> [--step <f> --verbose ]"<<std::endl;
+}
+
+int main(int argc,char **argv)
+{
+  int verbose=0;
+  int csv=0;
+  int kappa=0,specificity=0,sensitivity=0;
+  int jaccard=0;
+  double _step=1.0;
+  static struct option long_options[] = {
+    {"verbose", no_argument,             &verbose, 1},
+    {"quiet",   no_argument,             &verbose, 0},
+    {"step",    required_argument,       0, 's'},
+    {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 0:
+        break;
+      case 'v':
+        std::cout << "Version: 0.1" << std::endl;
+        return 0;
+      case 's':
+        _step=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;
+  }
+  
+  try
+  {
+    minc_1_reader rdr;
+    rdr.open(argv[optind],true);
+    minc::minc_info new_info;
+
+    //rdr.setup_read_float();
+    simple_volume<float> input_vol,output_vol;
+    load_simple_volume<float>(rdr,input_vol);
+    new_info=rdr.info();
+    
+    minc::fixed_vec<3,float> old_step;
+    minc::fixed_vec<3,float> old_start,new_start;
+    minc::fixed_vec<3,int> new_len;
+    
+    for(int i=1;i<4;i++)
+    {
+      old_step[i-1]=new_info[rdr.map_space(i)].step;
+      old_start[i-1]=new_info[rdr.map_space(i)].start;
+      
+      float len=(new_info[rdr.map_space(i)].length)*old_step[i-1];
+      
+      new_info[rdr.map_space(i)].start-=old_step[i-1];
+      new_info[rdr.map_space(i)].step=_step;
+      new_info[rdr.map_space(i)].start+=_step/2;
+      new_start[i-1]=new_info[rdr.map_space(i)].start;
+      
+      new_len[i-1]=new_info[rdr.map_space(i)].length=ceil(fabs(len/_step));
+    }
+    
+    output_vol.resize(new_len);
+    
+    for(int z=0;z<new_len[2];z++)
+      for(int y=0;y<new_len[1];y++)
+        for(int x=0;x<new_len[0];x++)
+        {
+          minc::fixed_vec<3,float> cc=IDX<float>(x*_step,y*_step,z*_step);
+          
+          cc+=new_start;
+          cc/=old_step;
+          cc-=old_start;
+          
+          output_vol.set(x,y,z,input_vol.interpolate(cc[0],cc[1],cc[2]));
+        }
+    minc_1_writer wrt;
+    wrt.open(argv[optind+1],new_info,2,NC_FLOAT);
+    save_simple_volume<float>(wrt,output_vol);
+  } catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  return 0;
+}
new file mode 100644
--- /dev/null
+++ b/ezminc/examples/volume_avg.cpp
@@ -0,0 +1,148 @@
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       : volume_avg
+@DESCRIPTION: an example of calculating volume average
+@COPYRIGHT  :
+              Copyright 2009 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 "minc_1_rw.h"
+#include <iostream>
+#include "minc_1_simple.h"
+#include "minc_1_simple_rw.h"
+#include <getopt.h>
+#include <math.h>
+
+using namespace minc;
+
+void show_usage(const char *name)
+{
+  std::cerr 
+	  << "Usage: "<<name<<" <input1> .... <inputn>  <output> " << std::endl
+    << "\tn should be more than 1"<< std::endl
+    << "Optional parameters:" << std::endl
+    << "\t--verbose be verbose" << std::endl
+    << "\t--clobber clobber the output files" << std::endl
+    << "\t--sd <sd.mnc> "<< std::endl;
+}
+
+int main(int argc,char **argv)
+{
+  int clobber=0;
+  int verbose=0;
+  std::string sd_f;
+  
+	static struct option long_options[] =
+  {
+    {"verbose", no_argument, &verbose, 1},
+    {"quiet", no_argument, &verbose, 0},
+    {"clobber", no_argument, &clobber, 1},
+    {"sd", required_argument, 0, 's'},
+    {0, 0, 0, 0}
+  };
+  
+	int c;
+	for (;;)
+	{
+		/* getopt_long stores the option index here. */
+		int option_index = 0;
+
+		c = getopt_long (argc, argv, "s", long_options, &option_index);
+
+		/* Detect the end of the options. */
+		if (c == -1)
+			break;
+
+		switch (c)
+		{
+		case 0:
+			break;
+    case 's':
+      sd_f=optarg;
+      break;
+		case '?':
+			/* getopt_long already printed an error message. */
+		default:
+			show_usage(argv[0]);
+			return 1;
+		}
+	}
+
+	if ((argc - optind) < 3)
+	{
+		show_usage(argv[0]);
+		return 1;
+	}
+  
+  std::string output=argv[argc-1]; //last argument is output file... maybe we should make it a parameter instead?
+  argc-=optind;
+  
+  if(!clobber && !access (output.c_str(), F_OK))
+  {
+    std::cerr << output.c_str () << " Exists!" << std::endl;
+    return 1;
+  }
+
+	try
+  {
+    
+    minc_1_reader rdr1;
+    rdr1.open(argv[optind]);
+    minc_float_volume _avg;
+    
+    load_simple_volume<float>(rdr1,_avg);
+    
+    minc_float_volume _sd(_avg);
+    minc_float_volume _tmp(_avg);
+    
+
+    for(int i=0;i<_avg.c_buf_size();i++)
+    {
+      _sd.c_buf()[i]=_avg.c_buf()[i]*_avg.c_buf()[i];
+    }
+    
+    for(int i=1;i<(argc-1);i++)
+    {
+      minc_1_reader rdr2;
+      rdr2.open(argv[optind+i]);
+      if(!is_same(rdr1,rdr2))
+      {
+        return 1;
+      }
+      
+      load_simple_volume<float>(rdr2,_tmp);
+      _avg+=_tmp;
+      _tmp*=_tmp;
+      _sd+=_tmp;
+    }
+    _avg/=(float)(argc-1);
+    for(int i=0;i<_avg.c_buf_size();i++)
+    {
+      _sd.c_buf()[i]=sqrt(_sd.c_buf()[i]/(argc-1) - _avg.c_buf()[i]*_avg.c_buf()[i] );
+    }
+    minc_1_writer wrt;
+    wrt.open(output.c_str(),rdr1.info(),2,NC_FLOAT);
+    save_simple_volume<float>(wrt,_avg);
+    
+    if(!sd_f.empty())
+    {
+      minc_1_writer wrt2;
+      wrt2.open(sd_f.c_str(),rdr1.info(),2,NC_FLOAT);
+      save_simple_volume<float>(wrt2,_sd);
+    }
+    
+	} catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  
+  return 0;
+}
+
new file mode 100644
--- /dev/null
+++ b/ezminc/examples/volume_msq_dist.cpp
@@ -0,0 +1,98 @@
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       :  volume_msq_dist
+@DESCRIPTION:  an example of calculating mean squared dissimilarity
+@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 "minc_1_rw.h"
+#include <iostream>
+#include "minc_1_simple.h"
+
+using namespace minc;
+
+int main(int argc,char **argv)
+{
+  
+	try
+  {
+		if(argc<4) {
+      std::cerr<<"Usage: "<<argv[0]<<" <input1.mnc> <input2.mnc> <mask.mnc>"<<std::endl;
+      return 1;
+    }
+    minc_1_reader rdr1;
+    rdr1.open(argv[1]);
+    
+    minc_1_reader rdr2;
+    rdr2.open(argv[2]);
+    
+    minc_1_reader rdr_m;
+    rdr_m.open(argv[3]);
+   
+    if(rdr1.dim_no()!=rdr2.dim_no() || rdr1.dim_no()!=rdr_m.dim_no())
+    {
+      std::cerr<<"Different number of dimensions!"<<std::endl;
+      return 1;
+    }
+    unsigned long size=1;
+    for(int i=0;i<5;i++)
+    {
+      if(rdr1.ndim(i)!=rdr2.ndim(i) || rdr1.ndim(i)!=rdr_m.ndim(i))
+      {
+        std::cerr<<"Different dimensions length! "<<std::endl;
+      }
+      if(rdr1.ndim(i)>0) size*=rdr1.ndim(i);
+    }
+    
+    for(int i=0;i<5;i++)
+    {
+      if(rdr1.nspacing(i)!=rdr2.nspacing(i) ||  rdr1.nspacing(i)!=rdr_m.nspacing(i))
+      {
+        std::cerr<<"Different step size! "<<std::endl;
+      }
+    }
+    
+    //std::cout<<size<<std::endl;
+    rdr1.setup_read_float();
+    rdr2.setup_read_float();
+    rdr_m.setup_read_byte();
+    
+    std::vector<float> buffer1(size),buffer2(size);
+    std::vector<unsigned char> mask(size);
+    
+    load_standard_volume<float>(rdr1,&buffer1[0]);
+    load_standard_volume<float>(rdr2,&buffer2[0]);
+    load_standard_volume<unsigned char>(rdr_m,&mask[0]);
+    
+    double avg=0;
+    int cnt=0;
+    for(int i=0;i<size;i++)
+    {
+      if(!mask[i]) continue;
+      double d=buffer1[i]-buffer2[i];
+      avg+=d*d;
+      cnt++;
+    }
+    if(cnt)
+      avg/=cnt;
+    if( rdr1.ndim(0)==3 )  //this is a grid file
+      avg*=3;
+    std::cout.precision(10);
+    std::cout<<avg<<std::endl;
+    
+	} catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  
+  return 0;
+}
+
new file mode 100644
--- /dev/null
+++ b/ezminc/examples/volume_similarity.cpp
@@ -0,0 +1,259 @@
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       :  volume_similarity
+@DESCRIPTION:  an example of calculating volume similarity metrics 
+@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 "minc_1_rw.h"
+#include <iostream>
+#include "minc_1_simple.h"
+#include <getopt.h>
+#include <algorithm>
+
+using namespace minc;
+
+void show_usage (const char * prog)
+{
+  std::cout<<"Program calculates volume similarity metrics"<<std::endl
+           <<"if multiple labels are present, similarity metrics will be calculated for each label separately"<<std::endl
+          <<"based on : M. Feuerman and A. R. Miller "
+          <<"\"Relationships between statistical measures of agreement: sensitivity, specificity and kappa\""
+          <<" Journal of Evaluation in Clinical Practice vol. 14 no. 5 pp 930-933 2008"<<std::endl
+          <<"http://dx.doi.org/10.1111/j.1365-2753.2008.00984.x"<<std::endl<<std::endl
+          <<"Usage: "<<prog<<" <input1.mnc> <input2.mnc> [--kappa --sensitivity --specificity --verbose --jaccard --csv]"<<std::endl;
+  
+}
+
+template<class T>class find_min_max
+{
+  public:
+  T _min;
+  T _max;
+  bool _initialized;
+  int _count;
+  
+  find_min_max():_initialized(false),_min(0),_max(0),_count(0)
+  {
+  }
+  
+  void operator()(const T& v)
+  {
+    if(_initialized)
+    {
+      if(v<_min) _min=v;
+      if(v>_max) _max=v;
+    } else {
+      _initialized=true;
+      _min=v;
+      _max=v;
+    }
+    _count++;
+  }
+  
+  const T& min(void) const {
+    return _min;
+  }
+  
+  const T& max(void) const{
+    return _max;
+  }
+  
+  int count(void) const {
+    return _count;
+  }
+};
+
+int main(int argc,char **argv)
+{
+  int verbose=0;
+  int csv=0;
+  int kappa=0,specificity=0,sensitivity=0;
+  int jaccard=0;
+  static struct option long_options[] = {
+    {"verbose", no_argument,             &verbose, 1},
+    {"quiet",   no_argument,             &verbose, 0},
+    {"kappa",   no_argument,             &kappa, 1},
+    {"sensitivity",   no_argument,       &sensitivity, 1},
+    {"specificity",   no_argument,       &specificity, 1},
+    {"jaccard",       no_argument,       &jaccard, 1},
+    {"csv",           no_argument,       &csv, 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 0:
+        break;
+      case 'v':
+        std::cout << "Version: 0.1" << std::endl;
+        return 0;
+      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(! (sensitivity||specificity||kappa||jaccard) ) //no options selected, do all 
+  {
+    verbose=1;
+    sensitivity=1;
+    specificity=1;
+    kappa=1;
+    jaccard=1;
+  }
+  
+  try
+  {
+    minc_1_reader rdr1;
+    rdr1.open(argv[optind]);
+    
+    minc_1_reader rdr2;
+    rdr2.open(argv[optind+1]);
+    
+    if(rdr1.dim_no()!=rdr2.dim_no() )
+    {
+      std::cerr<<"Different number of dimensions!"<<std::endl;
+      return 1;
+    }
+    unsigned long size=1;
+    
+    for(int i=0;i<5;i++)
+    {
+      if(rdr1.ndim(i)!=rdr2.ndim(i))
+        std::cerr<<"Different dimensions length! "<<std::endl;
+      
+      if(rdr1.ndim(i)>0) size*=rdr1.ndim(i);
+    }
+    
+    for(int i=0;i<5;i++)
+    {
+      if(rdr1.nspacing(i)!=rdr2.nspacing(i) )
+        std::cerr<<"Different step size! "<<std::endl;
+    }
+    
+    rdr1.setup_read_byte();
+    rdr2.setup_read_byte();
+    
+    std::vector<unsigned char> buffer1(size),buffer2(size);
+    
+    load_standard_volume<unsigned char>(rdr1,&buffer1[0]);
+    load_standard_volume<unsigned char>(rdr2,&buffer2[0]);
+    
+    find_min_max<unsigned char> f1,f2;
+    
+    //get min and max
+    f1=std::for_each(buffer1.begin(), buffer1.end(), f1);
+    f2=std::for_each(buffer2.begin(), buffer2.end(), f2);
+    
+    if( verbose && !csv)
+    {
+      //std::cout<<buffer1.size() << std::endl;
+      //std::cout<<buffer1.begin()<<" "<<buffer1.end()<<std::endl;
+
+      std::cout<<"Volume 1 min="<<(int)f1.min()<<" max="<<(int)f1.max()<<" count="<<f1.count()<<std::endl;
+      std::cout<<"Volume 2 min="<<(int)f2.min()<<" max="<<(int)f2.max()<<" count="<<f2.count()<<std::endl;
+    }
+
+    unsigned char low= std::min<unsigned char>(f1.min(), f2.min());
+    unsigned char hi = std::max<unsigned char>(f1.max(), f2.max());
+
+    if(low==0) low=1; //assume 0 is background
+
+
+    for(unsigned char label=low; label<=hi; label++)
+    {
+      int v1=0,v2=0;
+      double a=0.0,b=0.0,c=0.0,d=0.0;
+      for(int i=0; i<size ; i++ )
+      {
+        if( buffer1[i]==label ) v1++;
+        if( buffer2[i]==label ) v2++;
+        if( buffer1[i]==label && buffer2[i]==label ) a++;
+        if( buffer1[i]==label && buffer2[i]!=label ) c++;
+        if( buffer1[i]!=label && buffer2[i]==label ) b++;
+        if( buffer1[i]!=label && buffer2[i]!=label ) d++;
+      }
+      double _kappa=0.0;
+      double _sensitivity=0.0;
+      double _specificity=0.0;
+      double _jaccard=0.0;
+
+      if(v1>0 && v2>0)
+      {
+        //_kappa=v1v2*2.0/(v1+v2);
+        _kappa=a*2.0/(v1+v2);
+
+        _sensitivity=a/(a+c);
+        _specificity=d/(b+d);
+
+        _jaccard=a/(a+c+b);
+      }
+
+      std::cout.precision(10);
+
+      if( hi!=low && verbose && !csv)
+        std::cout<<"Label: "<<(int)label<<std::endl;
+
+      if( csv )
+      {
+        if(low!=hi)
+          std::cout<<(int)label<<",";
+
+        std::cout<<_kappa<<",";
+        std::cout<<_sensitivity<<",";
+        std::cout<<_specificity<<",";
+        std::cout<<_jaccard<<std::endl;
+      } else {
+        if( kappa ){
+          if(verbose) std::cout<<"Kappa ";
+          std::cout<<_kappa<<std::endl;
+        }
+        if( sensitivity )
+        {
+          if(verbose) std::cout<<"Sensitivity ";
+          std::cout<<_sensitivity<<std::endl;
+        }
+        if( specificity )
+        {
+          if(verbose) std::cout<<"Specificity ";
+          std::cout<<_specificity<<std::endl;
+        }
+        if( jaccard )
+        {
+          if(verbose) std::cout<<"Jaccard similarity ";
+          std::cout<<_jaccard<<std::endl;
+        }
+      }
+    }
+    
+  } catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  return 0;
+}
new file mode 100644
--- /dev/null
+++ b/ezminc/tests/minc_com.cpp
@@ -0,0 +1,48 @@
+#include "minc_1_rw.h"
+#include <iostream>
+#include "minc_1_simple.h"
+#include "minc_1_simple_rw.h"
+
+using namespace minc;
+
+int main(int argc,char **argv)
+{
+  try
+  {
+    if(argc<2) {
+      std::cerr<<"Usage: "<<argv[0]<<" <input.mnc> "<<std::endl;
+      return 1;
+    }
+    minc_1_reader rdr;
+    rdr.open(argv[1],true);
+    simple_volume<float> vol;
+    load_simple_volume<float>(rdr,vol);
+    //calculate COM
+    fixed_vec<3,int> i;
+    fixed_vec<3,double> w;
+    double total=0;
+    double tx=0.0,ty=0.0,tz=0.0;
+    
+    for(i[2]=0;i[2]<vol.dim(2);i[2]++)
+      for(i[1]=0;i[1]<vol.dim(1);i[1]++)
+        for(i[0]=0;i[0]<vol.dim(0);i[0]++)
+        {
+          w=vol.voxel_to_world(i);
+          tx+=w[0]*vol.get(i);
+          ty+=w[1]*vol.get(i);
+          tz+=w[2]*vol.get(i);
+          total+=vol.get(i);
+        }
+    tx/=total;
+    ty/=total;
+    tz/=total;
+    std::cout<<tx<<","<<ty<<","<<tz<<std::endl;
+  } catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  
+  return 0;
+}
+
new file mode 100644
--- /dev/null
+++ b/ezminc/tests/minc_com_4d.cpp
@@ -0,0 +1,70 @@
+#include "minc_1_rw.h"
+#include <iostream>
+#include "minc_1_simple_rw.h"
+
+using namespace minc;
+
+int main(int argc,char **argv)
+{
+  try
+  {
+    if(argc<2) {
+      std::cerr<<"Usage: "<<argv[0]<<" <input.mnc> "<<std::endl;
+      return 1;
+    }
+    minc_1_reader rdr;
+    rdr.open(argv[1],true);
+    simple_4d_volume<float> vol;
+    load_4d_volume<float>(rdr,vol);
+    //calculate COM
+    fixed_vec<3,int> i;
+    fixed_vec<3,double> w;
+    
+    double total=0.0;
+    double tx=0.0,ty=0.0,tz=0.0,tT=0.0;
+    
+    //std::cout<<"T-start:"<<vol.t_start()<<std::endl;
+    //std::cout<<"T-step:"<<vol.t_step()<<std::endl;
+    
+    for(int t=0;t<vol.frames();t++)
+    {
+      double rt=t*vol.t_step()+vol.t_start();
+      
+      for(i[2]=0;i[2]<vol.dim(2);i[2]++)
+        for(i[1]=0;i[1]<vol.dim(1);i[1]++)
+          for(i[0]=0;i[0]<vol.dim(0);i[0]++)
+      {
+        w=vol.voxel_to_world(i);
+        
+        double v=vol.get(i,t);
+        
+        tx+=w[0]*v;
+        ty+=w[1]*v;
+        tz+=w[2]*v;
+        tT+=rt  *v;
+        
+        total+=v;
+        //std::cout<<v<<"\t";
+      }
+      //std::cout<<std::endl;
+      //std::cout<<rt<<"\t"<<tT<<"\t"<<std::flush;
+      std::cout<<"."<<std::flush;
+    }
+    std::cout<<std::endl;
+    
+    tx/=total;
+    ty/=total;
+    tz/=total;
+    tT/=total;
+    
+    std::cout<<tx<<","<<ty<<","<<tz<<","<<tT<<std::endl;
+    
+  } catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  
+  return 0;
+}
+
new file mode 100644
--- /dev/null
+++ b/ezminc/tests/minc_rw_test.cpp
@@ -0,0 +1,129 @@
+#include "minc_1_rw.h"
+#include <iostream>
+#include "minc_1_simple.h"
+
+using namespace minc;
+
+template<class T>class volume_3d
+{
+  protected:
+    T* _data;
+    int _x,_y,_z;
+  public:
+    
+    volume_3d(T* data,int x,int y,int z):_data(data),_x(x),_y(y),_z(z)
+    {
+    }
+    
+    volume_3d(std::vector<T> data,int x,int y,int z):_data(&data[0]),_x(x),_y(y),_z(z)
+    {
+    }
+    
+    T& operator()(int x,int y,int z)
+    {
+      return _data[x+y*_x+z*_x*_y];
+    }
+  
+};
+
+int main(int argc,char **argv)
+{
+  try
+  {
+    if(argc<4) {
+      std::cerr<<"Usage: "<<argv[0]<<" <input.mnc> <output.mnc> <output2.mnc>"<<std::endl;
+      return 1;
+    }
+    minc_1_reader rdr;
+    rdr.open(argv[1]);
+    int i;
+    for(i=0;i<rdr.dim_no();i++)
+      std::cout<<rdr.dim(i).name<<" "<<rdr.dim(i).length<<std::endl;
+    std::cout<<"Slice len="<<rdr.slice_len()<<std::endl;
+    std::cout<<"History:"<<rdr.history().c_str()<<std::endl;
+    for(int v=0;v<rdr.var_number();v++)
+    {
+      std::string var=rdr.var_name(v);
+      for(int a=0;a<rdr.att_number(v);a++)
+      {
+        std::string aname=rdr.att_name(v,a);
+        nc_type dt=rdr.att_type(v,aname.c_str());
+        
+        std::cout<<var.c_str()<<":"<<aname.c_str()<<" ";
+        if(dt==NC_CHAR) 
+          std::cout<<rdr.att_value_string(v,aname.c_str());
+        else if(dt=NC_DOUBLE) 
+        {
+          std::vector<double> val=rdr.att_value_double(v,aname.c_str());
+          for(int d=0;d<val.size();d++)
+            std::cout<<val[d]<<"\t";
+        } else {
+          std::cout<<"???";
+        }
+        std::cout<<std::endl;
+      }
+    }
+    
+    rdr.setup_read_float();
+    
+    minc_1_writer wrt;
+    wrt.open(argv[2],rdr.info(),2,NC_FLOAT,true);
+    wrt.setup_write_float();
+    std::vector<float> v(rdr.slice_len());
+
+    int c=0;
+    double s=0;
+    for(rdr.begin(),wrt.begin();!rdr.last();rdr.next_slice(),wrt.next_slice())
+    {
+      //std::cout<<wrt.current_slice()[0]<<std::endl;
+      rdr.read(&v[0]);
+      for(int i=0;i<rdr.slice_len();i++)
+      {
+        s+=v[i];
+        c++;
+      }
+      wrt.write(&v[0]);
+    }
+    s/=c;
+    //std::cout<<wrt._image_range[0]<<" "<<wrt._image_range[1]<<std::endl;
+    std::cout<<s<<std::endl;
+    wrt.copy_headers(rdr);
+    wrt.append_history("minc_rw_test test1\n");
+    
+    //second test
+    unsigned long size=1;
+    for(i=0;i<rdr.dim_no();i++)
+      size*=rdr.dim(i).length;
+    
+    minc_1_writer wrt2;
+    wrt2.open(argv[3],rdr.info(),3,NC_SHORT,false);
+    wrt2.setup_write_float();
+    
+    std::vector<float> buffer(size);
+    load_standard_volume<float>(rdr,&buffer[0]);
+    double avg2=0;
+    //for(int i=0;i<size;i++)
+    //  avg2+=buffer[i];
+    volume_3d<float> volume(buffer,rdr.ndim(dim_info::DIM_X),rdr.ndim(dim_info::DIM_Y),rdr.ndim(dim_info::DIM_Z));
+    for(int z=0;z<rdr.ndim(dim_info::DIM_Z);z++)
+      for(int y=0;y<rdr.ndim(dim_info::DIM_Y);y++)
+        for(int x=0;x<rdr.ndim(dim_info::DIM_X);x++)
+        {
+          avg2+=volume(x,y,z);
+        }
+    
+    avg2/=size;
+    std::cout<<"avg2 "<<avg2<<std::endl;
+    
+    save_standard_volume<float>(wrt2,&buffer[0]);
+    wrt2.copy_headers(rdr);
+    wrt2.append_history("minc_rw_test test2\n");
+  } catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  
+  return 0;
+}
+
new file mode 100644
--- /dev/null
+++ b/ezminc/tests/minc_rw_test2.cpp
@@ -0,0 +1,42 @@
+#include "minc_1_rw.h"
+#include <iostream>
+#include "minc_1_simple.h"
+#include "minc_io_simple_volume.h"
+
+using namespace minc;
+
+int main(int argc,char **argv)
+{
+  try
+  {
+    if(argc<3) {
+      std::cerr<<"Usage: "<<argv[0]<<" <input.mnc> <output.mnc>"<<std::endl;
+      return 1;
+    }
+    minc_1_reader rdr;
+    rdr.open(argv[1],true);
+    rdr.setup_read_float();
+    simple_volume<float> vol;
+    vol.resize(rdr.ndim(1),rdr.ndim(2),rdr.ndim(3));
+    load_non_standard_volume<float>(rdr,vol.c_buf());
+    //rdr.close();
+    for(int z=0;z<vol.dim(2);z++)
+      for(int y=0;y<vol.dim(1);y++)
+        for(int x=0;x<vol.dim(0);x++)
+        {
+          vol(x,y,z)=x;
+        }
+    minc_1_writer wrt;
+    wrt.open(argv[2],rdr.info(),3,NC_FLOAT,false);
+    wrt.setup_write_float();
+    save_non_standard_volume<float>(wrt,vol.c_buf());
+    
+  } catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  
+  return 0;
+}
+
new file mode 100644
--- /dev/null
+++ b/ezminc/tests/minc_rw_test_4d.cpp
@@ -0,0 +1,51 @@
+#include "minc_1_simple_rw.h"
+#include <iostream>
+
+using namespace minc;
+
+int main(int argc,char **argv)
+{
+  try
+  {
+    if(argc<3) {
+      std::cerr<<"Usage: "<<argv[0]<<" <input.mnc> <output.mnc>"<<std::endl;
+      return 1;
+    }
+    minc_1_reader rdr;
+    rdr.open(argv[1]);
+    
+    if(rdr.dim_no()==3|| (rdr.dim_no()==4 && rdr.ndim(4)>0 ))
+    {
+      if(rdr.datatype()==NC_FLOAT || rdr.datatype()==NC_SHORT)
+      {
+        std::cout<<"Reading float volume"<<std::endl;
+        simple_4d_volume<float> vol;
+        
+        load_4d_volume(rdr,vol);
+        save_minc_file(argv[2],vol,"test",&rdr,rdr.datatype(),rdr.is_signed());
+      } else if(rdr.datatype()==NC_BYTE) {
+        std::cout<<"Reading byte volume"<<std::endl;
+        simple_4d_volume<unsigned char> vol;
+        
+        load_4d_volume(rdr,vol);
+        save_minc_file(argv[2],vol,"test",&rdr,rdr.datatype(),rdr.is_signed());
+      }
+    } else if((rdr.dim_no()==4|| rdr.dim_no()==5)&&rdr.ndim(0)==3) { //we are dealing with vectors
+        std::cout<<"Reading vector volume"<<std::endl;
+        simple_4d_volume< fixed_vec<3,float> > vol;
+        
+        load_4d_volume(rdr,vol);
+        save_minc_file(argv[2],vol,"test",&rdr,rdr.datatype(),rdr.is_signed());
+    }
+    
+  } catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  
+  return 0;
+}
+
+
+
new file mode 100644
--- /dev/null
+++ b/ezminc/tests/minc_rw_test_simple.cpp
@@ -0,0 +1,43 @@
+#include "minc_1_rw.h"
+#include <iostream>
+#include "minc_1_simple.h"
+
+using namespace minc;
+
+int main(int argc,char **argv)
+{
+  try
+  {
+    if(argc<3) {
+      std::cerr<<"Usage: "<<argv[0]<<" <input.mnc> <output.mnc>"<<std::endl;
+      return 1;
+    }
+    minc_1_reader rdr;
+    rdr.open(argv[1],true);
+    rdr.setup_read_float();
+    
+    minc_1_writer wrt;
+    for(int i=0;i<3;i++)
+    {
+      std::cout<<rdr.info()[i].step<<" "<<rdr.info()[i].start<<std::endl;
+    }
+    wrt.open(argv[2],rdr.info(),3,NC_FLOAT,false);
+    wrt.setup_write_float();
+    minc_input_iterator<float> in(rdr);
+    minc_output_iterator<float> out(wrt);
+    float c=0.0;
+    for(in.begin(),out.begin();!in.last();in.next(),out.next())
+    {
+      out.value(in.value()+c);
+      c++;
+    }
+    
+	} catch (const minc::generic_error & err) {
+    std::cerr << "Got an error at:" << err.file () << ":" << err.line () << std::endl;
+    std::cerr << err.msg()<<std::endl;
+    return 1;
+  }
+  
+  return 0;
+}
+
new file mode 100644
--- /dev/null
+++ b/ezminc/tests/transformpoint.cpp
@@ -0,0 +1,29 @@
+#include <volume_io.h>
+#include <iostream>
+
+
+int main(int argc, char **argv)
+{
+  if(argc<5)
+  {
+    std::cerr<<"Usage:"<<argv[0]<<"<XFM file> X Y Z"<<std::endl;
+    return 1;
+  }
+  double x,y,z;
+  double _x,_y,_z;
+  
+  x=atof(argv[2]);
+  y=atof(argv[3]);
+  z=atof(argv[4]);
+  
+  General_transform _xfm;
+  if(input_transform_file((char*)argv[1], &_xfm)!=OK)
+  {
+    std::cerr<<"Error reading:"<<argv[1]<<std::endl;
+    return 1;
+  }
+  general_transform_point( &_xfm, x, y, z,&_x, &_y, &_z);  
+  delete_general_transform(&_xfm);
+  std::cout<<_x<<" "<<_y<<" "<<_z<<std::endl;
+  return 0;
+};
\ No newline at end of file