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