docs: def trimim( fims, 
              affine=None,
              scale=2,
              divdim = 8**2,
              int_order=0,
              fmax = 0.05,
              outpath='',
              fname='',
              fcomment='',
              store_avg=False,
              store_img_intrmd=False,
              store_img=False,
              imdtype=np.float32,
              memlim=False,
              verbose=False):
      '''
      Trim and upsample PET image(s), e.g., for GPU execution,
      PVC correction, ROI sampling, etc.
      The input images 'fims' can be passed in multiple ways:
      1. as a string of the folder containing NIfTI files
      2. as a string of a NIfTI file path (this way a 4D image can be loaded).
      3. as a list of NIfTI file paths.
      4. as a 3D or 4D image
      Parameters:
      -----------
      affine: affine matrix, 4x4, for all the input images in case when images
              passed as a matrix (all have to have the same shape and data type)
      scale: the scaling factor for upsampling has to be greater than 1
      divdim: image divisor, i.e., the dimensions are a multiple of this number (default 64)
      int_order: interpolation order (0-nearest neighbour, 1-linear, as in scipy)
      fmax: fraction of the max image value used for finding image borders for trimming
      outpath: output folder path
      fname: file name when image given as a numpy matrix
      fcomment: part of the name of the output file, left for the user as a comment
      store_img_intrmd: stores intermediate images with suffix '_i'
      store_avg: stores the average image (if multiple images are given)
      imdtype: data type for output images
      memlim: Ture for cases when memory is limited and takes more processing time instead.
      verbose: verbose mode [True/False]
      '''
  
      using_multiple_files = False
  
      # case when input folder is given
      if isinstance(fims, basestring) and os.path.isdir(fims):
          # list of input images (e.g., PET)
          fimlist = [os.path.join(fims,f) for f in os.listdir(fims) if f.endswith('.nii') or f.endswith('.nii.gz')]
          imdic = imio.niisort(fimlist, memlim=memlim)
          if not (imdic['N']>50 and memlim):
              imin = imdic['im']
          imshape = imdic['shape']
          affine = imdic['affine']
          fldrin = fims
          fnms = [os.path.basename(f).split('.nii')[0] for f in imdic['files'] if f!=None]
          # number of images/frames
          Nim = imdic['N']
          using_multiple_files = True
  
      # case when input file is a 3D or 4D NIfTI image
      elif isinstance(fims, basestring) and os.path.isfile(fims) and (fims.endswith('nii') or fims.endswith('nii.gz')):
          imdic = imio.getnii(fims, output='all')
          imin = imdic['im']
          if imin.ndim==3:
              imin.shape = (1, imin.shape[0], imin.shape[1], imin.shape[2])
          imdtype = imdic['dtype']
          imshape = imdic['shape'][-3:]
          affine = imdic['affine']
          fldrin = os.path.dirname(fims)
          fnms = imin.shape[0] * [ os.path.basename(fims).split('.nii')[0] ]
          # number of images/frames
          Nim = imin.shape[0]
  
      # case when a list of input files is given
      elif isinstance(fims, list) and all([os.path.isfile(k) for k in fims]):
          imdic = imio.niisort(fims, memlim=memlim)
          if not (imdic['N']>50 and memlim):
              imin = imdic['im']
          imshape = imdic['shape']
          affine = imdic['affine']
          imdtype = imdic['dtype']
          fldrin = os.path.dirname(fims[0])
          fnms = [os.path.basename(f).split('.nii')[0] for f in imdic['files']]
          # number of images/frames
          Nim = imdic['N']
          using_multiple_files = True
  
      # case when an array [#frames, zdim, ydim, xdim].  Can be 3D or 4D
      elif isinstance(fims, (np.ndarray, np.generic)) and (fims.ndim==4 or fims.ndim==3):
          # check image affine
          if affine.shape!=(4,4):
              raise ValueError('Affine should be a 4x4 array.')
          # create a copy to avoid mutation when only one image (3D)
          imin = np.copy(fims)
          if fims.ndim==3:
              imin.shape = (1, imin.shape[0], imin.shape[1], imin.shape[2])
          imshape = imin.shape[-3:]
          fldrin = os.path.join(os.path.expanduser('~'), 'NIMPA_output')
          if fname=='':
              fnms = imin.shape[0] * [ 'NIMPA' ]
          else:
              fnms = imin.shape[0] * [ fname ]
          # number of images/frames
          Nim = imin.shape[0]
  
      else:
          raise TypeError('Wrong data type input.')
  
      
  
      #-------------------------------------------------------
      # store images in this folder
      if outpath=='':
          petudir = os.path.join(fldrin, 'trimmed')
      else:
          petudir = os.path.join(outpath, 'trimmed')
      imio.create_dir( petudir )
      #-------------------------------------------------------
  
      #-------------------------------------------------------
      # scale is preferred to be integer
      try:
          scale = int(scale)
      except ValueError:
          raise ValueError('e> scale has to be an integer.')
      # scale factor as the inverse of scale
      sf = 1/float(scale)
      if verbose:
          print 'i> upsampling scale {}, giving resolution scale factor {} for {} images.'.format(scale, sf, Nim)
      #-------------------------------------------------------
  
      #-------------------------------------------------------
      # scaled input image and get a sum image as the base for trimming
      if scale>1:
          newshape = (scale*imshape[0], scale*imshape[1], scale*imshape[2])
          imsum = np.zeros(newshape, dtype=imdtype)
          if not memlim:
              imscl = np.zeros((Nim,)+newshape, dtype=imdtype)
              for i in range(Nim):
                  imscl[i,:,:,:] = ndi.interpolation.zoom(imin[i,:,:,:], (scale, scale, scale), order=int_order )
                  imsum += imscl[i,:,:,:]
          else:
              for i in range(Nim):
                  if Nim>50 and using_multiple_files:
                      imin_temp = imio.getnii(imdic['files'][i])
                      imsum += ndi.interpolation.zoom(imin_temp, (scale, scale, scale), order=0 )
                      if verbose:
                          print 'i> image sum: read', imdic['files'][i]
                  else:
                      imsum += ndi.interpolation.zoom(imin[i,:,:,:], (scale, scale, scale), order=0 )
      else:
          imscl = imin
          imsum = np.sum(imin, axis=0)
  
      # smooth the sum image for improving trimming (if any)
      #imsum = ndi.filters.gaussian_filter(imsum, imio.fwhm2sig(4.0, voxsize=abs(affine[0,0])), mode='mirror')
      #-------------------------------------------------------
     
      # find the object bounding indexes in x, y and z axes, e.g., ix0-ix1 for the x axis
      qx = np.sum(imsum, axis=(0,1))
      ix0 = np.argmax( qx>(fmax*np.nanmax(qx)) )
      ix1 = ix0+np.argmin( qx[ix0:]>(fmax*np.nanmax(qx)) )
  
      qy = np.sum(imsum, axis=(0,2))
      iy0 = np.argmax( qy>(fmax*np.nanmax(qy)) )
      iy1 = iy0+np.argmin( qy[iy0:]>(fmax*np.nanmax(qy)) )
  
      qz = np.sum(imsum, axis=(1,2))
      iz0 = np.argmax( qz>(fmax*np.nanmax(qz)) )
  
      # find the maximum voxel range for x and y axes
      IX = ix1-ix0+1
      IY = iy1-iy0+1
      tmp = max(IX, IY)
      #> get the range such that it is divisible by 
      #> divdim (64 by default) for GPU execution
      IXY = divdim * ((tmp+divdim-1)/divdim)
      div = (IXY-IX)/2
      # x
      ix0 -= div
      ix1 += (IXY-IX)-div
      # y
      div = (IXY-IY)/2
      iy0 -= div
      iy1 += (IXY-IY)-div
      # z
      tmp = (len(qz)-iz0+1)
      IZ = divdim * ((tmp+divdim-1)/divdim)
      iz0 -= IZ-tmp+1
      
      # save the trimming parameters in a dic
      trimpar = {'x':(ix0, ix1), 'y':(iy0, iy1), 'z':(iz0), 'fmax':fmax}
  
      # new dims (z,y,x)
      newdims = (imsum.shape[0]-iz0, iy1-iy0+1, ix1-ix0+1)
      imtrim = np.zeros((Nim,)+newdims, dtype=imdtype)
      imsumt = np.zeros(newdims, dtype=imdtype)
      # in case of needed padding (negative indx resulting above)
      # the absolute values are supposed to work like padding in case the indx are negative
      iz0s, iy0s, ix0s = iz0, iy0, ix0
      iz0t, iy0t, ix0t = 0,0,0
      if iz0<0: 
          iz0s=0; 
          iz0t = abs(iz0)
          print '-----------------------------------------------------------------'
          print 'w> Correcting for trimming outside the original image (z-axis)'
          print '-----------------------------------------------------------------'
      if iy0<0: 
          iy0s=0; 
          iy0t = abs(iy0)
          print '-----------------------------------------------------------------'
          print 'w> Correcting for trimming outside the original image (y-axis)'
          print '-----------------------------------------------------------------'
      if ix0<0: 
          ix0s=0; 
          ix0t = abs(ix0)
          print '-----------------------------------------------------------------'
          print 'w> Correcting for trimming outside the original image (x-axis)'
          print '-----------------------------------------------------------------'
      
  
      #> in case the upper index goes beyond the scaled but untrimmed image
      iy1t = imsumt.shape[1]
      if iy1>=imsum.shape[1]:
          iy1t -= iy1+1
  
      #> the same for x
      ix1t = imsumt.shape[2]
      if ix1>=imsum.shape[2]:
          ix1t -= ix1+1
  
  
      # first trim the sum image
      imsumt[iz0t:, iy0t:iy1t, ix0t:ix1t] = imsum[iz0s:, iy0s:iy1+1, ix0s:ix1+1]
  
      #> new affine matrix for the scaled and trimmed image
      A = np.diag(sf*np.diag(affine))
      #> note half of new voxel offset is used for the new centre of voxels
      A[0,3] = affine[0,3] + A[0,0]*(ix0-0.5)        
      A[1,3] = affine[1,3] + (affine[1,1]*(imshape[1]-1) - A[1,1]*(iy1-0.5))
      A[2,3] = affine[2,3] - A[1,1]*0.5
      A[3,3] = 1
  
  
      # output dictionary
      dctout = {'affine': A, 'trimpar':trimpar, 'imsum':imsumt}
  
      # NIfTI image description (to be stored in the header)
      niidescr = 'trimm(x,y,z):' \
                 + str(trimpar['x'])+',' \
                 + str(trimpar['y'])+',' \
                 + str((trimpar['z'],)) \
                 + ';scale='+str(scale) \
                 + ';fmx='+str(fmax)
      
      # store the sum image
      if store_avg:
          fsum = os.path.join(petudir, 'avg_trimmed-upsampled-scale-'+str(scale)+fcomment+'.nii.gz')
          imio.array2nii( imsumt[::-1,::-1,:], A, fsum, descrip=niidescr)
          if verbose:  print 'i> saved averaged image to:', fsum
          dctout['fsum'] = fsum
  
      # list of file names for the upsampled and trimmed images
      fpetu = []
      # perform the trimming and save the intermediate images if requested
      for i in range(Nim):
  
          # memory saving option, second time doing interpolation
          if memlim:
              if Nim>50 and using_multiple_files:
                      imin_temp = imio.getnii(imdic['files'][i])
                      im = ndi.interpolation.zoom(imin_temp, (scale, scale, scale), order=int_order )
                      if verbose:
                          print 'i> image scaling:', imdic['files'][i]
              else:
                  im = ndi.interpolation.zoom(imin[i,:,:,:], (scale, scale, scale), order=int_order )
          else:
              im = imscl[i,:,:,:]
  
          # trim the scaled image
          imtrim[i, iz0t:, iy0t:iy1t, ix0t:ix1t] = im[iz0s:, iy0s:iy1+1, ix0s:ix1+1]
          
          # save the up-sampled and trimmed PET images
          if store_img_intrmd:
              _frm = '_trmfrm'+str(i)
              _fstr = '_trimmed-upsampled-scale-'+str(scale) + _frm*(Nim>1) +fcomment
              fpetu.append( os.path.join(petudir, fnms[i]+_fstr+'_i.nii.gz') )
              imio.array2nii( imtrim[i,::-1,::-1,:], A, fpetu[i], descrip=niidescr)
              if verbose:  print 'i> saved upsampled PET image to:', fpetu[i]
  
      if store_img:
          _nfrm = '_nfrm'+str(Nim)
          fim = os.path.join(petudir, 'trimmed-upsampled-scale-'+str(scale))+_nfrm*(Nim>1)+fcomment+'.nii.gz'
          imio.array2nii( np.squeeze(imtrim[:,::-1,::-1,:]), A, fim, descrip=niidescr)
          dctout['fim'] = fim
  
      # file names (with paths) for the intermediate PET images
      dctout['fimi'] = fpetu
      dctout['im'] = np.squeeze(imtrim)
      dctout['N'] = Nim
      dctout['affine'] = A
          
      return dctout