1
# -*- coding: utf-8 -*-
3
Created on Wed Dec 5 19:07:44 2018
7
Purpose: convert an elastix/transformix deformation field containing displacements
8
to displacement gradients.
10
By default, three files are created: grad.u_x, grad.u_y, grad.u_z
11
These files are saved as .mhd/raw with 32 bit signed floats
14
from skimage import io
17
# input full path to Elastix deformation field here:
18
infile = r"C:\Users\gcg\Documents\launchpad\trunk\DVC\bspline\deformationField.mhd"
20
def compute_du_dx(u_file):
21
""" Compute gradient of displacement field.
22
u_file: absolute filename to .MHD file
24
Creates three different files, grad u_x, grad u_y, nably u_z
26
u = io.imread(u_file, plugin="simpleitk")
27
print "read deformation field, size is %f GB" % (1.0 * sys.getsizeof(u) / (1024*1024*1024))
28
nz, ny, nx, ndim = np.shape(u)
29
print "dimensions : %d * %d * %d * %d" % (nz, ny, nx, ndim)
30
Ne = nz * ny * nx * ndim
31
print "number of entries:", Ne
32
print "bytes per entry: ", sys.getsizeof(u) / Ne
35
print "computing grad u_x"
36
gradient = np.transpose(np.gradient(u[:,:,:,0]), (1,2,3,0))
37
print "shape of gradient is:", np.shape(gradient)
38
outfile = os.path.join(path, "grad_ux.mhd")
39
print "... saving grad u_x"
40
io.imsave(outfile, gradient, "simpleitk")
43
print "computing grad u_y"
44
gradient = np.transpose(np.gradient(u[:,:,:,1]), (1,2,3,0))
45
outfile = os.path.join(path, "grad_uy.mhd")
46
print "... saving grad u_y"
47
io.imsave(outfile, gradient, "simpleitk")
50
print "computing grad u_z"
51
gradient = np.transpose(np.gradient(u[:,:,:,2]), (1,2,3,0))
52
outfile = os.path.join(path, "grad_uz.mhd")
53
print "... saving grad u_z"
54
io.imsave(outfile, gradient, "simpleitk")
57
if __name__ == "__main__":
58
path = os.path.dirname(infile)
59
print "path is: ", path
added
removed
DVC/displacement_gradient_from_deformation_field.py
