new type MRI_USHRT, env variables FS_MGZIO_USEVOXELBUF & FS_MGZIO_TIM…

…ING (freesurfer#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 <[email protected]>
ayendiki · Dec 6, 2021 · b82c3a8 · b82c3a8
1 parent a0c821d
commit b82c3a8
Show file tree

Hide file tree

Showing 29 changed files with 940 additions and 128 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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,