Mercurial > hg > minc-tools
changeset 2552:345f8c960657
added ezminc library
| author | Vladimir S. FONOV <vladimir.fonov@gmail.com> |
|---|---|
| date | Thu, 08 Dec 2011 18:47:56 -0500 |
| parents | 3fe122bac48b |
| children | 0b52afea4248 |
| files | CMakeLists.txt ezminc/CMakeLists.txt ezminc/minc_1_rw.cpp ezminc/minc_1_rw.h ezminc/minc_1_simple.h ezminc/minc_1_simple_rw.cpp ezminc/minc_1_simple_rw.h ezminc/minc_io_4d_volume.h ezminc/minc_io_exceptions.h ezminc/minc_io_fixed_vector.h ezminc/minc_io_simple_volume.h |
| diffstat | 11 files changed, 3677 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,6 +31,7 @@ OPTION(BUILD_MINC2 "Support minc2 file format" ON) OPTION(BUILD_TOOLS "Build minc tools (mincreshape,mincresample, etc)" ON) OPTION(BUILD_CONVERTERS "Build minc conversion programs (mnc2nii, nii2mnc , dcm2mnc...)" ON) +OPTION(BUILD_EZMINC "Build C++ interface library EZminc" ON) ADD_DEFINITIONS(-DHAVE_CONFIG_H) @@ -237,3 +238,7 @@ # and then the conversion subdir ADD_SUBDIRECTORY( conversion ) ENDIF(BUILD_CONVERTERS) + +IF(BUILD_EZMINC) +ADD_SUBDIRECTORY( ezminc ) +ENDIF(BUILD_EZMINC) \ No newline at end of file
new file mode 100644 --- /dev/null +++ b/ezminc/CMakeLists.txt @@ -0,0 +1,39 @@ +LINK_LIBRARIES( + m + z ) + +IF(USE_MINC2) + ADD_DEFINITIONS( -DMINC2 ) + LINK_LIBRARIES(volume_io2 minc2 netcdf hdf5 z) +ELSE(USE_MINC2) + LINK_LIBRARIES(volume_io minc netcdf) +ENDIF(USE_MINC2) + + +INCLUDE_DIRECTORIES( + ${VOLUME_TOOLS_BINARY_DIR} + ${CMAKE_CURRENT_BINARY_DIR} + ${CMAKE_CURRENT_SOURCE_DIR} +) + + +SET( MINC_IO_HEADERS + minc_io_exceptions.h + minc_io_fixed_vector.h + minc_io_simple_volume.h + minc_1_rw.h + minc_1_simple.h + minc_1_simple_rw.h + minc_io_4d_volume.h + ) + +SET( MINC_IO_SRC + minc_1_rw.cpp + minc_1_simple_rw.cpp + ) + +ADD_LIBRARY( minc_io ${MINC_IO_HEADERS} ${MINC_IO_SRC}) + + +INSTALL(TARGETS minc_io ARCHIVE DESTINATION lib) +INSTALL(FILES ${MINC_IO_HEADERS} DESTINATION include)
new file mode 100644 --- /dev/null +++ b/ezminc/minc_1_rw.cpp @@ -0,0 +1,1421 @@ +/* ----------------------------- MNI Header ----------------------------------- +@NAME : +@DESCRIPTION: Primitive" interface to minc files, using minctoraw and rawtominc programs + does not require linking with minc2 library. + Created during the days when minc2 didn't compile on 64bit architecture +@COPYRIGHT : + Copyright 2007 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 <stdio.h> +#include <iostream> +#include <sstream> +#include <string.h> +//minc stuff +#include <math.h> +#include <limits.h> +#include "minc_1_rw.h" + +namespace minc +{ + dim_info::dim_info(int l, double sta, + double spa,dimensions d, + bool hd): + length(l),start(sta),step(spa),dim(d),have_dir_cos(hd) + { + switch(dim) + { + case dim_info::DIM_X:name=MIxspace;break; + case dim_info::DIM_Y:name=MIyspace;break; + case dim_info::DIM_Z:name=MIzspace;break; + case dim_info::DIM_TIME:name=MItime;break; + case dim_info::DIM_VEC:name=MIvector_dimension;break; + default: REPORT_ERROR("Unknown Dimension!"); + } + } + + minc_1_base::minc_1_base(): + _dims(3,0), + _map_to_std(5,-1), + _mincid(MI_ERROR), + _imgid(MI_ERROR), + _icvid(MI_ERROR), + _cur(MAX_VAR_DIMS,0), + _slab(MAX_VAR_DIMS,1), + _last(false), + _datatype(MI_ORIGINAL_TYPE), + _io_datatype(MI_ORIGINAL_TYPE), + _slab_len(0), + _positive_directions(true), + _minc2(false) + { + _icvid=miicv_create(); + } + + //! destructor, closes minc file + minc_1_base::~minc_1_base() + { + close(); + } + + void minc_1_base::close(void) + { + if(_icvid!=MI_ERROR) + { + miicv_free(_icvid); + _icvid=MI_ERROR; + } + if(_mincid!=MI_ERROR) miclose(_mincid); + _mincid=MI_ERROR; + } + + + std::string minc_1_base::history(void) const + { + nc_type datatype; + int att_length; +#ifndef WIN32 + int op=ncopts; +#endif //WIN32 + //ncopts=0; + if ((ncattinq(_mincid, NC_GLOBAL, MIhistory, &datatype,&att_length) == MI_ERROR) || + (datatype != NC_CHAR)) + { + //ncopts=op; + return ""; + } + char* str = new char[att_length+1]; + str[0] = '\0'; + miattgetstr(_mincid, NC_GLOBAL, (char*)MIhistory, att_length,str); + //ncopts=op; + std::string r(str); + delete [] str; + return r; + } + + //code from mincinfo + int minc_1_base::var_number(void) const + { + int nvars; + if(ncinquire(_mincid, NULL, &nvars, NULL, NULL)!=MI_ERROR) + return nvars; + return 0; + } + + std::string minc_1_base::var_name(int no) const + { + char name[MAX_NC_NAME]; + if(ncvarinq(_mincid, no, name, NULL, NULL, NULL, NULL)!=MI_ERROR) + return name; + } + + int minc_1_base::att_number(const char *var_name) const + { + int varid; + int natts; + if (*var_name=='\0') { + varid = NC_GLOBAL; + } else { + if((varid = ncvarid(_mincid, var_name))==MI_ERROR) + return 0; + } + return att_number(varid); + } + + int minc_1_base::var_id(const char *var_name) const + { + return ncvarid(_mincid, var_name); + } + + long minc_1_base::var_length(const char *var_name) const + { + int varid=var_id(var_name); + if(varid!=MI_ERROR) + return var_length(varid); + else + return 0; + } + + long minc_1_base::var_length(int var_id) const + { + int vardims; + + if(ncvarinq(_mincid, var_id, NULL, NULL, &vardims, NULL, NULL)!=MI_ERROR) + { + if(vardims==0) return 1; + int *dims=new int[vardims]; + if(ncvarinq(_mincid, var_id, NULL, NULL, NULL, dims, NULL)!=MI_ERROR) + { + long varlength=1; + if(ncdiminq(_mincid,dims[0],NULL,&varlength)!=MI_ERROR) + { + delete[] dims; + return varlength; + } + delete[] dims; + return 1; + } else { + delete[] dims; + return 1; + } + } + return 0; + } + + + //! get the number of attributes associated with variable + int minc_1_base::att_number(int var_no) const + { + int natts; + if(ncvarinq(_mincid, var_no, NULL, NULL, NULL, NULL, &natts)!=MI_ERROR) + return natts; + return 0; + } + + + std::string minc_1_base::att_name(const char *var_name,int no) const + { + int varid; + int attid; + char name[MAX_NC_NAME]; + if (*var_name=='\0') + varid = NC_GLOBAL; + else + { + if((varid = ncvarid(_mincid, var_name))==MI_ERROR) + return ""; + } + return att_name(varid,no); + } + + std::string minc_1_base::att_name(int varid,int no) const + { + int attid; + char name[MAX_NC_NAME]; + if(ncattname(_mincid, varid, no, name)==MI_ERROR) + return ""; + + return name; + } + + std::string minc_1_base::att_value_string(const char *var_name,const char *att_name) const + { + int varid; + int attid; + char name[MAX_NC_NAME]; + if (*var_name=='\0') + varid = NC_GLOBAL; + else + { + if((varid = ncvarid(_mincid, var_name))==MI_ERROR) + return ""; + } + return att_value_string(varid,att_name); + } + + std::string minc_1_base::att_value_string(int varid,const char *att_name) const + { + int att_length; + nc_type datatype; + + //TODO: make this handle other (double?) data types correctly + if ((ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) || + (datatype != NC_CHAR)) + { + //ncopts=op; + return ""; + } + char* str = new char[att_length+1]; + str[0] = '\0'; + miattgetstr(_mincid, varid, (char *)att_name, att_length, str); + //ncopts=op; + std::string r(str); + delete [] str; + return r; + } + + std::vector<double> minc_1_base::att_value_double(const char *var_name,const char *att_name) const + { + int varid; + if (*var_name=='\0') + varid = NC_GLOBAL; + else + { + if((varid = ncvarid(_mincid, var_name))==MI_ERROR) + return std::vector<double>(0); + } + return att_value_double(varid,att_name); + } + + std::vector<short> minc_1_base::att_value_short(const char *var_name,const char *att_name) const + { + int varid; + if (*var_name=='\0') + varid = NC_GLOBAL; + else + { + if((varid = ncvarid(_mincid, var_name))==MI_ERROR) + return std::vector<short>(0); + } + return att_value_short(varid,att_name); + } + + std::vector<unsigned char> minc_1_base::att_value_byte(const char *var_name,const char *att_name) const + { + int varid; + if (*var_name=='\0') + varid = NC_GLOBAL; + else + { + if((varid = ncvarid(_mincid, var_name))==MI_ERROR) + return std::vector<unsigned char>(0); + } + return att_value_byte(varid,att_name); + } + + + std::vector<int> minc_1_base::att_value_int(const char *var_name,const char *att_name) const + { + int varid; + if (*var_name=='\0') + varid = NC_GLOBAL; + else + { + if((varid = ncvarid(_mincid, var_name))==MI_ERROR) + return std::vector<int>(0); + } + return att_value_int(varid,att_name); + } + + std::vector<int> minc_1_base::att_value_int(int varid,const char *att_name) const + { + int att_length; + nc_type datatype; + + if ((ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) || + (datatype != NC_INT)) + { + //ncopts=op; + return std::vector<int>(0); + } + std::vector<int> r(att_length); + miattget(_mincid, varid, (char*)att_name, NC_INT, att_length,&r[0], NULL) ; + //ncopts=op; + return r; + } + + + std::vector<double> minc_1_base::att_value_double(int varid,const char *att_name) const + { + int att_length; + nc_type datatype; + + //TODO: make this handle other (double?) data types correctly + if ((ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) || + (datatype != NC_DOUBLE)) + { + //ncopts=op; + return std::vector<double>(0); + } + std::vector<double> r(att_length); + miattget(_mincid, varid, (char*)att_name, NC_DOUBLE, att_length,&r[0], NULL) ; + //ncopts=op; + return r; + } + + std::vector<short> minc_1_base::att_value_short(int varid,const char *att_name) const + { + int att_length; + nc_type datatype; + + //TODO: make this handle other (double?) data types correctly + if ((ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) || + (datatype != NC_SHORT)) + { + //ncopts=op; + return std::vector<short>(0); + } + std::vector<short> r(att_length); + miattget(_mincid, varid, (char*)att_name, NC_SHORT, att_length,&r[0], NULL) ; + //ncopts=op; + return r; + } + + std::vector<unsigned char> minc_1_base::att_value_byte(int varid,const char *att_name) const + { + int att_length; + nc_type datatype; + + //TODO: make this handle other (double?) data types correctly + if ((ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) || + (datatype != NC_BYTE)) + { + //ncopts=op; + return std::vector<unsigned char>(0); + } + std::vector<unsigned char> r(att_length); + miattget(_mincid, varid, (char*)att_name, NC_BYTE, att_length,&r[0], NULL) ; + //ncopts=op; + return r; + } + + + nc_type minc_1_base::att_type(const char *var_name,const char *att_name) const + { + int varid; + if (*var_name=='\0') + varid = NC_GLOBAL; + else + { + if((varid = ncvarid(_mincid, var_name))==MI_ERROR) + return MI_ORIGINAL_TYPE; + } + return att_type(varid,att_name); + } + + nc_type minc_1_base::att_type(int varid,const char *att_name) const + { + int att_length; + nc_type datatype; + + if(ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) + return MI_ORIGINAL_TYPE; + return datatype; + } + + + int minc_1_base::att_length(const char *var_name,const char *att_name) const + { + int varid; + if (*var_name=='\0') + varid = NC_GLOBAL; + else + { + if((varid = ncvarid(_mincid, var_name))==MI_ERROR) + return 0; + } + return att_length(varid,att_name); + } + + int minc_1_base::att_length(int varid,const char *att_name) const + { + int att_length; + nc_type datatype; + + if(ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) + return 0; + + return att_length; + } + + + minc_1_reader::minc_1_reader(const minc_1_reader&):_metadate_only(false),_have_temp_file(false),_read_prepared(false) + { + } + + minc_1_reader::minc_1_reader():_metadate_only(false),_have_temp_file(false),_read_prepared(false) + { + } + + + //based on the code from mincextract + void minc_1_reader::open(const char *path,bool positive_directions/*=true*/,bool metadate_only/*=false*/,bool rw/*=false*/) + { +#ifndef WIN32 + ncopts = 0; +#endif + _metadate_only=metadate_only; + _read_prepared=false; + int element_size; + int idim; + int nstart, ncount; + _positive_directions=positive_directions; + //ncopts = 0; + + // Open the file + + // Expand file + if(_metadate_only) + { + int created_tempfile; + char * tempfile = miexpand_file((char*)path, NULL, true, &created_tempfile); + if (tempfile == NULL) REPORT_ERROR("Error expanding minc file"); + _tempfile=tempfile; + path=_tempfile.c_str(); + _have_temp_file=created_tempfile; + } + _mincid = miopen((char*)path, rw?NC_WRITE:NC_NOWRITE); + + if(_mincid == MI_ERROR) REPORT_ERROR("Can't open minc file for reading!"); +#ifdef MINC2 + if (MI2_ISH5OBJ(_mincid)) { + _minc2 = true; + } +#endif + + /* Inquire about the image variable */ + _imgid = ncvarid(_mincid, MIimage); + if(_imgid == MI_ERROR) REPORT_ERROR("Can't get Image ID"); + ncvarinq(_mincid, _imgid, NULL, NULL, &_ndims, mdims, NULL); + //get image data type... not used for now + miget_datatype(_mincid, _imgid, &_datatype, &_is_signed); + //dir_cos.SetIdentity(); + int dim_cnt=0; + miget_image_range(_mincid, _image_range); + //go through dimensions , calculating parameters for reshaping into ZYX array if needed + _info.resize(_ndims); + _world_matrix.resize(_ndims*4,0.0); + _voxel_matrix.resize(_ndims*4,0); + //_dir_cos.resize(_ndims*_ndims,0.0); + + for(int i=_ndims-1;i>=0;i--) + { + //_world_matrix[i*(_ndims+1)]=1.0; + //_voxel_matrix[i*(_ndims+1)]=1.0; + char dimname[MAX_NC_NAME]; + long dimlength; + //get dimensions info + ncdiminq(_mincid, mdims[i], dimname, &dimlength); + _info[i].name=dimname; + _info[i].length=dimlength; + _info[i].have_dir_cos=false; + int axis=-1; + unsigned int sz=0; + + if(!strcmp(dimname,MIxspace)) + { + _dims[0]=dimlength;axis=0; + _info[i].dim=dim_info::DIM_X; + _map_to_std[1]=i; + } else if(!strcmp(dimname,MIyspace)) { + _dims[1]=dimlength;axis=1; + _info[i].dim=dim_info::DIM_Y; + _map_to_std[2]=i; + } else if(!strcmp(dimname,MIzspace)) { + _dims[2]=dimlength;axis=2; + _info[i].dim=dim_info::DIM_Z; + _map_to_std[3]=i; + } else if(!strcmp(dimname,MIvector_dimension)) { + axis=-1; + _info[i].dim=dim_info::DIM_VEC; + _map_to_std[0]=i; + } else if(!strcmp(dimname,MItime)) { + axis=3; + _info[i].dim=dim_info::DIM_TIME; + _map_to_std[4]=i; + } else { + REPORT_ERROR ("Unknown dimension"); + _info[i].dim=dim_info::DIM_UNKNOWN; + } + + if(_info[i].dim!=dim_info::DIM_VEC) + { + //ncopts = 0; + int dimid = ncvarid(_mincid, dimname); + //ncopts = NC_VERBOSE | NC_FATAL; + if (dimid == MI_ERROR) continue; + + // Get dimension attributes + //ncopts = 0; + miattget1(_mincid, dimid, (char*)MIstep, NC_DOUBLE, &_info[i].step); + if(_info[i].step == 0.0) + _info[i].step = 1.0; + miattget1(_mincid, dimid, (char*)MIstart, NC_DOUBLE, &_info[i].start); + + if(_positive_directions && _info[i].step<0.0) + { + _info[i].start+=_info[i].step*(dimlength-1); + _info[i].step=-_info[i].step; + } + + if(miattget(_mincid, dimid, (char*)MIdirection_cosines, NC_DOUBLE, 3, &_info[i].dir_cos[0], NULL)!= MI_ERROR) + { + _info[i].have_dir_cos=true; + + /* Normalize the direction cosine */ + double len=sqrt(_info[i].dir_cos[0]*_info[i].dir_cos[0]+ + _info[i].dir_cos[1]*_info[i].dir_cos[1]+ + _info[i].dir_cos[2]*_info[i].dir_cos[2]); + + if(len>1e-6 && fabs(len-1.0)>1e-6) //TODO: use some epsiolon here? + { + for(int a=0;a<3;a++) + _info[i].dir_cos[a]/=len; + } + + } + // fill voxel matrix + _voxel_matrix[i*4+axis]=1; + } else { //vectors don't have spatial component! + _info[i].start=0; + _info[i].step=0.0; + _info[i].dir_cos[0]=_info[i].dir_cos[1]=_info[i].dir_cos[2]=0.0; + _info[i].have_dir_cos=false; + } + + //fill world matrix + for(int a=0;a<3;a++) + _world_matrix[i*4+a]=_info[i].dir_cos[a]*_info[i].step; + if(axis==3) //time + _world_matrix[i*4+3]=_info[i].step; + else + _world_matrix[i*4+3]=0.0; + } + //ncopts = NC_VERBOSE | NC_FATAL; + + // now let's find out the slice dimensions + int idmax = ncvarid(_mincid, MIimagemax); + _slice_dimensions=0; + if(idmax != MI_ERROR) + { + int nmax_dims; + int mmax_dims[MAX_VAR_DIMS]; + ncvarinq(_mincid, _imgid, NULL, NULL, &nmax_dims, mmax_dims, NULL); + if(nmax_dims>0) + _slice_dimensions=_ndims-nmax_dims; + } + + if(_slice_dimensions<=0) + { + if(_info[_ndims-1].dim==dim_info::DIM_VEC || _info[_ndims-1].dim==dim_info::DIM_TIME) + _slice_dimensions=std::min(_ndims,3); + else + _slice_dimensions=std::min(_ndims,2); + } + std::fill(_slab.begin(),_slab.end(),1); + _slab_len=1; + for(int i=0;i<_slice_dimensions;i++) + { + _slab[_ndims-i-1]=_info[_ndims-i-1].length; + _slab_len*=_info[_ndims-i-1].length; + } + } + + minc_1_reader::~minc_1_reader() + { + if(_have_temp_file) + remove(_tempfile.c_str()); + } + + minc_1_writer::minc_1_writer(): + _set_image_range(false),_set_slice_range(false), + _calc_min_max(true),_write_prepared(false) + { + } + + minc_1_writer::minc_1_writer(const minc_1_writer&): + _set_image_range(false),_set_slice_range(false), + _calc_min_max(true),_write_prepared(false) + { + } + + + void minc_1_writer::open(const char *path,const minc_info& inf,int slice_dimensions,nc_type datatype,int _s) + { +#ifndef WIN32 + ncopts = 0; +#endif + _info=inf; + //int mdims[MAX_VAR_DIMS]; + double vrange[2]; + _write_prepared=false; + + _mincid = micreate((char*)path, NC_CLOBBER/*|MI2_CREATE_V2*/); //TODO: add environment variable checking +#ifdef MINC2 + if (MI2_ISH5OBJ(_mincid)) { //micreate might create MINC2 file if environment variable is set + _minc2 = true; + } +#endif + _ndims=_info.size(); + _datatype=datatype; + _slice_dimensions=slice_dimensions; + _is_signed=_s; + fill(_map_to_std.begin(),_map_to_std.end(),-1); + for(int i=_ndims-1;i>=0;i--) + { + //just a precaution + switch(_info[i].dim) + { + case dim_info::DIM_X:_info[i].name=MIxspace;_map_to_std[1]=i;break; + case dim_info::DIM_Y:_info[i].name=MIyspace;_map_to_std[2]=i;break; + case dim_info::DIM_Z:_info[i].name=MIzspace;_map_to_std[3]=i;break; + case dim_info::DIM_TIME:_info[i].name=MItime;_map_to_std[4]=i;break; + default: + case dim_info::DIM_VEC:_info[i].name=MIvector_dimension;_map_to_std[0]=i;break; + //default: REPORT_ERROR("Unknown Dimension!"); + } + mdims[i]=ncdimdef(_mincid, _info[i].name.c_str(), _info[i].length); + if(_info[i].dim!=dim_info::DIM_VEC) + { + int dimid=micreate_std_variable(_mincid,(char*)_info[i].name.c_str(),NC_INT, 0, NULL); + miattputdbl(_mincid, dimid, (char*)MIstep,_info[i].step); + miattputdbl(_mincid, dimid, (char*)MIstart,_info[i].start); + + if(_info[i].have_dir_cos) + ncattput(_mincid, dimid, (char*)MIdirection_cosines,NC_DOUBLE, 3, _info[i].dir_cos); + } + } + _slab_len=1; + for(int i=0;i<_slice_dimensions;i++) + { + _slab[_ndims-i-1]=_info[_ndims-i-1].length; + _slab_len*=_info[_ndims-i-1].length; + } + + _icmax=_icmin=MI_ERROR; + //ncopts = NC_OPTS_VAL; + _imgid=micreate_std_variable(_mincid, (char*)MIimage, _datatype, _ndims, mdims); + _image_range[0]=DBL_MAX;_image_range[1]=-DBL_MAX; + + switch(_datatype) + { + case NC_DOUBLE: + vrange[0]=-DBL_MAX;vrange[1]=DBL_MAX; + _is_signed=1; + break; + case NC_FLOAT: + vrange[0]=-FLT_MAX;vrange[1]=FLT_MAX; + _is_signed=1; + break; + case NC_SHORT: + if(_is_signed) + { + vrange[0]=SHRT_MIN; + vrange[1]=SHRT_MAX; + } else { + vrange[0]=0;vrange[1]=USHRT_MAX; + } + break; + case NC_BYTE: + if(_is_signed) + { + vrange[0]=-128;vrange[1]=127; + } else { + vrange[0]=0;vrange[1]=255; + } + break; + case NC_INT: + if(_is_signed) + { + vrange[0]=INT_MIN;vrange[1]=INT_MAX; + }else{ + vrange[0]=0;vrange[1]=UINT_MAX; + } + break; + default:break; + }; + miattputstr(_mincid, _imgid, (char*)MIcomplete, (char*)MI_FALSE); + miattputstr(_mincid, _imgid, (char*)MIsigntype, (char*)(_is_signed?MI_SIGNED:MI_UNSIGNED)); + ncattput(_mincid, _imgid, (char*)MIvalid_range, NC_DOUBLE, 2, vrange); + miset_valid_range(_mincid, _imgid, vrange); + } + + void minc_1_writer::open(const char *path,const minc_1_base& imitate) + { + + open(path,imitate.info(),imitate.slice_dimensions(), + imitate.datatype(),imitate.is_signed()); + + copy_headers(imitate); + } + + void minc_1_writer::open(const char *path,const char *imitate_file) + { + minc_1_reader rdr; + //open minc file in metadate mode + rdr.open(imitate_file,false,true); + open(path,rdr); + //copy_headers(rdr); + } + + void minc_1_writer::setup_write_float() + { + _image_range[0]=DBL_MAX;_image_range[1]=-DBL_MAX; + + switch(_datatype) + { + case NC_DOUBLE: + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL); + + _set_image_range=true; + _set_slice_range=false; + break; + + case NC_FLOAT: + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL); + + _set_image_range=true; + _set_slice_range=false; + break; + + case NC_SHORT: + + _set_image_range=false; + _set_slice_range=true; + + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + break; + + case NC_BYTE: + _set_image_range=false; + _set_slice_range=true; + + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + break; + case NC_INT: + _set_image_range=false; + _set_slice_range=true; + + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + break; + + default: + break; + }; + ncendef(_mincid); + + if(_datatype==NC_DOUBLE || _datatype==NC_FLOAT) + { + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED); + miicv_setint(_icvid, MI_ICV_TYPE, NC_FLOAT); + miicv_setint(_icvid, MI_ICV_DO_NORM, true); + miicv_setint(_icvid, MI_ICV_USER_NORM, true); + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, -FLT_MAX); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, FLT_MAX); + _calc_min_max=true; + + } else { //do something smart here? + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED); + miicv_setint(_icvid, MI_ICV_TYPE, NC_FLOAT); + miicv_setint(_icvid, MI_ICV_DO_NORM, false); + //miicv_setint(_icvid, MI_ICV_USER_NORM, false); + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, -FLT_MAX); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, FLT_MAX); + _calc_min_max=true; + + } + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_FLOAT; + _write_prepared=true; + } + + void minc_1_writer::setup_write_double() + { + _image_range[0]=DBL_MAX;_image_range[1]=-DBL_MAX; + + switch(_datatype) + { + case NC_DOUBLE: + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL); + + _set_image_range=true; + _set_slice_range=false; + break; + + case NC_FLOAT: + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL); + + _set_image_range=true; + _set_slice_range=false; + break; + + case NC_SHORT: + + _set_image_range=false; + _set_slice_range=true; + + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + break; + + case NC_BYTE: + _set_image_range=false; + _set_slice_range=true; + + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + break; + case NC_INT: + _set_image_range=false; + _set_slice_range=true; + + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims); + break; + + default: + break; + }; + ncendef(_mincid); + + if(_datatype==NC_DOUBLE) + { + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED); + miicv_setint(_icvid, MI_ICV_TYPE, NC_DOUBLE); + miicv_setint(_icvid, MI_ICV_DO_NORM, true); + miicv_setint(_icvid, MI_ICV_USER_NORM, true); + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, -DBL_MAX); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, DBL_MAX); + _calc_min_max=true; + + } else if(_datatype==NC_FLOAT) { + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED); + miicv_setint(_icvid, MI_ICV_TYPE, NC_DOUBLE); + miicv_setint(_icvid, MI_ICV_DO_NORM, true); + miicv_setint(_icvid, MI_ICV_USER_NORM, true); + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, -DBL_MAX); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, DBL_MAX); + _calc_min_max=true; + + } else { //do something smart here? + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED); + miicv_setint(_icvid, MI_ICV_TYPE, NC_DOUBLE); + miicv_setint(_icvid, MI_ICV_DO_NORM, false); + //miicv_setint(_icvid, MI_ICV_USER_NORM, false); + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, -DBL_MAX); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, DBL_MAX); + _calc_min_max=true; + + } + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_DOUBLE; + _write_prepared=true; + } + + void minc_1_writer::setup_write_short(bool n) + { + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL); + _set_image_range=true; + _set_slice_range=false; + + ncendef(_mincid); + + miicv_setint(_icvid, MI_ICV_TYPE, NC_SHORT); + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED); + //miicv_setstr(_icvid, MI_ICV_SIGN, true); + /* Set range of values */ //TODO: set this to something sensible? + miicv_setint(_icvid, MI_ICV_VALID_MIN, SHRT_MIN); + miicv_setint(_icvid, MI_ICV_VALID_MAX, SHRT_MAX); + + /* No normalization so that pixels are scaled to the slice */ + miicv_setint(_icvid, MI_ICV_DO_NORM, false); + miicv_setint(_icvid, MI_ICV_DO_RANGE, false); + + miicv_attach(_icvid, _mincid, _imgid); + + _io_datatype=NC_SHORT; + _write_prepared=true; + } + + void minc_1_writer::setup_write_ushort(bool n) + { + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL); + _set_image_range=true; + _set_slice_range=false; + + ncendef(_mincid); + miicv_setint(_icvid, MI_ICV_TYPE, NC_SHORT); + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED); + miicv_setstr(_icvid, MI_ICV_SIGN, false); + + /* Set range of values */ //TODO: set this to something sensible? + miicv_setint(_icvid, MI_ICV_VALID_MIN, 0); + miicv_setint(_icvid, MI_ICV_VALID_MAX, USHRT_MAX); + + /* No normalization so that pixels are scaled to the slice */ + miicv_setint(_icvid, MI_ICV_DO_NORM, false); + miicv_setint(_icvid, MI_ICV_DO_RANGE, false); + + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_SHORT; + _write_prepared=true; + } + + void minc_1_writer::setup_write_byte(bool n) + { + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL); + _set_image_range=true; + _set_slice_range=false; + + ncendef(_mincid); + miicv_setint(_icvid, MI_ICV_TYPE, NC_BYTE); + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED); + + /* Set range of values */ //TODO: set this to something sensible? + miicv_setint(_icvid, MI_ICV_VALID_MIN, 0); + miicv_setint(_icvid, MI_ICV_VALID_MAX, UCHAR_MAX); + + /* No normalization so that pixels are scaled to the slice */ + miicv_setint(_icvid, MI_ICV_DO_NORM, false); + miicv_setint(_icvid, MI_ICV_DO_RANGE, false); + + miicv_attach(_icvid, _mincid, _imgid); + + _io_datatype=NC_BYTE; + _write_prepared=true; + } + + void minc_1_writer::setup_write_int(bool n) + { + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL); + _set_image_range=true; + _set_slice_range=false; + + ncendef(_mincid); + miicv_setint(_icvid, MI_ICV_TYPE, NC_INT); + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED); + + /* Set range of values */ //TODO: set this to something sensible? + miicv_setint(_icvid, MI_ICV_VALID_MIN, INT_MIN); + miicv_setint(_icvid, MI_ICV_VALID_MAX, INT_MAX); + + /* No normalization so that pixels are scaled to the slice */ + miicv_setint(_icvid, MI_ICV_DO_NORM, false); + miicv_setint(_icvid, MI_ICV_DO_RANGE, false); + + miicv_attach(_icvid, _mincid, _imgid); + + + _io_datatype=NC_INT; + _write_prepared=true; + } + + void minc_1_writer::setup_write_uint(bool n) + { + _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL); + _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL); + _set_image_range=true; + _set_slice_range=false; + + ncendef(_mincid); + miicv_setint(_icvid, MI_ICV_TYPE, NC_INT); + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED); + + /* Set range of values */ //TODO: set this to something sensible? + miicv_setint(_icvid, MI_ICV_VALID_MIN, 0); + miicv_setint(_icvid, MI_ICV_VALID_MAX, UINT_MAX); + + /* No normalization so that pixels are scaled to the slice */ + miicv_setint(_icvid, MI_ICV_DO_NORM, false); + miicv_setint(_icvid, MI_ICV_DO_RANGE, false); + + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_INT; + _write_prepared=true; + } + + void minc_1_reader::_setup_dimensions(void) + { + if(_metadate_only) + REPORT_ERROR("Minc file in metadate only mode!"); + if(_positive_directions) + { + /* We want to ensure that images have X, Y and Z dimensions in the + positive direction, giving patient left on left and for drawing from + bottom up. If we wanted patient right on left and drawing from + top down, we would set to MI_ICV_NEGATIVE. */ + + miicv_setint(_icvid, MI_ICV_DO_DIM_CONV, true); + //TODO: make sure to change only x,y,z conversions here + //miicv_setint(_icvid, MI_ICV_XDIM_DIR, 3); + //we want to convert only X,Y,Z dimensions if they are present + int num=(_map_to_std[1]>=0?1:0)+(_map_to_std[2]>=0?1:0)+(_map_to_std[3]>=0?1:0); + miicv_setint(_icvid, MI_ICV_NUM_IMGDIMS, num); + + if(_map_to_std[1]>=0) + { + miicv_setint(_icvid, MI_ICV_DIM_SIZE+_map_to_std[1],-1); + miicv_setint(_icvid, MI_ICV_XDIM_DIR, MI_ICV_POSITIVE); + } + + if(_map_to_std[2]>=0) + { + miicv_setint(_icvid, MI_ICV_DIM_SIZE+_map_to_std[2],-1); + miicv_setint(_icvid, MI_ICV_YDIM_DIR, MI_ICV_POSITIVE); + } + + if(_map_to_std[3]>=0) + { + miicv_setint(_icvid, MI_ICV_DIM_SIZE+_map_to_std[3],-1); + miicv_setint(_icvid, MI_ICV_ZDIM_DIR, MI_ICV_POSITIVE); + } + } + miicv_setint(_icvid, MI_ICV_DO_SCALAR, false); + } + + void minc_1_reader::setup_read_float(void) + { + if(_metadate_only) + REPORT_ERROR("Minc file in metadate only mode!"); + miicv_setint(_icvid, MI_ICV_TYPE, NC_FLOAT); + miicv_setint(_icvid, MI_ICV_DO_NORM, true); + miicv_setint(_icvid, MI_ICV_USER_NORM, true); + /* Make sure that any out of range values are mapped to lowest value + of type (for input only) */ + miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true); + + _setup_dimensions(); + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_FLOAT; + _read_prepared=true; + } + + void minc_1_reader::setup_read_double(void) + { + if(_metadate_only) + REPORT_ERROR("Minc file in metadate only mode!"); + miicv_setint(_icvid, MI_ICV_TYPE, NC_DOUBLE); + miicv_setint(_icvid, MI_ICV_DO_NORM, true); + miicv_setint(_icvid, MI_ICV_USER_NORM, true); + /* Make sure that any out of range values are mapped to lowest value + of type (for input only) */ + miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true); + + _setup_dimensions(); + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_DOUBLE; + _read_prepared=true; + } + + void minc_1_reader::setup_read_short(bool n) + { + if(_metadate_only) + REPORT_ERROR("Minc file in metadate only mode!"); + + miicv_setint(_icvid, MI_ICV_TYPE, NC_SHORT); + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED); + /* Set range of values */ + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, _image_range[0]); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, _image_range[1]); + + /* Do normalization so that all pixels are on same scale */ + miicv_setint(_icvid, MI_ICV_DO_NORM, true); + //miicv_setint(_icvid, MI_ICV_USER_NORM, true); + /* Make sure that any out of range values are mapped to lowest value + of type (for input only) */ + miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true); + + _setup_dimensions(); + + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_SHORT; + _read_prepared=true; + } + + void minc_1_reader::setup_read_ushort(bool n) + { + if(_metadate_only) + REPORT_ERROR("Minc file in metadate only mode!"); + + miicv_setint(_icvid, MI_ICV_TYPE, NC_SHORT); + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED); + /* Set range of values */ + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, _image_range[0]); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, _image_range[1]); + + /* Do normalization so that all pixels are on same scale */ + miicv_setint(_icvid, MI_ICV_DO_NORM, false); + //miicv_setint(_icvid, MI_ICV_USER_NORM, true); + /* Make sure that any out of range values are mapped to lowest value + of type (for input only) */ + miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true); + + _setup_dimensions(); + + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_SHORT; + _read_prepared=true; + } + + void minc_1_reader::setup_read_byte(bool n) + { + if(_metadate_only) + REPORT_ERROR("Minc file in metadate only mode!"); + miicv_setint(_icvid, MI_ICV_TYPE, NC_BYTE); + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED); + /* Set range of values */ + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, _image_range[0]); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, _image_range[1]); + + /* Do normalization so that all pixels are on same scale */ + miicv_setint(_icvid, MI_ICV_DO_NORM, true); + /* Make sure that any out of range values are mapped to lowest value + of type (for input only) */ + miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true); + + _setup_dimensions(); + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_BYTE; + _read_prepared=true; + } + + void minc_1_reader::setup_read_int(bool n) + { + if(_metadate_only) + REPORT_ERROR("Minc file in metadate only mode!"); + miicv_setint(_icvid, MI_ICV_TYPE, NC_INT); + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED); + /* Set range of values */ + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, _image_range[0]); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, _image_range[1]); + + /* Do normalization so that all pixels are on same scale */ + miicv_setint(_icvid, MI_ICV_DO_NORM, true); + /* Make sure that any out of range values are mapped to lowest value + of type (for input only) */ + miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true); + + _setup_dimensions(); + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_INT; + _read_prepared=true; + } + + void minc_1_reader::setup_read_uint(bool n) + { + if(_metadate_only) + REPORT_ERROR("Minc file in metadate only mode!"); + miicv_setint(_icvid, MI_ICV_TYPE, NC_INT); + miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED); + /* Set range of values */ + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, _image_range[0]); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, _image_range[1]); + + /* Do normalization so that all pixels are on same scale */ + miicv_setint(_icvid, MI_ICV_DO_NORM, true); + /* Make sure that any out of range values are mapped to lowest value + of type (for input only) */ + miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true); + + _setup_dimensions(); + miicv_attach(_icvid, _mincid, _imgid); + _io_datatype=NC_INT; + _read_prepared=true; + } + + void minc_1_reader::read(void* buffer) + { + if(!_read_prepared) + REPORT_ERROR("Not ready to read, use setup_read_XXXX"); + + miicv_get(_icvid, &_cur[0], &_slab[0], buffer); + } + + void minc_1_writer::write(void* buffer) + { + if(!_write_prepared) + REPORT_ERROR("Not ready to write, use setup_write_XXXX"); + + double r_min= DBL_MAX; //slab minimal value + double r_max=-DBL_MAX; //slab maximal value + //int irmin=0,irmax=0; + if(_calc_min_max ) + { + if(_io_datatype==NC_FLOAT) + {// + float *tmp=(float*)buffer; + for(int i=0;i<_slab_len;i++) + { + if(r_min>tmp[i]) r_min=tmp[i];//irmin=i; + if(r_max<tmp[i]) r_max=tmp[i];//irmax=i; + } + } else if(_io_datatype==NC_DOUBLE) { + double *tmp=(double*)buffer; + for(int i=0;i<_slab_len;i++) + { + if(r_min>tmp[i]) r_min=tmp[i];//irmin=i; + if(r_max<tmp[i]) r_max=tmp[i];//irmax=i; + } + } else if(_io_datatype==NC_SHORT && !_is_signed) { + unsigned short *tmp=(unsigned short *)buffer; + for(int i=0;i<_slab_len;i++) + { + if(r_min>tmp[i]) r_min=tmp[i]; + if(r_max<tmp[i]) r_max=tmp[i]; + } + } else if(_io_datatype==NC_SHORT && _is_signed) { + short *tmp=(short *)buffer; + for(int i=0;i<_slab_len;i++) + { + if(r_min>tmp[i]) r_min=tmp[i]; + if(r_max<tmp[i]) r_max=tmp[i]; + } + } else if(_io_datatype==NC_BYTE) { + unsigned char *tmp=(unsigned char *)buffer; + for(int i=0;i<_slab_len;i++) + { + if(r_min>tmp[i]) r_min=tmp[i]; + if(r_max<tmp[i]) r_max=tmp[i]; + } + } else if(_io_datatype==NC_INT && _is_signed) { + int *tmp=(int *)buffer; + for(int i=0;i<_slab_len;i++) + { + if(r_min>tmp[i]) r_min=tmp[i]; + if(r_max<tmp[i]) r_max=tmp[i]; + } + } else if(_io_datatype==NC_INT && !_is_signed) { + unsigned int *tmp=(unsigned int *)buffer; + for(int i=0;i<_slab_len;i++) + { + if(r_min>tmp[i]) r_min=tmp[i]; + if(r_max<tmp[i]) r_max=tmp[i]; + } + } + /* + if(r_min<-10000||r_max>10000) + { + std::cerr<<r_min<<":"<<r_max<<" "; + for(int i=0;i<4;i++) + std::cerr<<_cur[i]<<","; + std::cerr<<" "<<r_min<<":"<<r_max; + std::cerr<<" "<<irmin<<":"<<irmax<<" "<<_slab_len; + std::cerr<<std::endl; + }*/ + + if(_set_slice_range) + { + miicv_detach(_icvid); + miicv_setdbl(_icvid, MI_ICV_VALID_MIN, r_min); + miicv_setdbl(_icvid, MI_ICV_VALID_MAX, r_max); + miicv_attach(_icvid, _mincid, _imgid); + } + + if(_set_slice_range) + { + mivarput1(_mincid, _icmin, &_cur[0], NC_DOUBLE, NULL, &r_min); + mivarput1(_mincid, _icmax, &_cur[0], NC_DOUBLE, NULL, &r_max); + } + + if(_image_range[0]>r_min) _image_range[0]=r_min; + if(_image_range[1]<r_max) _image_range[1]=r_max; + } + miicv_put(_icvid, &_cur[0], &_slab[0], buffer); + } + + minc_1_writer::~minc_1_writer() + { + if(_set_image_range) + { + mivarput1(_mincid, _icmin, 0, NC_DOUBLE, NULL, &_image_range[0]); + mivarput1(_mincid, _icmax, 0, NC_DOUBLE, NULL, &_image_range[1]); + miset_valid_range(_mincid, _imgid, _image_range); + } + } + + void minc_1_writer::copy_headers(const minc_1_base& src) + { + + //code copied from mincresample + int nexcluded, excluded_vars[10]; + int varid; + + /* Create the list of excluded variables */ + nexcluded = 0; + //ncopts = 0; + + if ((varid=ncvarid(src.mincid(), MIxspace)) != MI_ERROR) + excluded_vars[nexcluded++] = varid; + if ((varid=ncvarid(src.mincid(), MIyspace)) != MI_ERROR) + excluded_vars[nexcluded++] = varid; + if ((varid=ncvarid(src.mincid(), MIzspace)) != MI_ERROR) + excluded_vars[nexcluded++] = varid; + if ((varid=ncvarid(src.mincid(), MItime)) != MI_ERROR) + excluded_vars[nexcluded++] = varid; + if ((varid=ncvarid(src.mincid(), MIimage)) != MI_ERROR) + excluded_vars[nexcluded++] = varid; + if ((varid=ncvarid(src.mincid(), MIimagemax)) != MI_ERROR) + excluded_vars[nexcluded++] = varid; + if ((varid=ncvarid(src.mincid(), MIimagemin)) != MI_ERROR) + excluded_vars[nexcluded++] = varid; + if ((varid=ncvarid(src.mincid(), "rootvariable")) != MI_ERROR) + excluded_vars[nexcluded++] = varid; +//#if MINC2 + if ((varid=ncvarid(src.mincid(), MIvector_dimension)) != MI_ERROR) + excluded_vars[nexcluded++] = varid; +//#endif /* MINC2 */ + //ncopts = NC_VERBOSE | NC_FATAL; + /* Copy all other variable definitions */ + micopy_all_var_defs(src.mincid(), _mincid, nexcluded, excluded_vars); + } + + //! append a line into minc history + void minc_1_writer::append_history(const char *append_history) + { + nc_type datatype; + int att_length; + //ncopts=0; + if ((ncattinq(_mincid, NC_GLOBAL, MIhistory, &datatype,&att_length) == MI_ERROR) || + (datatype != NC_CHAR)) + att_length = 0; + att_length += strlen(append_history) + 1; + char* str = new char[att_length]; + str[0] = '\0'; + miattgetstr(_mincid, NC_GLOBAL, (char*)MIhistory, att_length,str); + //ncopts=NC_VERBOSE | NC_FATAL; + strcat(str, append_history); + miattputstr(_mincid, NC_GLOBAL, (char*)MIhistory, str); + delete [] str; + } + + int minc_1_base::create_var_id(const char *varname) + { + int old_ncopts = ncopts; ncopts = 0; + int res=var_id(varname); + if(res==MI_ERROR) //need to create a variable + res=micreate_group_variable(_mincid,(char*)varname);//ncvardef(_mincid,varname,NC_INT,0,0); + if(res==MI_ERROR) //need to create a variable + res=ncvardef(_mincid,varname,NC_INT,0,0); + ncopts = old_ncopts; + return res; + } + + void minc_1_base::insert(const char *varname,const char *attname,double val) + { + ncattput(_mincid, create_var_id(varname),attname, NC_DOUBLE, 1, (void *) &val); + } + + void minc_1_base::insert(const char *varname,const char *attname,const char* val) + { + ncattput(_mincid, create_var_id(varname),attname, NC_CHAR, strlen(val) + 1, (void *) val); + } + + void minc_1_base::insert(const char *varname,const char *attname,const std::vector<double> &val) + { + ncattput(_mincid, create_var_id(varname),attname, NC_DOUBLE, val.size(), (void *) &val[0]); + } + + void minc_1_base::insert(const char *varname,const char *attname,const std::vector<int> &val) + { + ncattput(_mincid, create_var_id(varname),attname, NC_INT, val.size(), (void *) &val[0]); + } + + void minc_1_base::insert(const char *varname,const char *attname,const std::vector<short> &val) + { + ncattput(_mincid, create_var_id(varname),attname, NC_SHORT, val.size(), (void *) &val[0]); + } + + void minc_1_base::insert(const char *varname,const char *attname,const std::vector<unsigned char> &val) + { + ncattput(_mincid, create_var_id(varname),attname, NC_BYTE, val.size(), (void *) &val[0]); + } + +};
new file mode 100644 --- /dev/null +++ b/ezminc/minc_1_rw.h @@ -0,0 +1,391 @@ +#ifndef __MINC_1_RW__ +#define __MINC_1_RW__ + +#include <vector> +#include <string> + +#include "minc_io_exceptions.h" +//#include "minc_io_fixed_vector.h" + +#ifdef USE_MINC2 +#define MINC2 1 +#endif + +extern "C" { +#include <minc.h> +} + +#include <typeinfo> +#include <float.h> +#include <iostream> + +namespace minc +{ + + //! class for storing dimension information + struct dim_info + { + enum dimensions {DIM_UNKNOWN=0,DIM_X,DIM_Y,DIM_Z,DIM_TIME,DIM_VEC} ; + dim_info():length(0),step(0),start(0),have_dir_cos(false) + { + dir_cos[0]=dir_cos[1]=dir_cos[2]=0.0; + } + dim_info(int l, double sta,double spa,dimensions d,bool hd=false); + size_t length; + double step,start; + bool have_dir_cos; + double dir_cos[3]; + std::string name; + dimensions dim; + }; + + //! collection of dimensions describing a minc file + typedef std::vector<dim_info> minc_info; + + //! minc file rw base class + class minc_1_base + { + protected: + int _slab_len; + int _icvid; + std::vector<long> _cur,_slab; + size_t _slice_dimensions; + bool _last; + bool _positive_directions; + nc_type _datatype; + nc_type _io_datatype; + char _dimension_names[MAX_VAR_DIMS][MAX_NC_NAME]; + std::vector<double> _dir_cos; + long vcount[MAX_VAR_DIMS]; + std::vector<double> _world_matrix; + std::vector<int> _voxel_matrix; + int _ndims, mdims[MAX_VAR_DIMS]; + int _is_signed; + int _mincid, _imgid; + int _icmax,_icmin; + double _image_range[2]; + std::vector<long> _dims; + std::vector<int> _map_to_std; + minc_info _info; + bool _minc2; + public: + + //! get the minc handle + int mincid(void) const //this is not really a const ? + { + return _mincid; + } + + //! get the data type id (NC_BYTE,NC_INT etc) + nc_type datatype(void) const + { + return _datatype; + } + + //! byte size of the volume elements + unsigned int element_size(void) const + { + switch(_io_datatype) + { + case NC_FLOAT: return sizeof(float); + case NC_DOUBLE: return sizeof(double); + case NC_SHORT: return sizeof(short); + case NC_BYTE: return sizeof(char); + default:return 0;//maybe throw exception here? + } + } + + //! is data stored in signed format + bool is_signed(void) const + { + return _is_signed; + } + + //! constructor + minc_1_base(); + + //! destructor, closes minc file + virtual ~minc_1_base(); + + //! close the minc file + virtual void close(); + + //! is last slice was read? + bool last(void) const + { + return _last; + } + + //! go to the beginning of file + void begin(void) + { + fill(_cur.begin(),_cur.end(),0); + _last=false; + } + + //! advance to next slice + bool next_slice(void) + { + if(_last) return !_last; + + for(int i=_ndims-_slice_dimensions-1;i>=0;i--) + { + _cur[i]++; + if(_cur[i]<static_cast<long>(_info[i].length)) + break; + if(!i) + _last=true; + else + _cur[i]=0; + } + return !_last; + } + + //! slice length in elements + int slice_len(void) const + { + return _slab_len; + } + + //! number of dimensions + int dim_no(void) const + { + return _ndims; + } + + //! get the dimension information + const dim_info& dim(unsigned int n) const + { + if(n>=static_cast<unsigned int>(_ndims)) + REPORT_ERROR("Dimension is not defined"); + return _info[n]; + } + + //! get the pointer to the dimension description array + const minc_info& info(void) const + { + return _info; + } + + //! get the number of dimensions in one slice + int slice_dimensions(void) const + { + return _slice_dimensions; + } + + //! get the current slice index + const std::vector<long> & current_slice(void) const + { + return _cur; + } + + //! get the normalized dimensions sizes + //! ( 0 - vector_dimension, 1 - x, 2- y , 3 -z , 4 - time) + int ndim(int i) const + { + int j=_map_to_std[i]; + if(j>=0) return _info[j].length; + return 0; + } + //! get normalized dimension start coordinate (see ndim) + double nstart(int i) const + { + int j=_map_to_std[i]; + if(j>=0) return _info[j].start; + return 0.0; + } + + //! get normalized dimension spacing (see ndim) + double nspacing(int i) const + { + int j=_map_to_std[i]; + if(j>=0) return _info[j].step; + return 0.0; + } + + //! get normalized dimension direction cosine component (see ndim) + double ndir_cos(int i,int j) const + { + int k=_map_to_std[i]; + if(k>=0) return _info[k].dir_cos[j]; + return 0.0; + } + + //! check if a normalized dimension has direction cosine information + bool have_dir_cos(int i) const + { + int k=_map_to_std[i]; + if(k>=0) return _info[k].have_dir_cos; + return false; + } + + //! map file dimensions into normalized dimensions + int map_space(int i) + { + return _map_to_std[i]; + } + + //metadate info handling function: + //! read the minc history (:history attribute) + std::string history(void) const; + + //! retrive var id, if it exists, otherwise return MI_ERROR + int var_id(const char *var_name) const; + + //! get variable length + long var_length(const char *var_name) const; + + //! get variable length + long var_length(int var_id) const; + + //! read the number of variables + int var_number(void) const; + + //! get the variable name number no + std::string var_name(int no) const; + + //! get the number of attributes associated with variable + int att_number(const char *var_name) const; + + //! get the number of attributes associated with variable + int att_number(int var_no) const; + + //! get the attribute name, given the number + std::string att_name(const char *var_name,int no) const; + //! get the attribute name, given the number + std::string att_name(int varid,int no) const; + + //! get the string attribute value , given the name + std::string att_value_string(const char *var_name,const char *att_name) const; + //! get the string attribute value , given variable id + std::string att_value_string(int varid,const char *att_name) const; + + //! get the double attribute value , given the name + std::vector<double> att_value_double(const char *var_name,const char *att_name) const; + //! get the int attribute value , given the name + std::vector<int> att_value_int(const char *var_name,const char *att_name) const; + //! get the short attribute value , given the variable id + std::vector<short> att_value_short(const char *var_name,const char *att_name) const; + //! get the byte attribute value , given the variable id + std::vector<unsigned char> att_value_byte(const char *var_name,const char *att_name) const; + + //! get the double attribute value , given the variable id + std::vector<double> att_value_double(int varid,const char *att_name) const; + //! get the int attribute value , given the variable id + std::vector<int> att_value_int(int varid,const char *att_name) const; + //! get the short attribute value , given the variable id + std::vector<short> att_value_short(int varid,const char *att_name) const; + //! get the byte attribute value , given the variable id + std::vector<unsigned char> att_value_byte(int varid,const char *att_name) const; + + //! enquire about attribute data type + nc_type att_type(const char *var_name,const char *att_name) const; + //! enquire about attribute data type + nc_type att_type(int varid,const char *att_name) const; + + //! enquire about attribute length + int att_length(const char *var_name,const char *att_name) const; + //! enquire about attribute length + int att_length(int varid,const char *att_name) const; + + + //! return var_id for the given name (create one, if it doesn't exists) + int create_var_id(const char *varname); + + void insert(const char *varname,const char *attname,double val); + void insert(const char *varname,const char *attname,const char* val); + void insert(const char *varname,const char *attname,const std::vector<double> &val); + void insert(const char *varname,const char *attname,const std::vector<int> &val); + void insert(const char *varname,const char *attname,const std::vector<short> &val); + void insert(const char *varname,const char *attname,const std::vector<unsigned char> &val); + + + //! check if the file in MINC2 format + bool is_minc2(void) const + { + return _minc2; + } + + }; + + //! minc file reader + class minc_1_reader:public minc_1_base + { + protected: + bool _metadate_only; + std::string _tempfile; + bool _have_temp_file; + bool _read_prepared; + void _setup_dimensions(void); + + public: + //! copy constructor + minc_1_reader(const minc_1_reader&); + + //! default constructor + minc_1_reader(); + + //! destructor + virtual ~minc_1_reader(); + //! open a minc file + void open(const char *path,bool positive_directions=false,bool metadate_only=false,bool rw=false); + + //! read single slice + void read(void* slice); + //! setup reading in float format + void setup_read_float(void); + //! setup reading in double format + void setup_read_double(void); + //! setup reading in signed short format + void setup_read_short(bool normalized=false); + //! setup reading in unsigned short format + void setup_read_ushort(bool normalized=false); + //! setup reading in byte format + void setup_read_byte(bool normalized=false); + //! setup reading in int format + void setup_read_int(bool normalized=false); + //! setup reading in unsigned int format + void setup_read_uint(bool normalized=false); + }; + + //! minc file writer + class minc_1_writer:public minc_1_base + { + protected: + bool _set_image_range; + bool _set_slice_range; + bool _calc_min_max; + bool _write_prepared; + public: + void open(const char *path,const minc_info& inf,int slice_dimensions,nc_type datatype,int __signed=0); + void open(const char *path,const minc_1_base& imitate); + void open(const char *path,const char *imitate_file); + + void setup_write_float(void); + void setup_write_double(void); + void setup_write_short(bool normalize=false); + void setup_write_ushort(bool normalize=false); + void setup_write_byte(bool normalize=false); + void setup_write_int(bool normalize=false); + void setup_write_uint(bool normalize=false); + + //! copy header from another minc file + void copy_headers(const minc_1_base& src); + + //! append a line into minc history + void append_history(const char *append_history); + + + //! constructor + minc_1_writer(); + + minc_1_writer(const minc_1_writer&); + + //! destructor + virtual ~minc_1_writer(); + + //!write a single slice, size of the buffer should be more or equall to slab_len + void write(void* slice); + }; +}; +#endif //__PRIMITIVE_MINC_IO__
new file mode 100644 --- /dev/null +++ b/ezminc/minc_1_simple.h @@ -0,0 +1,290 @@ +#ifndef __MINC_1_SIMPLE_H__ +#define __MINC_1_SIMPLE_H__ + +#include "minc_1_rw.h" + +namespace minc +{ + + template <class T> class minc_input_iterator + { + protected: + mutable minc_1_reader* _rw; + std::vector<T> _buf; + std::vector<long> _cur; + bool _last; + size_t _count; + public: + + const std::vector<long>& cur(void) const + { + return _cur; + } + + + minc_input_iterator(const minc_input_iterator<T>& a):_rw(a._rw),_cur(a._cur),_last(a._last),_count(a._count) + { + } + + minc_input_iterator(minc_1_reader& rw):_rw(&rw),_last(false),_count(0) + { + } + + minc_input_iterator():_rw(NULL),_last(false),_count(0) + { + } + + void attach(minc_1_reader& rw) + { + _rw=&rw; + _last=false; + _count=0; + } + + bool next(void) + { + if(_last) return false; + _count++; + for(size_t i=static_cast<size_t>(_rw->dim_no()-1); + i>static_cast<size_t>(_rw->dim_no()-_rw->slice_dimensions()-1);i--) + { + _cur[i]++; + if(_cur[i]<static_cast<long>(_rw->dim(i).length)) + break; + if(i>static_cast<size_t>(_rw->dim_no()-_rw->slice_dimensions())) + _cur[i]=0; + else + { + //move to next slice + if(i==0) // the case when slice_dimensions==dim_no + { + _last=true; + _count=0; + break; + } + if(!_rw->next_slice()) + { + _last=true; + break; + } + _rw->read(&_buf[0]); + _cur=_rw->current_slice(); + _count=0; + break; + } + } + return !_last; + } + + bool last(void) + { + return _last; + } + + void begin(void) + { + _cur.resize(MAX_VAR_DIMS,0); + _buf.resize(_rw->slice_len()); + _count=0; + _rw->begin(); + _rw->read(&_buf[0]); + _cur=_rw->current_slice(); + } + + const T& value(void) const + { + return _buf[_count]; + } + }; + + template <class T> class minc_output_iterator + { + protected: + mutable minc_1_writer* _rw; + std::vector<T> _buf; + std::vector<long> _cur; + bool _last; + size_t _count; + public: + const std::vector<long>& cur(void) const + { + return _cur; + } + + minc_output_iterator(const minc_output_iterator<T>& a):_rw(a._rw),_cur(a._cur),_last(a._last),_count(a._count) + { + } + + minc_output_iterator(minc_1_writer& rw):_rw(&rw),_last(false),_count(0) + { + _buf.resize(rw.slice_len()); + } + + minc_output_iterator():_rw(NULL),_last(false),_count(0) + { + } + + void attach(minc_1_writer& rw) + { + _rw=&rw; + _last=false; + _count=0; + } + + ~minc_output_iterator() + { + if(_count && !_last) + _rw->write(&_buf[0]); + } + + bool next(void) + { + if(_last) return false; + _count++; + for(int i=_rw->dim_no()-1;i>(_rw->dim_no()-_rw->slice_dimensions()-1);i--) + { + _cur[i]++; + if(_cur[i]<static_cast<long>(_rw->dim(i).length)) + break; + if(i>(_rw->dim_no()-_rw->slice_dimensions())) + _cur[i]=0; + else + { + //write slice into minc file + _rw->write(&_buf[0]); + _count=0; + //move to next slice + if(i==0) // the case when slice_dimensions==dim_no + { + _last=true; + return false; + } + if(!_rw->next_slice()) + { + _last=true; + break; + } + _cur=_rw->current_slice(); + break; + } + } + return !_last; + } + + bool last(void) + { + return _last; + } + + void begin(void) + { + _buf.resize(_rw->slice_len()); + _cur.resize(MAX_VAR_DIMS,0); + _count=0; + _rw->begin(); + _cur=_rw->current_slice(); + } + + void value(const T& v) + { + _buf[_count]=v; + } + }; + + //! will attempt to laod the whole volume in T Z Y X V order into buffer, file should be prepared (setup_read_XXXX) + template<class T> void load_standard_volume(minc_1_reader& rw, T* volume) + { + std::vector<size_t> strides(MAX_VAR_DIMS,0); + size_t str=1; + for(size_t i=0;i<5;i++) + { + if(rw.map_space(i)<0) continue; + strides[rw.map_space(i)]=str; + str*=rw.ndim(i); + } + + minc_input_iterator<T> in(rw); + for(in.begin();!in.last();in.next()) + { + size_t address=0; + for(size_t i=0;i<static_cast<size_t>(rw.dim_no());i++) + address+=in.cur()[i]*strides[i]; + + volume[address]=in.value(); + } + } + + //! will attempt to save the whole volume in T Z Y X V order from buffer, file should be prepared (setup_read_XXXX) + template<class T> void save_standard_volume(minc_1_writer& rw, const T* volume) + { + std::vector<size_t> strides(MAX_VAR_DIMS,0); + size_t str=1; + for(size_t i=0;i<5;i++) + { + if(rw.map_space(i)<0) continue; + strides[rw.map_space(i)]=str; + str*=rw.ndim(i); + } + + minc_output_iterator<T> out(rw); + for(out.begin();!out.last();out.next()) + { + size_t address=0; + for(size_t i=0;i<static_cast<size_t>(rw.dim_no());i++) + address+=out.cur()[i]*strides[i]; + + out.value(volume[address]); + } + } + + //! will attempt to load the whole volume in Z Y X T V order into buffer, file should be prepared (setup_read_XXXX) + template<class T> void load_non_standard_volume(minc_1_reader& rw, T* volume) + { + std::vector<size_t> strides(MAX_VAR_DIMS,0); + size_t str=1; + const size_t dimorder[]={0,4,1,2,3}; + for(size_t i=0;i<5;i++) + { + if(rw.map_space(dimorder[i])<0|| !rw.ndim(dimorder[i]) ) continue; + strides[rw.map_space(dimorder[i])]=str; + str*=rw.ndim(dimorder[i]); + } + + minc_input_iterator<T> in(rw); + for(in.begin();!in.last();in.next()) + { + size_t address=0; + for(size_t i=0;i<rw.dim_no();i++) + address+=in.cur()[i]*strides[i]; + + volume[address]=in.value(); + } + } + + //! will attempt to save the whole volume in V T Z Y X order from buffer, file should be prepared (setup_read_XXXX) + template<class T> void save_non_standard_volume(minc_1_writer& rw, const T* volume) + { + std::vector<size_t> strides(MAX_VAR_DIMS,0); + size_t str=1; + const size_t dimorder[]={0,4,1,2,3}; + for(size_t i=0;i<5;i++) + { + if(rw.map_space(dimorder[i])<0 || !rw.ndim(dimorder[i]) ) continue; + strides[rw.map_space(dimorder[i])]=str; + str*=rw.ndim(dimorder[i]); + } + minc_output_iterator<T> out(rw); + for(out.begin();!out.last();out.next()) + { + size_t address=0; + for(size_t i=0;i<rw.dim_no();i++) + address+=out.cur()[i]*strides[i]; + + out.value(volume[address]); + } + } + + +};//minc + +#endif //__MINC_1_SIMPLE_H__
new file mode 100644 --- /dev/null +++ b/ezminc/minc_1_simple_rw.cpp @@ -0,0 +1,53 @@ +#include <iostream> +#include <math.h> + +#include "minc_1_simple_rw.h" + +namespace minc +{ + + const double minc_eps=1e-5; + + bool is_same(minc_1_reader& one,minc_1_reader& two,bool verbose) + { + if(one.dim_no()!=two.dim_no()) + { + if(verbose) + std::cerr<<"Unequal number of dimensions !"<<std::endl; + return false; + } + for(int j=0;j<5;j++) + { + if(one.ndim(j)!=two.ndim(j)) + { + if(verbose) + std::cerr<<"Unequal dimension sizes"<<std::endl; + return false; + } + + if(fabs(one.nstart(j)-two.nstart(j))>minc_eps) + { + if(verbose) + std::cerr<<"Unequal dimension sarts"<<std::endl; + return false; + } + + if(fabs(one.nspacing(j)-two.nspacing(j))>minc_eps) + { + if(verbose) + std::cerr<<"Unequal dimension steps"<<std::endl; + return false; + } + + for(int i=0;i<3;i++) + if(fabs(one.ndir_cos(j,i)-two.ndir_cos(j,i))>minc_eps) + { + if(verbose) + std::cerr<<"Unequal direction cosines"<<std::endl; + return false; + } + } + return true; + } + +};
new file mode 100644 --- /dev/null +++ b/ezminc/minc_1_simple_rw.h @@ -0,0 +1,334 @@ +#ifndef __MINC_1_SIMPLE_RW_H__ +#define __MINC_1_SIMPLE_RW_H__ + +#include "minc_1_simple.h" +#include "minc_io_simple_volume.h" +#include "minc_io_fixed_vector.h" +#include "minc_io_4d_volume.h" + +namespace minc +{ + + template<class T> void load_simple_volume(minc_1_reader& rw,simple_volume<T>& vol) + { + if(rw.ndim(1)<=0||rw.ndim(2)<=0||rw.ndim(3)<=0||rw.ndim(4)>0) + REPORT_ERROR("Need 3D minc file"); + + vol.resize(rw.ndim(1),rw.ndim(2),rw.ndim(3)); + + if(typeid(T)==typeid(unsigned char)) + { + rw.setup_read_byte(); + load_standard_volume(rw,vol.c_buf()); + } + else if(typeid(T)==typeid(int)) + { + rw.setup_read_int(); + load_standard_volume(rw,vol.c_buf()); + } + else if(typeid(T)==typeid(fixed_vec<3,float>)) + { + rw.setup_read_float(); + load_standard_volume<float>(rw,(float*)vol.c_buf()); + } + else if(typeid(T)==typeid(float)) + { + rw.setup_read_float(); + load_standard_volume(rw,vol.c_buf()); + } + else if(typeid(T)==typeid(fixed_vec<3,double>)) + { + rw.setup_read_double(); + load_standard_volume<double>(rw,(double*)vol.c_buf()); + } + else if(typeid(T)==typeid(double)) + { + rw.setup_read_double(); + load_standard_volume(rw,vol.c_buf()); + } else + REPORT_ERROR("Data type not supported for minc io"); + + //set coordinate transfer parameters + for(int i=0;i<3;i++) + { + vol.step()[i]=rw.nspacing(i+1); + vol.start()[i]=rw.nstart(i+1); + + if(rw.have_dir_cos(i+1)) + { + for(int j=0;j<3;j++) + vol.direction_cosines(i)[j]=rw.ndir_cos(i+1,j); + } else { + for(int j=0;j<3;j++) + vol.direction_cosines(i)[j]=(i==j?1.0:0.0); //identity + } + } + } + + template<class T> void save_simple_volume(minc_1_writer& rw,const simple_volume<T>& vol) + { + if(typeid(T)==typeid(unsigned char)) + { + rw.setup_write_byte(); + save_standard_volume(rw,vol.c_buf()); + } + else if(typeid(T)==typeid(int)) + { + rw.setup_write_int(); + save_standard_volume(rw,vol.c_buf()); + } + else if(typeid(T)==typeid(fixed_vec<3,float>)) + { + rw.setup_write_float(); + save_standard_volume<float>(rw,(float*)vol.c_buf()); + } + else if(typeid(T)==typeid(float)) + { + rw.setup_write_float(); + save_standard_volume(rw,vol.c_buf()); + } + else if(typeid(T)==typeid(fixed_vec<3,double>)) + { + rw.setup_write_double(); + save_standard_volume<double>(rw,(double*)vol.c_buf()); + } + else if(typeid(T)==typeid(double)) + { + rw.setup_write_double(); + save_standard_volume(rw,vol.c_buf()); + } + else + REPORT_ERROR("Data type not supported for minc io"); + } + + + template<class T> void load_4d_volume(minc_1_reader& rw,simple_4d_volume<T>& vol) + { + //if(rw.ndim(1)<=0||rw.ndim(2)<=0||rw.ndim(3)<=0||rw.ndim(4)<=0) + // REPORT_ERROR("Need 4D minc file"); + + vol.resize(rw.ndim(1),rw.ndim(2),rw.ndim(3),rw.ndim(4)>0?rw.ndim(4):1); //always assume 4 dimensions + + if(typeid(T)==typeid(unsigned char)) + rw.setup_read_byte(); + else if(typeid(T)==typeid(int)) + rw.setup_read_int(); + else if(typeid(T)==typeid(fixed_vec<3,float>)) + rw.setup_read_float(); + else if(typeid(T)==typeid(float)) + rw.setup_read_float(); + else if(typeid(T)==typeid(fixed_vec<3,double>)) + rw.setup_read_double(); + else if(typeid(T)==typeid(double)) + rw.setup_read_double(); + else + REPORT_ERROR("Data type not supported for minc io"); + + std::vector<size_t> strides(MAX_VAR_DIMS,0); + size_t str=1; + + for(size_t i=0;i<5;i++) //T is a special case + { + if(rw.map_space(i)<0) continue; + strides[rw.map_space(i)]=str; + str*=rw.ndim(i); + } + + if(rw.map_space(4)>=0) + strides[rw.map_space(4)]=0; //t dimension + + minc_input_iterator<T> in(rw); + for(in.begin();!in.last();in.next()) + { + size_t address=0; + size_t slice=0; + for(size_t i=0;i<rw.dim_no();i++) + { + if(strides[i]>0) + address+=in.cur()[i]*strides[i]; + else // + slice=in.cur()[i]; + } + vol.frame(slice).c_buf()[address]=in.value(); + } + + //set coordinate transfer parameters + for(int i=0;i<3;i++) + { + vol.step()[i]=rw.nspacing(i+1); + vol.start()[i]=rw.nstart(i+1); + + if(rw.have_dir_cos(i+1)) + { + for(int j=0;j<3;j++) + vol.direction_cosines(i)[j]=rw.ndir_cos(i+1,j); + } else { + for(int j=0;j<3;j++) + vol.direction_cosines(i)[j]=(i==j?1.0:0.0); //identity + } + } + if(rw.ndim(4)>0) + { + vol.t_start()=rw.nstart(4);//T + vol.t_step()=rw.nspacing(4);//T + } else { + vol.t_start()=0;//T + vol.t_step()=0;//T + } + } + + template<class T> void save_4d_volume(minc_1_writer& rw,const simple_4d_volume<T>& vol) + { + if(typeid(T)==typeid(unsigned char)) + rw.setup_write_byte(); + else if(typeid(T)==typeid(int)) + rw.setup_write_int(); + else if(typeid(T)==typeid(fixed_vec<3,float>)) + rw.setup_write_float(); + else if(typeid(T)==typeid(float)) + rw.setup_write_float(); + else if(typeid(T)==typeid(fixed_vec<3,double>)) + rw.setup_write_double(); + else if(typeid(T)==typeid(double)) + rw.setup_write_double(); + else + REPORT_ERROR("Data type not supported for minc io"); + + std::vector<size_t> strides(MAX_VAR_DIMS,0); + size_t str=1; + for(size_t i=0;i<4;i++)//T is a special + { + if(rw.map_space(i)<0) continue; + strides[rw.map_space(i)]=str; + str*=rw.ndim(i); + } + + if(rw.map_space(4)>=0) + strides[rw.map_space(4)]=0; //t dimension + + minc_output_iterator<T> out(rw); + for(out.begin();!out.last();out.next()) + { + size_t address=0; + size_t slice=0; + for(size_t i=0;i<rw.dim_no();i++) + { + if(strides[i]>0) + address+=out.cur()[i]*strides[i]; + else // + slice=out.cur()[i]; + } + out.value(vol.frame(slice).c_buf()[address]); + } + } + + bool is_same(minc_1_reader& one,minc_1_reader& two,bool verbose=true); + + template<class T> void load_minc_file(const char *file,simple_4d_volume<T>& vol) + { + minc_1_reader rdr; + rdr.open(file); + load_4d_volume(rdr,vol); + } + + template<class T> void generate_info(const simple_4d_volume<T>& vol,minc_info& info) + { + bool have_time=vol.frames()>1||vol.t_step()!=0.0; //assume that it is 3D file otherwise + + bool is_vector=false; + + if(typeid(T)==typeid(fixed_vec<3,float>)) { + is_vector=true; + } + + info.resize(3+(is_vector?1:0)+(have_time?1:0)); + + if(is_vector) + { + info[0].dim=dim_info::DIM_VEC; + info[0].length=3; + info[0].step=1; + + } + + for(int i=0;i<3;i++) + { + int ii=i+(is_vector?1:0); + info[ii].dim=dim_info::dimensions( dim_info::DIM_X+i); + + info[ii].length=vol.dim(i); + info[ii].step =vol.step()[i]; + info[ii].start =vol.start()[i]; + info[ii].have_dir_cos=true; + + for(int j=0;j<3;j++) + info[ii].dir_cos[j]=vol.direction_cosines(i)[j]; + } + + if(have_time) + { + info[3+(is_vector?1:0)].dim=dim_info::DIM_TIME; + info[3+(is_vector?1:0)].step=vol.t_step(); + info[3+(is_vector?1:0)].start=vol.t_start(); + info[3+(is_vector?1:0)].length=vol.frames(); + } + + } + + template<class T> void save_minc_file(const char *file,const simple_4d_volume<T>& vol, + const char* history=NULL,const minc_1_reader* original=NULL, + nc_type datatype=NC_NAT,bool is_signed=false) + { + minc_1_writer wrt; + //convert parameters to info + + if(typeid(T)==typeid(unsigned char)) + { + if(datatype==NC_NAT) datatype=NC_BYTE; + + } else if(typeid(T)==typeid(int)) { + if(datatype==NC_NAT) datatype=NC_INT; + + is_signed=true; + } else if(typeid(T)==typeid(unsigned int)) { + if(datatype==NC_NAT) datatype=NC_INT; + + is_signed=false; + } else if(typeid(T)==typeid(float)) { + if(datatype==NC_NAT) datatype=NC_FLOAT; + + is_signed=true; + } else if(typeid(T)==typeid(fixed_vec<3,float>)) { + if(datatype==NC_NAT) datatype=NC_FLOAT; + + is_signed=true; + } else if(typeid(T)==typeid(double)) { + if(datatype==NC_NAT) datatype=NC_DOUBLE; + + is_signed=true; + } else if(typeid(T)==typeid(fixed_vec<3,double>)) { + if(datatype==NC_NAT) datatype=NC_DOUBLE; + + is_signed=true; + } else + REPORT_ERROR("Unsupported data type!"); + + minc_info info; + generate_info<T>(vol,info); + + wrt.open(file,info,2,datatype,is_signed); + + if(original) + { + wrt.copy_headers(*original); + } + + if(history) + wrt.append_history(history); + + save_4d_volume(wrt,vol); + } + +}; + +#endif //__MINC_1_SIMPLE_RW_H__
new file mode 100644 --- /dev/null +++ b/ezminc/minc_io_4d_volume.h @@ -0,0 +1,157 @@ +#ifndef __MINC_IO_4D_VOLUME_H_ + +#include "minc_io_simple_volume.h" +#include <vector> +#include <cstring> + +namespace minc +{ + + //! simple 4D volume - collection of 3D volumes + template<class T> class simple_4d_volume + { + protected: + enum {ndims=3}; + typedef fixed_vec<ndims,int> idx; + typedef fixed_vec<ndims,double> vect; + + double _start_t; + double _step_t; + + void allocate(int n,const idx& sz) + { + _volumes.resize(n); + for(int i=0;i<n;i++) + { + _volumes[i].resize(sz); + } + } + + public: + typedef simple_volume<T> volume; + typedef std::vector<volume> volume_list; + + int dim(int i) const + { + return _volumes[0].dim(i); + } + + vect voxel_to_world(const idx& iii) const + { + return _volumes[0].voxel_to_world(iii); + } + + idx world_to_voxel(const vect& iii) const + { + return _volumes[0].world_to_voxel(iii); + } + + void resize(int x,int y,int z,int t) + { + allocate(t,IDX<int>(x,y,z)); + } + + //! number of temporal frames + size_t frames(void) const + { + return _volumes.size(); + } + + T& at(int x,int y,int z,int t) + { + return _volumes[t].at(x,y,z); + } + + const T& get(int x,int y,int z,int t) const + { + return _volumes[t].get(x,y,z); + } + + void set(int x,int y,int z,int t,const T& v) + { + _volumes[t].set(x,y,z,v); + } + + T& at(const idx& i,int t) + { + return _volumes[t].at(i); + } + + const T& get(const idx& i,int t) const + { + return _volumes[t].get(i); + } + + void set(const idx& i,int t,const T& v) + { + _volumes[t].set(i,v); + } + + volume& frame(int t) + { + return _volumes[t]; + } + + const volume& frame(int t) const + { + return _volumes[t]; + } + + vect& start(void) + { + return _volumes[0].start(); + } + + const vect& start(void) const + { + return _volumes[0].start(); + } + + vect& step(void) + { + return _volumes[0].step(); + } + + const vect& step(void) const + { + return _volumes[0].step(); + } + + vect& direction_cosines(int i) + { + return _volumes[0].direction_cosines(i); + } + + const vect& direction_cosines(int i) const + { + return _volumes[0].direction_cosines(i); + } + + double & t_step(void) + { + return _step_t; + } + + double t_step(void) const + { + return _step_t; + } + + double & t_start(void) + { + return _start_t; + } + + double t_start(void) const + { + return _start_t; + } + + protected: + + volume_list _volumes; + }; + +}; + +#endif //__MINC_IO_4D_VOLUME_H_
new file mode 100644 --- /dev/null +++ b/ezminc/minc_io_exceptions.h @@ -0,0 +1,58 @@ +/* ----------------------------- 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. +---------------------------------------------------------------------------- */ +#ifndef __EXCPETIONS_H__ +#define __EXCPETIONS_H__ + +#define REPORT_ERROR(MSG) throw minc::generic_error(__FILE__,__LINE__,MSG) + +namespace minc +{ + class generic_error + { + public: + const char *_file; + int _line; + const char *_msg; + int _code; + public: + + generic_error (const char *file, int line, const char *msg = "Error", int code = 0): + _file (file), _line (line), _msg (msg), _code (code) + { + // std::cerr<<"Exception created: "<<_file<<":"<<_line<<" "<<_msg<<std::endl; + } + + const char *file (void) const + { + return _file; + } + + const char *msg (void) const + { + return _msg; + } + + int line (void) const + { + return _line; + } + + int code (void) const + { + return _code; + } + }; +}; //minc +#endif //__EXCPETIONS_H__
new file mode 100644 --- /dev/null +++ b/ezminc/minc_io_fixed_vector.h @@ -0,0 +1,375 @@ +/* ----------------------------- 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. +---------------------------------------------------------------------------- */ +#ifndef __FIXED_VECTOR_H__ +#define __FIXED_VECTOR_H__ + +#include <limits> + +namespace minc +{ + //! fixed size array, which support arithmetic operations + template<int dim,class I=int> class fixed_vec + { + protected: + I c[dim]; + public: + //! default constructor, does nothing (i.e data is uninitilized) + fixed_vec() {} + + //! constructor which sets all the elements to the same value + explicit fixed_vec(I v) + { + for(unsigned int i=0;i<dim;i++) + c[i]=v; + } + + //! constructor which sets all the elements to be a copy of C-array + explicit fixed_vec(const I* v) + { + for(unsigned int i=0;i<dim;i++) + c[i]=v[i]; + } + + //! conversion to the C array + I* c_buf() + { + return c; + } + + //! conversion to const C array + const I* c_buf() const + { + return c; + } + + //! element access operator + I& operator[](int i) + { +#ifdef _INDEX_CHECK + if(i<0 || i>=dim) REPORT_ERROR("Index out of bounds"); +#endif //_INDEX_CHECK + return c[i]; + } + //! const element access operator + I operator[](int i) const + { +#ifdef _INDEX_CHECK + if(i<0 || i>=dim) REPORT_ERROR("Index out of bounds"); +#endif //_INDEX_CHECK + return c[i]; + } + + //! const element access operator + I get(int i) + { + return (*this)[i]; + } + //! element writing operator + void set(int i,I v) + { + (*this)[i]=v; + } + + //! find a maximum of elements + I max(void) const + { + I s=std::numeric_limits < I >::min ();; + for(unsigned int i=0;i<dim;i++) + if(c[i]>s) s=c[i]; + return s; + } + + //! find a minimum of elements + I min(void) const + { + I s=std::numeric_limits < I >::max ();; + for(unsigned int i=0;i<dim;i++) + if(c[i]<s) s=c[i]; + return s; + } + + //! calculate sum of all elements + I sum(void) const + { + I s=0; + for(unsigned int i=0;i<dim;i++) + s+=c[i]; + return s; + } + + //! modulus squared + I mod2(void) const + { + I s=0; + for(unsigned int i=0;i<dim;i++) + s+=c[i]*c[i]; + return s; + } + + //! volume (product of all elements) + I vol(void) const + { + I s=1; + for(unsigned int i=0;i<dim;i++) + s*=c[i]; + return s; + } + + //! \name fixed_vec arithmentic operations + //@{ + fixed_vec<dim,I>& operator *=(const fixed_vec<dim,I>& b) + { + for(unsigned int i=0;i<dim;i++) + c[i]*=b[i]; + return *this; + } + + fixed_vec<dim,I>& operator *=(const I b) + { + for(unsigned int i=0;i<dim;i++) + c[i]*=b; + return *this; + } + + fixed_vec<dim,I>& operator +=(const fixed_vec<dim,I>& b) + { + for(unsigned int i=0;i<dim;i++) + c[i]+=b[i]; + return *this; + } + + fixed_vec<dim,I>& operator -=(const fixed_vec<dim,I>& b) + { + for(unsigned int i=0;i<dim;i++) + c[i]-=b[i]; + return *this; + } + + fixed_vec<dim,I>& operator /=(const fixed_vec<dim,I>& b) + { + for(unsigned int i=0;i<dim;i++) + c[i]/=b[i]; + return *this; + } + + fixed_vec<dim,I>& operator /=(const I b) + { + for(unsigned int i=0;i<dim;i++) + c[i]/=b; + return *this; + } + + fixed_vec<dim,I> operator /(const I b) + { + fixed_vec<dim,I> tmp; + for(unsigned int i=0;i<dim;i++) + tmp[i]=c[i]/b; + return tmp; + } + + fixed_vec<dim,I> operator *(const I b) + { + fixed_vec<dim,I> tmp; + for(unsigned int i=0;i<dim;i++) + tmp[i]=c[i]*b; + return tmp; + } + + fixed_vec<dim,I> operator -(const fixed_vec<dim,I>& b) + { + fixed_vec<dim,I> tmp; + for(unsigned int i=0;i<dim;i++) + tmp[i]=c[i]-b[i]; + return tmp; + } + + fixed_vec<dim,I> operator +(const fixed_vec<dim,I>& b) + { + fixed_vec<dim,I> tmp; + for(unsigned int i=0;i<dim;i++) + tmp[i]=c[i]+b[i]; + return tmp; + } + + //@} + + //! assignement operator, copies contents + fixed_vec<dim,I>& operator=(const fixed_vec<dim,I>& b) + { + for(unsigned int i=0;i<dim;i++) + c[i]=b[i]; + return *this; + } + + //! assignement operator, copies contents, assumes that b have at least this size + fixed_vec<dim,I>& operator=(const I* b) + { + for(unsigned int i=0;i<dim;i++) + c[i]=b[i]; + return *this; + } + + //! assignement operator, sets all elements to the same value + fixed_vec<dim,I>& operator=(const I b) + { + for(unsigned int i=0;i<dim;i++) c[i]=b; + return *this; + } + + //! inequality operator, does element wise equality check + bool operator!=(const fixed_vec<dim,I>& b) const + { + for(int i=0;i<dim;i++) if(c[i]!=b[i]) return true; + return false; + } + + //! equality operator, does element wise equality check + bool operator==(const fixed_vec<dim,I>& b) const + { + for(int i=0;i<dim;i++) if(c[i]!=b[i]) return false; + return true; + } + + //! reverse the order of elements + void reverse(void) + { + for(int i=0;i<dim/2;i++) + { + I tmp=c[i]; + c[i]=c[dim-i-1]; + c[dim-i-1]=tmp; + } + } + }; + + //! element wise division + template<int dim,class I> fixed_vec<dim,I> operator/(const fixed_vec<dim,I> &l,const fixed_vec<dim,I> &r) + { + fixed_vec<dim,I> out=l; + out/=r; + return out; //this is not effecient - no return value optimisation + } + + //! element wise multiplication + template<int dim,class I> fixed_vec<dim,I> operator*(const fixed_vec<dim,I> &l,const fixed_vec<dim,I> &r) + { + fixed_vec<dim,I> out=l; + out*=r; + return out; //this is not effecient - no return value optimisation + } + + //! element wise addition + template<int dim,class I> fixed_vec<dim,I> operator+(const fixed_vec<dim,I> &l,const fixed_vec<dim,I> &r) + { + fixed_vec<dim,I> out=l; + out+=r; + return out; //this is not effecient - no return value optimisation + } + + //! element wise subtraction + template<int dim,class I> fixed_vec<dim,I> operator-(const fixed_vec<dim,I> &l,const fixed_vec<dim,I> &r) + { + fixed_vec<dim,I> out=l; + out-=r; + return out; //this is not effecient - no return value optimisation + } + + //! devide all elements by a value + template<int dim,class I> fixed_vec<dim,I> operator/(const fixed_vec<dim,I> &l,I r) + { + fixed_vec<dim,I> out=l; + out/=r; + return out; //this is not effecient - no return value optimisation + } + + //! multiply all elements by a value + template<int dim,class I> fixed_vec<dim,I> operator*(const fixed_vec<dim,I> &l,I r) + { + fixed_vec<dim,I> out=l; + out*=r; + return out; //this is not effecient - no return value optimisation + } + + //! add a value to all elements + template<int dim,class I> fixed_vec<dim,I> operator+(const fixed_vec<dim,I> &l,I r) + { + fixed_vec<dim,I> out=l; + out+=r; + return out; //this is not effecient - no return value optimisation + } + + //! subtract a value from all elements + template<int dim,class I> fixed_vec<dim,I> operator-(const fixed_vec<dim,I> &l,I r) + { + fixed_vec<dim,I> out=l; + out-=r; + return out; //this is not effecient - no return value optimisation + } + + //! create 1d fixed_vec + template<class I> fixed_vec<1,I> IDX(I i) + { + fixed_vec<1,I> d; + d[0]=i; + return d; + } + + //! create 2d fixed_vec + template<class I> fixed_vec<2,I> IDX(I i,I j) + { + fixed_vec<2,I> d; + d[0]=i; + d[1]=j; + return d; + } + + //! create 3d fixed_vec + template<class I> fixed_vec<3,I> IDX(I i,I j,I k) + { + fixed_vec<3,I> d; + d[0]=i; + d[1]=j; + d[2]=k; + return d; + } + + //! create 4d fixed_vec + template<class I> fixed_vec<4,I> IDX(I i,I j,I k,I l) + { + fixed_vec<3,I> d; + d[0]=i; + d[1]=j; + d[2]=k; + d[3]=l; + return d; + } + + //!average value of a vector + template<class T,int d>T AVG(const fixed_vec<d,T> &v) + { + return v.sum()/d; + } + + //!dot product of two vectors + template<class T,int d>T dot(const fixed_vec<d,T> &v1,const fixed_vec<d,T> &v2) + { + T val=0; + for(int i=0;i<d;i++) val+=v1[i]*v2[i]; + return val; + } + +}; + +#endif //__FIXED_VECTOR_H__
new file mode 100644 --- /dev/null +++ b/ezminc/minc_io_simple_volume.h @@ -0,0 +1,554 @@ +#ifndef __SIMPLE_VOLUME_H__ +#define __SIMPLE_VOLUME_H__ + +#include "minc_io_exceptions.h" +#include "minc_io_fixed_vector.h" +#include <string.h> +#include <math.h> + +namespace minc +{ + //! very simple 3D volume, initially created as an example but became usable + template<class T> class simple_volume + { + public: + + enum {ndims=3}; + typedef fixed_vec<ndims,int> idx; + typedef fixed_vec<ndims,double> vect; + + protected: + + T * _vol; //! the volume itself + idx _size; //! dimension sizes + idx _stride; //! used internally + int _count; //! total number of voxels + bool _free_memory; //! should the array be freed + + vect _step,_start; //! conversion to wold coordinates + vect _direction_cosines[3]; + + + void _allocate(T* data=NULL) + { + _stride[0]=1; + int total=_size[0]; + for(int i=1;i<ndims;i++) + { + _stride[i]=_size[i-1]*_stride[i-1]; + total*=_size[i]; + } + _count=total; + if(data) + { + _vol=data; + _free_memory=false; + } else { + _vol=new T[total]; + _free_memory=true; + } + + _step=IDX<double>(1.0,1.0,1.0); + _start=IDX<double>(0.0,0.0,0.0); + + _direction_cosines[0]=IDX<double>(1.0,0.0,0.0); + _direction_cosines[1]=IDX<double>(0.0,1.0,0.0); + _direction_cosines[2]=IDX<double>(0.0,0.0,1.0); + + } + + public: + vect& start(void) + { + return _start; + } + + const vect& start(void) const + { + return _start; + } + + vect& step(void) + { + return _step; + } + + + const vect& step(void) const + { + return _step; + } + + vect& direction_cosines(int i) + { + return _direction_cosines[i]; + } + + const vect& direction_cosines(int i) const + { + return _direction_cosines[i]; + } + + operator T*() + { + return _vol; + } + + T* c_buf() + { + return _vol; + } + + const T* c_buf() const + { + return _vol; + } + + int c_buf_size() const + { + return _count; + } + + explicit simple_volume(const int* dims):_vol(0),_size(dims) + { + _allocate(); + } + + simple_volume(const simple_volume<T>& a,bool copy_data=true):_vol(0) + { + for(int i=0;i<ndims;i++) + _size[i]=a._size[i]; + _allocate(); + + if(copy_data) + { + memmove(_vol,a._vol,_size[0]*_size[1]*_size[2]*sizeof(T)); + } + + _step=a._step; + _start=a._start; + for(int i=0;i<ndims;i++) + _direction_cosines[i]=a._direction_cosines[i]; + } + + simple_volume(int sx,int sy,int sz):_vol(0) + { + _size=IDX(sx,sy,sz); + _allocate(); + } + + explicit simple_volume(const idx& s):_vol(0) + { + _size=s; + _allocate(); + } + + simple_volume():_vol(0) + { + for(int i=0;i<ndims;i++) + { + _size[i]=0; + _step[i]=0.0; + _start[i]=0.0; + + _direction_cosines[i]=IDX<double>(0.0,0.0,0.0); + _direction_cosines[i][i]=1.0; + } + } + + bool empty(void) const + { + return !_size[0]||!_vol; + } + + void resize(int sx,int sy,int sz) + { + if( _size[0]==sx && _size[1]==sy && _size[2]==sz ) + return; + + if(_vol && _free_memory) + delete [] _vol; + _size=IDX(sx,sy,sz); + _allocate(); + } + + void resize(const idx& s) + { + if(_size==s) return; + + if(_vol&&_free_memory) + delete [] _vol; + _size=s; + _allocate(); + } + + virtual ~simple_volume() + { + if(_vol && _free_memory) + delete [] _vol; + } + + T& operator()(int x,int y,int z) + { + return _vol[x+y*_stride[1]+z*_stride[2]]; + } + + T& operator()(const idx& i) + { + return _vol[dot(i,_stride)]; + } + + const T& operator()(int x,int y,int z) const + { + return get(x,y,z); + } + + const T& operator()(const idx& i) const + { + return get(i); + } + + const T& get(int x,int y,int z) const + { + return _vol[x+y*_stride[1]+z*_stride[2]]; + } + + const T& get(const idx& i) const + { + return _vol[dot(i,_stride)]; + } + + const T& safe_get(int x,int y,int z) const + { + check_index(x,y,z); + return _vol[x+y*_stride[1]+z*_stride[2]]; + } + + const T& safe_get(idx i) const + { + check_index(i); + return get(i); + } + + //trilinear intrpolation + double interpolate(float _x,float _y,float _z) const + { + int x=floor(_x); + int y=floor(_y); + int z=floor(_z); + + float dx=_x-x; + float dy=_y-y; + float dz=_z-z; + + if(x<0) x=-x; + if(y<0) y=-y; + if(z<0) z=-z; + + if(x>=(_size[0]-1)) x=_size[0]*2-3-x; + if(y>=(_size[1]-1)) y=_size[1]*2-3-y; + if(z>=(_size[2]-1)) z=_size[2]*2-3-z; + + //trilinear intrpolation + return (1.0-dx)*(1.0-dy)*(1.0-dz)*get(x,y,z)+ + + dx*(1.0-dy)*(1.0-dz)*get(x+1,y,z)+ + (1.0-dx)*dy*(1.0-dz)*get(x,y+1,z)+ + (1.0-dx)*(1.0-dy)*dz*get(x,y,z+1)+ + + dx*(1.0-dy)*dz*get(x+1,y,z+1)+ + (1.0-dx)*dy*dz*get(x,y+1,z+1)+ + dx*dy*(1.0-dz)*get(x+1,y+1,z)+ + + dx*dy*dz*get(x+1,y+1,z+1); + } + + T set(int x,int y,int z, const T&v) + { + return _vol[x+y*_stride[1]+z*_stride[2]]=v; + } + + T set(const idx& i, const T&v) + { + return _vol[dot(i,_stride)]=v; + } + + int dim(int i) const + { + return _size[i]; + } + + const int* dims() const + { + return _size.c_buf(); + } + + const idx& size() const + { + return _size; + } + + void extract_subvolume(simple_volume<T>& dst,const idx& s, const idx& f) const + { + for(int k=s[2];k<f[2];k++) + for(int j=s[1];j<f[1];j++) + for(int i=s[0];i<f[0];i++) + dst(i,j,k)=get(i,j,k); + } + + void check_index(int &ii,int &jj,int &kk) const + { + if(ii<0) ii=-ii; + if(jj<0) jj=-jj; + if(kk<0) kk=-kk; + + if(ii>=dim(0)) ii=2*dim(0)-ii-1; + if(jj>=dim(1)) jj=2*dim(1)-jj-1; + if(kk>=dim(2)) kk=2*dim(2)-kk-1; + } + + void check_index(idx& iii) + { + for(int i=0;i<3;i++) + { + if(iii[i]<0) + iii[i]=-iii[i]; + + if(iii[i]>=dim(i)) + iii[i]=2*dim(i)-iii[i]-1; + + } + } + + bool hit(int ii,int jj,int kk) const + { + if(ii<0) return false; + if(jj<0) return false; + if(kk<0) return false; + + if(ii>=dim(0)) return false; + if(jj>=dim(1)) return false; + if(kk>=dim(2)) return false; + return true; + } + + bool hit(const idx iii) const + { + for(int i=0;i<3;i++) + { + if(iii[i]<0) return false; + if(iii[i]>=dim(i)) return false; + } + return true; + } + + simple_volume<T>& operator+=(const simple_volume<T>& a) + { + for(int i=0;i<ndims;i++) + if(_size[i]!=a._size[i]) + REPORT_ERROR("Dimensions are different"); + + for(int i=0;i<_count;i++) + _vol[i]+=a._vol[i]; + return *this; + } + + simple_volume<T>& operator+=(const T& a) + { + for(int i=0;i<_count;i++) + _vol[i]+=a; + return *this; + } + + simple_volume<T>& operator-=(const simple_volume<T>& a) + { + if(_size!=a._size) + REPORT_ERROR("Dimensions are different"); + + for(int i=0;i<_count;i++) + _vol[i]-=a._vol[i]; + return *this; + } + + simple_volume<T>& operator-=(const T& a) + { + for(int i=0;i<_count;i++) + _vol[i]-=a; + return *this; + } + + simple_volume<T>& operator*=(const simple_volume<T>& a) + { + for(int i=0;i<ndims;i++) + if(_size[i]!=a._size[i]) + REPORT_ERROR("Dimensions are different"); + + for(int i=0;i<_count;i++) + _vol[i]*=a._vol[i]; + return *this; + } + + simple_volume<T>& operator*=(const T& a) + { + for(int i=0;i<_count;i++) + _vol[i]*=a; + return *this; + } + + simple_volume<T>& operator/=(const simple_volume<T>& a) + { + if(_size!=a._size()) + REPORT_ERROR("Dimensions are different"); + + for(int i=0;i<_count;i++) + _vol[i]/=a._vol[i]; + return *this; + } + + simple_volume<T>& operator/=(const T& a) + { + for(int i=0;i<_count;i++) + _vol[i]/=a; + return *this; + } + + simple_volume& operator=(const simple_volume<T>&a) + { + resize(a.dim(0),a.dim(1),a.dim(2)); + + memmove(_vol, a._vol, _count*sizeof(T)); + + _step=a._step; + _start=a._start; + + for(int i=0;i<ndims;i++) + _direction_cosines[i]=a._direction_cosines[i]; + + return *this; + } + + simple_volume& operator=(const T&a) + { + for(int i=0;i<_count;i++) + _vol[i]=a; + return *this; + } + + void weighted_add(const simple_volume<T>&a, double w) + { + if(_size!=a._size) + REPORT_ERROR("Dimensions are different"); + + for(int i=0;i<_count;i++) + _vol[i]+=a._vol[i]*w; + } + + vect voxel_to_world(const idx& iii) const + { + vect ret=IDX<double>(0.0,0.0,0.0); + for(int i=0;i<ndims;i++) + { + for(int j=0;j<ndims;j++) + ret[i]+=(_step[j]*iii[j]+_start[j])*_direction_cosines[j][i]; + } + return ret; + } + + vect world_to_voxel_c(const vect& iii) const + { + vect ret=IDX<double>(0.0,0.0,0.0); + for(int i=0;i<ndims;i++) + { + for(int j=0;j<ndims;j++) + ret[i]+=((iii[j]-_start[j])/_step[j])*_direction_cosines[i][j]; //transpose! + } + return ret; + } + + idx world_to_voxel(const vect& iii) const + { + vect ret=world_to_voxel_c(iii); + + idx r; + + for(int i=0;i<ndims;i++) + r[i]=floor(ret[i]+0.5); + + return r; + } + + //!use provided buffer for storage + void assign(const idx& s,T* array) + { + _size=s; + allocate(array); + } + + }; + + //! remove (unpad) or add padding as needed, volume will be centered + template<class T>void pad_volume(const simple_volume<T> &src,simple_volume<T> &dst, const T& fill) + { + fixed_vec<3,int> sz1=src.size(); + fixed_vec<3,int> sz2=dst.size(); + fixed_vec<3,int> d=sz2-sz1; + fixed_vec<3,int> i; + + d/=2;//offset + + for( i[2]=0;i[2]<sz2[2];i[2]++) + for( i[1]=0;i[1]<sz2[1];i[1]++) + for( i[0]=0;i[0]<sz2[0];i[0]++) + { + fixed_vec<3,int> j= i-d; + + if( src.hit(j)) + dst.set(i,src.get(j)); + else + dst.set(i,fill); + } + } + + template<class T>void volume_min_max(const simple_volume<T>& v,T &_min,T &_max) + { + _min=_max=v.c_buf()[0]; + + for(int i=0;i<v.c_buf_size();i++) + { + if(isnan(v.c_buf()[i]) || isinf(v.c_buf()[i])) + continue; + + if(v.c_buf()[i]>_max) _max=v.c_buf()[i]; + else if(v.c_buf()[i]<_min) _min=v.c_buf()[i]; + } + } + + template<class T> void volume_min_max(const simple_volume<T>& v,const simple_volume<unsigned char>& mask,T &_min,T &_max) + { + if(v.size()!=mask.size()) + REPORT_ERROR("Volume size mismatch"); + + _min=1e10; + _max=-1e10; + + for(int i=0;i<v.c_buf_size();i++) + { + if(mask.c_buf()[i]) + { + if(isnan(v.c_buf()[i]) || isinf(v.c_buf()[i])) + continue; + if(v.c_buf()[i]>_max) _max=v.c_buf()[i]; + else if(v.c_buf()[i]<_min) _min=v.c_buf()[i]; + } + } + } + + + typedef simple_volume<float> minc_float_volume; + typedef simple_volume<fixed_vec<3,float> > minc_grid_volume; + typedef simple_volume<unsigned char> minc_byte_volume; + +}; //minc + + +#endif // __SIMPLE_VOLUME_H__