00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include "cmtkXformIO.h"
00034
00035 #include <Base/cmtkDeformationField.h>
00036
00037 #ifdef CMTK_BUILD_NRRD_TEEM
00038 # include <teem/nrrd.h>
00039 #else
00040 # ifdef CMTK_BUILD_NRRD
00041 # include <NrrdIO.h>
00042 # endif
00043 #endif
00044
00045 #ifdef CMTK_BUILD_NRRD
00046
00047 namespace
00048 cmtk
00049 {
00050
00053
00054 Xform::SmartPtr
00055 XformIO::ReadNrrd( const char* path, const bool )
00056 {
00057 DeformationField::SmartPtr dfield( NULL );
00058 try
00059 {
00060 Nrrd *nrrd = nrrdNew();
00061 if ( nrrdLoad( nrrd, path, NULL ) )
00062 throw biffGetDone(NRRD);
00063
00064 if ( nrrd->dim != 4 )
00065 {
00066 StdErr << "ERROR: deformation field must be stored as 4-dimensional Nrrd.\n";
00067 return dfield;
00068 }
00069
00070 if ( nrrd->axis[0].kind != nrrdKindVector )
00071 {
00072 StdErr << "ERROR: deformation field vectors in Nrrd must be stored together.\n";
00073 return dfield;
00074 }
00075
00076 if ( nrrd->axis[0].size != 3 )
00077 {
00078 StdErr << "ERROR: deformation field vectors in Nrrd must be three dimensional.\n";
00079 return dfield;
00080 }
00081
00082 NrrdAxisInfo* nrrdSpaceAxes = nrrd->axis+1;
00083 const int dims[3] = { nrrdSpaceAxes[0].size, nrrdSpaceAxes[1].size, nrrdSpaceAxes[2].size };
00084
00085
00086 double spacing[3] = { 1, 1, 1 };
00087 for ( size_t ax = 0; ax < 3; ++ax )
00088 {
00089 switch ( nrrdSpacingCalculate( nrrd, ax+1, spacing+ax, nrrd->axis[ax+1].spaceDirection ) )
00090 {
00091 case nrrdSpacingStatusScalarNoSpace:
00092 break;
00093 case nrrdSpacingStatusDirection:
00094 break;
00095 case nrrdSpacingStatusScalarWithSpace:
00096 StdErr << "WARNING: nrrdSpacingCalculate returned nrrdSpacingStatusScalarWithSpace\n";
00097 spacing[ax] = nrrdSpaceAxes[ax].spacing;
00098 break;
00099 case nrrdSpacingStatusNone:
00100 default:
00101 StdErr << "WARNING: no pixel spacings in Nrrd for axis " << ax << "; setting to 1.0\n";
00102 spacing[ax] = 1.0;
00103 break;
00104 }
00105 }
00106
00107 const Types::Coordinate size[3] = { (dims[0]-1) * spacing[0], (dims[1]-1) * spacing[1], (dims[2]-1) * spacing[2] };
00108 const Types::Coordinate origin[3] = { nrrd->spaceOrigin[0], nrrd->spaceOrigin[1], nrrd->spaceOrigin[2] };
00109 dfield = DeformationField::SmartPtr( new DeformationField( FixedVector<3,Types::Coordinate>( size ), DeformationField::IndexType( dims ), origin ) );
00110
00111 ScalarDataType type = TYPE_NONE;
00112 switch ( nrrd->type )
00113 {
00114 case nrrdTypeUChar: type = TYPE_BYTE; break;
00115 case nrrdTypeChar: type = TYPE_CHAR; break;
00116 case nrrdTypeUShort: type = TYPE_USHORT; break;
00117 case nrrdTypeShort: type = TYPE_SHORT; break;
00118 case nrrdTypeInt: type = TYPE_INT; break;
00119 case nrrdTypeFloat: type = TYPE_FLOAT; break;
00120 case nrrdTypeDouble: type = TYPE_DOUBLE; break;
00121 default: break;
00122 }
00123
00124 if ( type != TYPE_NONE )
00125 {
00126 TypedArray::SmartPtr data( TypedArray::Create( type, nrrd->data, 3 * dims[0] * dims[1] * dims[2] ) );
00127 data->ConvertSubArray( dfield->m_Parameters, TYPE_COORDINATE, 0, data->GetDataSize() );
00128 }
00129 else
00130 {
00131 StdErr << "ERROR: unsupported data type in nrrd file.\n";
00132 return dfield;
00133 }
00134
00135 const char* orientationSpace = NULL;
00136 switch ( nrrd->space )
00137 {
00138 case nrrdSpaceRightAnteriorSuperior:
00139 case nrrdSpaceRightAnteriorSuperiorTime:
00140 orientationSpace = "RAS";
00141 break;
00142 case nrrdSpaceLeftAnteriorSuperior:
00143 case nrrdSpaceLeftAnteriorSuperiorTime:
00144 orientationSpace = "LAS";
00145 break;
00146 case nrrdSpaceLeftPosteriorSuperior:
00147 case nrrdSpaceLeftPosteriorSuperiorTime:
00148 orientationSpace = "LPS";
00149 break;
00150 default:
00151 break;
00152 }
00153
00154 if ( orientationSpace )
00155 {
00156 dfield->m_MetaInformation[META_SPACE] = dfield->m_MetaInformation[META_SPACE_ORIGINAL] = orientationSpace;
00157
00158 const Types::Coordinate directions[3][3] =
00159 {
00160 { nrrdSpaceAxes[0].spaceDirection[0] * spacing[0],
00161 nrrdSpaceAxes[0].spaceDirection[1] * spacing[0],
00162 nrrdSpaceAxes[0].spaceDirection[2] * spacing[0] },
00163 { nrrdSpaceAxes[1].spaceDirection[0] * spacing[1],
00164 nrrdSpaceAxes[1].spaceDirection[1] * spacing[1],
00165 nrrdSpaceAxes[1].spaceDirection[2] * spacing[1] },
00166 { nrrdSpaceAxes[2].spaceDirection[0] * spacing[2],
00167 nrrdSpaceAxes[2].spaceDirection[1] * spacing[2],
00168 nrrdSpaceAxes[2].spaceDirection[2] * spacing[2] }
00169 };
00170
00171 const Matrix3x3<Types::Coordinate> m3( directions );
00172 Matrix4x4<Types::Coordinate> m4( m3 );
00173 for ( int i = 0; i < 3; ++i )
00174 m4[3][i] = nrrd->spaceOrigin[i];
00175
00176 AffineXform::SmartPtr xform( new AffineXform( m4 ) ) ;
00177 dfield->SetInitialAffineXform( xform );
00178
00179 char orientationImage[4];
00180 AnatomicalOrientation::GetOrientationFromDirections( orientationImage, m4, orientationSpace );
00181 dfield->m_MetaInformation[META_IMAGE_ORIENTATION] = orientationImage;
00182 dfield->m_MetaInformation[META_IMAGE_ORIENTATION_ORIGINAL] = orientationImage;
00183 }
00184
00185 nrrdNix( nrrd );
00186 }
00187 catch ( char* err )
00188 {
00189 StdErr << "ERROR: nrrd library returned error '" << err << "'\n";
00190 free( err );
00191 }
00192
00193 return dfield;
00194 }
00195
00196 void
00197 XformIO::WriteNrrd
00198 ( const Xform* xform, const char *path, const bool )
00199 {
00200 const DeformationField* dfield = dynamic_cast<const DeformationField*>( xform );
00201 if ( ! dfield )
00202 {
00203 StdErr << "ERROR: XformIO::WriteNrrd can only write DeformationField objects so far.\n"
00204 << " No data was written.\n";
00205 return;
00206 }
00207
00208 void* val = static_cast<void*>( dfield->m_Parameters );
00209 const int type = (sizeof(Types::Coordinate) == sizeof(float)) ? nrrdTypeFloat : nrrdTypeDouble;
00210
00211 Nrrd *nval = nrrdNew();
00212 NrrdIoState *nios = nrrdIoStateNew();
00213
00214 if ( nrrdEncodingGzip->available() )
00215 {
00216 nrrdIoStateEncodingSet( nios, nrrdEncodingGzip );
00217 nrrdIoStateSet( nios, nrrdIoStateZlibLevel, 9 );
00218 }
00219 else
00220 {
00221 StdErr << "WARNING: Nrrd library does not support Gzip compression encoding.\n"
00222 << " Please add -DTEEM_ZLIB to compiler options when building Nrrd library.\n";
00223 }
00224
00225 try
00226 {
00227 if ( nrrdWrap_va( nval, val, type, 4, (size_t)3, (size_t)dfield->m_Dims[0], (size_t)dfield->m_Dims[1], (size_t)dfield->m_Dims[2] ) )
00228 {
00229 throw( biffGetDone(NRRD) );
00230 }
00231
00232 nrrdSpaceDimensionSet( nval, 3 );
00233
00234 if ( dfield->MetaKeyExists(META_SPACE_UNITS_STRING) )
00235 {
00236 nval->spaceUnits[0] = strdup( dfield->m_MetaInformation[META_SPACE_UNITS_STRING].c_str() );
00237 }
00238
00239 int kind[NRRD_DIM_MAX] = { nrrdKindVector, nrrdKindDomain, nrrdKindDomain, nrrdKindDomain };
00240 nrrdAxisInfoSet_nva( nval, nrrdAxisInfoKind, kind );
00241 nrrdAxisInfoSet_va( nval, nrrdAxisInfoLabel, "Vx;Vy;Vz", "x", "y", "z" );
00242
00243 double origin[NRRD_DIM_MAX] = { dfield->m_Offset[0], dfield->m_Offset[1], dfield->m_Offset[2] };
00244 if ( nrrdSpaceOriginSet( nval, origin ) )
00245 {
00246 throw( biffGetDone(NRRD) );
00247 }
00248
00249 nval->space = nrrdSpaceRightAnteriorSuperior;
00250 double spaceDir[NRRD_DIM_MAX][NRRD_SPACE_DIM_MAX];
00251 for ( int i = 0; i < 4; ++i )
00252 {
00253 for ( int j = 0; j < 3; ++j )
00254 {
00255 if ( i )
00256 {
00257 if ( i-1 == j )
00258 spaceDir[i][j] = dfield->Spacing[ i-1 ];
00259 else
00260 spaceDir[i][j] = 0.0;
00261 }
00262 else
00263 {
00264 spaceDir[i][j] = AIR_NAN;
00265 }
00266 }
00267 }
00268 nrrdAxisInfoSet_nva( nval, nrrdAxisInfoSpaceDirection, spaceDir );
00269
00270 if ( nrrdSave( path, nval, nios ) )
00271 {
00272 throw( biffGetDone(NRRD) );
00273 }
00274 }
00275 catch ( char* err )
00276 {
00277 StdErr << "ERROR: NrrdIO library returned error '" << err << "'\n";
00278 free( err );
00279 }
00280
00281 nrrdIoStateNix( nios );
00282 nrrdNix(nval);
00283 }
00284
00286
00287 }
00288
00289 #endif // #ifdef CMTK_BUILD_NRRD
1.7.2
