-
Notifications
You must be signed in to change notification settings - Fork 2
/
fixto1809c_quick
97 lines (79 loc) · 3.29 KB
/
fixto1809c_quick
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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
#!/usr/bin/env python
import sys
print("For 2.3mm template extens differ by [-2.850006 , 0.1499939 , 0.55000305]")
print("no fast way to deal with this!")
sys.exit(1)
"""
one voxel shift
84 100 84 1 old 2.3.nii
83 99 83 1 2.3.nii
96 114 96 1 old 2.nii
97 115 97 1 2.nii
64 76 64 1 old 3.nii
65 77 65 1 3.nii
"""
import nibabel as nib
import numpy as np
badfile = "/opt/ni_tools/standard_old/mni_icbm152_nlin_asym_09c/mni_icbm152_t1_tal_nlin_asym_09c_2.3mm.nii"
#badfile = "/Volumes/Zeus/preproc/petrest_rac1/MHRest_FM_ica/11403_20151106/func_to_template.nii.gz"
#correct_ex = "/Volumes/Zeus/tmp/11403_20141106_func_to_template.1809c.nii.gz"
outfile = "/Volumes/Zeus/tmp/11403_20141106_func_to_template.nii.gz"
new_template = "/opt/ni_tools/standard/mni_icbm152_nlin_asym_09c/mni_icbm152_t1_tal_nlin_asym_09c_2.3mm.nii"
orig = nib.load(badfile)
template = nib.load(new_template)
orig_mat = orig.get_fdata()
o_dim = orig.shape
dim_diff = np.array(o_dim[0:3]) - np.array(template.shape)
# we only handle off by 1 voxel in all dims
assert np.all(np.abs(dim_diff) == 1)
if dim_diff[0] > 0: # 2.3mm
# dims: sag, cor, ax
# >>> template.shape
# (83, 99, 83)
# >>> orig.shape
# (84, 100, 84, 320)
# x=np.array([np.array([x.header['qoffset_'+k] for k in ['x','y','z']]) for x in [template, orig]])
# [[ -92.59999, -132. , -78. ], # template
# [ -95.45 , -131.85 , -77.45 ]] # bad
# np.diff(x,axis=0) # [-2.850006 , 0.1499939 , 0.55000305]
if orig.ndim == 3:
#adjusted = orig_mat[:(o_dim[0]-2),:(o_dim[1]-2),:(o_dim[2]-2)]
adjusted = orig_mat[2:, 2:, 2:]
elif orig.ndim == 4:
adjusted = orig_mat[:(o_dim[0]-1),:(o_dim[1]-1),:(o_dim[2]-1),:]
else:
raise Exception(f"can only handle 3 or 4 dims, not {orig.ndim}")
input_mat = adjusted
else:
input_mat = orig_mat
if orig.ndim == 3:
adjusted = np.zeros(template.shape)
adjusted[1:,1:,1:] = input_mat
if orig.ndim == 4:
adjusted = np.zeros([*template.shape, o_dim[3]])
adjusted[1:,1:,1:,:] = input_mat
# make sure we got the shape right
assert adjusted.shape[0:3] == template.shape
adjusted_nii = nib.Nifti1Image(adjusted, template.affine)
nib.save(adjusted_nii, outfile)
# from nipy.labs.viz import plot_map, mni_sform, coord_transform
#from nilearn.plotting import plot_anat
def edgeme(d):
edge = cv.Canny(np.uint8(d), 500, 1000)#.astype(bool)
return edge
def plotme(d1, d2, i=60, j=60, k=60):
# from matplotlib import pyplot as plt
# import cv2 as cv
d_diff = d1 - d2
plt.subplot(4,3,1); plt.imshow(d1[i,:,:])
plt.subplot(4,3,2); plt.imshow(d1[:,j,:])
plt.subplot(4,3,3); plt.imshow(d1[:,:,k])
plt.subplot(4,3,4); plt.imshow(d2[i,:,:])
plt.subplot(4,3,5); plt.imshow(d2[:,j,:])
plt.subplot(4,3,6); plt.imshow(d2[:,:,k])
plt.subplot(4,3,7); plt.imshow(d_diff[i,:,:])
plt.subplot(4,3,8); plt.imshow(d_diff[:,j,:])
plt.subplot(4,3,9); plt.imshow(d_diff[:,:,k])
plt.subplot(4,3,10); plt.imshow(edgeme(d1[i,:,:])*10+edgeme(d2[i,:,:]))
plt.subplot(4,3,11); plt.imshow(edgeme(d1[:,j,:])*10+edgeme(d2[:,j,:]))
plt.subplot(4,3,12); plt.imshow(edgeme(d1[:,:,k])*10+edgeme(d2[:,:,k]))