Statistical shape modelling with pyssam

This statistical shape modelling example visualises modes of left lower lung lobe shape correlations.

import pyssam

from copy import copy
from glob import glob
import matplotlib.pyplot as plt
import numpy as np

First, we source landmark data to use in our shape model

LANDMARK_DIR = "../../example_data/lung_landmarks/"
landmark_files = glob(LANDMARK_DIR + "/landmarks*.csv")
if len(landmark_files) == 0:
  raise AssertionError(
    "The directories you have declared are empty.",
    "\nPlease check your input arguments.",
  )

landmark_coordinates = np.array(
  [np.loadtxt(l, delimiter=",") for l in landmark_files]
)

Initialising the model

Here we convert to landmark coordinates into a parameterised shape model. We first initialise the class, which handles all pre-processing. Then, we can compute the shape model components and mean population shape

ssm_obj = pyssam.SSM(landmark_coordinates)
ssm_obj.create_pca_model(ssm_obj.landmarks_columns_scale)
mean_shape_columnvector = ssm_obj.compute_dataset_mean()
mean_shape = mean_shape_columnvector.reshape(-1, 3)
shape_model_components = ssm_obj.pca_model_components

Plotting and analysis

# Define some plotting functions

def plot_cumulative_variance(explained_variance, target_variance=-1):
    number_of_components = np.arange(0, len(explained_variance))+1
    fig, ax = plt.subplots(1,1)
    color = "blue"
    ax.plot(number_of_components, explained_variance*100.0, marker="o", ms=2, color=color, mec=color, mfc=color)
    if target_variance > 0.0:
        ax.axhline(target_variance*100.0)

    ax.set_ylabel("Variance [%]")
    ax.set_xlabel("Number of components")
    ax.grid(axis="x")
    plt.show()

def plot_shape_modes(
  mean_shape_columnvector,
  mean_shape,
  original_shape_parameter_vector,
  shape_model_components,
  mode_to_plot,
):
  weights = [-2, 0, 2]
  fig, ax = plt.subplots(1, 3)
  for j, weights_i in enumerate(weights):
    shape_parameter_vector = copy(original_shape_parameter_vector)
    shape_parameter_vector[mode_to_plot] = weights_i
    mode_i_coords = ssm_obj.morph_model(
        mean_shape_columnvector,
        shape_model_components,
        shape_parameter_vector
    ).reshape(-1, 3)

    offset_dist = pyssam.utils.euclidean_distance(
      mean_shape,
      mode_i_coords
    )
    # colour points blue if closer to point cloud centre than mean shape
    mean_shape_dist_from_centre = pyssam.utils.euclidean_distance(
      mean_shape,
      np.zeros(3),
    )
    mode_i_dist_from_centre = pyssam.utils.euclidean_distance(
      mode_i_coords,
      np.zeros(3),
    )
    offset_dist = np.where(
        mode_i_dist_from_centre<mean_shape_dist_from_centre,
        offset_dist*-1,
        offset_dist,
    )
    if weights_i == 0:
      ax[j].scatter(
        mode_i_coords[:, 0],
        mode_i_coords[:, 2],
        c="gray",
        s=1,
      )
      ax[j].set_title("mean shape")
    else:
      ax[j].scatter(
        mode_i_coords[:, 0],
        mode_i_coords[:, 2],
        c=offset_dist,
        cmap="seismic",
        vmin=-1,
        vmax=1,
        s=1,
      )
      ax[j].set_title(f"mode {mode_to_plot} \nweight {weights_i}")
    ax[j].axis('off')
    ax[j].margins(0,0)
    ax[j].xaxis.set_major_locator(plt.NullLocator())
    ax[j].yaxis.set_major_locator(plt.NullLocator())

  plt.show()

Generally, the first part in assessing the model once trained is visualising how the explained variance changes with an increasing number of modes. If few modes are required, it means it will be much simpler to fit the SSM to an image.

print(f"To obtain {ssm_obj.desired_variance*100}% variance, {ssm_obj.required_mode_number} modes are required")
plot_cumulative_variance(np.cumsum(ssm_obj.pca_object.explained_variance_ratio_), 0.9)

