diff --git a/+imutils/get_nii_coordinates.m b/+imutils/get_nii_coordinates.m
new file mode 100644
index 00000000..4028f246
--- /dev/null
+++ b/+imutils/get_nii_coordinates.m
@@ -0,0 +1,48 @@
+function [x, y, z] = get_nii_coordinates( nii )
+%GET_NII_COORDINATES Return voxel position "world-coordinates" in mm
+%     
+%     [x, y, z] = get_nii_coordinates( niiFile )
+%     [x, y, z] = get_nii_coordinates( niiInfo )
+% 
+% Returns `x, y, z`: three 3-D arrays of voxel positions in mm.
+%
+% The single input can be given either as:
+%
+% 1. `niiFile`—A path string to the NIfTI file of interest; or, 
+%
+% 2. `niiInfo`—A struct of the form `niiInfo = niftiinfo( niiFile )`
+%
+% __NOTE__
+% The function is implemented as
+% ```
+% [i, j, k] = ndgrid( [0:niiInfo.ImageSize(1)-1], [0:niiInfo.ImageSize(2)-1], [0:niiInfo.ImageSize(3)-1] ) ;
+% [x, y, z] = niiInfo.Transform.transformPointsForward( i, j, k ) ;  
+% ```
+% __TODO__
+% * function should be further tested.
+
+% -----------
+%% Check input
+narginchk(1,1);
+
+if [ isstring( nii ) | ischar( nii ) ] & isfile( nii )
+    
+    [~, info] = imutils.read_nii( nii );
+
+elseif isstruct( nii ) && isfield( nii, 'ImageSize' ) && ...
+        isfield( nii, 'Transform' ) && isa( nii.Transform, 'affine3d' ) 
+
+    info = nii ;
+
+else
+    error('Input must be a path to a .nii file, or a struct returned by `niftiinfo`')
+end
+
+% ---------------
+%% Get coordinates
+
+% voxel indices
+[i, j, k] = ndgrid( [0:info.ImageSize(1)-1], [0:info.ImageSize(2)-1], [0:info.ImageSize(3)-1] ) ;
+[x, y, z] = info.Transform.transformPointsForward( i, j, k ) ;
+
+end
diff --git a/+imutils/read_nii.m b/+imutils/read_nii.m
new file mode 100644
index 00000000..5445f12c
--- /dev/null
+++ b/+imutils/read_nii.m
@@ -0,0 +1,194 @@
+function [img,info,json] = read_nii( niiFile, options )
+%READ_NII Load NIfTI image, header, and (if present) the accompanying .json sidecar
+%     
+%     [img, info]       = read_nii( niiFile )
+%     [img, info, json] = read_nii( niiFile, options )
+% 
+% The input `niiFile` is the path to the NIfTI image as a string scalar or
+% character vector. When called with 2 output arguments, the function is
+% equivalent short-hand for
+%     
+%     info = niftiinfo( niiFile ); img = niftiread( info ) 
+%
+% The function checks the `niiFile` parent-folder for the presence of a
+% sidecar (an identically named file, but with a .json file extension).
+% When such a file exists, the 3rd output is returned as a struct via 
+% `json = jsondecode( fileread( jsonFile ) );` otherwise, `json = []`.
+% 
+% __OPTIONS__
+% 
+% The function accepts an `options` struct of parameters (for now, only one) 
+% as a 2nd argument for which the `.rescale` field can be assigned:
+%
+% | `options.rescale` |  Effect                                                |
+% | ------------------| -------------------------------------------------------|
+% | 'off'             | Rescaling disabled                                     |
+% | 'basic'           | Rescale according to nii header info                   | 
+% | 'auto' [default]  | Rescale and convert to physical units when possible    |
+%
+% __NOTE__  
+% For now, the sole effect of `'auto'` is to convert Siemens raw phase images
+% to physical units (radians), which requires converting from their original
+% integer type (between [0,4095]) to a 32-bit "single" float (between [-pi,pi)). 
+% The json sidecar must be available to verify the Manufacturer and ImageType
+% entries. Otherwise, and for all other image inputs `'auto'` reverts to
+% `'basic'`.
+
+DEFAULTS.rescale = 'auto';
+
+narginchk(1, 2);
+
+if nargin == 1 || ~isfield( options, 'rescale' ) || isempty( options.rescale )
+    options.rescale = DEFAULTS.rescale;
+else
+    assert( isstruct(options), 'Optional 2nd input must be a struct' );
+    assert( ismember( options.rescale, ["off" "basic" "auto"] ), ...
+        'Invalid assignment to `options.rescale`: Value must be "off", "basic", or "auto"]' );
+end
+
+info = niftiinfo( niiFile ); 
+img  = niftiread( info ); 
+
+% `extractBefore` should get the correct filename in both .nii and .nii.gz cases
+jsonFile = strcat( extractBefore( niiFile, '.nii' ), '.json' ); 
+
+if isfile( jsonFile )
+    if isOctave()
+      json = loadjson( jsonFile );
+    else % matlab
+      json = jsondecode( fileread( jsonFile ) );  
+    end
+else
+    json = [];
+end
+
+% -------------------
+%% Optional rescaling
+if strcmp(options.rescale, 'off')
+    return
+end
+
+if strcmp(options.rescale, 'basic')
+    [img, info] = rescale( img, info );
+    return
+end
+
+% NOTE: Other approaches to rescaling and/or converting from raw file values
+% could be added (including cases where the json sidecar is unavailable)
+if isfield( json, 'Manufacturer') && strcmp(json.Manufacturer, 'Siemens') ...
+  && strcmp( check_json_image_type(json), 'phase' )
+
+    [img, info] = convert_siemens_phase( img, info );
+else 
+    [img, info] = rescale( img, info );
+end
+
+end
+
+% -----------------------------------------------------------------------------
+%% Local functions
+% -----------------------------------------------------------------------------
+function [img, info] = convert_siemens_phase( img0, info0 )
+%CONVERT_SIEMENS_PHASE Return phase img in rad 
+%    
+%    [img, info] = convert_siemens_phase( img0, info0 )
+%
+% Converts from integer-type (between [0,4095]) to a 32-bit "single" float
+% (between [-pi,pi)). 
+
+PHASE_SCALING_SIEMENS = 4096;
+
+if (info0.AdditiveOffset == -PHASE_SCALING_SIEMENS) && (info0.MultiplicativeScaling == 2) 
+
+    [img, info] = rescale( img0, info0 );
+
+    img = single(img)*(pi/PHASE_SCALING_SIEMENS);
+
+    %% Update header: both Matlab simplified and original `raw` portions
+    info.Datatype     = 'single';
+    info.BitsPerPixel = 32;
+    % NB: 16 is the NIfTI code for float; bitpix = number of bits
+    info.raw.datatype = 16;
+    info.raw.bitpix   = 32;
+else
+    warning( 'The nii header differs from that expected of Siemens phase data.\n%s', ...
+              'Output values (units) are effectively unknown' );
+    [img, info] = rescale( img0, info0 );
+end
+
+end
+
+% -----------------------------------------------------------------------------
+function [imgType] = check_json_image_type( json )
+%CHECK_JSON_IMAGE_TYPE Check json; return 'magnitude', 'phase' or 'unknown'
+%    
+%    imgType = check_json_image_type( json )
+%
+% Checks the `ImageType` field of struct `json` (derived from the json sidecar) 
+% and returns a char vector `imgType` which may be:  
+%
+% - `magnitude` for magnitude images
+% - `phase` for phase images
+% - `unknown` otherwise 
+%
+% __NOTE__
+% The check is hardly exhaustive and has only been tested for Siemens MRI data.
+assert( nargin<2 )
+
+if isempty(json)
+    imgType = []; 
+    return
+end
+    
+assert( isstruct(json) );
+
+imgType = 'unknown'; % default
+
+if myisfield( json, 'ImageType' )
+
+    isPhase = any( ismember(json.ImageType, "P") );
+    isMag   = any( ismember(json.ImageType, "M") ); 
+
+    if isPhase && isMag 
+        % Both true: json file and/or DICOM issue (hopefully this doesn't occur?) 
+        warning('Ambiguous ImageType entry in json file: Indicates magnitude AND phase?');
+    elseif isPhase
+        imgType = 'phase';
+    elseif isMag
+        imgType = 'magnitude';
+    end
+
+    if strcmp( imgType, 'unknown' ) && ~isfield( json, 'Manufacturer' ) || ~strcmp( json.Manufacturer, 'Siemens' )
+        warning('Unknown image type. Possibly due to images sourced from non-Siemens MRI')
+    end
+end
+
+end
+% -----------------------------------------------------------------------------
+function [img, info] = rescale( img0, info0 )
+%RESCALE Rescale image according to NIfTI header info
+%      
+%     [img, info] = rescale( img, info )
+    
+    img1 = info0.AdditiveOffset + info0.MultiplicativeScaling .* img0;
+   
+    %% Check for possible integer overflow/type issues
+    if isequaln( img1, double(info0.AdditiveOffset) + double(info0.MultiplicativeScaling) * double(img0) )
+
+        img  = img1;
+        info = info0;
+        
+        %% Update header: both Matlab simplified and original `raw` portions
+        info.MultiplicativeScaling = 1;
+        info.AdditiveOffset        = 0;
+        info.raw.scl_slope         = 1;
+        info.raw.scl_inter         = 0;
+
+    else
+        warning( ['Aborting default image rescaling to avoid integer overflow.\n'
+                  'Original NIfTI header may contain errors.\n '] ) 
+    
+        img  = img0;
+        info = info0;
+    end
+end
diff --git a/Img/read_nii.m b/Img/read_nii.m
new file mode 100644
index 00000000..1926567d
--- /dev/null
+++ b/Img/read_nii.m
@@ -0,0 +1,37 @@
+function [img,info,json] = read_nii( niiFile )
+%READ_NII Load NIfTI image and header; and, when present, the accompanying .json sidecar
+%     
+%     [img, info]       = read_nii( niiFile )
+%     [img, info, json] = read_nii( niiFile )
+% 
+% The input `niiFile` is the path to the NIfTI image as a string scalar or
+% character vector. When called with 2 output arguments, the function is
+% equivalent short-hand for
+%     
+%     info = niftiinfo( niiFile ); img = niftiread( info ) 
+%
+% When called with the 3rd output argument, the function checks the parent
+% folder of `niiFile` for an identically named file but with a .json file extension.
+% When such a file is present, the 3rd output is returned as a struct via 
+% `json = jsondecode( fileread( jsonFile ) );` otherwise, `json = []`.
+arguments
+    niiFile(1,:) string {mustBeStringScalarOrCharVector, mustBeFile} ;
+end
+
+info = niftiinfo( niiFile ); 
+img  = niftiread( info ); 
+
+if nargout < 3
+    return;
+end
+
+[folder, name] = fileparts( niiFile );
+jsonFile       = fullfile( folder, name + ".json" ); 
+
+if isfile( jsonFile )
+    json = jsondecode( fileread( jsonFile ) );     
+else
+    json = [];
+end
+
+end
diff --git a/misc/isOctave.m b/misc/isOctave.m
new file mode 100644
index 00000000..5ec635de
--- /dev/null
+++ b/misc/isOctave.m
@@ -0,0 +1,8 @@
+function r = isOctave ()
+  % https://stackoverflow.com/a/9838357
+  persistent x;
+  if (isempty (x))
+    x = exist ('OCTAVE_VERSION', 'builtin');
+  end
+  r = x;
+end