"""
    @file  dtidroid-convert.py
    @brief Convert DTI-DROID deformation field to ITK format with
           LPS orientation of displacement vectors in spatial units.

    Copyright (c) 2012 University of Pennsylvania.

    This file is part of DTI-DROID.

    DTI-DROID is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    DTI-DROID is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with DTI-DROID.  If not, see <http://www.gnu.org/licenses/>.

    Contact: SBIA Group <sbia-software at uphs.upenn.edu>
"""

import sys
import numpy
import nifti

# ----------------------------------------------------------------------------
if __name__ == '__main__':
  # show usage if called with incorrect number of arguments
  if len(sys.argv) != 5:
    print "Usage: dtidroid-convert droidField templateImageFile movingImageFile outputField"
    sys.exit(1)
  # arguments
  droidFile    = sys.argv[1]
  templateFile = sys.argv[2]
  movingFile   = sys.argv[3]
  outFile      = sys.argv[4]
  # read template and moving image
  templateImage = nifti.NiftiImage(templateFile)
  movingImage   = nifti.NiftiImage(movingFile)
  # droid does not use the origen feild thus we need a global displacement vector
  # to overlay the orgins (convert it to LPS)
  globalVec = movingImage.getQForm()[0:3,3] - templateImage.getQForm()[0:3,3]
  globalVec[0] = globalVec[0] * -1
  globalVec[1] = globalVec[1] * -1
  # read DROID deformation field
  dims    = templateImage.extent[0:3]
  nVoxels = numpy.product(dims)
  f  = numpy.fromfile(droidFile,numpy.float32)
  f2 = numpy.reshape(f,[3,nVoxels],'F')
  f3 = f2[[1,0,2],:]
  # voxels -> RAS
  d = numpy.copy(templateImage.getQForm()[0:3,0:3])
  # voxels -> LPS
  d[0,:] = -1 * d[0,:]
  d[1,:] = -1 * d[1,:]
  # convert deformation field
  v1 = numpy.zeros(f3.shape,dtype=numpy.float32)
  for i in numpy.arange(0, nVoxels):
    v1[:,i] = numpy.dot(d,f3[:,i]) + globalVec
  # save converted deformation field
  v2 = numpy.transpose( numpy.reshape(v1,[3,1,dims[0],dims[1],dims[2]],'F') , (0,1,4,3,2))
  defIm = nifti.NiftiImage(v2,header=templateImage.header)
  defIm.raw_nimg.intent_code = 1007
  defIm.save(outFile)
