Skip to content
Snippets Groups Projects
Commit 1920c99b authored by s224362's avatar s224362
Browse files

Added some notebooks to try using gitlab

parent 67e3abe3
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Notebooks/Agamodon_anguliceps_processing.ipynb 0 → 100644
+ 1089
− 0
View file @ 1920c99b
%% Cell type:markdown id: tags:
### Imports
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from skimage import measure
from scipy.interpolate import RegularGridInterpolator
from scipy import signal
import open3d as o3d
import sympy as sp
from scipy.optimize import curve_fit
import sympy as sp
from mpl_toolkits.mplot3d import Axes3D
import nibabel as nib
import matplotlib.animation as animation
from IPython.display import HTML
import matplotlib
from ipywidgets import interactive, fixed, IntSlider
import math
from skimage.morphology import dilation, ball
from skimage.morphology import skeletonize_3d, skeletonize
from matplotlib.widgets import RectangleSelector
from scipy.ndimage import gaussian_filter
from mayavi import mlab
import nibabel as nib
from matplotlib.ticker import MaxNLocator
```
%% Cell type:markdown id: tags:
### Helper Functions
%% Cell type:code id: tags:
``` python
# Function to marching cube meshes
def plotMesh(ax, mesh, ax_x, ax_y, ax_z, azim, elev):
ax.add_collection3d(mesh)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_xlim(ax_x[0], ax_x[1])
ax.set_ylim(ax_y[0], ax_y[1])
ax.set_zlim(ax_z[0], ax_z[1])
ax.set_aspect('equal')
ax.azim = azim
ax.elev = elev
# Define the 3D polynomial function of 4th degree
def poly_3d(xy, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15):
x, y = xy
return (a0 + a1*x + a2*y + a3*x**2 + a4*y**2 + a5*x*y + a6*x**3 + a7*y**3 + a8*x**2*y + a9*x*y**2 +
a10*x**4 + a11*y**4 + a12*x**3*y + a13*x**2*y**2 + a14*x*y**3)
# Fits polynomial to set of datapoints and displays normalvectors (plot)
def fit_polynomial_surface(points, subset=10, normal_direction=0, num_points_quiver=10, num_points_surface=50, scale_factor=35):
'''
num_points_quiver = 10 # Adjust the number of points as needed for the quivers
num_points_surface = 50 # Adjust the number of points as needed for the surface
scale_factor = 35 # Adjust the scale of the normal vectors
normal_direction = 1 # Change this to 1 (or 0) to reverse the direction of the normals
WANT NORMALS FACING INTO BEANIE SKULL
extra_distance: The extra distance added around the points for sampling
step_size: The step size between each sampling along the normal vectors
steps: The amount of steps done for sampling along normal vectors
grid_size_x: Size of the grid of sample points along the x-axis (can be left out and calculated automatically)
grid_size_y: Size of the grid of sample points along the y-axis (can be left out and calculated automatically)
'''
# Extract x, y, z coordinates from the points (These are subsets defined earlier!)
x_num = points[::subset, 0]
y_num = points[::subset, 1]
z_num = points[::subset, 2]
# Initial guess for the parameters
p0 = np.ones(16) # Needs to match the amount of a's in poly_3d (helper function)
# Use curve_fit to fit the function to the data
popt, pcov = curve_fit(poly_3d, (x_num, y_num), z_num, p0)
# Create a grid of x, y values for the surface
x_grid_surface, y_grid_surface = np.meshgrid(np.linspace(min(x_num), max(x_num), num_points_surface), np.linspace(min(y_num), max(y_num), num_points_surface))
# Compute the corresponding z values for the surface
z_grid_surface = poly_3d((x_grid_surface, y_grid_surface), *popt)
# Create a grid of x, y values for the quivers
x_grid_quiver, y_grid_quiver = np.meshgrid(np.linspace(min(x_num), max(x_num), num_points_quiver), np.linspace(min(y_num), max(y_num), num_points_quiver))
# Compute the corresponding z values for the quivers
z_grid_quiver = poly_3d((x_grid_quiver, y_grid_quiver), *popt)
# Define the symbols
a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15 = sp.symbols('a0 a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15')
x, y = sp.symbols('x y')
# Define the polynomial function
poly_expr = (a0 + a1*x + a2*y + a3*x**2 + a4*y**2 + a5*x*y + a6*x**3 + a7*y**3 + a8*x**2*y + a9*x*y**2 +
a10*x**4 + a11*y**4 + a12*x**3*y + a13*x**2*y**2 + a14*x*y**3)
# Differentiate with respect to x and y
dz_dx = sp.diff(poly_expr, x)
dz_dy = sp.diff(poly_expr, y)
# Create lambdified functions for evaluating the derivatives
dz_dx_func = sp.lambdify((x, y, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15), dz_dx)
dz_dy_func = sp.lambdify((x, y, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15), dz_dy)
# Evaluate the derivatives at each point on the quiver surface
dz_dx_vals = dz_dx_func(x_grid_quiver, y_grid_quiver, *popt)
dz_dy_vals = dz_dy_func(x_grid_quiver, y_grid_quiver, *popt)
# Compute the normal vectors
normals = np.dstack((-dz_dx_vals, -dz_dy_vals, np.ones_like(dz_dx_vals)))
# To change direction of normals
normals *= (-1)**normal_direction
# Normalize the normals
normals /= np.linalg.norm(normals, axis=2, keepdims=True)
# Create a 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Plot the original points
ax.scatter(x_num, y_num, z_num, color='b')
# Plot the fitted surface using the surface grid
ax.plot_surface(x_grid_surface, y_grid_surface, z_grid_surface, color='r', alpha=0.5)
# Plot the normal vectors using the quiver grid
for i in range(len(x_grid_quiver)):
for j in range(len(y_grid_quiver)):
ax.quiver(x_grid_quiver[i, j], y_grid_quiver[i, j], z_grid_quiver[i, j], normals[i, j, 0], normals[i, j, 1], normals[i, j, 2],
length=scale_factor, linewidth=5, color='g')
ax.set_aspect('equal')
plt.show()
return popt, dz_dx_func, dz_dy_func
# Generates PCA-coordinates to sample
def generate_sample_points(points, popt, dz_dx_func, dz_dy_func, normal_direction=0, subset=10, extra_distance=15, step_size=0.1, min_steps=10, grid_size_x=None, grid_size_y=None):
'''
subset: amount of subsamples used
extra_distance: The extra distance added around the points for sampling
step_size: The step size between each sampling along the normal vectors
steps: The amount of steps done for sampling along normal vectors
grid_size_x: Size of the grid of sample points along the x-axis (can be left out and calculated automatically)
grid_size_y: Size of the grid of sample points along the y-axis (can be left out and calculated automatically)
STANDARD FOR AGAMODON:
extra_distance = 15
step_size = 0.1
min_steps: The number of steps in the least sampled direction! (The other direction is sampled 5x as much)
grid_size_x = 500
grid_size_y = 200
'''
# Extract x, y, z coordinates from the points (These are subsets defined earlier!)
x_num = points[::subset, 0]
y_num = points[::subset, 1]
z_num = points[::subset, 2]
# For calculating the span of the points to generate meshgrid
min_values = np.min(points, axis=0)
max_values = np.max(points, axis=0)
# If user has not chosen, make it the size of the spans
if grid_size_x is None:
grid_size_x = int(max_values[0] - min_values[0])
if grid_size_y is None:
grid_size_y = int(max_values[1] - min_values[1])
# Step 1: Create a grid in the x-y-plane
x_grid, y_grid = np.meshgrid(np.linspace(min(x_num) - extra_distance, max(x_num) + extra_distance, grid_size_x), np.linspace(min(y_num) - extra_distance, max(y_num) + extra_distance, grid_size_y))
# Step 2: Compute the height of the fitted polynomial surface at each point
z_grid = poly_3d((x_grid, y_grid), *popt)
# Step 3: Compute the normal vector at each point on the grid
dz_dx_vals = dz_dx_func(x_grid, y_grid, *popt) # defined priorly.
dz_dy_vals = dz_dy_func(x_grid, y_grid, *popt)
normals = np.dstack((-dz_dx_vals, -dz_dy_vals, np.ones_like(dz_dx_vals)))
normals *= (-1)**normal_direction # To change direction of normals
# Step 4: Normalize the normals for better visualization
normals /= np.linalg.norm(normals, axis=2, keepdims=True)
# Step 5: Go `steps` steps in the direction of the normal vector and `steps` steps in the opposite direction to create the layers
# Modify the number of steps in each direction
steps_positive = 5 * min_steps
steps_negative = min_steps
grid_3d = np.zeros((grid_size_y, grid_size_x, steps_positive + steps_negative + 1, 3)) # 3D array to store the x, y, and z coordinates of the points (each index in the 3D matrix contains a 3D-coordinate in PCA-space)
# Go steps_negative steps in the opposite direction
for i in range(-steps_negative, 0):
displacement = i * step_size * normals # displacement in x, y, and z directions
grid_3d[:, :, i+steps_negative, 0] = x_grid + displacement[:, :, 0] # x coordinate
grid_3d[:, :, i+steps_negative, 1] = y_grid + displacement[:, :, 1] # y coordinate
grid_3d[:, :, i+steps_negative, 2] = z_grid + displacement[:, :, 2] # z coordinate
# Go steps_positive steps in the direction of the normal vector
for i in range(steps_positive):
displacement = i * step_size * normals # displacement in x, y, and z directions
grid_3d[:, :, i+steps_negative+1, 0] = x_grid + displacement[:, :, 0] # x coordinate
grid_3d[:, :, i+steps_negative+1, 1] = y_grid + displacement[:, :, 1] # y coordinate
grid_3d[:, :, i+steps_negative+1, 2] = z_grid + displacement[:, :, 2] # z coordinate
# Flatten the 3D grid
points_flat = grid_3d.reshape(-1, 3)
# Plot a small subset of the found points
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points_flat[::2000, 0], points_flat[::2000, 1], points_flat[::2000, 2])
plt.show()
return points_flat, grid_3d.shape[0:3]
# Interpolate sample points in original volume file
def sample_in_original_volume(points_original, original_nii_path, grid_shape, intMethod='linear', expVal=0.0):
# Load in original volume
original = nib.load(original_nii_path)
imgDim = original.header['dim'][1:4]
# Set-up interpolators for moving image
x = np.arange(start=0, stop=imgDim[0], step=1)
y = np.arange(start=0, stop=imgDim[1], step=1)
z = np.arange(start=0, stop=imgDim[2], step=1)
F_moving = RegularGridInterpolator((x, y, z), original.get_fdata().astype('float16'), method=intMethod, bounds_error=False, fill_value=expVal)
# Evaluate transformed grid points in the moving image
fVal = F_moving(points_original)
# Reshape to voxel grid
volQ = np.reshape(fVal,newshape=grid_shape).astype('float32')
return volQ
# Animate moving through the layers of volume
def ani_through_volume(volume, indices=None):
'''
indices: List [a,b] of layer numbers you want animated
'''
def update_img(z):
ax.clear()
ax.imshow(volume[:,:,z])
ax.set_title('Layer ' + str(z))
# User can specify indices they want to see
if indices != None:
a = indices[0]
b = indices[1]
else:
a = 0
b = volume.shape[2] - 1
matplotlib.rcParams['animation.embed_limit'] = 2**128
fig, ax = plt.subplots()
ani = animation.FuncAnimation(fig, update_img, frames=range(a,b+1), interval=200)
return HTML(ani.to_jshtml())
# Create interactive plot for slicing
def interactive_plot_slicing(volume):
def get_slice(vol, x_start, x_end, y_start, y_end, z):
slice = vol[x_start:x_end, y_start:y_end, z]
plt.imshow(slice)
plt.show()
x_start_slider = IntSlider(min=0, max=volume.shape[0]-1, value=0)
x_end_slider = IntSlider(min=0, max=volume.shape[0]-1, value=volume.shape[0]-1)
y_start_slider = IntSlider(min=0, max=volume.shape[1]-1, value=0)
y_end_slider = IntSlider(min=0, max=volume.shape[1]-1, value=volume.shape[1]-1)
z_slider = IntSlider(min=0, max=volume.shape[2]-1, value=0)
def update_x_end_range(*args):
x_end_slider.min = x_start_slider.value
x_start_slider.observe(update_x_end_range, 'value')
def update_y_end_range(*args):
y_end_slider.min = y_start_slider.value
y_start_slider.observe(update_y_end_range, 'value')
interactive_plot = interactive(get_slice,
vol=fixed(volume),
x_start=x_start_slider,
x_end=x_end_slider,
y_start=y_start_slider,
y_end=y_end_slider,
z=z_slider)
return interactive_plot
# Allows user to remove unwanted sutures
def interactive_process_volume(volume, interactive_plot, z_vals):
# Define the slicing parameters
x_start = interactive_plot.kwargs['x_start']
x_end = interactive_plot.kwargs['x_end']
y_start = interactive_plot.kwargs['y_start']
y_end = interactive_plot.kwargs['y_end']
min_z, max_z = z_vals[0], z_vals[1]
# Crop the volume based on the provided slices
cropped_volume = volume[x_start:x_end, y_start:y_end, min_z:max_z+1]
# Calculate the median layer index
median_z = cropped_volume.shape[2] // 2
median_layer = cropped_volume[:, :, median_z]
# Create a new volume to apply changes
new_volume = cropped_volume.copy()
fig, ax = plt.subplots()
ax.imshow(median_layer, cmap='gray')
def line_select_callback(eclick, erelease):
nonlocal new_volume
x1, y1 = int(eclick.xdata), int(eclick.ydata)
x2, y2 = int(erelease.xdata), int(erelease.ydata)
if x1 > x2:
x1, x2 = x2, x1
if y1 > y2:
y1, y2 = y2, y1
# Find the maximum value within the selected area
max_value = np.max(new_volume[y1:y2+1, x1:x2+1, :])
# Assign the maximum value to all pixels within the selected area for each slice
for z in range(new_volume.shape[2]):
new_volume[y1:y2+1, x1:x2+1, z] = max_value
# Update the displayed image in real-time
ax.clear()
ax.imshow(new_volume[:, :, median_z], cmap='gray')
plt.draw()
rs = RectangleSelector(ax, line_select_callback,
useblit=True,
button=[1], # Only use left button
minspanx=5, minspany=5,
spancoords='pixels',
interactive=True)
plt.show(block=True) # Block until plot window is closed
# Return the new_volume after the plot is closed
return new_volume
# Interactive plot for deciding threshold value
def interactive_plot_threshold(volume):
# Define the function to generate the plot
def plot_threshold(threshold):
slice = volume[:,:, volume.shape[2]//2]
slice_bin = np.where(slice >= threshold, 1, 0)
slice_bin= 1 - slice_bin # Make suture the "foreground"
plt.imshow(slice_bin)
plt.show()
# Calculate the minimum value for the slider
min_val = (volume.max()+volume.min())//2
# Round up to the nearest hundred
min_val = math.ceil(min_val / 100.0) * 100
# Create a slider for the threshold value
threshold_slider = IntSlider(min=min_val, max=volume.max(), step=100, value=min_val)
# Create the interactive plot
interactive_plot_new = interactive(plot_threshold, threshold=threshold_slider)
return interactive_plot_new
# Extracts the suture depending on your former choices
def extract_suture(volume, points_orig, interactive_plot_slice, interactive_plot_thresh, z_vals, grid_shape):
'''
Input volume is sliced
'''
min_z = z_vals[0]
max_z = z_vals[1]
num_indices = max_z - min_z
# Extrac chosen indices
x_start = interactive_plot_slice.kwargs['x_start']
x_end = interactive_plot_slice.kwargs['x_end']
y_start = interactive_plot_slice.kwargs['y_start']
y_end = interactive_plot_slice.kwargs['y_end']
# Extract chosen threshold
threshold_value = interactive_plot_thresh.kwargs['threshold']
# Generate a slice
slice = volume[:,:, volume.shape[2]//2]
# Generate empty array to store suture points in
suture = np.zeros((slice.shape[0], slice.shape[1], volume.shape[2]))
# For animation
binary_slices = np.zeros((slice.shape[0], slice.shape[1], num_indices))
# Two first indicies were bad in segmentation
for index in range(num_indices):
slice = volume[:,:, index]
slice_bin = np.where(slice >= threshold_value, 1, 0)
slice_bin= 1 - slice_bin # Make suture the "foreground"
# Label the connected components in the binary image
labels = measure.label(slice_bin)
# Now, you can analyze each individual blob
properties = measure.regionprops(labels)
# Sort the regions by area_bbox in descending order
properties_sorted = sorted(properties, key=lambda x: x.bbox_area, reverse=True)
# Get the region with the largest area_bbox
largest_bbox_region = properties_sorted[0]
# Get the coordinates of the pixels in the blob with the largest area_bbox
coords_largest_bbox = largest_bbox_region.coords
# Extract x and y coordinates
x_coords = coords_largest_bbox[:, 1]
y_coords = coords_largest_bbox[:, 0]
suture[y_coords, x_coords, index] = 1
binary_slices[:,:, index] = slice_bin
### ANIMATION
# Create a figure with two subplots
fig, axs = plt.subplots(1, 2, figsize=(8,3))
# Initialize the images
im1 = axs[0].imshow(suture[:, :, 0], animated=True)
im2 = axs[1].imshow(volume[:,:, 0], animated=True)
# Initialize the text
txt = axs[0].text(0.5, 1.01, f'Layer: {0}', transform=axs[0].transAxes)
# Update function for the animation
def updatefig(i):
im1.set_array(suture[:, :, i])
im2.set_array(volume[:,:, i])
txt.set_text(f'Layer: {i}')
return im1, im2, txt
# Create the animation
ani = animation.FuncAnimation(fig, updatefig, frames=range(num_indices), interval=200, blit=True)
### GET ORIGINAL COORDINATES
new_shape = grid_shape + (3,)
coords_orig = points_orig.reshape(new_shape)
# cooresponding original coordinates to sliced area
coords_orig_sliced = coords_orig[x_start:x_end, y_start:y_end, min_z:min_z+num_indices]
### PLOT ORIGINAL COORDINATES OF SLICED AREA
# Get the indices where suture is 1
indices = np.where(suture == 1)
# Get the corresponding points in coords_orig_sliced
points = coords_orig_sliced[indices]
# Create a 3D scatter plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[::10, 0], points[::10, 1], points[::10, 2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
return ani, points, suture, coords_orig_sliced
```
%% Cell type:markdown id: tags:
## "unfolded" volume generation
%% Cell type:markdown id: tags:
### Plot point-cloud from input file
%% Cell type:code id: tags:
``` python
input_file = 'Agamodon_anguliceps_cut_PLY.ply'
pcd = o3d.io.read_point_cloud(input_file)
points = np.asarray(pcd.points)
# Plot a subset of the points
# Choose a subset of points and view point of display
subset = 10
viewPoint = [90, -45]
fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(1,1,1, projection='3d')
ax.scatter(points[::subset,0], points[::subset,1], points[::subset,2], c='k', marker='.')
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_aspect('auto')
ax.azim = viewPoint[0]
ax.elev = viewPoint[1]
ax.set_title('Mesh-points in original domain')
plt.tight_layout()
plt.show()
```
%% Cell type:markdown id: tags:
### Transform point-cloud into PCA-space
%% Cell type:code id: tags:
``` python
# Principal Component Analysis (PCA) of the point cloud
points_mu = np.mean(points, axis = 0) # Center the point cloud
cov = np.cov(np.transpose(points)) # Covariance matrix
# SVD or Eigendecomposition
U,S,V = np.linalg.svd(cov,full_matrices=False)
# Apply transform to rotated pointcloud
points_rot = (points - points_mu) @ U
```
%% Cell type:markdown id: tags:
### Display points in PCA-space
%% Cell type:code id: tags:
``` python
%matplotlib qt5
```
%% Cell type:code id: tags:
``` python
# Plot the rotated pointcloud
fig = plt.figure(figsize=(8,4))
# display the original point cloud with the principal axes
ax = fig.add_subplot(1,2,1, projection='3d')
ax.scatter(points[::subset,0], points[::subset,1], points[::subset,2], c='k', marker='.')
# Set the number of ticks on each axis
ax.xaxis.set_major_locator(MaxNLocator(nbins=5))
ax.yaxis.set_major_locator(MaxNLocator(nbins=5))
ax.zaxis.set_major_locator(MaxNLocator(nbins=5))
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_aspect('equal')
ax.azim = 40
ax.elev = -25
ax.set_title('Original point cloud')
# Add the principal axes to the plot using quiver
ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,0], U[1,0], U[2,0], color='r', length=100, normalize=True, linewidth=2.5)
ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,1], U[1,1], U[2,1], color='g', length=100, normalize=True, linewidth=2.5)
ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,2], U[1,2], U[2,2], color='b', length=100, normalize=True, linewidth=2.5)
# Display the rotated point cloud
ax = fig.add_subplot(1,2,2, projection='3d')
ax.scatter(points_rot[::subset,0], points_rot[::subset,1], points_rot[::subset,2], c='k', marker='.')
# Set the number of ticks on each axis
ax.xaxis.set_major_locator(MaxNLocator(nbins=5))
ax.yaxis.set_major_locator(MaxNLocator(nbins=5))
ax.zaxis.set_major_locator(MaxNLocator(nbins=2))
ax.set_xlabel("PC1")
ax.set_ylabel("PC2")
ax.set_zlabel("PC3")
ax.set_aspect('equal')
ax.azim = 90
ax.elev = -45
ax.set_title('Rotated point cloud')
plt.tight_layout()
```
%% Cell type:markdown id: tags:
### Fit polynomial surface in PCA-space
%% Cell type:code id: tags:
``` python
%matplotlib inline
normal_direction = 1
popt, dz_dx_func, dz_dy_func = fit_polynomial_surface(points_rot, normal_direction=normal_direction)
```
%% Cell type:markdown id: tags:
### Generate sample point coordinates in PCA-space (and plot the sample points in PCA-space)
%% Cell type:code id: tags:
``` python
points_flat, grid_shape = generate_sample_points(points_rot, popt, dz_dx_func, dz_dy_func, normal_direction=normal_direction, extra_distance=15, min_steps=20, step_size=0.4)
```
%% Cell type:markdown id: tags:
### Transform sample points into world coordinates
%% Cell type:code id: tags:
``` python
# From SVD earlier
points_orig = points_flat @ U.T + points_mu
```
%% Cell type:markdown id: tags:
#### Visualise the sample points and the original surface
%% Cell type:code id: tags:
``` python
%matplotlib qt5
```
%% Cell type:code id: tags:
``` python
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points_orig[::1000, 0], points_orig[::1000, 1], points_orig[::1000, 2], alpha=0.3)
ax.scatter(points[::30,0], points[::30,1], points[::30,2])
plt.show()
```
%% Cell type:markdown id: tags:
### Sample points_orig's values in the original volume
%% Cell type:code id: tags:
``` python
path_nii = 'G:/Mit drev/DTU/Fagprojekt/DATA/Agamodon_anguliceps_CAS153467/Agamodon_anguliceps_CAS153467000.nii'
volQ = sample_in_original_volume(points_orig, path_nii, grid_shape) # Can take some seconds
```
%% Cell type:code id: tags:
``` python
# Save the new sampled volume
niiPCA = nib.Nifti1Image(volQ, np.eye(4))
name = 'Agamodon_anguliceps_unfolded_volume_new.nii'
nib.save(niiPCA, name)
```
%% Cell type:markdown id: tags:
## Suture extraction
%% Cell type:markdown id: tags:
### Load the NIfTI-file back in
%% Cell type:code id: tags:
``` python
# Load the .nii file in
name = 'Agamodon_anguliceps_unfolded_volume_new.nii'
nii_img = nib.load(name)
# Get the data as a numpy array
volQ = nii_img.get_fdata()
# Print max and min values
print('Volume minimum: ', volQ.min())
print('Volume maximum: ', volQ.max())
```
%% Cell type:markdown id: tags:
### Animate moving through the layers of the "unfolded" volume
%% Cell type:code id: tags:
``` python
ani_through_volume(volQ)
```
%% Cell type:markdown id: tags:
### Crop slices
%% Cell type:code id: tags:
``` python
%matplotlib inline
interactive_plot_slice = interactive_plot_slicing(volQ)
interactive_plot_slice
```
%% Cell type:markdown id: tags:
### Remove unwanted sutures interactively
%% Cell type:code id: tags:
``` python
%matplotlib qt5
z_vals = [60,61] # Choose z_vals you want to work on
new_volume = interactive_process_volume(volQ, interactive_plot_slice, z_vals)
```
%% Cell type:markdown id: tags:
### Choose layers and thresholding value
%% Cell type:code id: tags:
``` python
%matplotlib inline
interactive_plot_thresh = interactive_plot_threshold(new_volume)
interactive_plot_thresh
```
%% Cell type:code id: tags:
``` python
%matplotlib qt5
ani, points, suture, coords_orig = extract_suture(new_volume, points_orig, interactive_plot_slice, interactive_plot_thresh, z_vals, grid_shape)
HTML(ani.to_jshtml())
```
%% Cell type:markdown id: tags:
### Skeletonize
%% Cell type:markdown id: tags:
- Chooses only ONE layer to skeletonize and represent the suture (curve)
%% Cell type:code id: tags:
``` python
# Choose layer
z_val = 0
# Create skeleton
skeleton = skeletonize(suture[:,:,z_val])
# Get the 2D coordinates of the skeleton points
y, x = np.where(skeleton)
# Create a 2D plot
fig, ax = plt.subplots()
ax.scatter(x, y, color='r')
# Set labels and title
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('2D Skeleton')
# Show the plot
plt.show()
```
%% Cell type:code id: tags:
``` python
# Get the 3D coordinates of the skeleton points
skeleton_indices = np.where(skeleton == 1)
skeleton_coords = coords_orig[skeleton_indices + (z_val,)]
skeleton_coords = np.squeeze(skeleton_coords)
# Separate the x, y, and z coordinates
x, y, z = skeleton_coords[:, 0], skeleton_coords[:, 1], skeleton_coords[:, 2]
# Create a 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, color='r')
# Set labels and title
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D Skeleton')
# Show the plot
plt.show()
```
%% Cell type:markdown id: tags:
Export the skeleton
%% Cell type:code id: tags:
``` python
np.save("Agamodon_anguliceps_skeleton.npy", skeleton_coords)
```
%% Cell type:markdown id: tags:
### Plot of original cranium and extracted suture
%% Cell type:code id: tags:
``` python
# Plot the original point cloud with the extracted skeleton line
input_file = 'Agamodon_anguliceps_cut_PLY.ply'
pcd = o3d.io.read_point_cloud(input_file)
points_from_input = np.asarray(pcd.points)
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points_from_input[::subset,0], points_from_input[::subset,1], points_from_input[::subset,2], c='k', marker='.', alpha=0.7) # Make points semi-transparent
ax.scatter(skeleton_coords[:, 0], skeleton_coords[:, 1], skeleton_coords[:, 2], color='r', s=50) # Add mean_points to the plot with larger size
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_aspect('equal')
ax.azim = 40
ax.elev = -25
ax.set_title('Original point cloud')
# Add the principal axes to the plot using quiver
ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,0], U[1,0], U[2,0], color='r', length=2, normalize=True, linewidth=2.5)
ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,1], U[1,1], U[2,1], color='g', length=2, normalize=True, linewidth=2.5)
ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,2], U[1,2], U[2,2], color='b', length=2, normalize=True, linewidth=2.5)
plt.tight_layout()
```
%% Cell type:markdown id: tags:
### Mayavi show cranium with extracted points
%% Cell type:code id: tags:
``` python
# Load the .nii file
img = nib.load('Agamodon_anguliceps_CAS153467000.nii')
data = img.get_fdata()
# Display the 3D volume
mlab.contour3d(data)
# WITH EXTRACTED POINTS
mlab.points3d(points[:, 0], points[:, 1], points[:, 2], color=(1, 0, 0), scale_factor=30)
mlab.show()
```
%% Cell type:markdown id: tags:
### Mayavi show cranium with skeleton
%% Cell type:code id: tags:
``` python
# Load the .nii file
img = nib.load('Agamodon_anguliceps_CAS153467000.nii')
data = img.get_fdata()
# Display the 3D volume
mlab.contour3d(data)
# WITH SKELETON-POINTS
mlab.points3d(skeleton_coords[:, 0], skeleton_coords[:, 1], skeleton_coords[:, 2], color=(1, 0, 0), scale_factor=32)
mlab.show()
```
%% Cell type:markdown id: tags:
### Display extracted points and skeleton in PCA-space
%% Cell type:code id: tags:
``` python
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
# Assuming your data is in a variable called data
pca = PCA(n_components=2)
pca_points = pca.fit_transform(points)
# Plot the transformed points
plt.scatter(pca_points[:, 0], pca_points[:, 1])
# Transform the skeleton_coords into the same PCA space
pca_skeleton_coords = pca.transform(skeleton_coords)
# Plot the skeleton_coords in the same PCA space
plt.scatter(pca_skeleton_coords[:, 0], pca_skeleton_coords[:, 1], color='r')
plt.show()
```
Notebooks/Beanies_processing.ipynb 0 → 100644
+ 1726
− 0
View file @ 1920c99b
This diff is collapsed.
Notebooks/Curve_Analysis_Script.ipynb 0 → 100644
+ 380
− 0
View file @ 1920c99b
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment