Now using N-1 instead of N as denominator for computing the STD. (#1441)

* OPT: maths: added option for setting output type

* OPT: maths: set ddof=1 for STD calculation (issue #1440)

* OPT: fmri_compute_tsnr: now using ddof=1 for computing STD (issue #1440)

* REF: fmri_compute_tsnr: cleaned up old comments

scripts/sct_fmri_compute_tsnr.py

@@ -17,10 +17,8 @@
 import sct_maths
 import sct_utils as sct
 from msct_parser import Parser
-
-
-# from sct_average_data_across_dimension import average_data_across_dimension
-
+from msct_image import Image
+import numpy as np
 
 class Param:
  def __init__(self):
@@ -41,51 +39,23 @@ def compute(self):
 
  fname_data = self.fmri
 
- # # create temporary folder
- # sct.printv('\nCreate temporary folder...', self.param.verbose)
- # path_tmp = 'tmp.'+time.strftime("%y%m%d%H%M%S/")
- # status, output = sct.run('mkdir '+path_tmp, self.param.verbose)
-
- # # motion correct the fmri data
- # # sct.printv('\nMotion correct the fMRI data...', self.param.verbose, 'normal')
- # path_fmri, fname_fmri, ext_fmri = sct.extract_fname(self.fmri)
- # fname_fmri_moco = fname_fmri
- # # print sct.slash_at_the_end(path_fmri) + fname_fmri
- # # sct.run('mcflirt -in ' + sct.slash_at_the_end(path_fmri, 1) + fname_fmri + ' -out ' + fname_fmri_moco)
+ # open data
+ nii_data = Image(fname_data)
+ data = nii_data.data
 
  # compute mean
- fname_data_mean = sct.add_suffix(fname_data, '_mean')
- sct_maths.main(args=[
- '-i', fname_data,
- '-o', fname_data_mean,
- '-mean', 't'
- ])
-
+ data_mean = np.mean(data, 3)
  # compute STD
- fname_data_std = sct.add_suffix(fname_data, '_std')
- sct_maths.main(args=[
- '-i', fname_data,
- '-o', fname_data_std,
- '-std', 't'
- ])
-
- # compute tSNR
- fname_tsnr = sct.add_suffix(fname_data, '_tsnr')
- from msct_image import Image
- nii_mean = Image(fname_data_mean)
- data_mean = nii_mean.data
- data_std = Image(fname_data_std).data
+ data_std = np.std(data, 3, ddof=1)
+ # compute TSNR
  data_tsnr = data_mean / data_std
- nii_tsnr = nii_mean
+
+ # save TSNR
+ fname_tsnr = sct.add_suffix(fname_data, '_tsnr')
+ nii_tsnr = nii_data
  nii_tsnr.data = data_tsnr
  nii_tsnr.setFileName(fname_tsnr)
- nii_tsnr.save()
-
- # Remove temp files
- sct.printv('\nRemove temporary files...', self.param.verbose, 'normal')
- import os
- os.remove(fname_data_mean)
- os.remove(fname_data_std)
+ nii_tsnr.save(type='float32')
 
  # to view results
  sct.printv('\nDone! To view results, type:', self.param.verbose, 'normal')

scripts/sct_maths.py

@@ -164,6 +164,11 @@ def get_parser():
  description='Symmetrize data along the specified dimension.',
  mandatory=False,
  example=['0', '1', '2'])
+ parser.add_option(name='-type',
+ type_value='multiple_choice',
+ description='Output type.',
+ mandatory=False,
+ example=['uint8', 'int16', 'int32', 'float32', 'complex64', 'float64', 'int8', 'uint16', 'uint32', 'int64', 'uint64'])
  parser.add_option(name="-v",
  type_value="multiple_choice",
  description="""Verbose. 0: nothing. 1: basic. 2: extended.""",
@@ -188,6 +193,10 @@ def main(args = None):
  fname_in = arguments["-i"]
  fname_out = arguments["-o"]
  verbose = int(arguments['-v'])
+ if '-type' in arguments:
+ output_type = arguments['-type']
+ else:
+ output_type = ''
 
  # Open file(s)
  im = Image(fname_in)
@@ -270,7 +279,7 @@ def main(args = None):
  dim = dim_list.index(arguments['-std'])
  if dim + 1 > len(np.shape(data)): # in case input volume is 3d and dim=t
  data = data[..., np.newaxis]
- data_out = std(data, dim)
+ data_out = std(data, dim, ddof=1)
 
  elif "-smooth" in arguments:
  sigmas = arguments["-smooth"]
@@ -332,7 +341,7 @@ def main(args = None):
  nii_out = Image(fname_in) # use header of input file
  nii_out.data = data_out
  nii_out.setFileName(fname_out)
- nii_out.save()
+ nii_out.save(type=output_type)
  # TODO: case of multiple outputs
  # assert len(data_out) == n_out
  # if n_in == n_out: