Nifti and compressed nifti, orientation support, load and save

External/mim/Gui/Controllers/MimChooseImagingFiles.m

@@ -27,7 +27,8 @@
         '*.xif', 'HDllab/ATL Ultrasound';
         '*.vtk', 'Visualization Toolkit (VTK)';
         '*.vff', 'MicroCT';
-        '*.par;*.rec', 'Philips PAR/REC'
+        '*.par;*.rec', 'Philips PAR/REC';
+        '*.nii.gz', 'Compressed NIFTI [beta]'
         };
 
     [image_path, filenames, filter_index] = CoreDiskUtilities.ChooseFiles('Select the file to import', image_path, true, filespec);
@@ -72,7 +73,7 @@
         image_info = PTKImageInfo(image_path, {filenames{1}}, image_type, [], [], []);
         return;
 
-    elseif (filter_index == 7)
+    elseif (filter_index == 7) || (filter_index == 14)
         image_type = MimImageFileFormat.Nifti;
         image_info = PTKImageInfo(image_path, {filenames{1}}, image_type, [], [], []);
         return;

External/mim/Gui/Controllers/MimSaveAs.m

@@ -63,6 +63,7 @@
                 '*.nii', 'NIFTI (*.nii)';
                 '*.mhd', '8-bit metaheader and raw data (*.mhd)';
                 '*.mhd', '16-bit metaheader and raw data (*.mhd)';
+                '*.nii.gz', 'Compressed NIFTI [beta] (*.nii.gz)';
                 };
 
     if path_name == 0
