|
@@ -6,7 +6,7 @@ from scipy.ndimage import rotate
|
|
|
# Define the file paths
|
|
|
templateFilePath = r"C:\Users\aswen\Desktop\Code\2024_Ruthe_SND\input\average_template_50_from_website.nii.gz"
|
|
|
densityFilePath = r"C:\Users\aswen\Desktop\Code\2024_Ruthe_SND\input\12_wks_coronal_100141780_50um_projection_density.nii.gz"
|
|
|
-maskFilePath = r"C:\Users\aswen\Desktop\Code\2024_Ruthe_SND\output\Nifti_Trakte_registered\CC_Mop_dilated_registered.nii.gz"
|
|
|
+maskFilePath = r"\\10.209.5.114\Projects\Student_projects\14_Aref_Kalantari_2021\Luca\CC_Mop_desnitymap_50um.nii.gz"
|
|
|
|
|
|
# Load the NIfTI files
|
|
|
templateNifti = nib.load(templateFilePath)
|
|
@@ -22,20 +22,64 @@ dpi = 400
|
|
|
# Get the number of slices along each dimension
|
|
|
num_slices_dim0, num_slices_dim1, num_slices_dim2 = templateNifti.shape
|
|
|
|
|
|
-# Create a 3D array to hold the overlay
|
|
|
-overlay_volume = np.zeros((num_slices_dim0, num_slices_dim1, num_slices_dim2, 4), dtype=np.uint8)
|
|
|
-
|
|
|
-# Loop through and overlay all slices
|
|
|
-for sliceIdx in range(num_slices_dim2):
|
|
|
- # Extract the desired slice from the template, density, and mask
|
|
|
- templateSlice = np.squeeze(templateNifti.get_fdata()[:, :, sliceIdx]).astype(np.uint16)
|
|
|
- densitySlice = np.squeeze(densityNifti.get_fdata()[:, :, sliceIdx]).astype(np.double)
|
|
|
-
|
|
|
- mask_data = maskNifti.get_fdata()
|
|
|
- mask_data = mask_data.swapaxes(0, 2)
|
|
|
- mask_data = np.flip(mask_data, axis=0)
|
|
|
-
|
|
|
- maskSlice = np.squeeze(mask_data[:, :, sliceIdx])
|
|
|
+# Set the slice step for display (every tenth slice)
|
|
|
+slice_step = 10
|
|
|
+
|
|
|
+# Choose the dimension for iteration (0, 1, or 2)
|
|
|
+iteration_dimension = 0 # Change this value to select the dimension
|
|
|
+
|
|
|
+# Calculate the number of rows and columns for the grid based on the selected dimension
|
|
|
+if iteration_dimension == 0:
|
|
|
+ num_slices = num_slices_dim0
|
|
|
+elif iteration_dimension == 1:
|
|
|
+ num_slices = num_slices_dim1
|
|
|
+else:
|
|
|
+ num_slices = num_slices_dim2
|
|
|
+
|
|
|
+grid_size = int(np.sqrt(num_slices // slice_step))
|
|
|
+if grid_size * grid_size < num_slices // slice_step:
|
|
|
+ grid_size += 1
|
|
|
+
|
|
|
+cm = 1/2.54
|
|
|
+# Adjust the figure size to be larger
|
|
|
+fig, axs = plt.subplots(grid_size, grid_size, figsize=(20, 20), dpi=dpi)
|
|
|
+
|
|
|
+# Hide axis numbers
|
|
|
+for ax_row in axs:
|
|
|
+ for ax in ax_row:
|
|
|
+ ax.axis('off')
|
|
|
+
|
|
|
+# Load the mask and perform necessary flips
|
|
|
+mask_data = np.squeeze(maskNifti.get_fdata())
|
|
|
+# Swap the mask over the x and y axes
|
|
|
+mask_data = mask_data.swapaxes(0, 2)
|
|
|
+mask_data = np.flip(mask_data, axis=0)
|
|
|
+
|
|
|
+# Create a new NIfTI image with the adjusted mask data
|
|
|
+adjusted_mask_nifti = nib.Nifti1Image(mask_data, maskNifti.affine)
|
|
|
+
|
|
|
+# Save the new NIfTI image with adjusted header
|
|
|
+output_mask_filepath = r"\\10.209.5.114\Projects\Student_projects\14_Aref_Kalantari_2021\Luca\CC_Mop_desnitymap_50um_reoriented.nii.gz" # Specify the desired output file path
|
|
|
+nib.save(adjusted_mask_nifti, output_mask_filepath)
|
|
|
+
|
|
|
+# Loop through and display every nth slice as subplots along the selected dimension
|
|
|
+for i, sliceIdx in enumerate(range(0, num_slices, slice_step)):
|
|
|
+ row = i // grid_size
|
|
|
+ col = i % grid_size
|
|
|
+
|
|
|
+ # Extract the desired slice from the template and density based on the selected dimension
|
|
|
+ if iteration_dimension == 0:
|
|
|
+ templateSlice = np.squeeze(templateNifti.get_fdata()[sliceIdx, :, :]).astype(np.uint16)
|
|
|
+ densitySlice = np.squeeze(densityNifti.get_fdata()[sliceIdx, :, :]).astype(np.double)
|
|
|
+ maskSlice = np.squeeze(mask_data[sliceIdx, :, :])
|
|
|
+ elif iteration_dimension == 1:
|
|
|
+ templateSlice = np.squeeze(templateNifti.get_fdata()[:, sliceIdx, :]).astype(np.uint16)
|
|
|
+ densitySlice = np.squeeze(densityNifti.get_fdata()[:, sliceIdx, :]).astype(np.double)
|
|
|
+ maskSlice = np.squeeze(mask_data[:, sliceIdx, :])
|
|
|
+ else:
|
|
|
+ templateSlice = np.squeeze(templateNifti.get_fdata()[:, :, sliceIdx]).astype(np.uint16)
|
|
|
+ densitySlice = np.squeeze(densityNifti.get_fdata()[:, :, sliceIdx]).astype(np.double)
|
|
|
+ maskSlice = np.squeeze(mask_data[:, :, sliceIdx])
|
|
|
|
|
|
# Normalize the templateSlice to [0, 1] for visualization
|
|
|
templateSlice = (templateSlice - np.min(templateSlice)) / (np.max(templateSlice) - np.min(templateSlice))
|
|
@@ -43,25 +87,18 @@ for sliceIdx in range(num_slices_dim2):
|
|
|
# Overlay the density and mask onto the template slice
|
|
|
overlayedImage = templateSlice + densitySlice + maskSlice
|
|
|
|
|
|
- # Normalize the overlay to [0, 1]
|
|
|
- overlayedImage = (overlayedImage - np.min(overlayedImage)) / (np.max(overlayedImage) - np.min(overlayedImage))
|
|
|
-
|
|
|
- # Convert the overlay image to RGBA format (4 channels for transparency)
|
|
|
- overlay_image_rgba = plt.cm.gray(overlayedImage)
|
|
|
- overlay_image_rgba[:, :, 3] = (overlayedImage * 255).astype(np.uint8) # Set alpha channel based on intensity
|
|
|
-
|
|
|
- # Update the 3D overlay volume
|
|
|
- overlay_volume[:, :, sliceIdx, :] = overlay_image_rgba
|
|
|
-
|
|
|
-# Create a figure for the 3D visualization
|
|
|
-fig = plt.figure(figsize=(10, 10), dpi=dpi)
|
|
|
-ax = fig.add_subplot(111, projection='3d')
|
|
|
+ # Display the overlaid image in the corresponding subplot
|
|
|
+ axs[row, col].imshow(templateSlice, cmap='gray')
|
|
|
+ axs[row, col].imshow(densitySlice, cmap='hot', alpha=0.4)
|
|
|
+ axs[row, col].imshow(maskSlice, cmap='hot', alpha=0.2) # Apply "hot" colormap to the mask
|
|
|
+ # Apply "cool" colormap to the mask
|
|
|
|
|
|
-# Plot the 3D overlay volume
|
|
|
-ax.voxels(overlay_volume[:, :, :, 0] > 0, facecolors=overlay_volume[:, :, :, :3] / 255., edgecolor='k')
|
|
|
+# Remove empty subplots if necessary
|
|
|
+for i in range(num_slices // slice_step, grid_size * grid_size):
|
|
|
+ fig.delaxes(axs.flatten()[i])
|
|
|
|
|
|
-# Set the aspect ratio to be equal for all dimensions
|
|
|
-ax.set_box_aspect([1, 1, 1])
|
|
|
+# Adjust spacing between subplots
|
|
|
+plt.tight_layout()
|
|
|
|
|
|
-# Show the 3D plot
|
|
|
+# Show the plot
|
|
|
plt.show()
|