<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1, minimum-scale=1" />
<meta name="generator" content="pdoc 0.10.0" />
<title>pymskt.statistics.main API documentation</title>
<meta name="description" content="" />
<link rel="preload stylesheet" as="style" href="https://cdnjs.cloudflare.com/ajax/libs/10up-sanitize.css/11.0.1/sanitize.min.css" integrity="sha256-PK9q560IAAa6WVRRh76LtCaI8pjTJ2z11v0miyNNjrs=" crossorigin>
<link rel="preload stylesheet" as="style" href="https://cdnjs.cloudflare.com/ajax/libs/10up-sanitize.css/11.0.1/typography.min.css" integrity="sha256-7l/o7C8jubJiy74VsKTidCy1yBkRtiUGbVkYBylBqUg=" crossorigin>
<link rel="stylesheet preload" as="style" href="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/10.1.1/styles/github.min.css" crossorigin>
<style>:root{--highlight-color:#fe9}.flex{display:flex !important}body{line-height:1.5em}#content{padding:20px}#sidebar{padding:30px;overflow:hidden}#sidebar > *:last-child{margin-bottom:2cm}.http-server-breadcrumbs{font-size:130%;margin:0 0 15px 0}#footer{font-size:.75em;padding:5px 30px;border-top:1px solid #ddd;text-align:right}#footer p{margin:0 0 0 1em;display:inline-block}#footer p:last-child{margin-right:30px}h1,h2,h3,h4,h5{font-weight:300}h1{font-size:2.5em;line-height:1.1em}h2{font-size:1.75em;margin:1em 0 .50em 0}h3{font-size:1.4em;margin:25px 0 10px 0}h4{margin:0;font-size:105%}h1:target,h2:target,h3:target,h4:target,h5:target,h6:target{background:var(--highlight-color);padding:.2em 0}a{color:#058;text-decoration:none;transition:color .3s ease-in-out}a:hover{color:#e82}.title code{font-weight:bold}h2[id^="header-"]{margin-top:2em}.ident{color:#900}pre code{background:#f8f8f8;font-size:.8em;line-height:1.4em}code{background:#f2f2f1;padding:1px 4px;overflow-wrap:break-word}h1 code{background:transparent}pre{background:#f8f8f8;border:0;border-top:1px solid #ccc;border-bottom:1px solid #ccc;margin:1em 0;padding:1ex}#http-server-module-list{display:flex;flex-flow:column}#http-server-module-list div{display:flex}#http-server-module-list dt{min-width:10%}#http-server-module-list p{margin-top:0}.toc ul,#index{list-style-type:none;margin:0;padding:0}#index code{background:transparent}#index h3{border-bottom:1px solid #ddd}#index ul{padding:0}#index h4{margin-top:.6em;font-weight:bold}@media (min-width:200ex){#index .two-column{column-count:2}}@media (min-width:300ex){#index .two-column{column-count:3}}dl{margin-bottom:2em}dl dl:last-child{margin-bottom:4em}dd{margin:0 0 1em 3em}#header-classes + dl > dd{margin-bottom:3em}dd dd{margin-left:2em}dd p{margin:10px 0}.name{background:#eee;font-weight:bold;font-size:.85em;padding:5px 10px;display:inline-block;min-width:40%}.name:hover{background:#e0e0e0}dt:target .name{background:var(--highlight-color)}.name > span:first-child{white-space:nowrap}.name.class > span:nth-child(2){margin-left:.4em}.inherited{color:#999;border-left:5px solid #eee;padding-left:1em}.inheritance em{font-style:normal;font-weight:bold}.desc h2{font-weight:400;font-size:1.25em}.desc h3{font-size:1em}.desc dt code{background:inherit}.source summary,.git-link-div{color:#666;text-align:right;font-weight:400;font-size:.8em;text-transform:uppercase}.source summary > *{white-space:nowrap;cursor:pointer}.git-link{color:inherit;margin-left:1em}.source pre{max-height:500px;overflow:auto;margin:0}.source pre code{font-size:12px;overflow:visible}.hlist{list-style:none}.hlist li{display:inline}.hlist li:after{content:',\2002'}.hlist li:last-child:after{content:none}.hlist .hlist{display:inline;padding-left:1em}img{max-width:100%}td{padding:0 .5em}.admonition{padding:.1em .5em;margin-bottom:1em}.admonition-title{font-weight:bold}.admonition.note,.admonition.info,.admonition.important{background:#aef}.admonition.todo,.admonition.versionadded,.admonition.tip,.admonition.hint{background:#dfd}.admonition.warning,.admonition.versionchanged,.admonition.deprecated{background:#fd4}.admonition.error,.admonition.danger,.admonition.caution{background:lightpink}</style>
<style media="screen and (min-width: 700px)">@media screen and (min-width:700px){#sidebar{width:30%;height:100vh;overflow:auto;position:sticky;top:0}#content{width:70%;max-width:100ch;padding:3em 4em;border-left:1px solid #ddd}pre code{font-size:1em}.item .name{font-size:1em}main{display:flex;flex-direction:row-reverse;justify-content:flex-end}.toc ul ul,#index ul{padding-left:1.5em}.toc > ul > li{margin-top:.5em}}</style>
<style media="print">@media print{#sidebar h1{page-break-before:always}.source{display:none}}@media print{*{background:transparent !important;color:#000 !important;box-shadow:none !important;text-shadow:none !important}a[href]:after{content:" (" attr(href) ")";font-size:90%}a[href][title]:after{content:none}abbr[title]:after{content:" (" attr(title) ")"}.ir a:after,a[href^="javascript:"]:after,a[href^="#"]:after{content:""}pre,blockquote{border:1px solid #999;page-break-inside:avoid}thead{display:table-header-group}tr,img{page-break-inside:avoid}img{max-width:100% !important}@page{margin:0.5cm}p,h2,h3{orphans:3;widows:3}h1,h2,h3,h4,h5,h6{page-break-after:avoid}}</style>
<script defer src="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/10.1.1/highlight.min.js" integrity="sha256-Uv3H6lx7dJmRfRvH8TH6kJD1TSK1aFcwgx+mdg3epi8=" crossorigin></script>
<script>window.addEventListener('DOMContentLoaded', () => hljs.initHighlighting())</script>
</head>
<body>
<main>
<article id="content">
<header>
<h1 class="title">Module <code>pymskt.statistics.main</code></h1>
</header>
<section id="section-intro">
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">from re import sub
import numpy as np
from datetime import date
import os
import pyvista as pv
import pyacvd
from pymskt.mesh.meshRegistration import get_icp_transform, non_rigidly_register
from pymskt.mesh.meshTools import get_mesh_physical_point_coords, set_mesh_physical_point_coords, resample_surface
from pymskt.mesh.meshTransform import apply_transform
from pymskt.mesh.utils import get_symmetric_surface_distance, vtk_deep_copy
from pymskt.mesh import io
today = date.today()
class FindReferenceMeshICP:
"""
For list of meshes perform all possible ICP registrations to identify mesh with smallest
surface error to all other meshes.
Parameters
----------
list_meshes : _type_
_description_
"""
def __init__(
self,
list_mesh_paths,
max_n_iter=1000,
n_landmarks=1000,
reg_mode='similarity',
verbose=True
):
"""
Perform ICP registration between all pairs of meshes. Calculate
symmetric surface distance for all registered meshes. Find target
mesh with smallest mean surface error to all other meshes.
This smallest error mesh is the refrence mesh for the next step of
SSM pipelines (procrustes using non-rigid registration)
Parameters
----------
list_mesh_paths : _type_
_description_
max_n_iter : int, optional
_description_, by default 1000
n_landmarks : int, optional
_description_, by default 1000
reg_mode : str, optional
_description_, by default 'similarity'
verbose : bool, optional
_description_, by default True
"""
self.list_mesh_paths = list_mesh_paths
self.n_meshes = len(list_mesh_paths)
self._symm_surface_distances = np.zeros((self.n_meshes, self.n_meshes), dtype=float)
self._mean_errors = None
self.max_n_iter = max_n_iter
self.n_landmarks = n_landmarks
self.reg_mode = reg_mode
self.verbose=verbose
self._ref_idx = None
self._ref_path = None
def register_meshes(self, idx1_target, idx2_source):
target = io.read_vtk(self.list_mesh_paths[idx1_target])
source = io.read_vtk(self.list_mesh_paths[idx2_source])
icp_transform = get_icp_transform(
source,
target,
max_n_iter=self.max_n_iter,
n_landmarks=self.n_landmarks,
reg_mode=self.reg_mode
)
transformed_source = apply_transform(source, icp_transform)
symmetric_surf_distance = get_symmetric_surface_distance(target, transformed_source)
self._symm_surface_distances[idx1_target, idx2_source] = symmetric_surf_distance
def get_template_idx(self):
self._mean_errors = np.mean(self._symm_surface_distances, axis=1)
self._ref_idx = np.argmin(self._mean_errors)
self._ref_path = self.list_mesh_paths[self._ref_idx]
def execute(self):
if self.verbose is True:
print(f'Starting registrations, there are {len(self.list_mesh_paths)} meshes')
for idx1_target, target_path in enumerate(self.list_mesh_paths):
if self.verbose is True:
print(f'\tStarting target mesh {idx1_target}')
for idx2_source, source_path in enumerate(self.list_mesh_paths):
if self.verbose is True:
print(f'\t\tStarting source mesh {idx2_source}')
# If the target & mesh are same skip, errors = 0
if idx1_target == idx2_source:
continue
else:
self.register_meshes(idx1_target, idx2_source)
if self.verbose is True:
print('Finished all registrations!')
self.get_template_idx()
@property
def ref_idx(self):
return self._ref_idx
@property
def ref_path(self):
return self._ref_path
@property
def symm_surface_distances(self):
return self._symm_surface_distances
@property
def mean_errors(self):
return self._mean_errors
class ProcrustesRegistration:
# https://en.wikipedia.org/wiki/Generalized_Procrustes_analysis
def __init__(
self,
path_ref_mesh,
list_mesh_paths,
tolerance1=2e-1,
tolerance2=1e-2,
max_n_registration_steps=10,
verbose=True,
remesh_each_step=False,
patience=2,
ref_mesh_eigenmap_as_reference=True,
registering_secondary_bone=False, # True if registering secondary bone of joint, after
# primary already used for initial registration. E.g.,
# already did femur for knee, now applying to tibia/patella
**kwargs
):
self.path_ref_mesh = path_ref_mesh
self.list_mesh_paths = list_mesh_paths
# Ensure that path_ref_mesh is in list & at index 0
if self.path_ref_mesh in self.list_mesh_paths:
path_ref_idx = self.list_mesh_paths.index(self.path_ref_mesh)
self.list_mesh_paths.pop(path_ref_idx)
self.list_mesh_paths.insert(0, self.path_ref_mesh)
self._ref_mesh = io.read_vtk(self.path_ref_mesh)
self.n_points = self._ref_mesh.GetNumberOfPoints()
self.ref_mesh_eigenmap_as_reference = ref_mesh_eigenmap_as_reference
self.mean_mesh = None
self.tolerance1 = tolerance1
self.tolerance2 = tolerance2
self.max_n_registration_steps = max_n_registration_steps
self.kwargs = kwargs
# ORIGINALLY THIS WAS THE LOGIC:
# Ensure that the source mesh (mean, or reference) is the base mesh
# We want all meshes aligned with this reference. Then we want
# to apply a "warp" of the ref/mean mesh to make it
# EXCETION - if we are registering a secondary bone in a joint model
# E.g., for registering tibia/patella in knee model.
self.kwargs['icp_register_first'] = True
if registering_secondary_bone is False:
self.kwargs['icp_reg_target_to_source'] = True
elif registering_secondary_bone is True:
self.kwargs['icp_reg_target_to_source'] = False
self._registered_pt_coords = np.zeros((len(list_mesh_paths), self.n_points, 3), dtype=float)
self._registered_pt_coords[0, :, :] = get_mesh_physical_point_coords(self._ref_mesh)
self.sym_error = 100
self.list_errors = []
self.list_ref_meshes = []
self.reg_idx = 0
self.patience = patience
self.patience_idx = 0
self._best_score = 100
self.verbose = verbose
self.remesh_each_step = remesh_each_step
self.error_2_error_change = 100
def register(self, ref_mesh_source, other_mesh_idx):
target_mesh = io.read_vtk(self.list_mesh_paths[other_mesh_idx])
registered_mesh = non_rigidly_register(
target_mesh=target_mesh,
source_mesh=ref_mesh_source,
target_eigenmap_as_reference=not self.ref_mesh_eigenmap_as_reference,
**self.kwargs
)
return get_mesh_physical_point_coords(registered_mesh)
def execute(self):
# create placeholder to store registered point clouds & update inherited one only if also storing
registered_pt_coords = np.zeros_like(self._registered_pt_coords)
# keep doing registrations until max# is hit, or the minimum error between registrations is hit.
while (
(self.reg_idx < self.max_n_registration_steps) &
(self.sym_error > self.tolerance1) &
(self.error_2_error_change > self.tolerance2)
):
if self.verbose is True:
print(f'Starting registration round {self.reg_idx}')
# If its not the very first iteration - check whether or not we want to re-mesh after every iteration.
if (self.reg_idx != 0) & (self.remesh_each_step is True):
n_points = self._ref_mesh.GetNumberOfPoints()
self._ref_mesh = resample_surface(self._ref_mesh, subdivisions=2, clusters=n_points)
if n_points != self.n_points:
print(f'Updating n_points for mesh from {self.n_points} to {self._ref_mesh.GetNumberOfPoints()}')
# re-create the array to store registered points as the # vertices might change after re-meshing.
# also update n_points.
self.n_points = n_points
registered_pt_coords = np.zeros((len(self.list_mesh_paths), self.n_points, 3), dtype=float)
# register the reference mesh to all other meshes
for idx, path in enumerate(self.list_mesh_paths):
if self.verbose is True:
print(f'\tRegistering to mesh # {idx}')
# skip the first mesh in the list if its the first round (its the reference)
if (self.reg_idx == 0) & (idx == 0):
# first iteration & ref mesh, just use points as they are.
registered_pt_coords[idx, :, :] = get_mesh_physical_point_coords(self._ref_mesh)
continue
# register & save registered coordinates in the pre-allocated array
registered_pt_coords[idx, :, :] = self.register(self._ref_mesh, idx)
# Calculate the mean bone shape & create new mean bone shape mesh
mean_shape = np.mean(registered_pt_coords, axis=0)
mean_mesh = vtk_deep_copy(self._ref_mesh)
set_mesh_physical_point_coords(mean_mesh, mean_shape)
# store in list of reference meshes
self.list_ref_meshes.append(mean_mesh)
# Get surface distance between previous reference mesh and the new mean
sym_error = get_symmetric_surface_distance(self._ref_mesh, mean_mesh)
self.error_2_error_change = np.abs(sym_error - self.sym_error)
self.sym_error = sym_error
if self.sym_error >= self._best_score:
# if the error isnt going down, then keep track of that and done save the
# new reference mesh.
self.patience_idx += 1
else:
self.patience_idx = 0
# ONLY UPDATE THE REF_MESH or the REGISTERED_PTS WHEN THE INTER-MESH (REF) ERROR GETS BETTER
# NOT SURE IF THIS IS A GOOD IDEA - MIGHT WANT A BETTER CRITERIA?
self._ref_mesh = mean_mesh
self._registered_pt_coords = registered_pt_coords
# Store the symmetric error values so they can be plotted later
self.list_errors.append(self.sym_error)
if self.verbose is True:
print(f'\t\tSymmetric surface error: {self.sym_error}')
if self.patience_idx >= self.patience:
print(f'Early stopping initiated - no improvment for {self.patience_idx} iterations, patience is: {self.patience}')
break
self.reg_idx += 1
def save_meshes(
self,
mesh_suffix=f'procrustes_registered_{today.strftime("%b")}_{today.day}_{today.year}',
folder=None
):
if folder is not None:
os.makedirs(folder, exist_ok=True)
mesh = vtk_deep_copy(self._ref_mesh)
for idx, path in enumerate(self.list_mesh_paths):
# parse folder / filename for saving
orig_folder = os.path.dirname(path)
orig_filename = os.path.basename(path)
base_filename = orig_filename[: orig_filename.rfind(".")]
filename = f'{base_filename}_{mesh_suffix}_{idx}.vtk'
if folder is None:
path_to_save = os.path.join(orig_folder, filename)
else:
path_to_save = os.path.join(os.path.abspath(folder), filename)
# Keep recycling the same base mesh, just move the x/y/z point coords around.
set_mesh_physical_point_coords(mesh, self._registered_pt_coords[idx, :, :])
# save mesh to disk
io.write_vtk(mesh, path_to_save)
def save_ref_mesh(self, path):
io.write_vtk(self._ref_mesh, path)
@property
def ref_mesh(self):
return self._ref_mesh
@property
def registered_pt_coords(self):
return self._registered_pt_coords
</code></pre>
</details>
</section>
<section>
</section>
<section>
</section>
<section>
</section>
<section>
<h2 class="section-title" id="header-classes">Classes</h2>
<dl>
<dt id="pymskt.statistics.main.FindReferenceMeshICP"><code class="flex name class">
<span>class <span class="ident">FindReferenceMeshICP</span></span>
<span>(</span><span>list_mesh_paths, max_n_iter=1000, n_landmarks=1000, reg_mode='similarity', verbose=True)</span>
</code></dt>
<dd>
<div class="desc"><p>For list of meshes perform all possible ICP registrations to identify mesh with smallest
surface error to all other meshes. </p>
<h2 id="parameters">Parameters</h2>
<dl>
<dt><strong><code>list_meshes</code></strong> : <code>_type_</code></dt>
<dd><em>description</em></dd>
</dl>
<p>Perform ICP registration between all pairs of meshes. Calculate
symmetric surface distance for all registered meshes. Find target
mesh with smallest mean surface error to all other meshes. </p>
<p>This smallest error mesh is the refrence mesh for the next step of
SSM pipelines (procrustes using non-rigid registration)</p>
<h2 id="parameters_1">Parameters</h2>
<dl>
<dt><strong><code>list_mesh_paths</code></strong> : <code>_type_</code></dt>
<dd><em>description</em></dd>
<dt><strong><code>max_n_iter</code></strong> : <code>int</code>, optional</dt>
<dd><em>description</em>, by default 1000</dd>
<dt><strong><code>n_landmarks</code></strong> : <code>int</code>, optional</dt>
<dd><em>description</em>, by default 1000</dd>
<dt><strong><code>reg_mode</code></strong> : <code>str</code>, optional</dt>
<dd><em>description</em>, by default 'similarity'</dd>
<dt><strong><code>verbose</code></strong> : <code>bool</code>, optional</dt>
<dd><em>description</em>, by default True</dd>
</dl></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">class FindReferenceMeshICP:
"""
For list of meshes perform all possible ICP registrations to identify mesh with smallest
surface error to all other meshes.
Parameters
----------
list_meshes : _type_
_description_
"""
def __init__(
self,
list_mesh_paths,
max_n_iter=1000,
n_landmarks=1000,
reg_mode='similarity',
verbose=True
):
"""
Perform ICP registration between all pairs of meshes. Calculate
symmetric surface distance for all registered meshes. Find target
mesh with smallest mean surface error to all other meshes.
This smallest error mesh is the refrence mesh for the next step of
SSM pipelines (procrustes using non-rigid registration)
Parameters
----------
list_mesh_paths : _type_
_description_
max_n_iter : int, optional
_description_, by default 1000
n_landmarks : int, optional
_description_, by default 1000
reg_mode : str, optional
_description_, by default 'similarity'
verbose : bool, optional
_description_, by default True
"""
self.list_mesh_paths = list_mesh_paths
self.n_meshes = len(list_mesh_paths)
self._symm_surface_distances = np.zeros((self.n_meshes, self.n_meshes), dtype=float)
self._mean_errors = None
self.max_n_iter = max_n_iter
self.n_landmarks = n_landmarks
self.reg_mode = reg_mode
self.verbose=verbose
self._ref_idx = None
self._ref_path = None
def register_meshes(self, idx1_target, idx2_source):
target = io.read_vtk(self.list_mesh_paths[idx1_target])
source = io.read_vtk(self.list_mesh_paths[idx2_source])
icp_transform = get_icp_transform(
source,
target,
max_n_iter=self.max_n_iter,
n_landmarks=self.n_landmarks,
reg_mode=self.reg_mode
)
transformed_source = apply_transform(source, icp_transform)
symmetric_surf_distance = get_symmetric_surface_distance(target, transformed_source)
self._symm_surface_distances[idx1_target, idx2_source] = symmetric_surf_distance
def get_template_idx(self):
self._mean_errors = np.mean(self._symm_surface_distances, axis=1)
self._ref_idx = np.argmin(self._mean_errors)
self._ref_path = self.list_mesh_paths[self._ref_idx]
def execute(self):
if self.verbose is True:
print(f'Starting registrations, there are {len(self.list_mesh_paths)} meshes')
for idx1_target, target_path in enumerate(self.list_mesh_paths):
if self.verbose is True:
print(f'\tStarting target mesh {idx1_target}')
for idx2_source, source_path in enumerate(self.list_mesh_paths):
if self.verbose is True:
print(f'\t\tStarting source mesh {idx2_source}')
# If the target & mesh are same skip, errors = 0
if idx1_target == idx2_source:
continue
else:
self.register_meshes(idx1_target, idx2_source)
if self.verbose is True:
print('Finished all registrations!')
self.get_template_idx()
@property
def ref_idx(self):
return self._ref_idx
@property
def ref_path(self):
return self._ref_path
@property
def symm_surface_distances(self):
return self._symm_surface_distances
@property
def mean_errors(self):
return self._mean_errors</code></pre>
</details>
<h3>Instance variables</h3>
<dl>
<dt id="pymskt.statistics.main.FindReferenceMeshICP.mean_errors"><code class="name">var <span class="ident">mean_errors</span></code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">@property
def mean_errors(self):
return self._mean_errors</code></pre>
</details>
</dd>
<dt id="pymskt.statistics.main.FindReferenceMeshICP.ref_idx"><code class="name">var <span class="ident">ref_idx</span></code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">@property
def ref_idx(self):
return self._ref_idx</code></pre>
</details>
</dd>
<dt id="pymskt.statistics.main.FindReferenceMeshICP.ref_path"><code class="name">var <span class="ident">ref_path</span></code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">@property
def ref_path(self):
return self._ref_path</code></pre>
</details>
</dd>
<dt id="pymskt.statistics.main.FindReferenceMeshICP.symm_surface_distances"><code class="name">var <span class="ident">symm_surface_distances</span></code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">@property
def symm_surface_distances(self):
return self._symm_surface_distances</code></pre>
</details>
</dd>
</dl>
<h3>Methods</h3>
<dl>
<dt id="pymskt.statistics.main.FindReferenceMeshICP.execute"><code class="name flex">
<span>def <span class="ident">execute</span></span>(<span>self)</span>
</code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">def execute(self):
if self.verbose is True:
print(f'Starting registrations, there are {len(self.list_mesh_paths)} meshes')
for idx1_target, target_path in enumerate(self.list_mesh_paths):
if self.verbose is True:
print(f'\tStarting target mesh {idx1_target}')
for idx2_source, source_path in enumerate(self.list_mesh_paths):
if self.verbose is True:
print(f'\t\tStarting source mesh {idx2_source}')
# If the target & mesh are same skip, errors = 0
if idx1_target == idx2_source:
continue
else:
self.register_meshes(idx1_target, idx2_source)
if self.verbose is True:
print('Finished all registrations!')
self.get_template_idx()</code></pre>
</details>
</dd>
<dt id="pymskt.statistics.main.FindReferenceMeshICP.get_template_idx"><code class="name flex">
<span>def <span class="ident">get_template_idx</span></span>(<span>self)</span>
</code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">def get_template_idx(self):
self._mean_errors = np.mean(self._symm_surface_distances, axis=1)
self._ref_idx = np.argmin(self._mean_errors)
self._ref_path = self.list_mesh_paths[self._ref_idx]</code></pre>
</details>
</dd>
<dt id="pymskt.statistics.main.FindReferenceMeshICP.register_meshes"><code class="name flex">
<span>def <span class="ident">register_meshes</span></span>(<span>self, idx1_target, idx2_source)</span>
</code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">def register_meshes(self, idx1_target, idx2_source):
target = io.read_vtk(self.list_mesh_paths[idx1_target])
source = io.read_vtk(self.list_mesh_paths[idx2_source])
icp_transform = get_icp_transform(
source,
target,
max_n_iter=self.max_n_iter,
n_landmarks=self.n_landmarks,
reg_mode=self.reg_mode
)
transformed_source = apply_transform(source, icp_transform)
symmetric_surf_distance = get_symmetric_surface_distance(target, transformed_source)
self._symm_surface_distances[idx1_target, idx2_source] = symmetric_surf_distance</code></pre>
</details>
</dd>
</dl>
</dd>
<dt id="pymskt.statistics.main.ProcrustesRegistration"><code class="flex name class">
<span>class <span class="ident">ProcrustesRegistration</span></span>
<span>(</span><span>path_ref_mesh, list_mesh_paths, tolerance1=0.2, tolerance2=0.01, max_n_registration_steps=10, verbose=True, remesh_each_step=False, patience=2, ref_mesh_eigenmap_as_reference=True, registering_secondary_bone=False, **kwargs)</span>
</code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">class ProcrustesRegistration:
# https://en.wikipedia.org/wiki/Generalized_Procrustes_analysis
def __init__(
self,
path_ref_mesh,
list_mesh_paths,
tolerance1=2e-1,
tolerance2=1e-2,
max_n_registration_steps=10,
verbose=True,
remesh_each_step=False,
patience=2,
ref_mesh_eigenmap_as_reference=True,
registering_secondary_bone=False, # True if registering secondary bone of joint, after
# primary already used for initial registration. E.g.,
# already did femur for knee, now applying to tibia/patella
**kwargs
):
self.path_ref_mesh = path_ref_mesh
self.list_mesh_paths = list_mesh_paths
# Ensure that path_ref_mesh is in list & at index 0
if self.path_ref_mesh in self.list_mesh_paths:
path_ref_idx = self.list_mesh_paths.index(self.path_ref_mesh)
self.list_mesh_paths.pop(path_ref_idx)
self.list_mesh_paths.insert(0, self.path_ref_mesh)
self._ref_mesh = io.read_vtk(self.path_ref_mesh)
self.n_points = self._ref_mesh.GetNumberOfPoints()
self.ref_mesh_eigenmap_as_reference = ref_mesh_eigenmap_as_reference
self.mean_mesh = None
self.tolerance1 = tolerance1
self.tolerance2 = tolerance2
self.max_n_registration_steps = max_n_registration_steps
self.kwargs = kwargs
# ORIGINALLY THIS WAS THE LOGIC:
# Ensure that the source mesh (mean, or reference) is the base mesh
# We want all meshes aligned with this reference. Then we want
# to apply a "warp" of the ref/mean mesh to make it
# EXCETION - if we are registering a secondary bone in a joint model
# E.g., for registering tibia/patella in knee model.
self.kwargs['icp_register_first'] = True
if registering_secondary_bone is False:
self.kwargs['icp_reg_target_to_source'] = True
elif registering_secondary_bone is True:
self.kwargs['icp_reg_target_to_source'] = False
self._registered_pt_coords = np.zeros((len(list_mesh_paths), self.n_points, 3), dtype=float)
self._registered_pt_coords[0, :, :] = get_mesh_physical_point_coords(self._ref_mesh)
self.sym_error = 100
self.list_errors = []
self.list_ref_meshes = []
self.reg_idx = 0
self.patience = patience
self.patience_idx = 0
self._best_score = 100
self.verbose = verbose
self.remesh_each_step = remesh_each_step
self.error_2_error_change = 100
def register(self, ref_mesh_source, other_mesh_idx):
target_mesh = io.read_vtk(self.list_mesh_paths[other_mesh_idx])
registered_mesh = non_rigidly_register(
target_mesh=target_mesh,
source_mesh=ref_mesh_source,
target_eigenmap_as_reference=not self.ref_mesh_eigenmap_as_reference,
**self.kwargs
)
return get_mesh_physical_point_coords(registered_mesh)
def execute(self):
# create placeholder to store registered point clouds & update inherited one only if also storing
registered_pt_coords = np.zeros_like(self._registered_pt_coords)
# keep doing registrations until max# is hit, or the minimum error between registrations is hit.
while (
(self.reg_idx < self.max_n_registration_steps) &
(self.sym_error > self.tolerance1) &
(self.error_2_error_change > self.tolerance2)
):
if self.verbose is True:
print(f'Starting registration round {self.reg_idx}')
# If its not the very first iteration - check whether or not we want to re-mesh after every iteration.
if (self.reg_idx != 0) & (self.remesh_each_step is True):
n_points = self._ref_mesh.GetNumberOfPoints()
self._ref_mesh = resample_surface(self._ref_mesh, subdivisions=2, clusters=n_points)
if n_points != self.n_points:
print(f'Updating n_points for mesh from {self.n_points} to {self._ref_mesh.GetNumberOfPoints()}')
# re-create the array to store registered points as the # vertices might change after re-meshing.
# also update n_points.
self.n_points = n_points
registered_pt_coords = np.zeros((len(self.list_mesh_paths), self.n_points, 3), dtype=float)
# register the reference mesh to all other meshes
for idx, path in enumerate(self.list_mesh_paths):
if self.verbose is True:
print(f'\tRegistering to mesh # {idx}')
# skip the first mesh in the list if its the first round (its the reference)
if (self.reg_idx == 0) & (idx == 0):
# first iteration & ref mesh, just use points as they are.
registered_pt_coords[idx, :, :] = get_mesh_physical_point_coords(self._ref_mesh)
continue
# register & save registered coordinates in the pre-allocated array
registered_pt_coords[idx, :, :] = self.register(self._ref_mesh, idx)
# Calculate the mean bone shape & create new mean bone shape mesh
mean_shape = np.mean(registered_pt_coords, axis=0)
mean_mesh = vtk_deep_copy(self._ref_mesh)
set_mesh_physical_point_coords(mean_mesh, mean_shape)
# store in list of reference meshes
self.list_ref_meshes.append(mean_mesh)
# Get surface distance between previous reference mesh and the new mean
sym_error = get_symmetric_surface_distance(self._ref_mesh, mean_mesh)
self.error_2_error_change = np.abs(sym_error - self.sym_error)
self.sym_error = sym_error
if self.sym_error >= self._best_score:
# if the error isnt going down, then keep track of that and done save the
# new reference mesh.
self.patience_idx += 1
else:
self.patience_idx = 0
# ONLY UPDATE THE REF_MESH or the REGISTERED_PTS WHEN THE INTER-MESH (REF) ERROR GETS BETTER
# NOT SURE IF THIS IS A GOOD IDEA - MIGHT WANT A BETTER CRITERIA?
self._ref_mesh = mean_mesh
self._registered_pt_coords = registered_pt_coords
# Store the symmetric error values so they can be plotted later
self.list_errors.append(self.sym_error)
if self.verbose is True:
print(f'\t\tSymmetric surface error: {self.sym_error}')
if self.patience_idx >= self.patience:
print(f'Early stopping initiated - no improvment for {self.patience_idx} iterations, patience is: {self.patience}')
break
self.reg_idx += 1
def save_meshes(
self,
mesh_suffix=f'procrustes_registered_{today.strftime("%b")}_{today.day}_{today.year}',
folder=None
):
if folder is not None:
os.makedirs(folder, exist_ok=True)
mesh = vtk_deep_copy(self._ref_mesh)
for idx, path in enumerate(self.list_mesh_paths):
# parse folder / filename for saving
orig_folder = os.path.dirname(path)
orig_filename = os.path.basename(path)
base_filename = orig_filename[: orig_filename.rfind(".")]
filename = f'{base_filename}_{mesh_suffix}_{idx}.vtk'
if folder is None:
path_to_save = os.path.join(orig_folder, filename)
else:
path_to_save = os.path.join(os.path.abspath(folder), filename)
# Keep recycling the same base mesh, just move the x/y/z point coords around.
set_mesh_physical_point_coords(mesh, self._registered_pt_coords[idx, :, :])
# save mesh to disk
io.write_vtk(mesh, path_to_save)
def save_ref_mesh(self, path):
io.write_vtk(self._ref_mesh, path)
@property
def ref_mesh(self):
return self._ref_mesh
@property
def registered_pt_coords(self):
return self._registered_pt_coords</code></pre>
</details>
<h3>Instance variables</h3>
<dl>
<dt id="pymskt.statistics.main.ProcrustesRegistration.ref_mesh"><code class="name">var <span class="ident">ref_mesh</span></code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">@property
def ref_mesh(self):
return self._ref_mesh</code></pre>
</details>
</dd>
<dt id="pymskt.statistics.main.ProcrustesRegistration.registered_pt_coords"><code class="name">var <span class="ident">registered_pt_coords</span></code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">@property
def registered_pt_coords(self):
return self._registered_pt_coords</code></pre>
</details>
</dd>
</dl>
<h3>Methods</h3>
<dl>
<dt id="pymskt.statistics.main.ProcrustesRegistration.execute"><code class="name flex">
<span>def <span class="ident">execute</span></span>(<span>self)</span>
</code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">def execute(self):
# create placeholder to store registered point clouds & update inherited one only if also storing
registered_pt_coords = np.zeros_like(self._registered_pt_coords)
# keep doing registrations until max# is hit, or the minimum error between registrations is hit.
while (
(self.reg_idx < self.max_n_registration_steps) &
(self.sym_error > self.tolerance1) &
(self.error_2_error_change > self.tolerance2)
):
if self.verbose is True:
print(f'Starting registration round {self.reg_idx}')
# If its not the very first iteration - check whether or not we want to re-mesh after every iteration.
if (self.reg_idx != 0) & (self.remesh_each_step is True):
n_points = self._ref_mesh.GetNumberOfPoints()
self._ref_mesh = resample_surface(self._ref_mesh, subdivisions=2, clusters=n_points)
if n_points != self.n_points:
print(f'Updating n_points for mesh from {self.n_points} to {self._ref_mesh.GetNumberOfPoints()}')
# re-create the array to store registered points as the # vertices might change after re-meshing.
# also update n_points.
self.n_points = n_points
registered_pt_coords = np.zeros((len(self.list_mesh_paths), self.n_points, 3), dtype=float)
# register the reference mesh to all other meshes
for idx, path in enumerate(self.list_mesh_paths):
if self.verbose is True:
print(f'\tRegistering to mesh # {idx}')
# skip the first mesh in the list if its the first round (its the reference)
if (self.reg_idx == 0) & (idx == 0):
# first iteration & ref mesh, just use points as they are.
registered_pt_coords[idx, :, :] = get_mesh_physical_point_coords(self._ref_mesh)
continue
# register & save registered coordinates in the pre-allocated array
registered_pt_coords[idx, :, :] = self.register(self._ref_mesh, idx)
# Calculate the mean bone shape & create new mean bone shape mesh
mean_shape = np.mean(registered_pt_coords, axis=0)
mean_mesh = vtk_deep_copy(self._ref_mesh)
set_mesh_physical_point_coords(mean_mesh, mean_shape)
# store in list of reference meshes
self.list_ref_meshes.append(mean_mesh)
# Get surface distance between previous reference mesh and the new mean
sym_error = get_symmetric_surface_distance(self._ref_mesh, mean_mesh)
self.error_2_error_change = np.abs(sym_error - self.sym_error)
self.sym_error = sym_error
if self.sym_error >= self._best_score:
# if the error isnt going down, then keep track of that and done save the
# new reference mesh.
self.patience_idx += 1
else:
self.patience_idx = 0
# ONLY UPDATE THE REF_MESH or the REGISTERED_PTS WHEN THE INTER-MESH (REF) ERROR GETS BETTER
# NOT SURE IF THIS IS A GOOD IDEA - MIGHT WANT A BETTER CRITERIA?
self._ref_mesh = mean_mesh
self._registered_pt_coords = registered_pt_coords
# Store the symmetric error values so they can be plotted later
self.list_errors.append(self.sym_error)
if self.verbose is True:
print(f'\t\tSymmetric surface error: {self.sym_error}')
if self.patience_idx >= self.patience:
print(f'Early stopping initiated - no improvment for {self.patience_idx} iterations, patience is: {self.patience}')
break
self.reg_idx += 1</code></pre>
</details>
</dd>
<dt id="pymskt.statistics.main.ProcrustesRegistration.register"><code class="name flex">
<span>def <span class="ident">register</span></span>(<span>self, ref_mesh_source, other_mesh_idx)</span>
</code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">def register(self, ref_mesh_source, other_mesh_idx):
target_mesh = io.read_vtk(self.list_mesh_paths[other_mesh_idx])
registered_mesh = non_rigidly_register(
target_mesh=target_mesh,
source_mesh=ref_mesh_source,
target_eigenmap_as_reference=not self.ref_mesh_eigenmap_as_reference,
**self.kwargs
)
return get_mesh_physical_point_coords(registered_mesh)</code></pre>
</details>
</dd>
<dt id="pymskt.statistics.main.ProcrustesRegistration.save_meshes"><code class="name flex">
<span>def <span class="ident">save_meshes</span></span>(<span>self, mesh_suffix='procrustes_registered_Jul_31_2022', folder=None)</span>
</code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">def save_meshes(
self,
mesh_suffix=f'procrustes_registered_{today.strftime("%b")}_{today.day}_{today.year}',
folder=None
):
if folder is not None:
os.makedirs(folder, exist_ok=True)
mesh = vtk_deep_copy(self._ref_mesh)
for idx, path in enumerate(self.list_mesh_paths):
# parse folder / filename for saving
orig_folder = os.path.dirname(path)
orig_filename = os.path.basename(path)
base_filename = orig_filename[: orig_filename.rfind(".")]
filename = f'{base_filename}_{mesh_suffix}_{idx}.vtk'
if folder is None:
path_to_save = os.path.join(orig_folder, filename)
else:
path_to_save = os.path.join(os.path.abspath(folder), filename)
# Keep recycling the same base mesh, just move the x/y/z point coords around.
set_mesh_physical_point_coords(mesh, self._registered_pt_coords[idx, :, :])
# save mesh to disk
io.write_vtk(mesh, path_to_save)</code></pre>
</details>
</dd>
<dt id="pymskt.statistics.main.ProcrustesRegistration.save_ref_mesh"><code class="name flex">
<span>def <span class="ident">save_ref_mesh</span></span>(<span>self, path)</span>
</code></dt>
<dd>
<div class="desc"></div>
<details class="source">
<summary>
<span>Expand source code</span>
</summary>
<pre><code class="python">def save_ref_mesh(self, path):
io.write_vtk(self._ref_mesh, path)</code></pre>
</details>
</dd>
</dl>
</dd>
</dl>
</section>
</article>
<nav id="sidebar">
<h1>Index</h1>
<div class="toc">
<ul></ul>
</div>
<ul id="index">
<li><h3>Super-module</h3>
<ul>
<li><code><a title="pymskt.statistics" href="index.html">pymskt.statistics</a></code></li>
</ul>
</li>
<li><h3><a href="#header-classes">Classes</a></h3>
<ul>
<li>
<h4><code><a title="pymskt.statistics.main.FindReferenceMeshICP" href="#pymskt.statistics.main.FindReferenceMeshICP">FindReferenceMeshICP</a></code></h4>
<ul class="">
<li><code><a title="pymskt.statistics.main.FindReferenceMeshICP.execute" href="#pymskt.statistics.main.FindReferenceMeshICP.execute">execute</a></code></li>
<li><code><a title="pymskt.statistics.main.FindReferenceMeshICP.get_template_idx" href="#pymskt.statistics.main.FindReferenceMeshICP.get_template_idx">get_template_idx</a></code></li>
<li><code><a title="pymskt.statistics.main.FindReferenceMeshICP.mean_errors" href="#pymskt.statistics.main.FindReferenceMeshICP.mean_errors">mean_errors</a></code></li>
<li><code><a title="pymskt.statistics.main.FindReferenceMeshICP.ref_idx" href="#pymskt.statistics.main.FindReferenceMeshICP.ref_idx">ref_idx</a></code></li>
<li><code><a title="pymskt.statistics.main.FindReferenceMeshICP.ref_path" href="#pymskt.statistics.main.FindReferenceMeshICP.ref_path">ref_path</a></code></li>
<li><code><a title="pymskt.statistics.main.FindReferenceMeshICP.register_meshes" href="#pymskt.statistics.main.FindReferenceMeshICP.register_meshes">register_meshes</a></code></li>
<li><code><a title="pymskt.statistics.main.FindReferenceMeshICP.symm_surface_distances" href="#pymskt.statistics.main.FindReferenceMeshICP.symm_surface_distances">symm_surface_distances</a></code></li>
</ul>
</li>
<li>
<h4><code><a title="pymskt.statistics.main.ProcrustesRegistration" href="#pymskt.statistics.main.ProcrustesRegistration">ProcrustesRegistration</a></code></h4>
<ul class="">
<li><code><a title="pymskt.statistics.main.ProcrustesRegistration.execute" href="#pymskt.statistics.main.ProcrustesRegistration.execute">execute</a></code></li>
<li><code><a title="pymskt.statistics.main.ProcrustesRegistration.ref_mesh" href="#pymskt.statistics.main.ProcrustesRegistration.ref_mesh">ref_mesh</a></code></li>
<li><code><a title="pymskt.statistics.main.ProcrustesRegistration.register" href="#pymskt.statistics.main.ProcrustesRegistration.register">register</a></code></li>
<li><code><a title="pymskt.statistics.main.ProcrustesRegistration.registered_pt_coords" href="#pymskt.statistics.main.ProcrustesRegistration.registered_pt_coords">registered_pt_coords</a></code></li>
<li><code><a title="pymskt.statistics.main.ProcrustesRegistration.save_meshes" href="#pymskt.statistics.main.ProcrustesRegistration.save_meshes">save_meshes</a></code></li>
<li><code><a title="pymskt.statistics.main.ProcrustesRegistration.save_ref_mesh" href="#pymskt.statistics.main.ProcrustesRegistration.save_ref_mesh">save_ref_mesh</a></code></li>
</ul>
</li>
</ul>
</li>
</ul>
</nav>
</main>
<footer id="footer">
<p>Generated by <a href="https://pdoc3.github.io/pdoc" title="pdoc: Python API documentation generator"><cite>pdoc</cite> 0.10.0</a>.</p>
</footer>
</body>
</html>
Datasets
Models
