Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Minor style updates and QOL changes #830

Closed
wants to merge 5 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -84,3 +84,7 @@ Thumbs.db
doc/source/reference
venv/
.buildbot.patch

# IDE Settings #
############
.vscode/
11 changes: 10 additions & 1 deletion doc/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -213,7 +213,16 @@
#html_use_smartypants = True

# Custom sidebar templates, maps document names to template names.
html_sidebars = {'index': ['localtoc.html', 'relations.html', 'sourcelink.html', 'indexsidebar.html', 'searchbox.html', 'reggie.html']}
html_sidebars = {
'index': [
'localtoc.html',
'relations.html',
'sourcelink.html',
'indexsidebar.html',
'searchbox.html',
'reggie.html',
]
}

# Additional templates that should be rendered to pages, maps page names to
# template names.
Expand Down
3 changes: 2 additions & 1 deletion doc/source/devel/register_me.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
OUR_META = pjoin(OUR_PATH, 'meta.ini')
DISCOVER_INIS = {'user': HOME_INI, 'system': SYS_INI}


def main():
# Get ini file to which to write
try:
Expand All @@ -23,7 +24,7 @@ def main():
reg_to = 'user'
if reg_to in ('user', 'system'):
ini_fname = DISCOVER_INIS[reg_to]
else: # it is an ini file name
else: # it is an ini file name
ini_fname = reg_to

# Read parameters for our distribution
Expand Down
9 changes: 5 additions & 4 deletions doc/source/dicom/derivations/dicom_mosaic.py
Original file line number Diff line number Diff line change
@@ -1,11 +1,12 @@
''' Just showing the mosaic simplification '''
import sympy
from sympy import Matrix, Symbol, symbols, zeros, ones, eye
from sympy import Matrix, Symbol, symbols


def numbered_matrix(nrows, ncols, symbol_prefix):
return Matrix(nrows, ncols, lambda i, j: Symbol(
symbol_prefix + '_{%d%d}' % (i+1, j+1)))


def numbered_vector(nrows, symbol_prefix):
return Matrix(nrows, 1, lambda i, j: Symbol(
symbol_prefix + '_{%d}' % (i+1)))
Expand All @@ -17,12 +18,12 @@ def numbered_vector(nrows, symbol_prefix):
'md_{cols}', 'md_{rows}', 'rd_{cols}', 'rd_{rows}')

md_adj = Matrix((mdc - 1, mdr - 1, 0)) / -2
rd_adj = Matrix((rdc - 1 , rdr - 1, 0)) / -2
rd_adj = Matrix((rdc - 1, rdr - 1, 0)) / -2

adj = -(RS * md_adj) + RS * rd_adj
adj.simplify()

