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 <cmtkconfig.h>
00034
00035 #ifdef CMTK_BUILD_NRRD
00036
00037 #include "cmtkVolumeFromFile.h"
00038
00039 #include <IO/cmtkVolumeIO.h>
00040
00041 #include <Base/cmtkTypes.h>
00042 #include <Base/cmtkUniformVolume.h>
00043 #include <Base/cmtkAnatomicalOrientation.h>
00044
00045 #include <System/cmtkConsole.h>
00046
00047 #ifdef CMTK_BUILD_NRRD_TEEM
00048 # include <teem/nrrd.h>
00049 #else
00050 # include <NrrdIO.h>
00051 #endif
00052
00053 namespace
00054 cmtk
00055 {
00056
00059
00060 const UniformVolume::SmartPtr
00061 VolumeFromFile::ReadNRRD( const char* pathHdr )
00062 {
00063 UniformVolume::SmartPtr volume( NULL );
00064 try
00065 {
00066 Nrrd *nrrd = nrrdNew();
00067 if ( nrrdLoad( nrrd, pathHdr, NULL ) )
00068 throw biffGetDone(NRRD);
00069
00070 if ( nrrd->dim > 3 )
00071 {
00072 StdErr << "ERROR: for now, nrrd input can only handle data with dimension 3 or less.\n";
00073 return UniformVolume::SmartPtr( NULL );
00074 }
00075
00076 const int dims[3] =
00077 {
00078 (nrrd->dim > 0) ? nrrd->axis[0].size : 1,
00079 (nrrd->dim > 1) ? nrrd->axis[1].size : 1,
00080 (nrrd->dim > 2) ? nrrd->axis[2].size : 1
00081 };
00082
00083
00084 double spacing[3] = { 1, 1, 1 };
00085 for ( size_t ax = 0; ax < nrrd->dim; ++ax )
00086 {
00087 switch ( nrrdSpacingCalculate( nrrd, ax, spacing+ax, nrrd->axis[ax].spaceDirection ) )
00088 {
00089 case nrrdSpacingStatusScalarNoSpace:
00090 break;
00091 case nrrdSpacingStatusDirection:
00092 break;
00093 case nrrdSpacingStatusScalarWithSpace:
00094 StdErr << "WARNING: nrrdSpacingCalculate returned nrrdSpacingStatusScalarWithSpace\n";
00095 spacing[ax] = nrrd->axis[ax].spacing;
00096 break;
00097 case nrrdSpacingStatusNone:
00098 default:
00099 StdErr << "WARNING: no pixel spacings in Nrrd for axis " << ax << "; setting to 1.0\n";
00100 spacing[ax] = 1.0;
00101 break;
00102 }
00103 }
00104 const Types::Coordinate size[3] = { (dims[0]-1) * spacing[0], (dims[1]-1) * spacing[1], (dims[2]-1) * spacing[2] };
00105 volume = UniformVolume::SmartPtr( new UniformVolume( DataGrid::IndexType( dims ),UniformVolume::CoordinateVectorType( size ) ) );
00106
00107 ScalarDataType type = TYPE_NONE;
00108 switch ( nrrd->type )
00109 {
00110 case nrrdTypeUChar: type = TYPE_BYTE; break;
00111 case nrrdTypeChar: type = TYPE_CHAR; break;
00112 case nrrdTypeUShort: type = TYPE_USHORT; break;
00113 case nrrdTypeShort: type = TYPE_SHORT; break;
00114 case nrrdTypeInt: type = TYPE_INT; break;
00115 case nrrdTypeFloat: type = TYPE_FLOAT; break;
00116 case nrrdTypeDouble: type = TYPE_DOUBLE; break;
00117 default: break;
00118 }
00119
00120 if ( type != TYPE_NONE )
00121 {
00122 TypedArray::SmartPtr data( TypedArray::Create( type, nrrd->data, volume->GetNumberOfPixels() ) );
00123 volume->SetData( data );
00124 }
00125 else
00126 {
00127 StdErr << "ERROR: unsupported data type in nrrd file.\n";
00128 return volume;
00129 }
00130
00131 const char* orientationSpaceAnatomical = NULL;
00132 switch ( nrrd->space )
00133 {
00134 case nrrdSpaceRightAnteriorSuperior:
00135 case nrrdSpaceRightAnteriorSuperiorTime:
00136 orientationSpaceAnatomical = "RAS";
00137 volume->m_MetaInformation[META_SPACE] = orientationSpaceAnatomical;
00138 break;
00139 case nrrdSpaceLeftAnteriorSuperior:
00140 case nrrdSpaceLeftAnteriorSuperiorTime:
00141 orientationSpaceAnatomical = "LAS";
00142 volume->m_MetaInformation[META_SPACE] = orientationSpaceAnatomical;
00143 break;
00144 case nrrdSpaceLeftPosteriorSuperior:
00145 case nrrdSpaceLeftPosteriorSuperiorTime:
00146 orientationSpaceAnatomical = "LPS";
00147 volume->m_MetaInformation[META_SPACE] = orientationSpaceAnatomical;
00148 break;
00149 case nrrdSpace3DRightHanded:
00150 volume->m_MetaInformation[META_SPACE] = "3DRH";
00151 break;
00152 case nrrdSpace3DLeftHanded:
00153 volume->m_MetaInformation[META_SPACE] = "3DLH";
00154 break;
00155 case nrrdSpace3DRightHandedTime:
00156 volume->m_MetaInformation[META_SPACE] = "3DRHT";
00157 break;
00158 case nrrdSpace3DLeftHandedTime:
00159 volume->m_MetaInformation[META_SPACE] = "3DLHT";
00160 break;
00161 default:
00162 break;
00163 }
00164
00165 volume->m_MetaInformation[META_SPACE_ORIGINAL] = volume->m_MetaInformation[META_SPACE];
00166
00167 const Types::Coordinate directions[3][3] =
00168 {
00169 { nrrd->axis[0].spaceDirection[0] * spacing[0],
00170 nrrd->axis[0].spaceDirection[1] * spacing[0],
00171 nrrd->axis[0].spaceDirection[2] * spacing[0] },
00172 { nrrd->axis[1].spaceDirection[0] * spacing[1],
00173 nrrd->axis[1].spaceDirection[1] * spacing[1],
00174 nrrd->axis[1].spaceDirection[2] * spacing[1] },
00175 { nrrd->axis[2].spaceDirection[0] * spacing[2],
00176 nrrd->axis[2].spaceDirection[1] * spacing[2],
00177 nrrd->axis[2].spaceDirection[2] * spacing[2] }
00178 };
00179
00180 const Matrix3x3<Types::Coordinate> m3( directions );
00181 Matrix4x4<Types::Coordinate> m4( m3 );
00182 for ( int i = 0; i < 3; ++i )
00183 m4[3][i] = nrrd->spaceOrigin[i];
00184 volume->m_IndexToPhysicalMatrix = m4;
00185
00186 if ( orientationSpaceAnatomical )
00187 {
00188 char orientationImage[4];
00189 AnatomicalOrientation::GetOrientationFromDirections( orientationImage, m4, orientationSpaceAnatomical );
00190 volume->m_MetaInformation[META_IMAGE_ORIENTATION] = volume->m_MetaInformation[META_IMAGE_ORIENTATION_ORIGINAL] = orientationImage;
00191 volume->ChangeCoordinateSpace( AnatomicalOrientation::ORIENTATION_STANDARD );
00192 }
00193 else
00194 {
00195 }
00196
00197 if ( nrrd->spaceUnits[0] )
00198 volume->m_MetaInformation[META_SPACE_UNITS_STRING] = nrrd->spaceUnits[0];
00199
00200 nrrdNix( nrrd );
00201 }
00202 catch ( char* err )
00203 {
00204 StdErr << "ERROR: nrrd library returned error '" << err << "'\n";
00205 free( err );
00206 }
00207
00208 return volume;
00209 }
00210
00211
00212 void
00213 VolumeFromFile::WriteNRRD
00214 ( const char* pathHdr, const UniformVolume& volume, const bool )
00215 {
00216 UniformVolume::SmartPtr writeVolume( volume.Clone() );
00217 if ( writeVolume->MetaKeyExists( META_SPACE_ORIGINAL ) )
00218 writeVolume->ChangeCoordinateSpace( writeVolume->m_MetaInformation[META_SPACE_ORIGINAL] );
00219
00220 void* val = const_cast<void*>( writeVolume->GetData()->GetDataPtr() );
00221 int type = nrrdTypeUnknown;
00222 switch ( writeVolume->GetData()->GetType() )
00223 {
00224 case TYPE_BYTE: type = nrrdTypeUChar; break;
00225 case TYPE_CHAR: type = nrrdTypeChar; break;
00226 case TYPE_USHORT: type = nrrdTypeUShort; break;
00227 case TYPE_SHORT: type = nrrdTypeShort; break;
00228 case TYPE_INT: type = nrrdTypeInt; break;
00229 case TYPE_FLOAT: type = nrrdTypeFloat; break;
00230 case TYPE_DOUBLE: type = nrrdTypeDouble; break;
00231 default: break;
00232 }
00233
00234 Nrrd *nval = nrrdNew();
00235 NrrdIoState *nios = nrrdIoStateNew();
00236
00237 if ( VolumeIO::GetWriteCompressed() && nrrdEncodingGzip->available() )
00238 {
00239 nrrdIoStateEncodingSet( nios, nrrdEncodingGzip );
00240 nrrdIoStateSet( nios, nrrdIoStateZlibLevel, 9 );
00241 }
00242 else
00243 {
00244 StdErr << "WARNING: Nrrd library does not support Gzip compression encoding.\n"
00245 << " Please add -DTEEM_ZLIB to compiler options when building Nrrd library.\n";
00246 }
00247
00248 try
00249 {
00250 if ( nrrdWrap_va( nval, val, type, (size_t)3, (size_t)writeVolume->m_Dims[0], (size_t)writeVolume->m_Dims[1], (size_t)writeVolume->m_Dims[2] ) )
00251 {
00252 throw( biffGetDone(NRRD) );
00253 }
00254
00255 nrrdSpaceDimensionSet( nval, 3 );
00256
00257 if ( writeVolume->MetaKeyExists(META_SPACE_UNITS_STRING) )
00258 {
00259 nval->spaceUnits[0] = strdup( writeVolume->m_MetaInformation[META_SPACE_UNITS_STRING].c_str() );
00260 nval->spaceUnits[1] = strdup( writeVolume->m_MetaInformation[META_SPACE_UNITS_STRING].c_str() );
00261 nval->spaceUnits[2] = strdup( writeVolume->m_MetaInformation[META_SPACE_UNITS_STRING].c_str() );
00262 }
00263
00264 int kind[NRRD_DIM_MAX] = { nrrdKindDomain, nrrdKindDomain, nrrdKindDomain };
00265 nrrdAxisInfoSet_nva( nval, nrrdAxisInfoKind, kind );
00266
00267 const std::string space = writeVolume->m_MetaInformation[META_SPACE];
00268
00269
00270 if ( space == "RAS" )
00271 {
00272 nval->space = nrrdSpaceRightAnteriorSuperior;
00273 }
00274 else if ( space == "LAS" )
00275 {
00276 nval->space = nrrdSpaceLeftAnteriorSuperior;
00277 }
00278 else if ( space == "LPS" )
00279 {
00280 nval->space = nrrdSpaceLeftPosteriorSuperior;
00281 }
00282 else if ( space == "3DRH" )
00283 {
00284 nval->space = nrrdSpace3DRightHanded;
00285 }
00286 else if ( space == "3DLH" )
00287 {
00288 nval->space = nrrdSpace3DLeftHanded;
00289 }
00290 else if ( space == "3DRHT" )
00291 {
00292 nval->space = nrrdSpace3DRightHandedTime;
00293 }
00294 else if ( space == "3DLHT" )
00295 {
00296 nval->space = nrrdSpace3DLeftHandedTime;
00297 }
00298 else
00299 {
00300 if ( space.length() == 3 )
00301 {
00302 writeVolume->ChangeCoordinateSpace( "RAS" );
00303 nval->space = nrrdSpaceRightAnteriorSuperior;
00304 }
00305 }
00306
00307 const AffineXform::MatrixType& matrix = writeVolume->m_IndexToPhysicalMatrix;
00308 double spaceDir[NRRD_DIM_MAX][NRRD_SPACE_DIM_MAX];
00309 for ( int i = 0; i < 3; ++i )
00310 {
00311 for ( int j = 0; j < 3; ++j )
00312 spaceDir[i][j] = matrix[i][j];
00313 }
00314 nrrdAxisInfoSet_nva( nval, nrrdAxisInfoSpaceDirection, spaceDir );
00315
00316 double origin[NRRD_DIM_MAX] = { matrix[3][0], matrix[3][1], matrix[3][2] };
00317 if ( nrrdSpaceOriginSet( nval, origin ) )
00318 {
00319 throw( biffGetDone(NRRD) );
00320 }
00321
00322 nrrdAxisInfoSet_va( nval, nrrdAxisInfoLabel, "x", "y", "z" );
00323
00324 if ( nrrdSave( pathHdr, nval, nios ) )
00325 {
00326 throw( biffGetDone(NRRD) );
00327 }
00328 }
00329 catch ( char* err )
00330 {
00331 StdErr << "ERROR: NrrdIO library returned error '" << err << "'\n";
00332 free( err );
00333 }
00334
00335 nrrdIoStateNix( nios );
00336 nrrdNix(nval);
00337 }
00338
00340
00341 }
00342
00343 #endif // #ifdef CMTK_BUILD_NRRD
1.7.2