@@ -83,7 +84,7 @@ function SaveImage(image_data, filename, pathname, filter_index, patient_name, i
                 MimSaveImageAsDicom(image_data, pathname, filename, patient_name, is_secondary_capture, dicom_metadata, reporting)
             case 2
                 MimSaveAsMatlab(image_data, pathname, filename, reporting);
-            case 3
+            case {3,6}
                 MimSaveAsNifti(image_data, pathname, filename, reporting);
             case 4
                 MimSaveAsMetaheaderAndRaw(image_data, pathname, filename, 'char', reporting)
@@ -92,7 +93,7 @@ function SaveImage(image_data, filename, pathname, filter_index, patient_name, i
                     MimSaveAsMetaheaderAndRaw(image_data, pathname, filename, 'short', reporting)
                 else
                     MimSaveAsMetaheaderAndRaw(image_data, pathname, filename, 'ushort', reporting)
-                end
+                end                
         end
     end
 end

External/mim/Library/File/MimGuessFileType.m

@@ -37,6 +37,11 @@
         principal_filename = {image_filename};
         secondary_filenames = {};
         return;
+    elseif strcmp(ext, '.gz')
+        image_type = MimImageFileFormat.Nifti;
+        principal_filename = {image_filename};
+        secondary_filenames = {};
+        return;
 
     elseif strcmp(ext, '.img')
         hdr_filename = [name '.hdr'];

External/mim/Library/File/MimLoadOtherFormat.m

@@ -61,10 +61,32 @@
     switch image_file_format
 
         case MimImageFileFormat.Nifti
-            header_data = nii_read_header(header_filename);
+            [~,~,ext] = fileparts(header_filename);
+            % if compressed (*.nii.gz) extract to read
+            if strcmp(ext,'.gz')
+                compressed = 1;
+                if isfile(header_filename(1:end-3))
+                    reporting.Error('MimLoadNiiGz:extractedniiexists', ...
+                        ['Non-compressed *.nii of ', header_filename,...
+                        ' already exists at path. Move/delete and try again.']);
+                end
+                gunzip(header_filename)
+                readfile = header_filename(1:end-3);
+            else
+                readfile = header_filename(1:end-3);
+                compressed = 0;
+            end
+
+            header_data = nii_read_header(readfile);
             data = nii_read_volume(header_data);
+            if compressed == 1
+                % remove temp decompressed file
+                delete(readfile);
+                % change header data to point back to original file
+                header_data.Filename = header_data.Filename;
+            end
             [new_dimension_order, flip_orientation] = MimImageCoordinateUtilities.GetDimensionPermutationVectorFromNiiOrientation(header_data, reporting);
-        
+
         case MimImageFileFormat.Analyze % Experimental: assumes fixed orientation
             header_data = hdr_read_header(header_filename);
             data = hdr_read_volume(header_data);

External/mim/Library/File/MimSaveAsNifti.m

@@ -27,25 +27,62 @@ function MimSaveAsNifti(image_to_save, path, filename, reporting)
     if ~isa(image_to_save, 'PTKImage')
         reporting.Error('MimSaveAsNifti:InputMustBePTKImage', 'Requires a PTKImage as input');
     end
+    full_filename = fullfile(path, filename);
 
+    % check if .nii.gz given
+    [~,~,ext] = fileparts(filename);
+    if strcmp(ext,'.gz')
+        compressed = 1;
+        if isfile(fullfile(path, filename(1:end-3)))
+            reporting.Error('MimSaveAsNifti:rawniiexists', ...
+                        ['Non-compressed *.nii of ', full_filename,...
+                        ' already exists at path. Move/delete and try again.']);
+        end
+        full_filename = fullfile(path, filename(1:end-3));
+    else
+        compressed = 0;
+    end
+
     image_data = image_to_save.RawImage;
 
-    full_filename = fullfile(path, filename);
-
     resolution = image_to_save.VoxelSize([2, 1, 3]);
 
     offset = [0 0 0];
-    if isa(image_data, 'PTKDicomImage')
-        metadata = image_data. MetaHeader;
-        if isfield(metadata, 'Offset')
-            offset = metadata.Offset;
+    if isa(image_to_save, 'PTKDicomImage')
+        metadata = image_to_save. MetaHeader;
+        if isfield(metadata, 'ImagePositionPatient')
+            offset = metadata.ImagePositionPatient;
         end
     end
 
     image_data = permute(image_data, [2, 1, 3]);
-    image_data = flip(flip(flip(image_data, 3), 2), 1);
+    image_data = flip(image_data, 3); % only flip last dimension
+
 
     nii_data = make_nii(image_data, resolution, offset, [], image_to_save.Title);
     save_nii(nii_data, full_filename);
+
+    % set orientation
+    if isa(image_to_save, 'PTKDicomImage')
+        % get LPS form from dicom
+        if isfield(metadata, 'ImagePositionPatient')
+            dicomorientation = metadata.ImageOrientationPatient;
+            d1 = dicomorientation(1:3);
+            d2 = dicomorientation(4:6);
+            d3 = cross(d1,d2);
+            LPS_affine4 = [d1, d2, d3, offset];
+            LPS_affine4 = [LPS_affine4; [0,0,0,1]];
+            % convert to RAS for nifti
+            RAS_affine4 = diag([-1,-1,1,1])*LPS_affine4;
+            WriteNiftiOrientation(full_filename, RAS_affine4)
+        end
+    end
+
+    % if .nii.gz save as .nii.gz
+    if compressed == 1
+        gzip(full_filename)
+        delete(full_filename)
+    end
+
 end
 

External/mim/Library/File/WriteNiftiOrientation.m

@@ -0,0 +1,71 @@
+function WriteNiftiOrientation(inputfile, aff)
+% By Ashkan Pakzad May 2021 (ashkanpakzad.github.io) 
+% Copyright Ashkan Pakzad 2021. Distributed under MIT licence.
+
+% Give name of nifti file ending either .nii or .nii.gz to have orientation
+% written from given affine
+
+% aff is a 4x4 matrix that represents the orientation of the given image in
+% inputfile. As a nifti affine it maps the data coordinate directions ijk 
+% to image coordinates in RAS. The last row of aff must be [0 0 0 1].
+% The input aff matrix should not consider image spacing the top left 3x3
+% matrix should only be of [-1,0,1]. The image spacing saved in the
+% original image is written. 
+
+% For more information regarding medical image orientation consider
+% https://medium.com/@ashkanpakzad/understanding-3d-medical-image-orientation-for-programmers-fcf79c7beed0
+
+% decompress and load
+[~,~,ext] = fileparts(inputfile);
+if strcmp(ext,'.gz')
+    filename = string(gunzip(inputfile));
+else
+    filename = inputfile;
+end
+
+% Find Voxel dimensions; format = numdims,x,y,z,t
+vox_dims = fieldread(76, 8, 'float');
+vox_dims = vox_dims(2:4);
+
+% write given affine into file
+% only use sform, set qform to 0 and sform to 1
+fieldwrite(0, 252, 'short'); % q
+fieldwrite(1, 254, 'short'); % s
+
+% add spacing info to input affine
+aff(:,1) = aff(:,1)*vox_dims(1);
+aff(:,2) = aff(:,2)*vox_dims(2);
+aff(:,3) = aff(:,3)*vox_dims(3);
+
+Sx = aff(1,:);
+Sy = aff(2,:);
+Sz = aff(3,:);
+
+fieldwrite(Sx, 280, 'float')
+fieldwrite(Sy, 296, 'float')
+fieldwrite(Sz, 312, 'float')
+
+% re-gzip if originally gzipped
+if strcmp(ext,'.gz')
+    gzip(filename);
+    delete(filename)
+end
+
+% read and write functions
+
+function field = fieldread(offset, size, precision)
+    % open file and read binary field from given offset into size and type.
+    fid=fopen(filename);
+    fseek(fid,offset,'bof');
+    field = fread(fid,size,precision);
+    fclose(fid);
+end
+
+function fieldwrite(field, offset, precision)
+    % open file and read binary field from given offset into size and type.
+    fid=fopen(filename, 'r+');
+    fseek(fid,offset,'bof');
+    fwrite(fid,field,precision);
+    fclose(fid);
+end
+end