Q = RS[:,:2] * Matrix((
Q = RS[:, :2] * Matrix((
(mdc - rdc) / 2,
(mdr - rdr) / 2))
Q.simplify()
Expand Down
95 changes: 51 additions & 44 deletions doc/source/dicom/derivations/spm_dicom_orient.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,22 +12,25 @@
import sympy
from sympy import Matrix, Symbol, symbols, zeros, ones, eye


# The code below is general (independent of SPMs code)
def numbered_matrix(nrows, ncols, symbol_prefix):
return Matrix(nrows, ncols, lambda i, j: Symbol(
symbol_prefix + '_{%d%d}' % (i+1, j+1)))


def numbered_vector(nrows, symbol_prefix):
return Matrix(nrows, 1, lambda i, j: Symbol(
symbol_prefix + '_{%d}' % (i+1)))


# premultiplication matrix to go from 0 based to 1 based indexing
one_based = eye(4)
one_based[:3,3] = (1,1,1)
one_based[:3, 3] = (1, 1, 1)
# premult for swapping row and column indices
row_col_swap = eye(4)
row_col_swap[:,0] = eye(4)[:,1]
row_col_swap[:,1] = eye(4)[:,0]
row_col_swap[:, 0] = eye(4)[:, 1]
row_col_swap[:, 1] = eye(4)[:, 0]

# various worming matrices
orient_pat = numbered_matrix(3, 2, 'F')
Expand All @@ -40,44 +43,46 @@ def numbered_vector(nrows, symbol_prefix):
slice_thickness = Symbol('\Delta{s}')

R3 = orient_pat * np.diag(pixel_spacing)
R = zeros((4,2))
R[:3,:] = R3

# The following is specific to the SPM algorithm.
x1 = ones((4,1))
y1 = ones((4,1))
y1[:3,:] = pos_pat_0

to_inv = zeros((4,4))
to_inv[:,0] = x1
to_inv[:,1] = symbols('abcd')
to_inv[0,2] = 1
to_inv[1,3] = 1
inv_lhs = zeros((4,4))
inv_lhs[:,0] = y1
inv_lhs[:,1] = symbols('efgh')
inv_lhs[:,2:] = R
R = zeros((4, 2))
R[:3, :] = R3

# The following is specific to the SPM algorithm.
x1 = ones((4, 1))
y1 = ones((4, 1))
y1[:3, :] = pos_pat_0

to_inv = zeros((4, 4))
to_inv[:, 0] = x1
to_inv[:, 1] = symbols('abcd')
to_inv[0, 2] = 1
to_inv[1, 3] = 1
inv_lhs = zeros((4, 4))
inv_lhs[:, 0] = y1
inv_lhs[:, 1] = symbols('efgh')
inv_lhs[:, 2:] = R


def spm_full_matrix(x2, y2):
rhs = to_inv[:,:]
rhs[:,1] = x2
lhs = inv_lhs[:,:]
lhs[:,1] = y2
rhs = to_inv[:, :]
rhs[:, 1] = x2
lhs = inv_lhs[:, :]
lhs[:, 1] = y2
return lhs * rhs.inv()


# single slice case
orient = zeros((3,3))
orient[:3,:2] = orient_pat
orient[:,2] = orient_cross
x2_ss = Matrix((0,0,1,0))
y2_ss = zeros((4,1))
y2_ss[:3,:] = orient * Matrix((0,0,slice_thickness))
orient = zeros((3, 3))
orient[:3, :2] = orient_pat
orient[:, 2] = orient_cross
x2_ss = Matrix((0, 0, 1, 0))
y2_ss = zeros((4, 1))
y2_ss[:3, :] = orient * Matrix((0, 0, slice_thickness))
A_ss = spm_full_matrix(x2_ss, y2_ss)

# many slice case
x2_ms = Matrix((1,1,NZ,1))
y2_ms = ones((4,1))
y2_ms[:3,:] = pos_pat_N
x2_ms = Matrix((1, 1, NZ, 1))
y2_ms = ones((4, 1))
y2_ms[:3, :] = pos_pat_N
A_ms = spm_full_matrix(x2_ms, y2_ms)

# End of SPM algorithm
Expand All @@ -89,22 +94,22 @@ def spm_full_matrix(x2, y2):
single_aff = eye(4)
rot = orient
rot_scale = rot * np.diag(pixel_spacing[:] + [slice_thickness])
single_aff[:3,:3] = rot_scale
single_aff[:3,3] = pos_pat_0
single_aff[:3, :3] = rot_scale
single_aff[:3, 3] = pos_pat_0

# For multi-slice case, we have the start and the end slice position
# patient. This gives us the third column of the affine, because,
# ``pat_pos_N = aff * [[0,0,ZN-1,1]].T
multi_aff = eye(4)
multi_aff[:3,:2] = R3
trans_z_N = Matrix((0,0, NZ-1, 1))
multi_aff[:3, :2] = R3
trans_z_N = Matrix((0, 0, NZ-1, 1))
multi_aff[:3, 2] = missing_r_col
multi_aff[:3, 3] = pos_pat_0
est_pos_pat_N = multi_aff * trans_z_N
eqns = tuple(est_pos_pat_N[:3,0] - pos_pat_N)
solved = sympy.solve(eqns, tuple(missing_r_col))
multi_aff_solved = multi_aff[:,:]
multi_aff_solved[:3,2] = multi_aff_solved[:3,2].subs(solved)
eqns = tuple(est_pos_pat_N[:3, 0] - pos_pat_N)
solved = sympy.solve(eqns, tuple(missing_r_col))
multi_aff_solved = multi_aff[:, :]
multi_aff_solved[:3, 2] = multi_aff_solved[:3, 2].subs(solved)

# Check that SPM gave us the same result
A_ms_0based = A_ms * one_based
Expand All @@ -118,10 +123,10 @@ def spm_full_matrix(x2, y2):
A_i = single_aff
nz_trans = eye(4)
NZT = Symbol('d')
nz_trans[2,3] = NZT
nz_trans[2, 3] = NZT
A_j = A_i * nz_trans
IPP_i = A_i[:3,3]
IPP_j = A_j[:3,3]
IPP_i = A_i[:3, 3]
IPP_j = A_j[:3, 3]

# SPM does it with the inner product of the vectors
spm_z = IPP_j.T * orient_cross
Expand All @@ -132,11 +137,13 @@ def spm_full_matrix(x2, y2):
ipp_sum_div = sum(IPP_j) / sum(orient_cross)
ipp_sum_div = sympy.simplify(ipp_sum_div)


# Dump out the formulae here to latex for the RST docs
def my_latex(expr):
S = sympy.latex(expr)
return S[1:-1]


print 'Latex stuff'
print ' R = ' + my_latex(to_inv)
print ' '
Expand Down
45 changes: 24 additions & 21 deletions doc/source/scripts/make_coord_examples.py
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@
img = nipy.load_image(img_fname)
# Set affine as for FOV, not AC
RZS = img.affine[:3, :3]
vox_fov_center = -(np.array(img.shape) - 1) / 2.
vox_fov_center = -(np.array(img.shape) - 1) / 2.0
T = RZS.dot(vox_fov_center)
img.affine[:3, 3] = T
# Take stuff off the top of the full image, to emphasize FOV
Expand All @@ -63,18 +63,24 @@
epi_br = np.array((92, 70)) * 2
epi_tl = np.array((7, 63)) * 2
# Find lengths of sides
epi_y_len = np.sqrt((np.subtract(epi_bl, epi_tl)**2).sum())
epi_x_len = np.sqrt((np.subtract(epi_bl, epi_br)**2).sum())
epi_y_len = np.sqrt((np.subtract(epi_bl, epi_tl) ** 2).sum())
epi_x_len = np.sqrt((np.subtract(epi_bl, epi_br) ** 2).sum())
x, y = 0, 1
# Make a rectangular box with these sides


def make_ortho_box(bl, x_len, y_len):
""" Make a box with sides parallel to the axes
"""
return np.array((bl,
[bl[x] + x_len, bl[y]],
[bl[x], bl[y] + y_len],
[bl[x] + x_len, bl[y] + y_len]))
return np.array(
(
bl,
[bl[x] + x_len, bl[y]],
[bl[x], bl[y] + y_len],
[bl[x] + x_len, bl[y] + y_len],
)
)


orth_epi_box = make_ortho_box(epi_bl, epi_x_len, epi_y_len)

Expand All @@ -86,8 +92,7 @@ def make_ortho_box(bl, x_len, y_len):


def plot_line(pt1, pt2, fmt='r-', label=None):
plt.plot([pt1[0], pt2[0]], [pt1[1], pt2[1]], fmt,
label=label)
plt.plot([pt1[0], pt2[0]], [pt1[1], pt2[1]], fmt, label=label)


def plot_box(box_def, fmt='r-', label=None):
Expand All @@ -103,18 +108,14 @@ def rotate_box(box_def, angle, origin):
box_def_zeroed = box_def - origin
cost = math.cos(angle)
sint = math.sin(angle)
rot_array = np.array([[cost, -sint],
[sint, cost]])
rot_array = np.array([[cost, -sint], [sint, cost]])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was more readable before.

box_def_zeroed = np.dot(rot_array, box_def_zeroed.T).T
return box_def_zeroed + origin


def labeled_point(pt, marker, text, markersize=10, color='k'):
plt.plot(pt[0], pt[1], marker, markersize=markersize)
plt.text(pt[0] + markersize / 2,
pt[1] - markersize / 2,
text,
color=color)
plt.text(pt[0] + markersize / 2, pt[1] - markersize / 2, text, color=color)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here.



def plot_localizer():
Expand All @@ -126,8 +127,10 @@ def plot_localizer():
def save_plot():
# Plot using global variables
plot_localizer()

def vx2mm(pts):
return pts - iso_center

plot_box(vx2mm(rot_box), label='EPI bounding box')
plot_box(vx2mm(anat_box), 'b-', label='Structural bounding box')
labeled_point(vx2mm(epi_center), 'ro', 'EPI FOV center')
Expand All @@ -145,7 +148,7 @@ def vx2mm(pts):
anat_center = np.mean(anat_box, axis=0)
# y axis on the plot is first axis of image
sag_y, sag_x = sagittal.shape
iso_center = (np.array([sag_x, sag_y]) - 1) / 2.
iso_center = (np.array([sag_x, sag_y]) - 1) / 2.0
sag_extents = [-iso_center[0], iso_center[0], -iso_center[1], iso_center[1]]

# Back to image coordinates
Expand All @@ -155,7 +158,7 @@ def vx2mm(pts):
rot = np.eye(4)
rot[:3, :3] = euler.euler2mat(0, 0, -angle)
# downsample to make smaller output image
downsamp = 1/3
downsamp = 1 / 3
epi_scale = np.diag([downsamp, downsamp, downsamp, 1])
# template voxels to epi box image voxels
vox2epi_vox = epi_scale.dot(rot.dot(epi_trans))
Expand All @@ -165,8 +168,7 @@ def vx2mm(pts):
epi_vox_shape = np.array([data.shape[0], epi_x_len, epi_y_len]) * downsamp
# Make sure dimensions are odd by rounding up or down
# This makes the voxel center an integer index, which is convenient
epi_vox_shape = [np.floor(d) if np.floor(d) % 2 else np.ceil(d)
for d in epi_vox_shape]
epi_vox_shape = [np.floor(d) if np.floor(d) % 2 else np.ceil(d) for d in epi_vox_shape]
# resample, preserving affine
epi_cmap = nca.vox2mni(epi_vox2mm)
epi = rsm.resample(t2_img, epi_cmap, np.eye(4), epi_vox_shape)
Expand All @@ -178,8 +180,9 @@ def vx2mm(pts):
anat_trans[:3, 3] = -np.array([0, anat_box[0, 0], anat_box[0, 1]])
vox2anat_vox = anat_scale.dot(anat_trans)
anat_vox2mm = t1_img.affine.dot(npl.inv(vox2anat_vox))
anat_vox_shape = np.round(np.divide(
[data.shape[0], anat_x_len, anat_y_len], anat_vox_sizes))
anat_vox_shape = np.round(
np.divide([data.shape[0], anat_x_len, anat_y_len], anat_vox_sizes)
)
anat_cmap = nca.vox2mni(anat_vox2mm)
anat = rsm.resample(t1_img, anat_cmap, np.eye(4), anat_vox_shape)
anat_data = anat.get_fdata()
Expand Down
Loading