To obtain 90.0% variance, 24 modes are required

Now the interesting part. We visualise the first three principal components, where the points are coloured by their displacement.

mode_to_plot = 0
print(f"explained variance is {ssm_obj.pca_object.explained_variance_ratio_[mode_to_plot]}")

plot_shape_modes(
    mean_shape_columnvector,
    mean_shape,
    ssm_obj.model_parameters,
    ssm_obj.pca_model_components,
    mode_to_plot,
)

explained variance is 0.2862721047372064

mode_to_plot = 1
print(f"explained variance is {ssm_obj.pca_object.explained_variance_ratio_[mode_to_plot]}")

plot_shape_modes(
    mean_shape_columnvector,
    mean_shape,
    ssm_obj.model_parameters,
    ssm_obj.pca_model_components,
    mode_to_plot,
)

explained variance is 0.11562461259525639

mode_to_plot = 2
print(f"explained variance is {ssm_obj.pca_object.explained_variance_ratio_[mode_to_plot]}")

plot_shape_modes(
    mean_shape_columnvector,
    mean_shape,
    ssm_obj.model_parameters,
    ssm_obj.pca_model_components,
    mode_to_plot,
)

explained variance is 0.07577465318882828

When we look at the 15th mode, which accounts for around 1% for the total variance, we can see a very minor change in the lung structure (shown by very light red and blue patches).

mode_to_plot = 15
print(f"explained variance is {ssm_obj.pca_object.explained_variance_ratio_[mode_to_plot]}")

plot_shape_modes(
    mean_shape_columnvector,
    mean_shape,
    ssm_obj.model_parameters,
    ssm_obj.pca_model_components,
    mode_to_plot,
)

explained variance is 0.012525399253420878

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/Users/josh.williams/gitrepos/pyssam_base/docs/tutorial/ssm_example.ipynb Cell 19 line 4
      <a href='vscode-notebook-cell:/Users/josh.williams/gitrepos/pyssam_base/docs/tutorial/ssm_example.ipynb#X24sZmlsZQ%3D%3D?line=0'>1</a> mode_to_plot = 15
      <a href='vscode-notebook-cell:/Users/josh.williams/gitrepos/pyssam_base/docs/tutorial/ssm_example.ipynb#X24sZmlsZQ%3D%3D?line=1'>2</a> print(f"explained variance is {ssm_obj.pca_object.explained_variance_ratio_[mode_to_plot]}")
----> <a href='vscode-notebook-cell:/Users/josh.williams/gitrepos/pyssam_base/docs/tutorial/ssm_example.ipynb#X24sZmlsZQ%3D%3D?line=3'>4</a> plot_shape_modes(
      <a href='vscode-notebook-cell:/Users/josh.williams/gitrepos/pyssam_base/docs/tutorial/ssm_example.ipynb#X24sZmlsZQ%3D%3D?line=4'>5</a>     mean_shape_columnvector,
      <a href='vscode-notebook-cell:/Users/josh.williams/gitrepos/pyssam_base/docs/tutorial/ssm_example.ipynb#X24sZmlsZQ%3D%3D?line=5'>6</a>     mean_shape,
      <a href='vscode-notebook-cell:/Users/josh.williams/gitrepos/pyssam_base/docs/tutorial/ssm_example.ipynb#X24sZmlsZQ%3D%3D?line=6'>7</a>     ssm_obj.model_parameters,
      <a href='vscode-notebook-cell:/Users/josh.williams/gitrepos/pyssam_base/docs/tutorial/ssm_example.ipynb#X24sZmlsZQ%3D%3D?line=7'>8</a>     ssm_obj.pca_model_components,
      <a href='vscode-notebook-cell:/Users/josh.williams/gitrepos/pyssam_base/docs/tutorial/ssm_example.ipynb#X24sZmlsZQ%3D%3D?line=8'>9</a>     # mode_to_plot,
     <a href='vscode-notebook-cell:/Users/josh.williams/gitrepos/pyssam_base/docs/tutorial/ssm_example.ipynb#X24sZmlsZQ%3D%3D?line=9'>10</a> )

TypeError: plot_shape_modes() missing 1 required positional argument: 'mode_to_plot'