From b82c3a8dfeb283260ba6b1be6cb088194833116b Mon Sep 17 00:00:00 2001 From: yhuang43 <92936691+yhuang43@users.noreply.github.com> Date: Mon, 6 Dec 2021 15:06:41 -0500 Subject: [PATCH] new type MRI_USHRT, env variables FS_MGZIO_USEVOXELBUF & FS_MGZIO_TIMING (#899) * 1. check in package dcm2niix 2. add DICOMRead3/niiRead3 to convert dicom using dcm2niix 3. fix dcmGetDWIParams() to use stricmp() to compare the Manufacturer tag * set TI to -1 when there is no TI in dicom file for MGH_FREESURFER * remove unused codes * This is the correct version. * redo the check in * For some unknown reasons, dcm2niix and dcm2niixfsexe executables failed linking in nightly builds. They are excluded from the build until we figure out why. libdcm2niixfs.a is renamed to libdcm2niix.a * 1. rename dcm2fsWrapper.h => dcm2niix_fswrapper.h, dcm2fsWrapper.cpp => dcm2niix_fswrapper.cpp 2. add comments to dcm2niix changes * fix bvecs * 1. introduced new type MRI_USHRT (10) - unsigned short, convert DT_UINT16 to MRI_USHRT, MRI_USHRT to DT_UINT16 2. added macro MRIUSseq_vox & MRIUSvox to access voxel in pre-allocated memory 3. added handling of type MRI_USHRT in various functions 4. introduced environment variable FS_MGZIO_USEVOXELBUF to read/write (mgzRead/mgzWrite) voxels in one buffer access. FS_MGZIO_USEVOXELBUF is ignored if mri voxel buffer is not allocated as a big chunk. 5. introduced environment variable FS_MGZIO_TIMING to time voxel read/write in mgzRead/mgzWrite * ignore FS_MGZIO_USEVOXELBUF if mri voxel buffer is not allocated as a big chunk Co-authored-by: Yujing Huang --- freeview/FSVolume.cpp | 42 +++ freeview/LayerODF.cpp | 3 + freeview/LayerSurface.cpp | 7 + freeview/SurfaceAnnotation.cpp | 3 + freeview/VolumeFilter.cpp | 12 +- include/fio.h | 1 + include/machine.h | 2 + include/mri.h | 4 + matlab/HippoSF/myMRIread.m | 7 + matlab/load_mgh.m | 5 + matlab/load_mgh2.m | 5 +- mri_ca_normalize/mri_ca_normalize.cpp | 8 + mri_info/mri_info.cpp | 7 +- mri_ms_fitparms/mri_ms_fitparms.cpp | 19 + mri_rigid_register/mri_rigid_register.cpp | 17 + mri_robust_register/CostFunctions.h | 10 + mri_robust_register/MyMRI.cpp | 20 ++ mri_volcluster/mri_volcluster.cpp | 3 + pointset2label/pointset2label.cpp | 1 + utils/fio.cpp | 8 + utils/fsgdf.cpp | 4 + utils/gca.cpp | 6 + utils/gcamorph.cpp | 14 + utils/machine.cpp | 22 ++ utils/mri.cpp | 369 +++++++++++++++++++- utils/mri2.cpp | 6 + utils/mri_topology.cpp | 20 ++ utils/mrihisto.cpp | 40 +++ utils/mriio.cpp | 403 +++++++++++++++------- 29 files changed, 940 insertions(+), 128 deletions(-) diff --git a/freeview/FSVolume.cpp b/freeview/FSVolume.cpp index 6361e751cbc..2b7337e5c7b 100644 --- a/freeview/FSVolume.cpp +++ b/freeview/FSVolume.cpp @@ -583,6 +583,9 @@ bool FSVolume::Create( FSVolume* src_vol, bool bCopyVoxelData, int data_type ) case MRI_SHORT: m_imageData->AllocateScalars(VTK_SHORT, 1); break; + case MRI_USHRT: + m_imageData->AllocateScalars(VTK_UNSIGNED_SHORT, 1); + break; default: break; } @@ -605,6 +608,9 @@ bool FSVolume::Create( FSVolume* src_vol, bool bCopyVoxelData, int data_type ) case MRI_SHORT: m_imageData->SetScalarTypeToShort(); break; + case MRI_USHRT: + m_imageData->SetScalarTypeToUnsignedShort(); + break; default: break; } @@ -1195,6 +1201,9 @@ bool FSVolume::UpdateMRIFromImage( vtkImageData* rasImage, bool resampleToOrigin case MRI_SHORT: MRISseq_vox( mri, i, j, k, nFrame ) = (short)val; break; + case MRI_USHRT: + MRIUSseq_vox( mri, i, j, k, nFrame ) = (unsigned short)val; + break; default: break; } @@ -2042,6 +2051,9 @@ bool FSVolume::CreateImage( MRI* rasMRI ) case MRI_SHORT: imageData->AllocateScalars(VTK_SHORT, zFrames); break; + case MRI_USHRT: + imageData->AllocateScalars(VTK_UNSIGNED_SHORT, zFrames); + break; default: return false; } @@ -2065,6 +2077,9 @@ bool FSVolume::CreateImage( MRI* rasMRI ) case MRI_SHORT: imageData->SetScalarTypeToShort(); break; + case MRI_USHRT: + imageData->SetScalarTypeToUnsignedShort(); + break; default: return false; } @@ -2150,6 +2165,9 @@ bool FSVolume::ResizeRotatedImage( MRI* rasMRI, MRI* refTarget, vtkImageData* re case MRI_SHORT: imageData->AllocateScalars(VTK_SHORT, zFrames); break; + case MRI_USHRT: + imageData->AllocateScalars(VTK_UNSIGNED_SHORT, zFrames); + break; default: return false; } @@ -2173,6 +2191,9 @@ bool FSVolume::ResizeRotatedImage( MRI* rasMRI, MRI* refTarget, vtkImageData* re case MRI_SHORT: imageData->SetScalarTypeToShort(); break; + case MRI_USHRT: + imageData->SetScalarTypeToUnsignedShort(); + break; default: break ; } @@ -2406,6 +2427,9 @@ void FSVolume::CopyMRIDataToImage( MRI* mri, case MRI_SHORT: ((short*)ptr)[nTuple*zFrames+nFrame] = MRISseq_vox( mri, nX, nY, nZ, nFrame ); break; + case MRI_USHRT: + ((unsigned short*)ptr)[nTuple*zFrames+nFrame] = MRIUSseq_vox( mri, nX, nY, nZ, nFrame ); + break; default: break; } @@ -3163,6 +3187,24 @@ void FSVolume::GetFrameValueRange(int frame, double *range) } } break ; + case MRI_USHRT: + { + for (z = 0 ; z < depth ; z++) + { + for (y = 0 ; y < height ; y++) + { + for (x = 0 ; x < width ; x++) + { + val = (float)MRIUSseq_vox(mri, x, y, z, frame) ; + if (val < fmin) + fmin = val ; + if (val > fmax) + fmax = val ; + } + } + } + } + break ; case MRI_UCHAR: { for (z = 0 ; z < depth ; z++) diff --git a/freeview/LayerODF.cpp b/freeview/LayerODF.cpp index d39d290cd0b..e78d71e5683 100644 --- a/freeview/LayerODF.cpp +++ b/freeview/LayerODF.cpp @@ -197,6 +197,9 @@ bool LayerODF::Load(const QString &fn, const QString& vertex_fn, const QString& case MRI_SHORT: MRISseq_vox( mri2, j, k, n, i ) = MRISseq_vox( mri, i, j, k, n ); break; + case MRI_USHRT: + MRIUSseq_vox( mri2, j, k, n, i ) = MRIUSseq_vox( mri, i, j, k, n ); + break; default: break; } diff --git a/freeview/LayerSurface.cpp b/freeview/LayerSurface.cpp index b2401c15d4e..11c4524cae1 100644 --- a/freeview/LayerSurface.cpp +++ b/freeview/LayerSurface.cpp @@ -2402,6 +2402,13 @@ bool LayerSurface::LoadRGBFromFile(const QString &filename) for (int j = 0; j < 3; j++) map.data << MRISseq_vox( mri, i, j, 0, 0 ); break; + + case MRI_USHRT: + for (int i = 0; i < GetNumberOfVertices(); i++) + for (int j = 0; j < 3; j++) + map.data << MRIUSseq_vox( mri, i, j, 0, 0 ); + break; + default: MRIfree(&mri); return false; diff --git a/freeview/SurfaceAnnotation.cpp b/freeview/SurfaceAnnotation.cpp index ca48ec4a797..664ee035e0b 100644 --- a/freeview/SurfaceAnnotation.cpp +++ b/freeview/SurfaceAnnotation.cpp @@ -163,6 +163,9 @@ bool SurfaceAnnotation::LoadFromSegmentation(const QString &fn) case MRI_SHORT: n = MRIseq_vox( mri, i, 0, 0, 0); break; + case MRI_USHRT: + n = (int)MRIUSseq_vox( mri, i, 0, 0, 0); + break; default: break; } diff --git a/freeview/VolumeFilter.cpp b/freeview/VolumeFilter.cpp index 480be8709f9..496c5bdf20d 100644 --- a/freeview/VolumeFilter.cpp +++ b/freeview/VolumeFilter.cpp @@ -125,9 +125,11 @@ MRI* VolumeFilter::CreateMRIFromVolume( LayerMRI* layer ) mri_type = MRI_LONG; break; case VTK_SHORT: - case VTK_UNSIGNED_SHORT: mri_type = MRI_SHORT; break; + case VTK_UNSIGNED_SHORT: + mri_type = MRI_USHRT; + break; default: break; } @@ -180,6 +182,10 @@ MRI* VolumeFilter::CreateMRIFromVolume( LayerMRI* layer ) MRISseq_vox( mri, i, j, k, nFrame ) = (int)MyVTKUtils::GetImageDataComponent(ptr, dim, n_frames, i, j, k, nFrame, scalar_type); break; + case MRI_USHRT: + MRIUSseq_vox( mri, i, j, k, nFrame ) = + (unsigned short)MyVTKUtils::GetImageDataComponent(ptr, dim, n_frames, i, j, k, nFrame, scalar_type); + break; default: break; } @@ -234,6 +240,10 @@ void VolumeFilter::MapMRIToVolume( MRI* mri, LayerMRI* layer ) MyVTKUtils::SetImageDataComponent(ptr, dim, n_frames, nX, nY, nZ, nFrame, scalar_type, MRISseq_vox( mri, nX, nY, nZ, nFrame ) ); break; + case MRI_USHRT: + MyVTKUtils::SetImageDataComponent(ptr, dim, n_frames, nX, nY, nZ, nFrame, scalar_type, + MRIUSseq_vox( mri, nX, nY, nZ, nFrame ) ); + break; default: break; } diff --git a/include/fio.h b/include/fio.h index f8225811531..931f7d24720 100644 --- a/include/fio.h +++ b/include/fio.h @@ -69,6 +69,7 @@ int znzreadShortEx (short *ps, znzFile fp) ; int znzwriteDouble (double d, znzFile fp) ; int znzwriteFloat (float f, znzFile fp) ; int znzwriteShort (short s, znzFile fp) ; +int znzwriteUShort (unsigned short s, znzFile fp) ; int znzwriteInt (int v, znzFile fp) ; int znzwriteLong (long long v, znzFile fp) ; int znzwrite1 (int v, znzFile fp) ; diff --git a/include/machine.h b/include/machine.h index ebc550b7905..b6f192fdfd3 100644 --- a/include/machine.h +++ b/include/machine.h @@ -30,6 +30,7 @@ typedef int long32; typedef long long long64; +unsigned short swapUShort(unsigned short us); short swapShort(short s) ; long32 swapLong32(long32 l); long64 swapLong64(long64 l); @@ -47,6 +48,7 @@ int ByteSwap8(void *buf8, long int nitems); #define orderIntBytes(i) swapInt(i) #define orderShortBytes(i) swapShort(i) +#define orderUShortBytes(i) swapUShort(i) #define orderFloatBytes(i) swapFloat(i) #define orderDoubleBytes(i) swapDouble(i) #define orderLong32Bytes(i) swapLong32(i) diff --git a/include/mri.h b/include/mri.h index 661ae83c4a7..1e7e6cd7024 100644 --- a/include/mri.h +++ b/include/mri.h @@ -58,6 +58,7 @@ #define MRI_FLOAT_COMPLEX 7 #define MRI_DOUBLE_COMPLEX 8 #define MRI_RGB 9 +#define MRI_USHRT 10 #define NEAREST_NEIGHBOR_FACE 1 #define NEAREST_NEIGHBOR_EDGE 2 @@ -927,6 +928,7 @@ int MRIsurfaceRASToRAS(MRI *mri, double xsr, double ysr, double zsr, #define MRIclear_bit(mri,x,y,z) MRIvox(mri,(x)/8,y,z) &= ~(0x001 << ((x)%8)) #define MRISvox(mri,x,y,z) (((short *)mri->slices[z][y])[x]) +#define MRIUSvox(mri,x,y,z) (((unsigned short *)mri->slices[z][y])[x]) #define MRIFvox(mri,x,y,z) (((float *)(mri->slices[z][y]))[x]) #define MRIvox(mri,x,y,z) (((BUFTYPE *)mri->slices[z][y])[x]) #define MRISCvox(mri,x,y,z) (((signed char *)mri->slices[z][y])[x]) @@ -935,6 +937,8 @@ int MRIsurfaceRASToRAS(MRI *mri, double xsr, double ysr, double zsr, #define MRISseq_vox(mri,x,y,z,n) (((short*)\ mri->slices[z+(n)*mri->depth][y])[x]) +#define MRIUSseq_vox(mri,x,y,z,n) (((unsigned short*)\ +mri->slices[z+(n)*mri->depth][y])[x]) #define MRISCseq_vox(mri,x,y,z,n) (((signed char*)\ mri->slices[z+(n)*mri->depth][y])[x]) #define MRIFseq_vox(mri,x,y,z,n) (((float*)\ diff --git a/matlab/HippoSF/myMRIread.m b/matlab/HippoSF/myMRIread.m index 3410f763030..0c600ef9c30 100755 --- a/matlab/HippoSF/myMRIread.m +++ b/matlab/HippoSF/myMRIread.m @@ -537,6 +537,7 @@ MRI_FLOAT = 3 ; MRI_SHORT = 4 ; MRI_BITMAP = 5 ; +MRI_USHRT = 10 ; % Determine number of bytes per voxel switch type @@ -546,6 +547,8 @@ nbytespervox = 1; case MRI_SHORT, nbytespervox = 2; + case MRI_USHRT, + nbytespervox = 2; case MRI_INT, nbytespervox = 4; end @@ -572,6 +575,8 @@ vol = fread(fid, nv, 'uchar') ; case MRI_SHORT, vol = fread(fid, nv, 'short') ; + case MRI_USHRT, + vol = fread(fid, nv, 'uint16') ; case MRI_INT, vol = fread(fid, nv, 'int') ; end @@ -620,6 +625,8 @@ [tmpslice nread] = fread(fid, nvslice, 'uchar') ; case MRI_SHORT, [tmpslice nread] = fread(fid, nvslice, 'short') ; + case MRI_USHRT, + [tmpslice nread] = fread(fid, nvslice, 'uint16') ; case MRI_INT, [tmpslice nread] = fread(fid, nvslice, 'int') ; end diff --git a/matlab/load_mgh.m b/matlab/load_mgh.m index ab9609cb513..3299ef3179c 100644 --- a/matlab/load_mgh.m +++ b/matlab/load_mgh.m @@ -157,6 +157,7 @@ MRI_FLOAT = 3 ; MRI_SHORT = 4 ; MRI_BITMAP = 5 ; +MRI_USHRT = 10 ; % Determine number of bytes per voxel switch type @@ -166,6 +167,8 @@ nbytespervox = 1; case MRI_SHORT, nbytespervox = 2; + case MRI_USHRT, + nbytespervox = 2; case MRI_INT, nbytespervox = 4; end @@ -196,6 +199,8 @@ dtype = 'short' ; case MRI_INT, dtype = 'int' ; + case MRI_USHRT, + dtype = 'uint16' ; end % preserve volume datatype if env var is set to 1 diff --git a/matlab/load_mgh2.m b/matlab/load_mgh2.m index 64ea5bc6703..d9b628ccb1b 100644 --- a/matlab/load_mgh2.m +++ b/matlab/load_mgh2.m @@ -82,6 +82,7 @@ MRI_FLOAT = 3 ; MRI_SHORT = 4 ; MRI_BITMAP = 5 ; +MRI_USHRT = 10 ; fseek(fid, unused_space_size, 'cof') ; @@ -95,7 +96,9 @@ case MRI_SHORT, vol = fread(fid, nv, 'short') ; case MRI_INT, - vol = fread(fid, nv, 'int') ; + vol = fread(fid, nv, 'int') ; + case MRI_USHRT, + vol = fread(fid, nv, 'uint16') ; end if(~feof(fid)) [mr_parms count] = fread(fid,4,'float32'); diff --git a/mri_ca_normalize/mri_ca_normalize.cpp b/mri_ca_normalize/mri_ca_normalize.cpp index 9aa7b286695..c6d19729b13 100644 --- a/mri_ca_normalize/mri_ca_normalize.cpp +++ b/mri_ca_normalize/mri_ca_normalize.cpp @@ -1946,6 +1946,10 @@ normalizeFromLabel(MRI *mri_in, MRI *mri_dst, MRI *mri_seg, double *fas) MRISseq_vox(mri_dst, x, y, z, input) = (short)nint(val) ; break ; + case MRI_USHRT: + MRIUSseq_vox(mri_dst, x, y, z, input) = + (unsigned short)nint(val) ; + break ; case MRI_FLOAT: MRIFseq_vox(mri_dst, x, y, z, input) = val ; @@ -2146,6 +2150,10 @@ normalizeChannelFromLabel(MRI *mri_in, MRI *mri_dst, MRI *mri_seg, MRISseq_vox(mri_dst, x, y, z, input_index) = (short)nint(val) ; break ; + case MRI_USHRT: + MRIUSseq_vox(mri_dst, x, y, z, input_index) = + (unsigned short)nint(val) ; + break ; case MRI_FLOAT: MRIFseq_vox(mri_dst, x, y, z, input_index) = val ; diff --git a/mri_info/mri_info.cpp b/mri_info/mri_info.cpp index 3cf201697a3..205bf648c59 100644 --- a/mri_info/mri_info.cpp +++ b/mri_info/mri_info.cpp @@ -696,6 +696,7 @@ static void do_file(char *fname) case MRI_FLOAT: fprintf(fpout,"float\n") ; break ; case MRI_LONG: fprintf(fpout,"long\n") ; break ; case MRI_SHORT: fprintf(fpout,"short\n") ; break ; + case MRI_USHRT: fprintf(fpout,"ushrt\n") ; break ; case MRI_INT: fprintf(fpout,"int\n") ; break ; case MRI_TENSOR: fprintf(fpout,"tensor\n") ; break ; } @@ -761,7 +762,10 @@ static void do_file(char *fname) case MRI_SHORT: fprintf(fpout,"short\n") ; break ; - case MRI_INT: + case MRI_USHRT: + fprintf(fpout,"ushrt\n") ; + break ; + case MRI_INT: fprintf(fpout,"int\n") ; break ; case MRI_TENSOR: @@ -1153,6 +1157,7 @@ static void do_file(char *fname) printf(" type: %s (%d)\n", mri->type == MRI_UCHAR ? "UCHAR" : mri->type == MRI_SHORT ? "SHORT" : + mri->type == MRI_USHRT ? "USHRT" : mri->type == MRI_INT ? "INT" : mri->type == MRI_LONG ? "LONG" : mri->type == MRI_BITMAP ? "BITMAP" : diff --git a/mri_ms_fitparms/mri_ms_fitparms.cpp b/mri_ms_fitparms/mri_ms_fitparms.cpp index bccf0996335..f3aac38f3fb 100644 --- a/mri_ms_fitparms/mri_ms_fitparms.cpp +++ b/mri_ms_fitparms/mri_ms_fitparms.cpp @@ -2636,6 +2636,7 @@ MRIssqrt(MRI *mri_src, MRI *mri_dst) { int width, height, depth, x, y, z, frame ; short *psrc, *pdst ; + unsigned short *psrc2, *pdst2 ; width = mri_src->width ; height = mri_src->height ; @@ -2662,6 +2663,15 @@ MRIssqrt(MRI *mri_src, MRI *mri_dst) *pdst++ = sqrt(*psrc++) ; } break ; + case MRI_USHRT: + psrc2 = &MRIUSseq_vox(mri_src, 0, y, z, frame) ; + pdst2 = &MRIUSseq_vox(mri_dst, 0, y, z, frame) ; + for (x = 0 ; x < width ; x++) + { + check_finite(sqrt(*psrc2)) ; + *pdst2++ = sqrt(*psrc2++) ; + } + break ; default: ErrorReturn(NULL, (ERROR_UNSUPPORTED, @@ -2685,6 +2695,7 @@ MRIsscalarMul(MRI *mri_src, MRI *mri_dst, float scalar) { int width, height, depth, x, y, z, frame ; short *psrc, *pdst ; + unsigned short *psrc2, *pdst2 ; width = mri_src->width ; height = mri_src->height ; @@ -2710,6 +2721,14 @@ MRIsscalarMul(MRI *mri_src, MRI *mri_dst, float scalar) *pdst++ = *psrc++ * scalar ; } break ; + case MRI_USHRT: + psrc2 = &MRIUSseq_vox(mri_src, 0, y, z, frame) ; + pdst2 = &MRIUSseq_vox(mri_dst, 0, y, z, frame) ; + for (x = 0 ; x < width ; x++) + { + *pdst2++ = *psrc2++ * scalar ; + } + break ; default: ErrorReturn(NULL, (ERROR_UNSUPPORTED, diff --git a/mri_rigid_register/mri_rigid_register.cpp b/mri_rigid_register/mri_rigid_register.cpp index 4f0e2e6f359..0d9bc3fc333 100644 --- a/mri_rigid_register/mri_rigid_register.cpp +++ b/mri_rigid_register/mri_rigid_register.cpp @@ -1678,6 +1678,9 @@ apply_transform(MRI *mri_in, MRI *mri_target, MATRIX *M_reg, MRI *mri_xformed) { case MRI_SHORT: MRISvox(mri_xformed, x, y, z) = (short)nint(val) ; break ; + case MRI_USHRT: + MRIUSvox(mri_xformed, x, y, z) = (unsigned short)nint(val) ; + break ; case MRI_FLOAT: MRIFvox(mri_xformed, x, y, z) = (float)(val) ; break ; @@ -1797,6 +1800,7 @@ MRI * MRIssqrt(MRI *mri_src, MRI *mri_dst) { int width, height, depth, x, y, z, frame ; short *psrc, *pdst ; + unsigned short *psrc2, *pdst2 ; width = mri_src->width ; height = mri_src->height ; @@ -1814,6 +1818,12 @@ MRIssqrt(MRI *mri_src, MRI *mri_dst) { for (x = 0 ; x < width ; x++) *pdst++ = sqrt(*psrc++) ; break ; + case MRI_USHRT: + psrc2 = &MRIUSseq_vox(mri_src, 0, y, z, frame) ; + pdst2 = &MRIUSseq_vox(mri_dst, 0, y, z, frame) ; + for (x = 0 ; x < width ; x++) + *pdst2++ = sqrt(*psrc2++) ; + break ; default: ErrorReturn(NULL, (ERROR_UNSUPPORTED, @@ -1835,6 +1845,7 @@ MRI * MRIsscalarMul(MRI *mri_src, MRI *mri_dst, float scalar) { int width, height, depth, x, y, z, frame ; short *psrc, *pdst ; + unsigned short *psrc2, *pdst2 ; width = mri_src->width ; height = mri_src->height ; @@ -1852,6 +1863,12 @@ MRIsscalarMul(MRI *mri_src, MRI *mri_dst, float scalar) { for (x = 0 ; x < width ; x++) *pdst++ = *psrc++ * scalar ; break ; + case MRI_USHRT: + psrc2 = &MRIUSseq_vox(mri_src, 0, y, z, frame) ; + pdst2 = &MRIUSseq_vox(mri_dst, 0, y, z, frame) ; + for (x = 0 ; x < width ; x++) + *pdst2++ = *psrc2++ * scalar ; + break ; default: ErrorReturn(NULL, (ERROR_UNSUPPORTED, diff --git a/mri_robust_register/CostFunctions.h b/mri_robust_register/CostFunctions.h index 42a194a1f60..309532f4f0c 100644 --- a/mri_robust_register/CostFunctions.h +++ b/mri_robust_register/CostFunctions.h @@ -203,6 +203,7 @@ class MRIiterator float fromUCHAR(void); float fromSHORT(void); + float fromUSHRT(void); float fromINT(void); float fromLONG(void); float fromFLOAT(void); @@ -236,6 +237,10 @@ inline MRIiterator::MRIiterator(MRI * i) : getVal = &MRIiterator::fromSHORT; bytes_per_voxel = sizeof(short); break; + case MRI_USHRT: + getVal = &MRIiterator::fromUSHRT; + bytes_per_voxel = sizeof(unsigned short); + break; case MRI_INT: getVal = &MRIiterator::fromINT; bytes_per_voxel = sizeof(int); @@ -332,6 +337,11 @@ inline float MRIiterator::fromSHORT() return ((float) *(short *) pos); } +inline float MRIiterator::fromUSHRT() +{ + return ((float) *(unsigned short *) pos); +} + inline float MRIiterator::fromINT() { return ((float) *(int *) pos); diff --git a/mri_robust_register/MyMRI.cpp b/mri_robust_register/MyMRI.cpp index abe557dd904..e9f85a38080 100644 --- a/mri_robust_register/MyMRI.cpp +++ b/mri_robust_register/MyMRI.cpp @@ -102,6 +102,26 @@ MRI * MyMRI::MRIvalscale(MRI *mri_src, MRI *mri_dst, double s, double b) } } break; + case MRI_USHRT: + for (f = 0; f < nf ; f++ ) + for (z = 0; z < depth; z++) + { + for (y = 0; y < height; y++) + { + unsigned short *ps_src2 = &MRIUSseq_vox(mri_src, 0, y, z, f); + unsigned short *ps_dst2 = &MRIUSseq_vox(mri_dst, 0, y, z, f); + for (x = 0; x < width; x++) + { + val = (float)(*ps_src2++); + val -= b; + val *= s; + if (val < 0) val = 0; + if (val > USHRT_MAX) val = USHRT_MAX; + *ps_dst2++ = (unsigned short)nint(val); + } + } + } + break; case MRI_UCHAR: assert(s > 0); for (f = 0; f < nf ; f++ ) diff --git a/mri_volcluster/mri_volcluster.cpp b/mri_volcluster/mri_volcluster.cpp index e69f49531e3..6d98e3a2785 100644 --- a/mri_volcluster/mri_volcluster.cpp +++ b/mri_volcluster/mri_volcluster.cpp @@ -1860,6 +1860,9 @@ static MRI *MRIbinarize01(MRI *vol, float thmin, float thmax, case MRI_SHORT: val = (float)(MRISseq_vox(vol,col,row,slice,frame)); break; + case MRI_USHRT: + val = (float)(MRIUSseq_vox(vol,col,row,slice,frame)); + break; } if (ithsign == 0) val = fabs(val); diff --git a/pointset2label/pointset2label.cpp b/pointset2label/pointset2label.cpp index 83484a50d93..18f9e9992fd 100644 --- a/pointset2label/pointset2label.cpp +++ b/pointset2label/pointset2label.cpp @@ -21,6 +21,7 @@ void set_voxel_value(MRI* mri, double dx, double dy, double dz, double val) MRIIseq_vox( mri, x, y, z, 0 ) = (int)val; break; case MRI_SHORT: + case MRI_USHRT: MRIIseq_vox( mri, x, y, z, 0 ) = (int)val; break; default: diff --git a/utils/fio.cpp b/utils/fio.cpp index 1d0b8c5874d..59eaa873408 100644 --- a/utils/fio.cpp +++ b/utils/fio.cpp @@ -428,6 +428,14 @@ int znzwriteShort(short s, znzFile fp) return (znzwrite(&s, sizeof(short), 1, fp)); } +int znzwriteUShort(unsigned short s, znzFile fp) +{ +#if (BYTE_ORDER == LITTLE_ENDIAN) + s = swapUShort(s); +#endif + return (znzwrite(&s, sizeof(unsigned short), 1, fp)); +} + double znzreadDouble(znzFile fp) { double d; diff --git a/utils/fsgdf.cpp b/utils/fsgdf.cpp index 3cb855c72ef..c7b7648132b 100644 --- a/utils/fsgdf.cpp +++ b/utils/fsgdf.cpp @@ -1744,6 +1744,10 @@ int gdfGetNthSubjectMeasurement(FSGD *gd, int nsubject, v = MRISseq_vox(gd->data,x,y,z,nsubject); *value = v; break; + case MRI_USHRT: + v = MRIUSseq_vox(gd->data,x,y,z,nsubject); + *value = v; + break; default: break ; } diff --git a/utils/gca.cpp b/utils/gca.cpp index 2230ae8a7f7..49fa6ca4e1a 100644 --- a/utils/gca.cpp +++ b/utils/gca.cpp @@ -5416,6 +5416,9 @@ MRI *GCAlabelMri(GCA *gca, MRI *mri, int label, TRANSFORM *transform) case MRI_SHORT: MRISseq_vox(mri, x, y, z, frame) = (short)val; break; + case MRI_USHRT: + MRIUSseq_vox(mri, x, y, z, frame) = (unsigned short)val; + break; case MRI_FLOAT: MRIFseq_vox(mri, x, y, z, frame) = (float)val; break; @@ -8205,6 +8208,9 @@ MRI *GCAnormalizeSamples( case MRI_SHORT: MRISseq_vox(mri_dst, x, y, z, input) = (short)nint(val); break; + case MRI_USHRT: + MRIUSseq_vox(mri_dst, x, y, z, input) = (unsigned short)nint(val); + break; case MRI_FLOAT: MRIFseq_vox(mri_dst, x, y, z, input) = val; break; diff --git a/utils/gcamorph.cpp b/utils/gcamorph.cpp index 0a0c8b25af9..4a8cc60875a 100644 --- a/utils/gcamorph.cpp +++ b/utils/gcamorph.cpp @@ -4129,6 +4129,11 @@ MRI *GCAMmorphFromAtlas(MRI *mri_in, GCA_MORPH *gcam, MRI *mri_morphed, int samp MRISseq_vox(mri_morphed, x, y, z, f) = (short)MRIFseq_vox(mri_s_morphed, x, y, z, f); } break; + case MRI_USHRT: + for (f = 0; f < frames; f++) { + MRIUSseq_vox(mri_morphed, x, y, z, f) = (unsigned short)MRIFseq_vox(mri_s_morphed, x, y, z, f); + } + break; case MRI_FLOAT: for (f = 0; f < frames; f++) { MRIFseq_vox(mri_morphed, x, y, z, f) = MRIFseq_vox(mri_s_morphed, x, y, z, f); @@ -7259,6 +7264,9 @@ MRI *GCAMbuildMostLikelyVolume(GCA_MORPH *gcam, MRI *mri) case MRI_SHORT: MRISseq_vox(mri, x, y, z, n) = nint(val); break; + case MRI_USHRT: + MRIUSseq_vox(mri, x, y, z, n) = nint(val); + break; case MRI_UCHAR: MRIseq_vox(mri, x, y, z, n) = nint(val); break; @@ -7303,6 +7311,9 @@ MRI *GCAMbuildMostLikelyVolume(GCA_MORPH *gcam, MRI *mri) case MRI_SHORT: MRISseq_vox(mri, x, y, z, n) = nint(val); break; + case MRI_USHRT: + MRIUSseq_vox(mri, x, y, z, n) = nint(val); + break; case MRI_UCHAR: MRIseq_vox(mri, x, y, z, n) = nint(val); break; @@ -7367,6 +7378,9 @@ MRI *GCAMbuildVolume(GCA_MORPH *gcam, MRI *mri) case MRI_SHORT: MRISseq_vox(mri, x, y, z, n) = nint(val); break; + case MRI_USHRT: + MRIUSseq_vox(mri, x, y, z, n) = nint(val); + break; case MRI_UCHAR: MRIseq_vox(mri, x, y, z, n) = nint(val); break; diff --git a/utils/machine.cpp b/utils/machine.cpp index b99d4ae2552..54c85645768 100644 --- a/utils/machine.cpp +++ b/utils/machine.cpp @@ -120,6 +120,28 @@ short swapShort(short s) return (ss.s); } +/* addition for unsigned short */ +typedef union { + unsigned short s; + unsigned char buf[sizeof(unsigned short)]; +} SWAP_USHORT; + +unsigned short swapUShort(unsigned short s) +{ + SWAP_USHORT ss; + unsigned char c; + + /* first swap bytes in word */ + // ss.s and ss.buf share the same memory location + ss.s = s; + c = ss.buf[0]; + ss.buf[0] = ss.buf[1]; + ss.buf[1] = c; + + return (ss.s); +} +/* end of addition for unsigned short */ + typedef union { double d; long l[sizeof(double) / sizeof(long)]; diff --git a/utils/mri.cpp b/utils/mri.cpp index 29c8230f31c..6c2bd379f9a 100644 --- a/utils/mri.cpp +++ b/utils/mri.cpp @@ -219,6 +219,8 @@ void MRI::initSlices() ptr = (void *)((unsigned char*)ptr + vox_per_row); break; case MRI_SHORT: ptr = (void *)((short *)ptr + vox_per_row); break; + case MRI_USHRT: + ptr = (void *)((unsigned short *)ptr + vox_per_row); break; case MRI_RGB: case MRI_INT: ptr = (void *)((int *)ptr + vox_per_row); break; @@ -1339,6 +1341,7 @@ size_t MRIsizeof(int mritype) case MRI_FLOAT: return sizeof(float); case MRI_TENSOR: return sizeof(float); case MRI_SHORT: return sizeof(short); + case MRI_USHRT: return sizeof(unsigned short); } return (-1); // should never get here } @@ -1363,6 +1366,9 @@ double MRIptr2dbl(void *pmric, int mritype) case MRI_SHORT: v = (double)(*((short *)pmric)); break; + case MRI_USHRT: + v = (double)(*((unsigned short *)pmric)); + break; case MRI_INT: v = (double)(*((int *)pmric)); break; @@ -1395,6 +1401,9 @@ void MRIdbl2ptr(double v, void *pmric, int mritype) case MRI_SHORT: *((short *)pmric) = nint(v); break; + case MRI_USHRT: + *((unsigned short *)pmric) = nint(v); + break; case MRI_INT: *((int *)pmric) = nint(v); break; @@ -1492,6 +1501,9 @@ float MRIgetVoxVal(const MRI *mri, int c, int r, int s, int f) case MRI_SHORT: return (float)* ((short *)mri->chunk + c + r * mri->vox_per_row + s * mri->vox_per_slice + f * mri->vox_per_vol); break; + case MRI_USHRT: + return (float)* ((unsigned short *)mri->chunk + c + r * mri->vox_per_row + s * mri->vox_per_slice + f * mri->vox_per_vol); + break; case MRI_RGB: case MRI_INT: return (float)* ((int *)mri->chunk + c + r * mri->vox_per_row + s * mri->vox_per_slice + f * mri->vox_per_vol); @@ -1512,6 +1524,9 @@ float MRIgetVoxVal(const MRI *mri, int c, int r, int s, int f) case MRI_SHORT: return ((float)MRISseq_vox(mri, c, r, s, f)); break; + case MRI_USHRT: + return ((float)MRIUSseq_vox(mri, c, r, s, f)); + break; case MRI_RGB: case MRI_INT: return ((float)MRIIseq_vox(mri, c, r, s, f)); @@ -1549,6 +1564,10 @@ int MRIsetVoxVal(MRI *mri, int c, int r, int s, int f, float voxval) if (voxval < SHORT_MIN) voxval = SHORT_MIN; if (voxval > SHORT_MAX) voxval = SHORT_MAX; break; + case MRI_USHRT: + if (voxval < 0) voxval = 0; + if (voxval > USHRT_MAX) voxval = USHRT_MAX; + break; case MRI_INT: if (voxval < INT_MIN) voxval = INT_MIN; if (voxval > INT_MAX) voxval = INT_MAX; @@ -1567,6 +1586,9 @@ int MRIsetVoxVal(MRI *mri, int c, int r, int s, int f, float voxval) case MRI_SHORT: *((short *)mri->chunk + c + r * mri->vox_per_row + s * mri->vox_per_slice + f * mri->vox_per_vol) = nint(voxval); break; + case MRI_USHRT: + *((unsigned short *)mri->chunk + c + r * mri->vox_per_row + s * mri->vox_per_slice + f * mri->vox_per_vol) = nint(voxval); + break; case MRI_RGB: case MRI_INT: *((int *)mri->chunk + c + r * mri->vox_per_row + s * mri->vox_per_slice + f * mri->vox_per_vol) = nint(voxval); @@ -1587,6 +1609,9 @@ int MRIsetVoxVal(MRI *mri, int c, int r, int s, int f, float voxval) case MRI_SHORT: MRISseq_vox(mri, c, r, s, f) = nint(voxval); break; + case MRI_USHRT: + MRIUSseq_vox(mri, c, r, s, f) = nint(voxval); + break; case MRI_RGB: case MRI_INT: MRIIseq_vox(mri, c, r, s, f) = nint(voxval); @@ -1650,6 +1675,7 @@ int MRIprecisionCode(const char *PrecisionString) { if (!strcasecmp(PrecisionString, "uchar")) return (MRI_UCHAR); if (!strcasecmp(PrecisionString, "short")) return (MRI_SHORT); + if (!strcasecmp(PrecisionString, "ushrt")) return (MRI_USHRT); if (!strcasecmp(PrecisionString, "int")) return (MRI_INT); if (!strcasecmp(PrecisionString, "long")) return (MRI_LONG); if (!strcasecmp(PrecisionString, "float")) return (MRI_FLOAT); @@ -1671,6 +1697,9 @@ const char *MRIprecisionString(int PrecisionCode) case MRI_SHORT: return ("short"); break; + case MRI_USHRT: + return ("ushrt"); + break; case MRI_INT: return ("int"); break; @@ -1764,6 +1793,7 @@ MRI *MRIscalarMul(MRI *mri_src, MRI *mri_dst, float scalar) BUFTYPE *psrc, *pdst; float *pfsrc, *pfdst, dval; short *pssrc, *psdst; + unsigned short *pussrc, *pusdst; width = mri_src->width; height = mri_src->height; @@ -1779,7 +1809,7 @@ MRI *MRIscalarMul(MRI *mri_src, MRI *mri_dst, float scalar) MRIsetVoxVal(mri_dst, x, y, z, frame, dval * scalar); } } - else + else // not used switch (mri_src->type) { case MRI_UCHAR: psrc = &MRIseq_vox(mri_src, 0, y, z, frame); @@ -1801,6 +1831,11 @@ MRI *MRIscalarMul(MRI *mri_src, MRI *mri_dst, float scalar) psdst = &MRISseq_vox(mri_dst, 0, y, z, frame); for (x = 0; x < width; x++) *psdst++ = (short)nint((float)*pssrc++ * scalar); break; + case MRI_USHRT: + pussrc = &MRIUSseq_vox(mri_src, 0, y, z, frame); + pusdst = &MRIUSseq_vox(mri_dst, 0, y, z, frame); + for (x = 0; x < width; x++) *pusdst++ = (unsigned short)nint((float)*pussrc++ * scalar); + break; default: ErrorReturn(NULL, (ERROR_UNSUPPORTED, "MRIscalarMul: unsupported type %d", mri_src->type)); } @@ -1817,6 +1852,7 @@ MRI *MRIscalarMulFrame(MRI *mri_src, MRI *mri_dst, float scalar, int frame) BUFTYPE *psrc, *pdst; float *pfsrc, *pfdst, dval; short *pssrc, *psdst; + unsigned short *pussrc, *pusdst; width = mri_src->width; height = mri_src->height; @@ -1846,6 +1882,11 @@ MRI *MRIscalarMulFrame(MRI *mri_src, MRI *mri_dst, float scalar, int frame) psdst = &MRISseq_vox(mri_dst, 0, y, z, frame); for (x = 0; x < width; x++) *psdst++ = (short)nint((float)*pssrc++ * scalar); break; + case MRI_USHRT: + pussrc = &MRIUSseq_vox(mri_src, 0, y, z, frame); + pusdst = &MRIUSseq_vox(mri_dst, 0, y, z, frame); + for (x = 0; x < width; x++) *pusdst++ = (unsigned short)nint((float)*pussrc++ * scalar); + break; default: ErrorReturn(NULL, (ERROR_UNSUPPORTED, "MRIscalarMulFrame: unsupported type %d", mri_src->type)); } @@ -1940,6 +1981,19 @@ int MRIvalRange(MRI *mri, float *pmin, float *pmax) } } break; + case MRI_USHRT: + for (frame = 0; frame < mri->nframes; frame++) { + for (z = 0; z < depth; z++) { + for (y = 0; y < height; y++) { + for (x = 0; x < width; x++) { + val = (float)MRIUSseq_vox(mri, x, y, z, frame); + if (val < fmin) fmin = val; + if (val > fmax) fmax = val; + } + } + } + } + break; case MRI_UCHAR: for (frame = 0; frame < mri->nframes; frame++) { for (z = 0; z < depth; z++) { @@ -2050,6 +2104,17 @@ int MRIvalRangeFrame(MRI *mri, float *pmin, float *pmax, int frame) } } break; + case MRI_USHRT: + for (z = 0; z < depth; z++) { + for (y = 0; y < height; y++) { + for (x = 0; x < width; x++) { + val = (float)MRIUSseq_vox(mri, x, y, z, frame); + if (val < fmin) fmin = val; + if (val > fmax) fmax = val; + } + } + } + break; case MRI_UCHAR: for (z = 0; z < depth; z++) { for (y = 0; y < height; y++) { @@ -2210,6 +2275,18 @@ MRI *MRIvalScale(MRI *mri_src, MRI *mri_dst, float flo, float fhi) } } break; + case MRI_USHRT: + for (z = 0; z < depth; z++) { + for (y = 0; y < height; y++) { + unsigned short *pus_src = &MRIUSvox(mri_src, 0, y, z); + unsigned short *pus_dst = &MRIUSvox(mri_dst, 0, y, z); + for (x = 0; x < width; x++) { + val = (float)(*pus_src++); + *pus_dst++ = (unsigned short)nint((val - fmin) * scale + flo); + } + } + } + break; case MRI_UCHAR: for (z = 0; z < depth; z++) { for (y = 0; y < height; y++) { @@ -2343,6 +2420,36 @@ int MRIboundingBoxNbhd(MRI *mri, int thresh, int wsize, MRI_REGION *box) } } break; + case MRI_USHRT: + for (z = 0; z < depth; z++) { + for (y = 0; y < height; y++) { + unsigned short *pussrc = &MRIUSvox(mri, 0, y, z); + for (x = 0; x < width; x++) { + if (*pussrc++ > thresh) { + in_brain = 1; + for (zk = -whalf; in_brain && zk <= whalf; zk++) { + zi = mri->zi[z + zk]; + for (yk = -whalf; in_brain && yk <= whalf; yk++) { + yi = mri->yi[y + yk]; + for (xk = -whalf; in_brain && xk <= whalf; xk++) { + xi = mri->xi[x + xk]; + if (MRIUSvox(mri, xi, yi, zi) < thresh) in_brain = 0; + } + } + } + if (in_brain) { + if (x < box->x) box->x = x; + if (y < box->y) box->y = y; + if (z < box->z) box->z = z; + if (x > x1) x1 = x; + if (y > y1) y1 = y; + if (z > z1) z1 = z; + } + } + } + } + } + break; case MRI_FLOAT: for (z = 0; z < depth; z++) { for (y = 0; y < height; y++) { @@ -2725,6 +2832,23 @@ int MRIboundingBox(MRI *mri, int thresh, MRI_REGION *box) } } break; + case MRI_USHRT: + for (z = 0; z < depth; z++) { + for (y = 0; y < height; y++) { + unsigned short *pussrc = &MRIUSvox(mri, 0, y, z); + for (x = 0; x < width; x++) { + if (*pussrc++ > thresh) { + if (x < box->x) box->x = x; + if (y < box->y) box->y = y; + if (z < box->z) box->z = z; + if (x > x1) x1 = x; + if (y > y1) y1 = y; + if (z > z1) z1 = z; + } + } + } + } + break; case MRI_INT: for (z = 0; z < depth; z++) { for (y = 0; y < height; y++) { @@ -3662,6 +3786,9 @@ MRI *MRIextractInto(MRI *mri_src, MRI *mri_dst, int x0, int y0, int z0, int dx, case MRI_SHORT: bytes *= sizeof(short); break; + case MRI_USHRT: + bytes *= sizeof(unsigned short); + break; case MRI_UCHAR: break; } @@ -3682,6 +3809,9 @@ MRI *MRIextractInto(MRI *mri_src, MRI *mri_dst, int x0, int y0, int z0, int dx, case MRI_SHORT: memmove(&MRISseq_vox(mri_dst, x1, yd, zd, frame), &MRISseq_vox(mri_src, x0, ys, zs, frame), bytes); break; + case MRI_USHRT: + memmove(&MRIUSseq_vox(mri_dst, x1, yd, zd, frame), &MRIUSseq_vox(mri_src, x0, ys, zs, frame), bytes); + break; case MRI_LONG: memmove(&MRILseq_vox(mri_dst, x1, yd, zd, frame), &MRILseq_vox(mri_src, x0, ys, zs, frame), bytes); break; @@ -3964,6 +4094,9 @@ int MRIclear(MRI *mri) case MRI_SHORT: bytes *= sizeof(short); break; + case MRI_USHRT: + bytes *= sizeof(unsigned short); + break; default: ErrorReturn(ERROR_UNSUPPORTED, (ERROR_UNSUPPORTED, "MRIclear: unsupported input type %d", mri->type)); break; @@ -4656,6 +4789,7 @@ MRI *MRIsubtract(MRI *mri1, MRI *mri2, MRI *mri_dst) float v1, v2, v = 0.0; BUFTYPE *p1 = NULL, *p2 = NULL, *pdst = NULL; short *ps1 = NULL, *ps2 = NULL, *psdst = NULL; + unsigned short *pus1 = NULL, *pus2 = NULL, *pusdst = NULL; int *pi1 = NULL, *pi2 = NULL, *pidst = NULL; long *pl1 = NULL, *pl2 = NULL, *pldst = NULL; float *pf1 = NULL, *pf2 = NULL, *pfdst = NULL; @@ -4699,6 +4833,9 @@ MRI *MRIsubtract(MRI *mri1, MRI *mri2, MRI *mri_dst) case MRI_SHORT: psdst = (short *)mri_dst->slices[s][y]; break; + case MRI_USHRT: + pusdst = (unsigned short *)mri_dst->slices[s][y]; + break; case MRI_INT: pidst = (int *)mri_dst->slices[s][y]; break; @@ -4718,6 +4855,10 @@ MRI *MRIsubtract(MRI *mri1, MRI *mri2, MRI *mri_dst) ps1 = (short *)mri1->slices[s][y]; ps2 = (short *)mri2->slices[s][y]; break; + case MRI_USHRT: + pus1 = (unsigned short *)mri1->slices[s][y]; + pus2 = (unsigned short *)mri2->slices[s][y]; + break; case MRI_INT: pi1 = (int *)mri1->slices[s][y]; pi2 = (int *)mri2->slices[s][y]; @@ -4740,6 +4881,9 @@ MRI *MRIsubtract(MRI *mri1, MRI *mri2, MRI *mri_dst) case MRI_SHORT: v = (float)(*ps1++) - (float)(*ps2++); break; + case MRI_USHRT: + v = (float)(*pus1++) - (float)(*pus2++); + break; case MRI_INT: v = (float)(*pi1++) - (float)(*pi2++); break; @@ -4758,6 +4902,9 @@ MRI *MRIsubtract(MRI *mri1, MRI *mri2, MRI *mri_dst) case MRI_SHORT: (*psdst++) = (short)nint(v); break; + case MRI_USHRT: + (*pusdst++) = (unsigned short)nint(v); + break; case MRI_INT: (*pidst++) = (int)nint(v); break; @@ -4925,6 +5072,7 @@ MRI *MRIadd(MRI *mri1, MRI *mri2, MRI *mri_dst) float v1, v2, v = 0.0; BUFTYPE *p1 = NULL, *p2 = NULL, *pdst = NULL; short *ps1 = NULL, *ps2 = NULL, *psdst = NULL; + unsigned short *pus1 = NULL, *pus2 = NULL, *pusdst = NULL; int *pi1 = NULL, *pi2 = NULL, *pidst = NULL; long *pl1 = NULL, *pl2 = NULL, *pldst = NULL; float *pf1 = NULL, *pf2 = NULL, *pfdst = NULL; @@ -4970,6 +5118,9 @@ MRI *MRIadd(MRI *mri1, MRI *mri2, MRI *mri_dst) case MRI_SHORT: psdst = (short *)mri_dst->slices[s][y]; break; + case MRI_USHRT: + pusdst = (unsigned short *)mri_dst->slices[s][y]; + break; case MRI_INT: pidst = (int *)mri_dst->slices[s][y]; break; @@ -4990,6 +5141,10 @@ MRI *MRIadd(MRI *mri1, MRI *mri2, MRI *mri_dst) ps1 = (short *)mri1->slices[s][y]; ps2 = (short *)mri2->slices[s][y]; break; + case MRI_USHRT: + pus1 = (unsigned short *)mri1->slices[s][y]; + pus2 = (unsigned short *)mri2->slices[s][y]; + break; case MRI_INT: pi1 = (int *)mri1->slices[s][y]; pi2 = (int *)mri2->slices[s][y]; @@ -5015,6 +5170,9 @@ MRI *MRIadd(MRI *mri1, MRI *mri2, MRI *mri_dst) case MRI_SHORT: (*pdst++) = (BUFTYPE)((*ps1++) + (*ps2++)); break; + case MRI_USHRT: + (*pdst++) = (BUFTYPE)((*pus1++) + (*pus2++)); + break; case MRI_INT: (*pdst++) = (BUFTYPE)((*pi1++) + (*pi2++)); break; @@ -5034,6 +5192,9 @@ MRI *MRIadd(MRI *mri1, MRI *mri2, MRI *mri_dst) case MRI_SHORT: (*psdst++) = (short)((*ps1++) + (*ps2++)); break; + case MRI_USHRT: + (*psdst++) = (short)((*pus1++) + (*pus2++)); + break; case MRI_INT: (*psdst++) = (short)((*pi1++) + (*pi2++)); break; @@ -5045,6 +5206,28 @@ MRI *MRIadd(MRI *mri1, MRI *mri2, MRI *mri_dst) break; } break; + case MRI_USHRT: + switch (mri1->type) { + case MRI_UCHAR: + (*pusdst++) = ((unsigned short)(*p1++) + (*p2++)); + break; + case MRI_SHORT: + (*pusdst++) = (unsigned short)((*ps1++) + (*ps2++)); + break; + case MRI_USHRT: + (*pusdst++) = (unsigned short)((*pus1++) + (*pus2++)); + break; + case MRI_INT: + (*pusdst++) = (unsigned short)((*pi1++) + (*pi2++)); + break; + case MRI_LONG: + (*pusdst++) = (unsigned short)((*pl1++) + (*pl2++)); + break; + case MRI_FLOAT: + (*pusdst++) = (unsigned short)nint((*pf1++) + (*pf2++)); + break; + } + break; case MRI_INT: switch (mri1->type) { case MRI_UCHAR: @@ -5053,6 +5236,9 @@ MRI *MRIadd(MRI *mri1, MRI *mri2, MRI *mri_dst) case MRI_SHORT: (*pidst++) = ((int)(*ps1++) + (*ps2++)); break; + case MRI_USHRT: + (*pidst++) = ((int)(*pus1++) + (*pus2++)); + break; case MRI_INT: (*pidst++) = (int)((*pi1++) + (*pi2++)); break; @@ -5072,6 +5258,9 @@ MRI *MRIadd(MRI *mri1, MRI *mri2, MRI *mri_dst) case MRI_SHORT: (*pldst++) = ((long)(*ps1++) + (*ps2++)); break; + case MRI_USHRT: + (*pldst++) = ((long)(*pus1++) + (*pus2++)); + break; case MRI_INT: (*pldst++) = ((long)(*pi1++) + (*pi2++)); break; @@ -5091,6 +5280,9 @@ MRI *MRIadd(MRI *mri1, MRI *mri2, MRI *mri_dst) case MRI_SHORT: (*pfdst++) = ((float)(*ps1++) + (*ps2++)); break; + case MRI_USHRT: + (*pfdst++) = ((float)(*pus1++) + (*pus2++)); + break; case MRI_INT: (*pfdst++) = ((float)(*pi1++) + (*pi2++)); break; @@ -5376,6 +5568,9 @@ MRI *MRIcopy(MRI *mri_src, MRI *mri_dst) case MRI_SHORT: bytes *= sizeof(short); break; + case MRI_USHRT: + bytes *= sizeof(unsigned short); + break; case MRI_FLOAT: bytes *= sizeof(float); break; @@ -5413,6 +5608,20 @@ MRI *MRIcopy(MRI *mri_src, MRI *mri_dst) } } break; + case MRI_USHRT: /* float --> unsigned short */ + for (frame = 0; frame < mri_src->nframes; frame++) { + for (z = 0; z < depth; z++) { + for (y = 0; y < height; y++) { + unsigned short *usdst = &MRIUSseq_vox(mri_dst, 0, y, z, frame); + fsrc = &MRIFseq_vox(mri_src, 0, y, z, frame); + for (x = 0; x < width; x++) { + val = nint(*fsrc++); + *usdst++ = (unsigned short)val; + } + } + } + } + break; case MRI_UCHAR: /* float --> unsigned char */ for (frame = 0; frame < mri_src->nframes; frame++) { for (z = 0; z < depth; z++) { @@ -6797,6 +7006,9 @@ MRI *ImageToMRI(IMAGE *I) case MRI_SHORT: MRISseq_vox(mri, x, y, 0, frames) = *(IMAGESpix(I, x, yp) + frame_offset); break; + case MRI_USHRT: + MRIUSseq_vox(mri, x, y, 0, frames) = *(IMAGESpix(I, x, yp) + frame_offset); + break; case MRI_INT: { // int val; @@ -7202,6 +7414,7 @@ MRI *MRIdownsample2(MRI *mri_src, MRI *mri_dst) int width, depth, height, x, y, z, x1, y1, z1; BUFTYPE *psrc; short *pssrc; + unsigned short *pussrc; float *pfsrc; float val; @@ -7232,6 +7445,10 @@ MRI *MRIdownsample2(MRI *mri_src, MRI *mri_dst) pssrc = &MRISvox(mri_src, 2 * x, y1, z1); for (x1 = 2 * x; x1 <= 2 * x + 1; x1++) val += *pssrc++; break; + case MRI_USHRT: + pussrc = &MRIUSvox(mri_src, 2 * x, y1, z1); + for (x1 = 2 * x; x1 <= 2 * x + 1; x1++) val += *pussrc++; + break; case MRI_FLOAT: pfsrc = &MRIFvox(mri_src, 2 * x, y1, z1); for (x1 = 2 * x; x1 <= 2 * x + 1; x1++) val += *pfsrc++; @@ -7251,6 +7468,9 @@ MRI *MRIdownsample2(MRI *mri_src, MRI *mri_dst) case MRI_SHORT: MRISvox(mri_dst, x, y, z) = (short)nint(val / 8.0f); break; + case MRI_USHRT: + MRIUSvox(mri_dst, x, y, z) = (unsigned short)nint(val / 8.0f); + break; } } } @@ -7388,6 +7608,7 @@ MRI *MRIdownsampleNOld(MRI *mri_src, MRI *mri_dst, int N) int width, depth, height, x, y, z, x1, y1, z1; BUFTYPE *psrc; short *pssrc; + unsigned short *pussrc; float *pfsrc; float val; @@ -7418,6 +7639,10 @@ MRI *MRIdownsampleNOld(MRI *mri_src, MRI *mri_dst, int N) pssrc = &MRISvox(mri_src, N * x, y1, z1); for (x1 = N * x; x1 <= N * x + (N - 1); x1++) val += *pssrc++; break; + case MRI_USHRT: + pussrc = &MRIUSvox(mri_src, N * x, y1, z1); + for (x1 = N * x; x1 <= N * x + (N - 1); x1++) val += *pussrc++; + break; case MRI_FLOAT: pfsrc = &MRIFvox(mri_src, N * x, y1, z1); for (x1 = N * x; x1 <= N * x + (N - 1); x1++) val += *pfsrc++; @@ -7437,6 +7662,9 @@ MRI *MRIdownsampleNOld(MRI *mri_src, MRI *mri_dst, int N) case MRI_SHORT: MRISvox(mri_dst, x, y, z) = (short)nint(val / (N * N * N)); break; + case MRI_USHRT: + MRIUSvox(mri_dst, x, y, z) = (unsigned short)nint(val / (N * N * N)); + break; } } } @@ -7492,6 +7720,9 @@ int MRIvalueFill(MRI *mri, float value) case MRI_SHORT: MRISseq_vox(mri, c, r, s, f) = (short)(nint(value)); break; + case MRI_USHRT: + MRIUSseq_vox(mri, c, r, s, f) = (unsigned short)(nint(value)); + break; case MRI_INT: MRIIseq_vox(mri, c, r, s, f) = (int)(nint(value)); break; @@ -8601,6 +8832,16 @@ int MRIsampleVolumeFrame(const MRI *mri, double x, double y, double z, const int xmd * ymd * zpd * (double)MRISseq_vox(mri, xp, yp, zm, frame) + xmd * ymd * zmd * (double)MRISseq_vox(mri, xp, yp, zp, frame); break; + case MRI_USHRT: + *pval = val = xpd * ypd * zpd * (double)MRIUSseq_vox(mri, xm, ym, zm, frame) + + xpd * ypd * zmd * (double)MRIUSseq_vox(mri, xm, ym, zp, frame) + + xpd * ymd * zpd * (double)MRIUSseq_vox(mri, xm, yp, zm, frame) + + xpd * ymd * zmd * (double)MRIUSseq_vox(mri, xm, yp, zp, frame) + + xmd * ypd * zpd * (double)MRIUSseq_vox(mri, xp, ym, zm, frame) + + xmd * ypd * zmd * (double)MRIUSseq_vox(mri, xp, ym, zp, frame) + + xmd * ymd * zpd * (double)MRIUSseq_vox(mri, xp, yp, zm, frame) + + xmd * ymd * zmd * (double)MRIUSseq_vox(mri, xp, yp, zp, frame); + break; case MRI_INT: *pval = val = xpd * ypd * zpd * (double)MRIIseq_vox(mri, xm, ym, zm, frame) + xpd * ypd * zmd * (double)MRIIseq_vox(mri, xm, ym, zp, frame) + @@ -8894,6 +9135,9 @@ int MRIsampleVolumeType(const MRI *mri, double x, double y, double z, double *pv case MRI_SHORT: *pval = (float)MRISvox(mri, xv, yv, zv); break; + case MRI_USHRT: + *pval = (float)MRIUSvox(mri, xv, yv, zv); + break; case MRI_INT: *pval = (float)MRIIvox(mri, xv, yv, zv); break; @@ -8966,6 +9210,9 @@ int MRIsampleVolumeFrameType( case MRI_SHORT: *pval = (float)MRISseq_vox(mri, xv, yv, zv, frame); break; + case MRI_USHRT: + *pval = (float)MRIUSseq_vox(mri, xv, yv, zv, frame); + break; case MRI_INT: *pval = (float)MRIIseq_vox(mri, xv, yv, zv, frame); break; @@ -9026,6 +9273,10 @@ int MRIsampleVolumeFrameType_xyzInt_nRange_SAMPLE_NEAREST_floats(const MRI *mr for (frame = frameBegin; frame < frameEnd; frame++) valForEachFrame[frame - frameBegin] = (float)MRISseq_vox(mri, xv, yv, zv, frame); break; + case MRI_USHRT: + for (frame = frameBegin; frame < frameEnd; frame++) + valForEachFrame[frame - frameBegin] = (float)MRIUSseq_vox(mri, xv, yv, zv, frame); + break; case MRI_INT: for (frame = frameBegin; frame < frameEnd; frame++) valForEachFrame[frame - frameBegin] = (float)MRIIseq_vox(mri, xv, yv, zv, frame); @@ -9118,6 +9369,16 @@ int MRIinterpolateIntoVolumeFrame(MRI *mri, double x, double y, double z, int fr MRISseq_vox(mri, xp, yp, zm, frame) += nint(xmd * ymd * zpd * val); MRISseq_vox(mri, xp, yp, zp, frame) += nint(xmd * ymd * zmd * val); break; + case MRI_USHRT: + MRIUSseq_vox(mri, xm, ym, zm, frame) += nint(xpd * ypd * zpd * val); + MRIUSseq_vox(mri, xm, ym, zp, frame) += nint(xpd * ypd * zmd * val); + MRIUSseq_vox(mri, xm, yp, zm, frame) += nint(xpd * ymd * zpd * val); + MRIUSseq_vox(mri, xm, yp, zp, frame) += nint(xpd * ymd * zmd * val); + MRIUSseq_vox(mri, xp, ym, zm, frame) += nint(xmd * ypd * zpd * val); + MRIUSseq_vox(mri, xp, ym, zp, frame) += nint(xmd * ypd * zmd * val); + MRIUSseq_vox(mri, xp, yp, zm, frame) += nint(xmd * ymd * zpd * val); + MRIUSseq_vox(mri, xp, yp, zp, frame) += nint(xmd * ymd * zmd * val); + break; case MRI_INT: MRIIseq_vox(mri, xm, ym, zm, frame) += nint(xpd * ypd * zpd * val); MRIIseq_vox(mri, xm, ym, zp, frame) += nint(xpd * ypd * zmd * val); @@ -9213,6 +9474,13 @@ int MRIsampleVolume(const MRI *mri, double x, double y, double z, double *pval) xmd * ypd * zpd * (double)MRISvox(mri, xp, ym, zm) + xmd * ypd * zmd * (double)MRISvox(mri, xp, ym, zp) + xmd * ymd * zpd * (double)MRISvox(mri, xp, yp, zm) + xmd * ymd * zmd * (double)MRISvox(mri, xp, yp, zp); break; + case MRI_USHRT: + *pval = val = + xpd * ypd * zpd * (double)MRIUSvox(mri, xm, ym, zm) + xpd * ypd * zmd * (double)MRIUSvox(mri, xm, ym, zp) + + xpd * ymd * zpd * (double)MRIUSvox(mri, xm, yp, zm) + xpd * ymd * zmd * (double)MRIUSvox(mri, xm, yp, zp) + + xmd * ypd * zpd * (double)MRIUSvox(mri, xp, ym, zm) + xmd * ypd * zmd * (double)MRIUSvox(mri, xp, ym, zp) + + xmd * ymd * zpd * (double)MRIUSvox(mri, xp, yp, zm) + xmd * ymd * zmd * (double)MRIUSvox(mri, xp, yp, zp); + break; case MRI_INT: *pval = val = xpd * ypd * zpd * (double)MRIIvox(mri, xm, ym, zm) + xpd * ypd * zmd * (double)MRIIvox(mri, xm, ym, zp) + @@ -9391,6 +9659,16 @@ int MRIsampleSeqVolume(const MRI *mri, double x, double y, double z, float *valv xmd * ymd * zpd * (double)MRISseq_vox(mri, xp, yp, zm, f) + xmd * ymd * zmd * (double)MRISseq_vox(mri, xp, yp, zp, f); break; + case MRI_USHRT: + valvect[f] = xpd * ypd * zpd * (double)MRIUSseq_vox(mri, xm, ym, zm, f) + + xpd * ypd * zmd * (double)MRIUSseq_vox(mri, xm, ym, zp, f) + + xpd * ymd * zpd * (double)MRIUSseq_vox(mri, xm, yp, zm, f) + + xpd * ymd * zmd * (double)MRIUSseq_vox(mri, xm, yp, zp, f) + + xmd * ypd * zpd * (double)MRIUSseq_vox(mri, xp, ym, zm, f) + + xmd * ypd * zmd * (double)MRIUSseq_vox(mri, xp, ym, zp, f) + + xmd * ymd * zpd * (double)MRIUSseq_vox(mri, xp, yp, zm, f) + + xmd * ymd * zmd * (double)MRIUSseq_vox(mri, xp, yp, zp, f); + break; case MRI_INT: valvect[f] = xpd * ypd * zpd * (double)MRIIseq_vox(mri, xm, ym, zm, f) + xpd * ypd * zmd * (double)MRIIseq_vox(mri, xm, ym, zp, f) + @@ -9528,6 +9806,9 @@ int MRIsampleSeqVolumeType( xmd * ymd * zpd * (double)MRISseq_vox(mri, xp, yp, zm, f) + xmd * ymd * zmd * (double)MRISseq_vox(mri, xp, yp, zp, f) ;*/ break; + case MRI_USHRT: + valvect[f] = (double)MRIUSseq_vox(mri, xv, yv, zv, f); + break; case MRI_INT: valvect[f] = (double)MRIIseq_vox(mri, xv, yv, zv, f); /*xpd * ypd * zpd * (double)MRIIseq_vox(mri, xm, ym, zm, f) + @@ -9772,6 +10053,9 @@ int MRIcubicSampleVolume(const MRI *mri, double x, double y, double z, double *p case MRI_SHORT: vv[ix][iy][iz] = (double)MRISvox(mri, ix_low - 1 + ix, iy_low - 1 + iy, iz_low - 1 + iz); break; + case MRI_USHRT: + vv[ix][iy][iz] = (double)MRIUSvox(mri, ix_low - 1 + ix, iy_low - 1 + iy, iz_low - 1 + iz); + break; case MRI_INT: vv[ix][iy][iz] = (double)MRIIvox(mri, ix_low - 1 + ix, iy_low - 1 + iy, iz_low - 1 + iz); break; @@ -9976,6 +10260,9 @@ int MRIcubicSampleVolumeFrame(MRI *mri, double x, double y, double z, int frame, case MRI_SHORT: vv[ix][iy][iz] = (double)MRISseq_vox(mri, ix_low - 1 + ix, iy_low - 1 + iy, iz_low - 1 + iz, frame); break; + case MRI_USHRT: + vv[ix][iy][iz] = (double)MRIUSseq_vox(mri, ix_low - 1 + ix, iy_low - 1 + iy, iz_low - 1 + iz, frame); + break; case MRI_INT: vv[ix][iy][iz] = (double)MRIIseq_vox(mri, ix_low - 1 + ix, iy_low - 1 + iy, iz_low - 1 + iz, frame); break; @@ -10099,6 +10386,9 @@ int MRIsincSampleVolume(const MRI *mri, double x, double y, double z, int hw, do case MRI_SHORT: sum_x += (coeff_x[jx_rel] / coeff_x_sum) * (double)MRISvox(mri, jx1, jy1, jz1); break; + case MRI_USHRT: + sum_x += (coeff_x[jx_rel] / coeff_x_sum) * (double)MRIUSvox(mri, jx1, jy1, jz1); + break; case MRI_INT: sum_x += (coeff_x[jx_rel] / coeff_x_sum) * (double)MRIIvox(mri, jx1, jy1, jz1); break; @@ -10201,6 +10491,9 @@ int MRIsincSampleVolumeFrame(MRI *mri, double x, double y, double z, int frame, case MRI_SHORT: sum_x += (coeff_x[jx_rel] / coeff_x_sum) * (double)MRISseq_vox(mri, jx1, jy1, jz1, frame); break; + case MRI_USHRT: + sum_x += (coeff_x[jx_rel] / coeff_x_sum) * (double)MRIUSseq_vox(mri, jx1, jy1, jz1, frame); + break; case MRI_INT: sum_x += (coeff_x[jx_rel] / coeff_x_sum) * (double)MRIIseq_vox(mri, jx1, jy1, jz1, frame); break; @@ -11363,6 +11656,10 @@ MRI *MRIchangeType(MRI *src, int dest_type, float f_low, float f_high, int no_sc dest_min = SHORT_MIN; dest_max = SHORT_MAX; } + if (dest_type == MRI_USHRT) { + dest_min = 0; + dest_max = USHRT_MAX; + } if (dest_type == MRI_INT) { dest_min = INT_MIN; dest_max = INT_MAX; @@ -11729,6 +12026,32 @@ MRI *MRIresampleFill(MRI *src, MRI *template_vol, int resample_type, float fill_ : (float)MRISseq_vox(src, si + 1, sj + 1, sk + 1, nframe)); } + if (src->type == MRI_USHRT) { + val000 = + (!i_good_flag || !j_good_flag || !k_good_flag ? 0.0 : (float)MRIUSseq_vox(src, si, sj, sk, nframe)); + val001 = + (!i_good_flag || !j_good_flag || !k1_good_flag ? 0.0 + : (float)MRIUSseq_vox(src, si, sj, sk + 1, nframe)); + val010 = + (!i_good_flag || !j1_good_flag || !k_good_flag ? 0.0 + : (float)MRIUSseq_vox(src, si, sj + 1, sk, nframe)); + val011 = (!i_good_flag || !j1_good_flag || !k1_good_flag + ? 0.0 + : (float)MRIUSseq_vox(src, si, sj + 1, sk + 1, nframe)); + val100 = + (!i1_good_flag || !j_good_flag || !k_good_flag ? 0.0 + : (float)MRIUSseq_vox(src, si + 1, sj, sk, nframe)); + val101 = (!i1_good_flag || !j_good_flag || !k1_good_flag + ? 0.0 + : (float)MRIUSseq_vox(src, si + 1, sj, sk + 1, nframe)); + val110 = (!i1_good_flag || !j1_good_flag || !k_good_flag + ? 0.0 + : (float)MRIUSseq_vox(src, si + 1, sj + 1, sk, nframe)); + val111 = (!i1_good_flag || !j1_good_flag || !k1_good_flag + ? 0.0 + : (float)MRIUSseq_vox(src, si + 1, sj + 1, sk + 1, nframe)); + } + if (src->type == MRI_INT) { val000 = (!i_good_flag || !j_good_flag || !k_good_flag ? 0.0 : (float)MRIIseq_vox(src, si, sj, sk, nframe)); @@ -11955,6 +12278,11 @@ MRI *MRIresampleFill(MRI *src, MRI *template_vol, int resample_type, float fill_ if (val > SHORT_MAX) val = SHORT_MAX; MRISseq_vox(dest, di, dj, dk, nframe) = (short)nint(val); } + if (dest->type == MRI_USHRT) { + if (val < 0) val = 0; + if (val > USHRT_MAX) val = USHRT_MAX; + MRIUSseq_vox(dest, di, dj, dk, nframe) = (unsigned short)nint(val); + } if (dest->type == MRI_INT) { if (val < INT_MIN) val = INT_MIN; if (val > INT_MAX) val = INT_MAX; @@ -12141,6 +12469,23 @@ int MRIstats(MRI *mri, float *min, float *max, int *n_voxels, float *mean, float } } + if (mri->type == MRI_USHRT) { + *min = *max = (float)MRIUSseq_vox(mri, 0, 0, 0, 0); + + for (i = 0; i < mri->width; i++) + for (j = 0; j < mri->height; j++) + for (k = 0; k < mri->depth; k++) + for (t = 0; t < mri->nframes; t++) { + val = (float)MRIUSseq_vox(mri, i, j, k, t); + + if (val < *min) *min = val; + if (val > *max) *max = val; + + sum += val; + sq_sum += val * val; + } + } + if (mri->type == MRI_INT) { *min = *max = (float)MRILseq_vox(mri, 0, 0, 0, 0); @@ -12948,6 +13293,9 @@ MRI *MRIlog(MRI *in, MRI *mask, double a, double b, MRI *out) case MRI_SHORT: v = (double) *((short *)pin + nvox); break; + case MRI_USHRT: + v = (double) *((unsigned short *)pin + nvox); + break; case MRI_INT: v = (double) *((int *)pin + nvox); break; @@ -13140,6 +13488,7 @@ MRI *MRIdrand48(int ncols, int nrows, int nslices, int nframes, float min, float float range, v; BUFTYPE *pmri = NULL; short *psmri = NULL; + unsigned short *pusmri = NULL; int *pimri = NULL; long *plmri = NULL; float *pfmri = NULL; @@ -13170,6 +13519,9 @@ MRI *MRIdrand48(int ncols, int nrows, int nslices, int nframes, float min, float case MRI_SHORT: psmri = (short *)mri->slices[n][r]; break; + case MRI_USHRT: + pusmri = (unsigned short *)mri->slices[n][r]; + break; case MRI_INT: pimri = (int *)mri->slices[n][r]; break; @@ -13189,6 +13541,9 @@ MRI *MRIdrand48(int ncols, int nrows, int nslices, int nframes, float min, float case MRI_SHORT: *psmri++ = (short)nint(v); break; + case MRI_USHRT: + *pusmri++ = (unsigned short)nint(v); + break; case MRI_INT: *pimri++ = (int)nint(v); break; @@ -13256,6 +13611,7 @@ MRI *MRIconst(int ncols, int nrows, int nslices, int nframes, float val, MRI *mr int c, r, s, f, n; BUFTYPE *pmri = NULL; short *psmri = NULL; + unsigned short *pusmri = NULL; int *pimri = NULL; long *plmri = NULL; float *pfmri = NULL; @@ -13285,6 +13641,9 @@ MRI *MRIconst(int ncols, int nrows, int nslices, int nframes, float val, MRI *mr case MRI_SHORT: psmri = (short *)mri->slices[n][r]; break; + case MRI_USHRT: + pusmri = (unsigned short *)mri->slices[n][r]; + break; case MRI_INT: pimri = (int *)mri->slices[n][r]; break; @@ -13303,6 +13662,9 @@ MRI *MRIconst(int ncols, int nrows, int nslices, int nframes, float val, MRI *mr case MRI_SHORT: *psmri++ = (short)nint(val); break; + case MRI_USHRT: + *pusmri++ = (unsigned short)nint(val); + break; case MRI_INT: *pimri++ = (int)nint(val); break; @@ -13347,6 +13709,9 @@ int MRInormalizeSequence(MRI *mri, float target) case MRI_SHORT: MRISseq_vox(mri, x, y, z, frame) = MRISseq_vox(mri, x, y, z, frame) / norm; break; + case MRI_USHRT: + MRIUSseq_vox(mri, x, y, z, frame) = MRIUSseq_vox(mri, x, y, z, frame) / norm; + break; case MRI_FLOAT: MRIFseq_vox(mri, x, y, z, frame) /= norm; break; @@ -15802,6 +16167,8 @@ const char *MRItype2str(int type) return ("uchar"); case MRI_SHORT: return ("short"); + case MRI_USHRT: + return ("ushrt"); case MRI_INT: return ("int"); case MRI_LONG: diff --git a/utils/mri2.cpp b/utils/mri2.cpp index c07e4cceb51..26cd075ba3f 100644 --- a/utils/mri2.cpp +++ b/utils/mri2.cpp @@ -537,6 +537,9 @@ size_t mri_sizeof(MRI *vol) case MRI_SHORT: bytes = sizeof(short); break; + case MRI_USHRT: + bytes = sizeof(unsigned short); + break; case MRI_FLOAT: bytes = sizeof(float); break; @@ -596,6 +599,9 @@ MRI *mri_reshape(MRI *vol, int ncols, int nrows, int nslices, int nframes) case MRI_SHORT: MRISseq_vox(outvol, c2, r2, s2, f2) = MRISseq_vox(vol, c, r, s, f); break; + case MRI_USHRT: + MRIUSseq_vox(outvol, c2, r2, s2, f2) = MRIUSseq_vox(vol, c, r, s, f); + break; case MRI_FLOAT: MRIFseq_vox(outvol, c2, r2, s2, f2) = MRIFseq_vox(vol, c, r, s, f); break; diff --git a/utils/mri_topology.cpp b/utils/mri_topology.cpp index 92d69f411c0..b7effe097d6 100644 --- a/utils/mri_topology.cpp +++ b/utils/mri_topology.cpp @@ -1408,6 +1408,26 @@ static void setBorderValue(MRI *mri, int val, int dst) for (k = dst; k < depth - dst; k++) for (j = dst; j < height - dst; j++) MRISvox(mri, i, j, k) = val; break; + case MRI_USHRT: + k = dst; + for (i = dst; i < mri->width - dst; i++) + for (j = dst; j < height - dst; j++) MRIUSvox(mri, i, j, k) = val; + k = depth - 1 - dst; + for (i = dst; i < width - dst; i++) + for (j = dst; j < height - dst; j++) MRIUSvox(mri, i, j, k) = val; + j = dst; + for (i = dst; i < width - dst; i++) + for (k = dst; k < depth - dst; k++) MRIUSvox(mri, i, j, k) = val; + j = height - 1 - dst; + for (i = dst; i < width - dst; i++) + for (k = dst; k < depth - dst; k++) MRIUSvox(mri, i, j, k) = val; + i = dst; + for (k = dst; k < depth - dst; k++) + for (j = dst; j < height - dst; j++) MRIUSvox(mri, i, j, k) = val; + i = width - 1 - dst; + for (k = dst; k < depth - dst; k++) + for (j = dst; j < height - dst; j++) MRIUSvox(mri, i, j, k) = val; + break; case MRI_FLOAT: k = dst; for (i = dst; i < mri->width - dst; i++) diff --git a/utils/mrihisto.cpp b/utils/mrihisto.cpp index c99185eec96..41dd6f8ca9b 100644 --- a/utils/mrihisto.cpp +++ b/utils/mrihisto.cpp @@ -278,6 +278,19 @@ static HISTOGRAM *mriHistogramLabel(MRI *mri, int nbins, HISTOGRAM *histo, LABEL histo->counts[bin_no]++; } break; + case MRI_USHRT: + for (i = 0; i < label->n_points; i++) { + // MRIworldToVoxel(mri, label->lv[i].x, label->lv[i].y, label->lv[i].z, + // &xv, &yv, &zv) ; + MRIsurfaceRASToVoxel(mri, label->lv[i].x, label->lv[i].y, label->lv[i].z, &xv, &yv, &zv); + x = nint(xv); + y = nint(yv); + z = nint(zv); + val = (float)MRIUSvox(mri, x, y, z); + bin_no = nint((float)(val - fmin) / (float)bin_size); + histo->counts[bin_no]++; + } + break; case MRI_FLOAT: for (i = 0; i < label->n_points; i++) { // MRIworldToVoxel(mri, label->lv[i].x, label->lv[i].y, label->lv[i].z, @@ -392,6 +405,19 @@ static HISTOGRAM *mriHistogramRegion(MRI *mri, int nbins, HISTOGRAM *histo, MRI_ } } break; + case MRI_USHRT: + for (z = z0; z < depth; z++) { + for (y = y0; y < height; y++) { + unsigned short* uspsrc = &MRIUSvox(mri, x0, y, z); + for (x = x0; x < width; x++) { + bin_no = nint((float)(*uspsrc++ - fmin) / (float)histo->bin_size); + if (bin_no < 0) bin_no = 0; + if (bin_no >= histo->nbins) bin_no = histo->nbins - 1; + histo->counts[bin_no]++; + } + } + } + break; case MRI_FLOAT: for (z = z0; z < depth; z++) { for (y = y0; y < height; y++) { @@ -775,6 +801,11 @@ HISTOGRAM *MRIhistogramLabelRegion(MRI *mri, MRI *mri_labeled, MRI_REGION *regio bin_no = nint((float)(val - bmin) / (float)histo->bin_size); histo->counts[bin_no]++; break; + case MRI_USHRT: + val = MRIUSvox(mri, x, y, z); + bin_no = nint((float)(val - bmin) / (float)histo->bin_size); + histo->counts[bin_no]++; + break; case MRI_FLOAT: fval = MRIFvox(mri, x, y, z); bin_no = nint((fval - fmin) / (float)histo->bin_size); @@ -843,6 +874,15 @@ HISTOGRAM *MRIhistogramLabel(MRI *mri, MRI *mri_labeled, int label, int nbins) bin_no = histo->nbins - 1; histo->counts[bin_no]++; break; + case MRI_USHRT: + val = MRIUSvox(mri, x, y, z); + bin_no = nint((float)(val - bmin) / (float)histo->bin_size); + if (bin_no < 0) + bin_no = 0; + else if (bin_no >= histo->nbins) + bin_no = histo->nbins - 1; + histo->counts[bin_no]++; + break; case MRI_FLOAT: fval = MRIFvox(mri, x, y, z); bin_no = nint((fval - fmin) / (float)histo->bin_size); diff --git a/utils/mriio.cpp b/utils/mriio.cpp index 34f5b377c2a..b5701e30027 100644 --- a/utils/mriio.cpp +++ b/utils/mriio.cpp @@ -831,35 +831,57 @@ MRI *mri_read(const char *fname, int type, int volume_flag, int start_frame, int mri2->imnr0 = 1; mri2->imnr0 = mri2->nframes; - if (mri2->type == MRI_UCHAR) - for (t = 0; t < mri2->nframes; t++) - for (i = 0; i < mri2->width; i++) - for (j = 0; j < mri2->height; j++) - for (k = 0; k < mri2->depth; k++) MRIseq_vox(mri2, i, j, k, t) = MRIseq_vox(mri, i, j, k, t + start_frame); - - if (mri2->type == MRI_SHORT) - for (t = 0; t < mri2->nframes; t++) - for (i = 0; i < mri2->width; i++) - for (j = 0; j < mri2->height; j++) - for (k = 0; k < mri2->depth; k++) MRISseq_vox(mri2, i, j, k, t) = MRISseq_vox(mri, i, j, k, t + start_frame); - - if (mri2->type == MRI_INT) - for (t = 0; t < mri2->nframes; t++) - for (i = 0; i < mri2->width; i++) - for (j = 0; j < mri2->height; j++) - for (k = 0; k < mri2->depth; k++) MRIIseq_vox(mri2, i, j, k, t) = MRIIseq_vox(mri, i, j, k, t + start_frame); - - if (mri2->type == MRI_LONG) - for (t = 0; t < mri2->nframes; t++) - for (i = 0; i < mri2->width; i++) - for (j = 0; j < mri2->height; j++) - for (k = 0; k < mri2->depth; k++) MRILseq_vox(mri2, i, j, k, t) = MRILseq_vox(mri, i, j, k, t + start_frame); - - if (mri2->type == MRI_FLOAT) + if (0) // disable it for now, don't know a case to execute this code path // (getenv("FS_MGZIO_USEVOXELBUF")) + { + printf("mri_read(): copy voxel by row ...\n"); for (t = 0; t < mri2->nframes; t++) - for (i = 0; i < mri2->width; i++) - for (j = 0; j < mri2->height; j++) - for (k = 0; k < mri2->depth; k++) MRIFseq_vox(mri2, i, j, k, t) = MRIFseq_vox(mri, i, j, k, t + start_frame); + for (k = 0; k < mri2->depth; k++) { + //for (j = 0; j < mri2->height; j++) { + int bytes_to_copy = MRIsizeof(mri2->type) * mri2->width * mri2->height; + void *pbufsrc = &MRIseq_vox(mri, 0, 0, k, t + start_frame); + void *pbufdst = &MRIseq_vox(mri2, 0, 0, k, t); + memcpy(pbufdst, pbufsrc, bytes_to_copy); + //} + } + } + else + { + if (mri2->type == MRI_UCHAR) + for (t = 0; t < mri2->nframes; t++) + for (i = 0; i < mri2->width; i++) + for (j = 0; j < mri2->height; j++) + for (k = 0; k < mri2->depth; k++) MRIseq_vox(mri2, i, j, k, t) = MRIseq_vox(mri, i, j, k, t + start_frame); + + if (mri2->type == MRI_SHORT) + for (t = 0; t < mri2->nframes; t++) + for (i = 0; i < mri2->width; i++) + for (j = 0; j < mri2->height; j++) + for (k = 0; k < mri2->depth; k++) MRISseq_vox(mri2, i, j, k, t) = MRISseq_vox(mri, i, j, k, t + start_frame); + + if (mri2->type == MRI_USHRT) + for (t = 0; t < mri2->nframes; t++) + for (i = 0; i < mri2->width; i++) + for (j = 0; j < mri2->height; j++) + for (k = 0; k < mri2->depth; k++) MRIUSseq_vox(mri2, i, j, k, t) = MRIUSseq_vox(mri, i, j, k, t + start_frame); + + if (mri2->type == MRI_INT) + for (t = 0; t < mri2->nframes; t++) + for (i = 0; i < mri2->width; i++) + for (j = 0; j < mri2->height; j++) + for (k = 0; k < mri2->depth; k++) MRIIseq_vox(mri2, i, j, k, t) = MRIIseq_vox(mri, i, j, k, t + start_frame); + + if (mri2->type == MRI_LONG) + for (t = 0; t < mri2->nframes; t++) + for (i = 0; i < mri2->width; i++) + for (j = 0; j < mri2->height; j++) + for (k = 0; k < mri2->depth; k++) MRILseq_vox(mri2, i, j, k, t) = MRILseq_vox(mri, i, j, k, t + start_frame); + + if (mri2->type == MRI_FLOAT) + for (t = 0; t < mri2->nframes; t++) + for (i = 0; i < mri2->width; i++) + for (j = 0; j < mri2->height; j++) + for (k = 0; k < mri2->depth; k++) MRIFseq_vox(mri2, i, j, k, t) = MRIFseq_vox(mri, i, j, k, t + start_frame); + } if (nan_inf_check(mri) != NO_ERROR) { MRIfree(&mri); @@ -4787,10 +4809,10 @@ static dsr *ReadAnalyzeHeader(const char *hdrfile, int *swap, int *mritype, int } else if (hdr->dime.datatype == DT_UINT16) { // Can happen if this is nifti - printf( - "Unsigned short not supported, but trying to read it \n" - "in as a signed short. Will be ok if no vals >= 32k.\n"); - *mritype = MRI_SHORT; + printf("INFO: This is an unsigned short.\n"); + // "Unsigned short not supported, but trying to read it \n" + // "in as a signed short. Will be ok if no vals >= 32k.\n"); + *mritype = MRI_USHRT; *bytes_per_voxel = 2; } else if (hdr->dime.datatype == DT_SIGNED_INT) { @@ -5382,6 +5404,10 @@ static int analyzeWriteFrame(MRI *mri, const char *fname, int frame) hdr.dime.datatype = DT_SIGNED_SHORT; bytes_per_voxel = 2; } + else if (mri->type == MRI_USHRT) { + hdr.dime.datatype = DT_UINT16; + bytes_per_voxel = 2; + } /* --- assuming long and int are identical --- */ else if (mri->type == MRI_INT || mri->type == MRI_LONG) { hdr.dime.datatype = DT_SIGNED_INT; @@ -5654,6 +5680,10 @@ static int analyzeWrite4D(MRI *mri, const char *fname) hdr.dime.datatype = DT_SIGNED_SHORT; bytes_per_voxel = 2; } + else if (mri->type == MRI_USHRT) { + hdr.dime.datatype = DT_UINT16; + bytes_per_voxel = 2; + } /* --- assuming long and int are identical --- */ else if (mri->type == MRI_INT || mri->type == MRI_LONG) { hdr.dime.datatype = DT_SIGNED_INT; @@ -7668,9 +7698,9 @@ static MRI *nifti1Read(const char *fname, int read_volume) } else if (hdr.datatype == DT_UINT16) { // This will not always work ... - printf("INFO: this is an unsiged short. I'll try to read it, but\n"); - printf(" it might not work if there are values over 32k\n"); - fs_type = MRI_SHORT; + printf("INFO: This is an unsigned short.\n"); // I'll try to read it, but\n"); + //printf(" it might not work if there are values over 32k\n"); + fs_type = MRI_USHRT; bytes_per_voxel = 2; } else if (hdr.datatype == DT_SIGNED_INT) { @@ -8135,6 +8165,10 @@ static int nifti1Write(MRI *mri0, const char *fname) hdr.datatype = DT_SIGNED_SHORT; hdr.bitpix = 16; } + else if (mri->type == MRI_USHRT) { + hdr.datatype = DT_UINT16; + hdr.bitpix = 16; + } else if (mri->type == MRI_BITMAP) { ErrorReturn(ERROR_UNSUPPORTED, (ERROR_UNSUPPORTED, "nifti1Write(): data type MRI_BITMAP unsupported")); } @@ -8353,9 +8387,9 @@ static MRI *niiRead(const char *fname, int read_volume) } else if (hdr.datatype == DT_UINT16) { // This will not always work ... - printf("INFO: this is an unsiged short. I'll try to read it, but\n"); - printf(" it might not work if there are values over 32k\n"); - fs_type = MRI_SHORT; + printf("INFO: This is an unsigined short.\n"); // I'll try to read it, but\n"); + //printf(" it might not work if there are values over 32k\n"); + fs_type = MRI_USHRT; bytes_per_voxel = 2; } else if (hdr.datatype == DT_SIGNED_INT) { @@ -8873,9 +8907,9 @@ static MRI *niiRead3(MRIFSSTRUCT *mrifsStruct) } else if (hdr->datatype == DT_UINT16) { // This will not always work ... - printf("INFO: this is an unsiged short. I'll try to read it, but\n"); - printf(" it might not work if there are values over 32k\n"); - fs_type = MRI_SHORT; + printf("INFO: This is an unsigned short.\n");// I'll try to read it, but\n"); + //printf(" it might not work if there are values over 32k\n"); + fs_type = MRI_USHRT; bytes_per_voxel = 2; } else if (hdr->datatype == DT_SIGNED_INT) { @@ -9290,6 +9324,10 @@ static int niiWrite(MRI *mri0, const char *fname) hdr.datatype = DT_SIGNED_SHORT; hdr.bitpix = 16; } + else if (mri->type == MRI_USHRT) { + hdr.datatype = DT_UINT16; + hdr.bitpix = 16; + } else if (mri->type == MRI_BITMAP) { ErrorReturn(ERROR_UNSUPPORTED, (ERROR_UNSUPPORTED, "niiWrite(): data type MRI_BITMAP unsupported")); } @@ -10695,6 +10733,9 @@ static MRI *mghRead(const char *fname, int read_volume, int frame) case MRI_SHORT: bpv = sizeof(short); break; + case MRI_USHRT: + bpv = sizeof(unsigned short); + break; case MRI_INT: bpv = sizeof(int); break; @@ -10743,52 +10784,108 @@ static MRI *mghRead(const char *fname, int read_volume, int frame) buf = (BUFTYPE *)calloc(bytes, sizeof(BUFTYPE)); mri = MRIallocSequence(width, height, depth, type, nframes); mri->dof = dof; - for (frame = start_frame; frame <= end_frame; frame++) { - for (z = 0; z < depth; z++) { - if ((int)znzread(buf, sizeof(char), bytes, fp) != bytes) { - // fclose(fp) ; - znzclose(fp); - free(buf); - ErrorReturn(NULL, (ERROR_BADFILE, "mghRead(%s): could not read %d bytes at slice %d", fname, bytes, z)); - } - switch (type) { - case MRI_INT: - for (i = y = 0; y < height; y++) { - for (x = 0; x < width; x++, i++) { - ival = orderIntBytes(((int *)buf)[i]); - MRIIseq_vox(mri, x, y, z, frame - start_frame) = ival; + + struct timespec begin, end; + if (getenv("FS_MGZIO_TIMING")) + { + printf("INFO: Environment variable FS_MGZIO_TIMING set\n"); + clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &begin); + } + + int USEVOXELBUF = 0; + if (mri->ischunked && getenv("FS_MGZIO_USEVOXELBUF")) + { + USEVOXELBUF = 1; + printf("INFO: Environment variable FS_MGZIO_USEVOXELBUF set\n"); + printf("mghRead(): read voxel buffer ...\n"); + if (buf) free(buf); + + int bytes_to_read = bytes * depth * (end_frame-start_frame+1); + //BUFTYPE *bufsrc = (BUFTYPE*)calloc(bytes_to_read, sizeof(BUFTYPE)); + BUFTYPE *bufsrc = &MRIseq_vox(mri, 0, 0, 0, start_frame); + if ((int)znzread(bufsrc, sizeof(BUFTYPE), bytes_to_read, fp) != bytes_to_read) { + znzclose(fp); + free(bufsrc); + ErrorReturn(NULL, (ERROR_BADFILE, "mghRead (read voxel buffer): could not read %d bytes", bytes_to_read)); + } + +#if (BYTE_ORDER == LITTLE_ENDIAN) + if (bpv == 2) byteswapbufshort(bufsrc, bytes_to_read); + if (bpv == 4) byteswapbuffloat(bufsrc, bytes_to_read); + if (bpv == 8) byteswapbuffloat(bufsrc, bytes_to_read); +#endif + + //void *voxelbuf = &MRIseq_vox(mri, 0, 0, 0, start_frame); + //memcpy(voxelbuf, bufsrc, bytes_to_read); + + //if (bufsrc) free(bufsrc); + } + else // copy voxel point by point + { + for (frame = start_frame; frame <= end_frame; frame++) { + for (z = 0; z < depth; z++) { + if ((int)znzread(buf, sizeof(char), bytes, fp) != bytes) { + // fclose(fp) ; + znzclose(fp); + free(buf); + ErrorReturn(NULL, (ERROR_BADFILE, "mghRead(%s): could not read %d bytes at slice %d", fname, bytes, z)); + } + switch (type) { + case MRI_INT: + for (i = y = 0; y < height; y++) { + for (x = 0; x < width; x++, i++) { + ival = orderIntBytes(((int *)buf)[i]); + MRIIseq_vox(mri, x, y, z, frame - start_frame) = ival; + } } - } - break; - case MRI_SHORT: - for (i = y = 0; y < height; y++) { - for (x = 0; x < width; x++, i++) { - sval = orderShortBytes(((short *)buf)[i]); - MRISseq_vox(mri, x, y, z, frame - start_frame) = sval; + break; + case MRI_SHORT: + for (i = y = 0; y < height; y++) { + for (x = 0; x < width; x++, i++) { + sval = orderShortBytes(((short *)buf)[i]); + MRISseq_vox(mri, x, y, z, frame - start_frame) = sval; + } } - } - break; - case MRI_TENSOR: - case MRI_FLOAT: - for (i = y = 0; y < height; y++) { - for (x = 0; x < width; x++, i++) { - fval = orderFloatBytes(((float *)buf)[i]); - MRIFseq_vox(mri, x, y, z, frame - start_frame) = fval; + break; + case MRI_USHRT: + for (i = y = 0; y < height; y++) { + for (x = 0; x < width; x++, i++) { + unsigned short usval = orderUShortBytes(((unsigned short *)buf)[i]); + MRIUSseq_vox(mri, x, y, z, frame - start_frame) = usval; + } } - } - break; - case MRI_UCHAR: - local_buffer_to_image(buf, mri, z, frame - start_frame); - break; - default: - errno = 0; - ErrorReturn(NULL, (ERROR_UNSUPPORTED, "mghRead: unsupported type %d", mri->type)); - break; - } - exec_progress_callback(z, depth, frame - start_frame, end_frame - start_frame + 1); - } + break; + case MRI_TENSOR: + case MRI_FLOAT: + for (i = y = 0; y < height; y++) { + for (x = 0; x < width; x++, i++) { + fval = orderFloatBytes(((float *)buf)[i]); + MRIFseq_vox(mri, x, y, z, frame - start_frame) = fval; + } + } + break; + case MRI_UCHAR: + local_buffer_to_image(buf, mri, z, frame - start_frame); + break; + default: + errno = 0; + ErrorReturn(NULL, (ERROR_UNSUPPORTED, "mghRead: unsupported type %d", mri->type)); + break; + } // switch + exec_progress_callback(z, depth, frame - start_frame, end_frame - start_frame + 1); + } // depth + } // frame + if (buf) free(buf); + } // end of copying voxel point by point + + if (getenv("FS_MGZIO_TIMING")) + { + clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &end); + printf("Total time (mghRead) = %ld.%09ld seconds%s\n", + (end.tv_nsec < begin.tv_nsec) ? (end.tv_sec - 1 - begin.tv_sec) : (end.tv_sec - begin.tv_sec), + (end.tv_nsec < begin.tv_nsec) ? (1000000000 + end.tv_nsec - begin.tv_nsec) : (end.tv_nsec - begin.tv_nsec), + (USEVOXELBUF) ? " (USEVOXELBUF)" : ""); } - if (buf) free(buf); } if (good_ras_flag > 0) { @@ -11035,49 +11132,104 @@ static int mghWrite(MRI *mri, const char *fname, int frame) memset(buf, 0, UNUSED_SPACE_SIZE * sizeof(char)); znzwrite(buf, sizeof(char), unused_space_size, fp); - for (frame = start_frame; frame <= end_frame; frame++) { - for (z = 0; z < depth; z++) { - for (y = 0; y < height; y++) { - switch (mri->type) { - case MRI_SHORT: - for (x = 0; x < width; x++) { - if (z == 74 && y == 16 && x == 53) DiagBreak(); - sval = MRISseq_vox(mri, x, y, z, frame); - znzwriteShort(sval, fp); - } - break; - case MRI_INT: - for (x = 0; x < width; x++) { - if (z == 74 && y == 16 && x == 53) DiagBreak(); - ival = MRIIseq_vox(mri, x, y, z, frame); - znzwriteInt(ival, fp); - } - break; - case MRI_FLOAT: - for (x = 0; x < width; x++) { - if (z == 74 && y == 16 && x == 53) DiagBreak(); - // printf("mghWrite: MRI_FLOAT: curr (x, y, z, frame) = (%d, %d, %d, %d)\n", x, y, z, frame); - fval = MRIFseq_vox(mri, x, y, z, frame); - // if(x==10 && y == 0 && z == 0 && frame == 67) - // printf("MRIIO: %g\n",fval); - znzwriteFloat(fval, fp); - } - break; - case MRI_UCHAR: - if ((int)znzwrite(&MRIseq_vox(mri, 0, y, z, frame), sizeof(BUFTYPE), width, fp) != width) { - errno = 0; - ErrorReturn(ERROR_BADFILE, (ERROR_BADFILE, "mghWrite: could not write %d bytes to %s", width, fname)); - } - break; - default: - errno = 0; - ErrorReturn(ERROR_UNSUPPORTED, (ERROR_UNSUPPORTED, "mghWrite: unsupported type %d", mri->type)); - break; - } - } - exec_progress_callback(z, depth, frame - start_frame, end_frame - start_frame + 1); + struct timespec begin, end; + if (getenv("FS_MGZIO_TIMING")) + { + printf("INFO: Environment variable FS_MGZIO_TIMING set\n"); + clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &begin); + } + + int USEVOXELBUF = 0; + if (mri->ischunked && getenv("FS_MGZIO_USEVOXELBUF")) + { + USEVOXELBUF = 1; + printf("INFO: Environment variable FS_MGZIO_USEVOXELBUF set\n"); + printf("mghWrite(): output voxel bufer ...\n"); + int bytes_per_voxel = MRIsizeof(mri->type); + BUFTYPE *voxelbuf = &MRIseq_vox(mri, 0, 0, 0, start_frame); + int bytes_to_write = bytes_per_voxel * width * height * depth * (end_frame-start_frame+1); + +#if 0 + printf(" voxelbuf = %p - %p (%p - %p) %d (frame=%d, depth=%d, height=%d) (%d %d)\n", + voxelbuf, voxelbuf+(MRIsizeof(mri->type)*(width-1)), + &MRISseq_vox(mri, 0, y, z, frame), &MRISseq_vox(mri, width-1, z, y, frame), + bytes_to_write, frame, y, z, (int)mri->bytes_per_vox, width); +#endif + +#if (BYTE_ORDER == LITTLE_ENDIAN) + if (bytes_per_voxel == 2) byteswapbufshort(voxelbuf, bytes_to_write); + if (bytes_per_voxel == 4) byteswapbuffloat(voxelbuf, bytes_to_write); + if (bytes_per_voxel == 8) byteswapbuffloat(voxelbuf, bytes_to_write); +#endif + + if ((int)znzwrite(voxelbuf, sizeof(BUFTYPE), bytes_to_write, fp) != bytes_to_write) + { + errno = 0; + ErrorReturn(ERROR_BADFILE, (ERROR_BADFILE, "mghWrite: (output voxel buffer) could not write %d bytes", bytes_to_write)); } } + else + { + for (frame = start_frame; frame <= end_frame; frame++) { + for (z = 0; z < depth; z++) { + for (y = 0; y < height; y++) { + switch (mri->type) { + case MRI_SHORT: + for (x = 0; x < width; x++) { + if (z == 74 && y == 16 && x == 53) DiagBreak(); + sval = MRISseq_vox(mri, x, y, z, frame); + znzwriteShort(sval, fp); + } + break; + case MRI_USHRT: + for (x = 0; x < width; x++) { + if (z == 74 && y == 16 && x == 53) DiagBreak(); + unsigned short usval = MRIUSseq_vox(mri, x, y, z, frame); + znzwriteUShort(usval, fp); + } + break; + case MRI_INT: + for (x = 0; x < width; x++) { + if (z == 74 && y == 16 && x == 53) DiagBreak(); + ival = MRIIseq_vox(mri, x, y, z, frame); + znzwriteInt(ival, fp); + } + break; + case MRI_FLOAT: + for (x = 0; x < width; x++) { + if (z == 74 && y == 16 && x == 53) DiagBreak(); + // printf("mghWrite: MRI_FLOAT: curr (x, y, z, frame) = (%d, %d, %d, %d)\n", x, y, z, frame); + fval = MRIFseq_vox(mri, x, y, z, frame); + // if(x==10 && y == 0 && z == 0 && frame == 67) + // printf("MRIIO: %g\n",fval); + znzwriteFloat(fval, fp); + } + break; + case MRI_UCHAR: + if ((int)znzwrite(&MRIseq_vox(mri, 0, y, z, frame), sizeof(BUFTYPE), width, fp) != width) { + errno = 0; + ErrorReturn(ERROR_BADFILE, (ERROR_BADFILE, "mghWrite: could not write %d bytes to %s", width, fname)); + } + break; + default: + errno = 0; + ErrorReturn(ERROR_UNSUPPORTED, (ERROR_UNSUPPORTED, "mghWrite: unsupported type %d", mri->type)); + break; + } // switch + } // height + exec_progress_callback(z, depth, frame - start_frame, end_frame - start_frame + 1); + } // depth + } // frame + } // end of writing voxel point by point + + if (getenv("FS_MGZIO_TIMING")) + { + clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &end); + printf("Total time (mghWrite) = %ld.%09ld seconds%s\n", + (end.tv_nsec < begin.tv_nsec) ? (end.tv_sec - 1 - begin.tv_sec) : (end.tv_sec - begin.tv_sec), + (end.tv_nsec < begin.tv_nsec) ? (1000000000 + end.tv_nsec - begin.tv_nsec) : (end.tv_nsec - begin.tv_nsec), + (USEVOXELBUF) ? " (USEVOXELBUF)" : ""); + } znzwriteFloat(mri->tr, fp); znzwriteFloat(mri->flip_angle, fp); @@ -11424,6 +11576,9 @@ MRI *MRIreorder(MRI *mri_src, MRI *mri_dst, int xdim, int ydim, int zdim) case MRI_SHORT: MRISseq_vox(mri_dst, xd, yd, zd, f) = MRISseq_vox(mri_src, xs, ys, zs, f); break; + case MRI_USHRT: + MRIUSseq_vox(mri_dst, xd, yd, zd, f) = MRIUSseq_vox(mri_src, xs, ys, zs, f); + break; case MRI_FLOAT: MRIFseq_vox(mri_dst, xd, yd, zd, f) = MRIFseq_vox(mri_src, xs, ys, zs, f); break;