diff --git a/utilities/fiber_generation/DOCUMENTATION.md b/utilities/fiber_generation/DOCUMENTATION.md new file mode 100644 index 000000000..550c0b129 --- /dev/null +++ b/utilities/fiber_generation/DOCUMENTATION.md @@ -0,0 +1,542 @@ +# Cardiac Fiber Generation Documentation + +This document describes the mathematical framework for generating myocardial fiber orientations in cardiac geometries using Laplace-Dirichlet rule-based methods. + +# General Overview of the Fiber Generation Methods +To explain the main steps and functions of the code, in the following we explain the main concepts and methods to generate fibers. In the next section these concepts are used to describe the Bayer and Doste method. + +## 1. Laplace Problem + +The foundation of the rule-based fiber generation is the solution of Laplace-Dirichlet boundary value problems. These provide scalar fields $\phi$ (surface1 → surface2) that satisfy: + +$$ +\begin{align} +\nabla^2 \phi = 0 &\quad \text{in } \Omega \\ +\phi_t = 0 & \quad\text{on } S_{\text{surface1}} \\ +\phi_t = 1 & \quad\text{on } S_{\text{surface2}} +\end{align} +$$ + +### 1.1 Transmural Direction + +The transmural field $\phi_t$ (endo → epi) characterizes the wall thickness direction, varying from endocardium to epicardium: + +$$ +\begin{cases} +\nabla^2 \phi_t = 0 & \text{in } \Omega \\ +\phi_t = 0 & \text{on } S_{\text{endo}} \\ +\phi_t = 1 & \text{on } S_{\text{epi}} +\end{cases} +$$ + +This field is normalized to $[0, 1]$ range where $\phi_t = 0$ at the endocardium and $\phi_t = 1$ at the epicardium. + +### 1.2 Longitudinal Direction + +The longitudinal field $\phi_\ell$ (apex → base) characterizes the apex-to-base direction: + +$$ +\begin{cases} +\nabla^2 \phi_\ell = 0 & \text{in } \Omega \\ +\phi_\ell = 0 & \text{on } S_{\text{apex}} \\ +\phi_\ell = 1 & \text{on } S_{\text{base}} +\end{cases} +$$ + +This field is also normalized to $[0, 1]$ where $\phi_\ell = 0$ at the apex and $\phi_\ell = 1$ at the base. + +**Implementation Note**: The Laplace equations are solved using the SVMultiphysics solver configured to solve a steady-state heat equation (which is equivalent to the Laplace equation). The solver is configured with: +- `Conductivity = 1.0`, `Source_term = 0.0`, `Density = 0.0` +- `Spectral_radius_of_infinite_time_step = 0.0` +- Single time step to obtain the steady-state solution directly + +## 2. Definition of Basis + +A local orthonormal basis $\{\mathbf{e}_c, \mathbf{e}_\ell, \mathbf{e}_t\}$ is constructed at each point in the myocardium, where: +- $\mathbf{e}_c$: circumferential direction +- $\mathbf{e}_\ell$: longitudinal direction +- $\mathbf{e}_t$: transmural direction + +### 2.1 Obtain Gradients from Laplace Solutions + +The gradients of the Laplace fields provide natural directional vectors. The gradients are computed at mesh nodes and then averaged to cell centers for smoother results: + +$$ +\mathbf{g}_t = \nabla \phi_t, \quad \mathbf{g}_\ell = \nabla \phi_\ell +$$ + +These gradients are normalized to unit vectors: + +$$ +\hat{\mathbf{g}}_t = \frac{\mathbf{g}_t}{\|\mathbf{g}_t\|}, \quad \hat{\mathbf{g}}_\ell = \frac{\mathbf{g}_\ell}{\|\mathbf{g}_\ell\|} +$$ + +### 2.2 Calculate Circumferential Direction + +The local orthonormal basis is constructed using the `axis` function. Given the longitudinal direction $\hat{\mathbf{g}}_\ell$ and the transmural direction $\hat{\mathbf{g}}_t$: + +1. **Longitudinal basis vector**: + $$\mathbf{e}_\ell = \frac{\hat{\mathbf{g}}_\ell}{\|\hat{\mathbf{g}}_\ell\|}$$ + +2. **Transmural basis vector** (orthogonalized to $\mathbf{e}_\ell$): + $$\mathbf{e}_t' = \hat{\mathbf{g}}_t - (\hat{\mathbf{g}}_t \cdot \mathbf{e}_\ell)\mathbf{e}_\ell$$ + $$\mathbf{e}_t = \frac{\mathbf{e}_t'}{\|\mathbf{e}_t'\|}$$ + +3. **Circumferential basis vector** (orthogonal to both): + $$\mathbf{e}_c = \mathbf{e}_\ell \times \mathbf{e}_t$$ + +This ensures a right-handed orthonormal coordinate system at each element. We implement this in the **axis** function, $\mathbf Q=(\mathbf{e}_c, \mathbf{e}_\ell, \mathbf{e}_t)=\text{axis}(\hat{\mathbf{g}}_\ell, \hat{\mathbf{g}}_t)$. + +## 3. Definition of Angles over the Geometry + +Two angles define the fiber orientation relative to the local basis $\mathbf Q=(\mathbf e_c, \mathbf e_\ell, \mathbf e_t)$: + +- **$\alpha$ (helix angle)**: Rotation angle using as axis of rotation $\mathbf{e}_t$. Returns a rotated basis $\mathbf Q^{(\alpha)}=(\mathbf e_c', \mathbf e_\ell', \mathbf e_t)$ +- **$\beta$ (transverse angle)**: Rotation angle using as axis of rotation $\mathbf{e}_\ell'$. Returns the fiber ($\mathbf f$), sheetnormal ($\mathbf n$), and sheet ($\mathbf s$) vectors $\mathbf Q^{(\alpha, \beta)} = (\mathbf f, \mathbf n, \mathbf s)$. + +The angles vary linearly across the wall thickness based on the transmural coordinate $\phi_t$. For the single ventricle this is: + +$$ +\alpha(\phi_t) = \alpha_{\text{endo}}(1 - \phi_t) + \alpha_{\text{epi}}\phi_t +$$ + +$$ +\beta(\phi_t) = \beta_{\text{endo}}(1 - \phi_t) + \beta_{\text{epi}}\phi_t +$$ + +where: +- $\alpha_{\text{endo}}, \alpha_{\text{epi}}$: helix angles at endocardium and epicardium +- $\beta_{\text{endo}}, \beta_{\text{epi}}$: transverse angles at endocardium and epicardium + + +## 4. Rotation of the Basis + +The fiber direction is obtained by applying two successive rotations to the local basis. + +### 4.1 Rotation Using Matrices (Bayer) + +The two-step rotation can be expressed with standard rotation matrices applied to the local basis $\mathbf{Q} = [\mathbf{e}_c, \mathbf{e}_\ell, \mathbf{e}_t]$: + +- Helix rotation by $\alpha$ about the transmural axis $\mathbf{e}_t$: + +$$ +\mathbf{R}_\alpha = \begin{bmatrix} +\cos\alpha & -\sin\alpha & 0 \\ +\sin\alpha & \cos\alpha & 0 \\ +0 & 0 & 1 +\end{bmatrix}, \quad +\mathbf{Q}^{(\alpha)} = \mathbf{Q}\,\mathbf{R}_\alpha. +$$ + +- Transverse rotation by $\beta$ about the rotated longitudinal axis $\mathbf{e}_\ell^{(\alpha)}$ (the second column of $\mathbf{Q}^{(\alpha)}$): + +$$ +\mathbf{R}_\beta = \begin{bmatrix} +\cos\beta & 0 & \sin\beta \\ +0 & 1 & 0 \\ +-\sin\beta & 0 & \cos\beta +\end{bmatrix}, \quad +\mathbf{Q}^{(\alpha,\beta)} = \mathbf{Q}^{(\alpha)}\,\mathbf{R}_\beta. +$$ + +We implement this in the **orient\_matrix** function, $\mathbf Q^{(\alpha, \beta)}=(\mathbf{f}, \mathbf{n}, \mathbf{s})=\text{orient\_matrix}(\mathbf Q, \alpha, \beta)$. + +Note: In the original Bayer paper the second matrix was written as, +$$ +\mathbf{R}_\beta = \begin{bmatrix} +1 & 0 & 0 \\ +0 & \cos\beta & \sin\beta \\ +0 & -\sin\beta & \cos\beta +\end{bmatrix}. +$$ +Given how the orthogonal basis are ordered $(\mathbf{e}_c', \mathbf{e}_\ell', \mathbf{e}_t)$, this is equivalent to rotate over the first vector $\mathbf{e}_c'$ which is not what we want and keeps the fiber orientation unchanged and unaffected by $\beta$ angles. + + +### 4.2 Rotation Using Rodrigues' Formula (Doste) + +These rotations can also be achieved using the Rodrigues rotation formula. For a unit axis $\mathbf{n}$ and angle $\theta$: + +$$ +\mathbf{R}(\theta, \mathbf{n}) = \mathbf{I}\cos\theta + [\mathbf{n}]_\times \sin\theta + \mathbf{n}\,\mathbf{n}^T (1 - \cos\theta), +$$ + +where the skew-symmetric matrix $[\mathbf{n}]_\times$ is + +$$ +[\mathbf{n}]_\times = \begin{bmatrix} +0 & -n_z & n_y \\ +n_z & 0 & -n_x \\ +-n_y & n_x & 0 +\end{bmatrix}, \quad \mathbf{n} = (n_x, n_y, n_z)^T. +$$ + +Applying the two-step rotation to the local basis: + +- Helix rotation by $\alpha$ about $\mathbf{e}_t$: + +$$ +\mathbf{Q}^{(\alpha)} = \mathbf{Q}\,\mathbf{R}(\alpha, \mathbf{e}_t). +$$ + +- Transverse rotation by $\beta$ about the rotated longitudinal axis $\mathbf{e}_\ell^{(\alpha)}$: + +$$ +\mathbf{Q}^{(\alpha,\beta)} = \mathbf{Q}^{(\alpha)}\,\mathbf{R}(\beta, \mathbf{e}_\ell^{(\alpha)}), \quad \mathbf{e}_\ell^{(\alpha)} = \big(\mathbf{Q}^{(\alpha)}\big)[:, 1]. +$$ + +We implement this in the **orient\_rodrigues** function, $\mathbf Q^{(\alpha, \beta)}=(\mathbf{f}, \mathbf{n}, \mathbf{s})=\text{orient\_rodrigues}(\mathbf Q, \alpha, \beta)$. + +## 5. Basis Interpolation + +When working with biventricular geometries, different orthogonal basis are computed for the left ventricle (LV) and right ventricle (RV). For this, the orthogonal basis are represented as quaternions, which are then interpolated using spherical linear interpolation (SLERP). + +### Spherical Linear Interpolation (SLERP) + +Simple linear interpolation of rotation matrices can produce non-orthogonal results. Instead, **interpolate_basis** (bilinear spherical interpolation) is used, which operates via quaternion representation: + +Given two rotation matrices $\mathbf{Q}_1$ and $\mathbf{Q}_2$, and interpolation parameter $t \in [0, 1]$: + +1. **Convert to quaternions**: + $$\mathbf{q}_1 = \text{rotm2quat}(\mathbf{Q}_1), \quad \mathbf{q}_2 = \text{rotm2quat}(\mathbf{Q}_2)$$ + +2. **Ensure shortest path** (quaternion double cover): + $$\text{if } \mathbf{q}_1 \cdot \mathbf{q}_2 < 0: \quad \mathbf{q}_2 \leftarrow -\mathbf{q}_2$$ + +3. **SLERP formula**: + $$\theta_0 = \arccos(\mathbf{q}_1 \cdot \mathbf{q}_2)$$ + $$\mathbf{q}(t) = \frac{\sin((1-t)\theta_0)}{\sin\theta_0}\mathbf{q}_1 + \frac{\sin(t\theta_0)}{\sin\theta_0}\mathbf{q}_2$$ + +4. **Convert back to rotation matrix**: + $$\mathbf{Q}(t) = \text{quat2rotm}(\mathbf{q}(t))$$ + +For nearly parallel quaternions ($\sin\theta_0 < 10^{-6}$), linear interpolation is used instead to avoid numerical issues. + +We implement this in the **interpolate_basis** function, $\mathbf{Q}(t) = \text{interpolate\_basis}(\mathbf{Q}_1, \mathbf{Q}_2, t)$. + +# Bayer Method + +The Bayer et al. (2012) method is designed for truncated biventricular geometries without outflow tracts. + +## Required Boundaries + +The fiber generation algorithm requires specific boundary surfaces to be defined on the cardiac mesh: + +- **Epicardium**: The outer surface of the biventricle +- **LV Endocardium**: The inner surface of the left ventricle +- **RV Endocardium**: The inner surface of the right ventricle +- **Base**: The basal (top) boundary of the geometry +- **Apex**: The epicardial apex region + +The apex surface is automatically generated from the epicardium by identifying the point furthest from the base. Specifically, the apex point $\mathbf{p}_{\text{apex}}$ is found as: + +$$ +\mathbf{p}_{\text{apex}} = \arg\min_{\mathbf{p} \in S_{\text{epi}} \setminus S_{\text{base}}} \|\mathbf{p} - \mathbf{c}_{\text{base}}\| +$$ + +where $S_{\text{epi}}$ is the epicardial surface, $S_{\text{base}}$ is the base surface, and $\mathbf{c}_{\text{base}}$ is the centroid of the base. The apex surface consists of all elements in the epicardium that contain this apex point. + +### Required Laplace Fields + +Four Laplace problems are solved: + +1. **Epi transmural**: $\phi_{\text{epi}}$ (LV endo and RV endo → epi) +2. **LV transmural**: $\phi_{\text{LV}}$ (RV endo and epi → LV endo) +2. **RV transmural**: $\phi_{\text{RV}}$ (LV endo and epi → RV endo) +4. **Apex-to-base**: $\phi_{\text{AB}}$ (apex → base) + +## Input Angles + +The Bayer method requires four input angle parameters (typically specified in degrees): + +- **$\alpha_{\text{endo}}$**: Endocardial helix angle, typically $60°$ +- **$\alpha_{\text{epi}}$**: Epicardial helix angle, typically $-60°$ +- **$\beta_{\text{endo}}$**: Endocardial transverse angle, typically $-20°$ +- **$\beta_{\text{epi}}$**: Epicardial transverse angle, typically $20°$ + +These angles define the fiber architecture that varies smoothly from endocardium to epicardium across both ventricles. + +### Angle Definition + +The angles vary based on the interventricular coordinate $d$ and transmural coordinate $\phi_{\text{epi}}$: + +$$ +d = \frac{\phi_{\text{RV}}}{\phi_{\text{LV}} + \phi_{\text{RV}}} +$$ + +**Septum angles** (interpolated between LV and RV): +$$ +\alpha_s = \alpha_{\text{endo}}(1 - d) - \alpha_{\text{endo}} d +$$ +$$ +\beta_s = \beta_{\text{endo}}(1 - d) - \beta_{\text{endo}} d +$$ + +**Wall angles** (transmural variation): +$$ +\alpha_w = \alpha_{\text{endo}}(1 - \phi_{\text{epi}}) + \alpha_{\text{epi}}\phi_{\text{epi}} +$$ +$$ +\beta_w = \beta_{\text{endo}}(1 - \phi_{\text{epi}}) + \beta_{\text{epi}}\phi_{\text{epi}} +$$ + +### Algorithm Steps (as in Bayer's paper) + +1. **Construct LV and RV basis**: + - $\mathbf{Q}_{\text{LV}}^0 = \text{axis}(\nabla\phi_{\text{AB}}, -\nabla\phi_{\text{LV}})$ + - $\mathbf{Q}_{\text{RV}}^0 = \text{axis}(\nabla\phi_{\text{AB}}, \nabla\phi_{\text{RV}})$ + +2. **Rotate by septum angles**: + - $\mathbf{Q}_{\text{LV}} = \text{orient\_matrix}(\mathbf{Q}_{\text{LV}}^0, \alpha_s, \beta_s)$ + - $\mathbf{Q}_{\text{RV}} = \text{orient\_matrix}(\mathbf{Q}_{\text{RV}}^0, \alpha_s, \beta_s)$ + +3. **Interpolate endocardial basis**: + - $\mathbf{Q}_{\text{endo}} = \text{bislerp}(\mathbf{Q}_{\text{LV}}, \mathbf{Q}_{\text{RV}}, d)$ + +4. **Construct epicardial basis**: + - $\mathbf{Q}_{\text{epi}}^0 = \text{axis}(\nabla\phi_{\text{AB}}, \nabla\phi_{\text{epi}})$ + +5. **Rotate by wall angles**: + - $\mathbf{Q}_{\text{epi}} = \text{orient\_matrix}(\mathbf{Q}_{\text{epi}}^0, \alpha_w, \beta_w)$ + +6. **Interpolate final basis**: + - $\mathbf{Q} = \text{bislerp}(\mathbf{Q}_{\text{endo}}, \mathbf{Q}_{\text{epi}}, \phi_{\text{epi}})$ + - Extract: $\mathbf{f} = \mathbf{Q}[:, 0]$, $\mathbf{n} = \mathbf{Q}[:, 1]$, $\mathbf{s} = \mathbf{Q}[:, 2]$ + + +Note. The bislerp function performs the same SLERP operation as interpolate_basis, but includes an additional correction to account for the fact that fiber directions are equivalent up to sign; that is, for the physical applications considered, using $\mathbf{f}$ or $-\mathbf{f}$ is equivalent. This correction flips each vector of the basis $\mathbf{Q}_1$ and selects the flipped configuration that is closest to the target basis $\mathbf{Q}_2$. However, this approach can introduce issues (see the next section), particularly when $\mathbf{Q}_1$ and $\mathbf{Q}_2$ are nearly orthogonal. In such cases, small perturbations in the basis can lead to drastically different interpolation results. + +### Modified algorithm steps +When running the original implementation, we observed the resulting fibers showed discontinuities that arise due to the **bislerp** function. To solve these issues, we modify the algorithm as follows: + + + - In step 2, we do $\mathbf{Q}_{\text{LV}} = \text{orient\_matrix}(\mathbf{Q}_{\text{LV}}^0, \alpha_s, \text{abs}(\beta_s))$ and $\mathbf{Q}_{\text{RV}} = \text{orient\_matrix}(\mathbf{Q}_{\text{LV}}^0, \alpha_s, \text{abs}(\beta_s))$. + + Note that $\mathbf{Q}_{\mathrm{LV}}^{0}$ and $\mathbf{Q}_{\mathrm{RV}}^{0}$ share equivalent circumferential, longitudinal, and transmural directions \textbf{within the septum}. By definition, $\beta_s > 0$ on the LV side (assuming $\beta_{\mathrm{endo}} > 0$), causing the fiber vector to rotate outward from the septum. On the RV side, $\beta_s < 0$ which also causes the fiber vector to rotate away from the septum. However, this is a negative angle at the RV endocardium, which is not what we want. Taking the absolute value $|\beta_s|$ yields the correct fiber angles while preserving the transmural variation of $\beta_{\mathrm{endo}}$ (positive at both side of the septum and 0 at the center of the septum). + + + ![Illustration of beta angle effect at endocardium](example/biv_truncated/betaendo.png) + + - After step 3, for elements where $d > 0.5$, flip the first and third basis vectors (fiber and sheet) of $\mathbf{Q}_{\text{endo}}$. + + Note that $\mathbf{Q}_{\text{LV}}^{0}$ and $\mathbf{Q}_{\text{RV}}^{0}$ are constructed with opposite signs for the transmural direction. As a result, the LV basis rotates counterclockwise, whereas the RV basis rotates clockwise. At the septum, the circumferential vectors of both bases point in the same direction, which allows for a straightforward SLERP interpolation along the shortest path to obtain $\mathbf{Q}_{\text{endo}}$. However, on the RV side this construction causes the $\mathbf{Q}_{\text{endo}}$ basis to point exactly opposite to $\mathbf{Q}_{\text{epi}}$, leading to issues with the SLERP interpolation. Flipping the vectors on the RV side resolves this problem and ensures that the second interpolation remains smooth. + + + ![Illustration of beta angle effect at endocardium](example/biv_truncated/flipping.png) + + +# Doste Method + +The Doste et al. (2019) method extends the fiber generation to biventricular geometries with **outflow tracts**, requiring valve surfaces to be explicitly defined. + +## Required Boundaries + +The fiber generation algorithm requires the following boundary surfaces: + +- **Epicardium**: The outer surface of the biventricle +- **LV Endocardium**: The inner surface of the left ventricle +- **RV Endocardium**: The inner surface of the right ventricle +- **Mitral valve**: LV inflow boundary +- **Aortic valve**: LV outflow boundary +- **Tricuspid valve**: RV inflow boundary +- **Pulmonary valve**: RV outflow boundary +- **Apex**: The epicardial apex region +- **Top**: Surface that includes all valves. Only needed if the Apex needs to be generated. + +The apex surface is automatically generated from the epicardium by identifying the point furthest from the base/top boundary. The apex point $\mathbf{p}_{\text{apex}}$ is found as: + +$$ +\mathbf{p}_{\text{apex}} = \arg\min_{\mathbf{p} \in S_{\text{epi}} \setminus S_{\text{top}}} \|\mathbf{p} - \mathbf{c}_{\text{top}}\| +$$ + +where $S_{\text{epi}}$ is the epicardial surface, $S_{\text{top}}$ is the top surface, and $\mathbf{c}_{\text{top}}$ is the centroid of the top boundary. + +### Required Laplace Fields + +Eleven Laplace problems are solved: + +**Interventricular**: +1. $\phi_{\text{BiV}}$: RV endocardium → LV endocardium (interventricular field) +2. $\phi_{\text{Trans}}$: Combined transmural field for septal angle calculation + +**Left Ventricle**: +1. $\phi_{\text{LV,trans}}$: Epicardium → LV endocardium (LV transmural for basis construction) +2. $\phi_{\text{LV,av}}$: Aortic valve → apex (LV longitudinal from AV) +3. $\phi_{\text{LV,mv}}$: Mitral valve → apex (LV longitudinal from MV) +4. $\phi_{\text{LV,weight}}$: Aortic valve → mitral valve (LV valve weight) + +**Right Ventricle**: +1. $\phi_{\text{RV,trans}}$: Epicardium → RV endocardium (RV transmural for basis construction) +2. $\phi_{\text{RV,pv}}$: Pulmonary valve → apex (RV longitudinal from PV) +3. $\phi_{\text{RV,tv}}$: Tricuspid valve → apex (RV longitudinal from TV) +4. $\phi_{\text{RV,weight}}$: Pulmonary valve → tricuspid valve (RV valve weight) + +**Global Transmural**: +1. $\phi_{\text{epi,trans}}$: LV and RV endocardium → epicardium (global transmural for angle interpolation) + +**Note**: All Laplace fields are automatically normalized to the $[0, 1]$ range **except** for $\phi_{\text{BiV}}$ and $\phi_{\text{Trans}}$, which are kept in their original range for proper septal field calculation and basis construction. The $\phi_{\text{BiV}}$ field typically has negative values on the RV side and positive values on the LV side. + +## Input Angles + +The Doste method requires twelve input angle parameters (typically specified in degrees) to handle both ventricles and outflow tracts: + +- **$\alpha_{\text{endo,LV}}$**: LV endocardial helix angle, typically $60°$ +- **$\alpha_{\text{epi,LV}}$**: LV epicardial helix angle, typically $-60°$ +- **$\beta_{\text{endo,LV}}$**: LV endocardial transverse angle, typically $-20°$ +- **$\beta_{\text{epi,LV}}$**: LV epicardial transverse angle, typically $20°$ +- **$\alpha_{\text{endo,RV}}$**: RV endocardial helix angle, typically $90°$ +- **$\alpha_{\text{epi,RV}}$**: RV epicardial helix angle, typically $-25°$ +- **$\beta_{\text{endo,RV}}$**: RV endocardial transverse angle, typically $-20°$ +- **$\beta_{\text{epi,RV}}$**: RV epicardial transverse angle, typically $20°$ +- **$\alpha_{\text{OT,endo,LV}}$**: LV outflow tract endocardial helix angle, typically $90°$ +- **$\alpha_{\text{OT,epi,LV}}$**: LV outflow tract epicardial helix angle, typically $0°$ +- **$\alpha_{\text{OT,endo,RV}}$**: RV outflow tract endocardial helix angle, typically $90°$ +- **$\alpha_{\text{OT,epi,RV}}$**: RV outflow tract epicardial helix angle, typically $0°$ + +The outflow tract angles are blended with ventricular angles using the valve weight functions to create smooth transitions between the ventricular and outflow tract regions. + +### Valve Weight Functions + +To localize the influence of each valve, weight functions are computed and redistributed: + +$$ +w_{\text{LV}} = \text{redistribute}(\phi_{\text{LV,mv}}, q_{\text{up}} = 0.7, q_{\text{low}} = 0.01) +$$ +$$ +w_{\text{RV}} = \text{redistribute}(\phi_{\text{RV,tv}}, q_{\text{up}} = 0.1, q_{\text{low}} = 0.001) +$$ + +The redistribution clips values at quantiles $q_{\text{low}}$ and $q_{\text{up}}$, then renormalizes to $[0, 1]$. + +### Angle Definition + +The angles are region-specific with outflow tract blending: + +**LV wall angles**: +$$ +\alpha_{\text{LV,endo}} = \alpha_{\text{endo,LV}} w_{\text{LV}} + \alpha_{\text{OT,endo,LV}}(1 - w_{\text{LV}}) +$$ +$$ +\alpha_{\text{LV,epi}} = \alpha_{\text{epi,LV}} w_{\text{LV}} + \alpha_{\text{OT,epi,LV}}(1 - w_{\text{LV}}) +$$ +$$ +\alpha_{\text{wall,LV}} = \alpha_{\text{LV,endo}}(1 - \phi_{\text{epi,trans}}) + \alpha_{\text{LV,epi}}\phi_{\text{epi,trans}} +$$ +$$ +\beta_{\text{wall,LV}} = \left[\beta_{\text{endo,LV}}(1 - \phi_{\text{epi,trans}}) + \beta_{\text{epi,LV}}\phi_{\text{epi,trans}}\right] w_{\text{LV}} +$$ + +**RV wall angles**: +$$ +\alpha_{\text{RV,endo}} = \alpha_{\text{endo,RV}} w_{\text{RV}} + \alpha_{\text{OT,endo,RV}}(1 - w_{\text{RV}}) +$$ +$$ +\alpha_{\text{RV,epi}} = \alpha_{\text{epi,RV}} w_{\text{RV}} + \alpha_{\text{OT,epi,RV}}(1 - w_{\text{RV}}) +$$ +$$ +\alpha_{\text{wall,RV}} = \alpha_{\text{RV,endo}}(1 - \phi_{\text{epi,trans}}) + \alpha_{\text{RV,epi}}\phi_{\text{epi,trans}} +$$ +$$ +\beta_{\text{wall,RV}} = \left[\beta_{\text{endo,RV}}(1 - \phi_{\text{epi,trans}}) + \beta_{\text{epi,RV}}\phi_{\text{epi,trans}}\right] w_{\text{RV}} +$$ + +**Septum angles**: Computed by blending LV and RV contributions weighted by septal position. + +First, compute a septal field $s$ that is 1 at both endocardia but assigns 2/3 of the septum to the LV: +$$ +s = \begin{cases} +\phi_{\text{Trans}} / 2 & \text{if } \phi_{\text{Trans}} < 0 \\ +\phi_{\text{Trans}} & \text{otherwise} +\end{cases} +$$ +$$ +s = |s| +$$ + +Then compute septum angles based on which ventricle the element belongs to: +$$ +\alpha_{\text{septum}} = \begin{cases} +\alpha_{\text{LV,endo}} \cdot s & \text{if } \phi_{\text{BiV}} < 0 \\ +\alpha_{\text{RV,endo}} \cdot s & \text{if } \phi_{\text{BiV}} > 0 +\end{cases} +$$ +$$ +\beta_{\text{septum}} = \begin{cases} +\beta_{\text{endo,LV}} \cdot s \cdot w_{\text{LV}} & \text{if } \phi_{\text{BiV}} < 0 \\ +\beta_{\text{endo,RV}} \cdot s \cdot w_{\text{RV}} & \text{if } \phi_{\text{BiV}} > 0 +\end{cases} +$$ + +### Algorithm Steps + +1. **Construct LV and RV longitudinal directions**: + - Blend valve gradients: $\mathbf{g}_{\ell,\text{LV}} = w_{\text{LV}} \nabla\phi_{\text{LV,mv}} + (1 - w_{\text{LV}})\nabla\phi_{\text{LV,av}}$ + - Blend valve gradients: $\mathbf{g}_{\ell,\text{RV}} = w_{\text{RV}} \nabla\phi_{\text{RV,tv}} + (1 - w_{\text{RV}})\nabla\phi_{\text{RV,pv}}$ + +2. **Construct LV and RV basis** using chamber-specific transmural fields: + - $\mathbf{Q}_{\text{LV}} = \text{axis}(-\mathbf{g}_{\ell,\text{LV}}, -\nabla\phi_{\text{LV,trans}})$ + - $\mathbf{Q}_{\text{RV}} = \text{axis}(-\mathbf{g}_{\ell,\text{RV}}, -\nabla\phi_{\text{RV,trans}})$ + + Note: The minus signs ensure the basis vectors follow the FibGen convention (apex → base, endo → epi). + +3. **Rotate by septum angles to create septal bases**: + - $\mathbf{Q}_{\text{LV,septum}} = \text{orient\_rodrigues}(\mathbf{Q}_{\text{LV}}, \alpha_{\text{septum}}, \beta_{\text{septum}})$ + - $\mathbf{Q}_{\text{RV,septum}} = \text{orient\_rodrigues}(\mathbf{Q}_{\text{RV}}, \alpha_{\text{septum}}, \beta_{\text{septum}})$ + +4. **Rotate by wall angles to create wall bases**: + - $\mathbf{Q}_{\text{LV,wall}} = \text{orient\_rodrigues}(\mathbf{Q}_{\text{LV}}, \alpha_{\text{wall,LV}}, \beta_{\text{wall,LV}})$ + - $\mathbf{Q}_{\text{RV,wall}} = \text{orient\_rodrigues}(\mathbf{Q}_{\text{RV}}, \alpha_{\text{wall,RV}}, \beta_{\text{wall,RV}})$ + +5. **Create discontinuous septal basis**: + + Normalize $\phi_{\text{BiV}}$ for interpolation: + $$ + \phi_{\text{BiV}}^{\text{norm}} = \begin{cases} + \phi_{\text{BiV}} / 2 & \text{if } \phi_{\text{BiV}} < 0 \\ + \phi_{\text{BiV}} & \text{otherwise} + \end{cases} + $$ + $$ + \phi_{\text{BiV}}^{\text{interp}} = \frac{\phi_{\text{BiV}}^{\text{norm}} + 1}{2} + $$ + + Then assign septal basis discontinuously: + $$\mathbf{Q}_{\text{sep}} = \begin{cases} + \mathbf{Q}_{\text{LV,septum}} & \text{if } \phi_{\text{BiV}}^{\text{interp}} < 0.5 \\ + \mathbf{Q}_{\text{RV,septum}} & \text{if } \phi_{\text{BiV}}^{\text{interp}} \geq 0.5 + \end{cases}$$ + +6. **Interpolate epicardial basis**: + - $\mathbf{Q}_{\text{epi}} = \text{interpolate\_basis}(\mathbf{Q}_{\text{LV,wall}}, \mathbf{Q}_{\text{RV,wall}}, \phi_{\text{BiV}}^{\text{interp}})$ + +7. **Interpolate from endocardium to epicardium**: + - $\mathbf{Q} = \text{interpolate\_basis}(\mathbf{Q}_{\text{sep}}, \mathbf{Q}_{\text{epi}}, \phi_{\text{epi,trans}})$ + - Extract: $\mathbf{f} = \mathbf{Q}[:, 0]$, $\mathbf{n} = \mathbf{Q}[:, 1]$, $\mathbf{s} = \mathbf{Q}[:, 2]$ + + + + +### Modified algorithm steps +In the original Doste paper, only one transmural field is mentioned ($\phi_{\text{Trans}}$) which is used to define the transmural gradient $\mathbf g_t$. Piersanti et al (2021) adds a second transmural field (see Fig. 4) ($\phi_{\text{BiV}}$) to differentiate between LV and RV. We use these two plus the individual transmural fields ($\phi_{\text{LV,trans}}$) and ($\phi_{\text{RV,trans}}$) which provides a smoother transition between LV and RV. + + +# On the convention of the orthogonal basis and angles + +Different papers use different conventions to define the orthogonal circumferential, longitudinal, and transmural basis (see Table below). + +| Method | Longitudinal | Transmural | sign(alpha) endo/epi | sign(beta) endo/epi | +|--------|--------------|------------|----------------------|---------------------| +| Bayer | apex -> base | endo -> epi | +/- | -/+ | +| Doste | base -> apex | epi -> endo (rv) endo -> epi (lv) | -/+ (lv); +/- (rv) | -/+ | +| Piersanti (Doste) | apex -> base | endo -> epi (rv) epi -> endo (lv) | -/+ (lv); +/- (rv) in text section 2.6; +/- (lv); +/- (rv) in eq. 8 | -/+ | + +For coherency, for all methods and for all chambers, we consider the transmural vector $\mathbf e_t$ pointing outwards (endo to epi), the longitudinal vector $\mathbf e_\ell$ pointint towards the base (apex to base), and the circumferential vector $\mathbf e_c=\mathbf e_\ell \times \mathbf e_t$. This way, when using literature angle values, we consider, + +| Method | Longitudinal | Transmural | sign(alpha) endo/epi | sign(beta) endo/epi | +|--------|--------------|------------|----------------------|---------------------| +| FibGen | apex -> base | endo -> epi | +/- | -/+ | + + +# References +1. Bayer, J. D., Blake, R. C., Plank, G., & Trayanova, N. A. (2012). A Novel Rule-Based Algorithm for Assigning Myocardial Fiber Orientation to Computational Heart Models. Annals of Biomedical Engineering, 40(10), 2243–2254. https://doi.org/10.1007/s10439-012-0593-5 +2. Doste, R., Soto‐Iglesias, D., Bernardino, G., Alcaine, A., Sebastian, R., Giffard‐Roisin, S., Sermesant, M., Berruezo, A., Sanchez‐Quintana, D., & Camara, O. (2019). A rule‐based method to model myocardial fiber orientation in cardiac biventricular geometries with outflow tracts. International Journal for Numerical Methods in Biomedical Engineering, 35(4). https://doi.org/10.1002/cnm.3185 +3. Piersanti, R., Africa, P. C., Fedele, M., Vergara, C., Dedè, L., Corno, A. F., & Quarteroni, A. (2021). Modeling cardiac muscle fibers in ventricular and atrial electrophysiology simulations. Computer Methods in Applied Mechanics and Engineering, 373, 113468. https://doi.org/10.1016/j.cma.2020.113468 \ No newline at end of file diff --git a/utilities/fiber_generation/README.md b/utilities/fiber_generation/README.md new file mode 100644 index 000000000..d2659239d --- /dev/null +++ b/utilities/fiber_generation/README.md @@ -0,0 +1,28 @@ +# svMultiPhysics fiber generation codes +Python + svMultiPhysics codes for fiber generation. Two methods are implemented: +* Bayer et al. (2012). [link](https://doi.org/10.1007/s10439-012-0593-5) +* Doste et al. (2018). [link](https://doi.org/10.1002/cnm.3185) + +## Installation + +### Installing as a Python Package +You can install the `fiber_generation` codes as a package using pip: + +```bash +pip install -e . +``` +This will install all the required packages and will allow you to call the functions in these packages from any directory. + +## Examples +The `main_bayer.py` and `main_doste.py` are scripts to run both methods in the geometry described in the `example/biv_truncated` and `example/biv_with_outflow_tracts` folders respectively. + +Results for truncated BiV (Bayer) +Results for BiV w/ outflow tracts (Doste) + +Note that the Doste methods needs a geometry with outflow tracts to be run (each valve needs to be defined as a separated surface). Bayer can be run in any biventricular geometry. + + +### Documentation +For details ont the implementation and an study of the results, see: +- [DOCUMENTATION.md](DOCUMENTATION.md) - Comprehensive documentation on methods and usage +- [VALIDATION.md](VALIDATION.md) - Validation studies \ No newline at end of file diff --git a/utilities/fiber_generation/VALIDATION.md b/utilities/fiber_generation/VALIDATION.md new file mode 100644 index 000000000..a8cec0ca7 --- /dev/null +++ b/utilities/fiber_generation/VALIDATION.md @@ -0,0 +1,114 @@ +# Results and validation + +This document presents validation results for fiber generation using two methods: the Bayer method and the Doste method. + +The `examples` folder contains two geometries: a truncated biventricle (`examples/truncated`) for the Bayer method and a biventricle with outflow tracts (`examples/ot`) for the Doste method. + +Validation scripts (`validation_bayer.py` and `validation_doste.py`) can be run to generate the correlation plots shown below. ParaView Python scripts (`paraview_bayer.py` and `paraview_doste.py`) can be used to generate the different visualization views (fiber, sheet, sheet-normal orientations) for each case. + +The validation codes first calculate the $\alpha$ and $\beta$ angles using scalar interpolations. Then, the $\alpha$ and $\beta$ angles are calculated from the fiber direction $\mathbf f$ and the orthogonal basis $(\mathbf e_c, \mathbf e_\ell, \mathbf e_t)$. + +We do this for three combinations of parameters: +1. Setting all $\alpha$ to default values, and all $\beta$ to 0. +2. Setting all $\alpha$ to 0, and all $\beta$ to default values. +3. Setting all $\alpha$ and $\beta$ to default values. + +This allows to isolate and identify issues in the $\alpha$ and $\beta$ rotations. + +Ideally both scalar interpolated and fiber derived angles should match exactly, but given that the orthogonal basis $(\mathbf e_c, \mathbf e_\ell, \mathbf e_t)$ also needs to be interpolated some differences arise. + + +Notes: +- To run the `validation*.py` scripts the file with the Laplace results must be created. This can be done using the `main*.py` scripts. +- To run the `paraview*.py` you must run the `validation*.py` codes first. +- `paraview*.py` codes are run within the Paraview GUI. To do so, + 1. Open Paraview, go to the Python Shell (if not visible, go to View, and click so it appears, usually at the bottom panel). + 2. Click on Run Script and select the desire script. + +--- + +## Bayer Method + +The Bayer method results are demonstrated on a truncated biventricular geometry. + +### Fiber Orientation + +![Bayer Fiber Full View](example/biv_truncated/bayer_fiber.png) + +*Figure 1: Fiber orientation field generated using the Bayer method - full view* + +![Bayer Fiber Slice View](example/biv_truncated/bayer_fiber_slice.png) + +*Figure 2: Fiber orientation field generated using the Bayer method - slice view* + +### Sheet Orientation + +![Bayer Sheet Full View](example/biv_truncated/bayer_sheet.png) + +*Figure 3: Sheet orientation field generated using the Bayer method - full view* + +![Bayer Sheet Slice View](example/biv_truncated/bayer_sheet_slice.png) + +*Figure 4: Sheet orientation field generated using the Bayer method - slice view* + +### Sheet-Normal Orientation + +![Bayer Sheet-Normal Full View](example/biv_truncated/bayer_sheet-normal.png) + +*Figure 5: Sheet-normal orientation field generated using the Bayer method - full view* + +![Bayer Sheet-Normal Slice View](example/biv_truncated/bayer_sheet-normal_slice.png) + +*Figure 6: Sheet-normal orientation field generated using the Bayer method - slice view* + +### Angle Correlations + +To check the code, we first calculate the $\alpha$ and $\beta$ angles using scalar interpolations. Then, we calculate the $\alpha$ and $\beta$ angles using the fiber direction $\mathbf f$ and the orthogonal basis $\mathbf e_c$, $\mathbf e_\ell$, $\mathbf e_t$. + +![Bayer Angle Correlations](example/biv_truncated/bayer_angle_correlations.png) + +*Figure 7: Correlation plots comparing scalar interpolation angles with fiber derived angles for the Bayer method. Blue and red dots show the $\alpha$ and $\beta$ angles. For reference, the original Bayer method with no modifications is shown.* + +--- + +## Doste Method + +The Doste method results are demonstrated on a complete biventricular geometry with outflow tracts. + +### Fiber Orientation + +![Doste Fiber Full View](example/biv_with_outflow_tracts/doste_fiber.png) + +*Figure 8: Fiber orientation field generated using the Doste method - full view* + +![Doste Fiber Slice View](example/biv_with_outflow_tracts/doste_fiber_slice.png) + +*Figure 9: Fiber orientation field generated using the Doste method - slice view* + +### Sheet Orientation + +![Doste Sheet Full View](example/biv_with_outflow_tracts/doste_sheet.png) + +*Figure 10: Sheet orientation field generated using the Doste method - full view* + +![Doste Sheet Slice View](example/biv_with_outflow_tracts/doste_sheet_slice.png) + +*Figure 11: Sheet orientation field generated using the Doste method - slice view* + +### Sheet-Normal Orientation + +![Doste Sheet-Normal Full View](example/biv_with_outflow_tracts/doste_sheet-normal.png) + +*Figure 12: Sheet-normal orientation field generated using the Doste method - full view* + +![Doste Sheet-Normal Slice View](example/biv_with_outflow_tracts/doste_sheet-normal_slice.png) + +*Figure 13: Sheet-normal orientation field generated using the Doste method - slice view* + +### Angle Correlations + +![Doste Angle Correlations](example/biv_with_outflow_tracts/doste_angle_correlations.png) + +*Figure 14: Correlation plots comparing scalar interpolation angles with fiber derived angles for the Doste method. Blue and red dots show the $\alpha$ and $\beta$ angles** + +--- \ No newline at end of file diff --git a/utilities/fiber_generation/example/biv_truncated/VOLUME.vtu b/utilities/fiber_generation/example/biv_truncated/VOLUME.vtu new file mode 100644 index 000000000..bec495b6e --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/VOLUME.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bcd011125435852298177578498ebcf2b5aac1080b824163d1facdefc48c2a27 +size 24524913 diff --git a/utilities/fiber_generation/example/biv_truncated/bayer_angle_correlations.png b/utilities/fiber_generation/example/biv_truncated/bayer_angle_correlations.png new file mode 100644 index 000000000..361db6b3a --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/bayer_angle_correlations.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d04384b8f0304f368cd2db3635be78d4f60c717f9fe0c49106d095c0fac8b1f9 +size 167822 diff --git a/utilities/fiber_generation/example/biv_truncated/bayer_fiber.png b/utilities/fiber_generation/example/biv_truncated/bayer_fiber.png new file mode 100644 index 000000000..25d1d20ed --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/bayer_fiber.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8172518297fffd88e9edb10a9cb91abe8279247ea39873b52d306c4aaf86b89c +size 721193 diff --git a/utilities/fiber_generation/example/biv_truncated/bayer_fiber_slice.png b/utilities/fiber_generation/example/biv_truncated/bayer_fiber_slice.png new file mode 100644 index 000000000..1c7c85324 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/bayer_fiber_slice.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:297d9e63397879871fae4a58b36e8be8b0b3c7dae89480dc1276c6ca346ebd7f +size 613428 diff --git a/utilities/fiber_generation/example/biv_truncated/bayer_sheet-normal.png b/utilities/fiber_generation/example/biv_truncated/bayer_sheet-normal.png new file mode 100644 index 000000000..67adce301 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/bayer_sheet-normal.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7c77fc04580afb81e80db91c2b7f9ab3cfa140122cc2603dfd5ce75341a18287 +size 628260 diff --git a/utilities/fiber_generation/example/biv_truncated/bayer_sheet-normal_slice.png b/utilities/fiber_generation/example/biv_truncated/bayer_sheet-normal_slice.png new file mode 100644 index 000000000..b51b84c39 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/bayer_sheet-normal_slice.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:37515ace4857f6b356fa74b3e0624e5e84ae5ed61e68b4dcc222c0c42b8728d3 +size 550082 diff --git a/utilities/fiber_generation/example/biv_truncated/bayer_sheet.png b/utilities/fiber_generation/example/biv_truncated/bayer_sheet.png new file mode 100644 index 000000000..0f9015b15 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/bayer_sheet.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:36b62daabc388c793ce07ca1c91818110215827a7579cca12eba6f7778d5ce94 +size 742911 diff --git a/utilities/fiber_generation/example/biv_truncated/bayer_sheet_slice.png b/utilities/fiber_generation/example/biv_truncated/bayer_sheet_slice.png new file mode 100644 index 000000000..4ceea9423 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/bayer_sheet_slice.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:077bdb97704b0ca90dc6971c0c422fa5cb5a44d0cc3fa8f01092c5c10a57f22f +size 652930 diff --git a/utilities/fiber_generation/example/biv_truncated/betaendo.png b/utilities/fiber_generation/example/biv_truncated/betaendo.png new file mode 100644 index 000000000..7dc749ec5 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/betaendo.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:41c9b81bf0f974e17be22d36e2bbc57d544087e352d6e6c7dede1cbc179e55aa +size 32445 diff --git a/utilities/fiber_generation/example/biv_truncated/flipping.png b/utilities/fiber_generation/example/biv_truncated/flipping.png new file mode 100644 index 000000000..4ad255935 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/flipping.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e4f6dc5aafb8ee191ec68b5817bc13a011103a0f40d89b43b486c145accd9994 +size 58679 diff --git a/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/BASE.vtp b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/BASE.vtp new file mode 100644 index 000000000..6268fb43c --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/BASE.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b56eb592388c9cc2226c6372ffc8e0199952698f84ffa81330672d272503d7aa +size 191256 diff --git a/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/EPI.vtp b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/EPI.vtp new file mode 100644 index 000000000..1aa8cca64 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/EPI.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2c5ecc05750d086d04dfe9902e7a7e200831ed5cb11ffe6339186112d8bac0ce +size 526632 diff --git a/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/EPI_APEX.vtp b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/EPI_APEX.vtp new file mode 100644 index 000000000..dc75156db --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/EPI_APEX.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d251d13ec1f1245cf1328ddbfd5b93db9573f406548c8cae004d591317364fa0 +size 2857 diff --git a/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/EPI_MID.vtp b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/EPI_MID.vtp new file mode 100644 index 000000000..f0b984ad4 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/EPI_MID.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:cc6a73358eff58939c32e73dfd9326b9a53cf25afa409e5c27b7a87ef3ec08ad +size 525056 diff --git a/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/LV.vtp b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/LV.vtp new file mode 100644 index 000000000..9e9dcbb00 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/LV.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e16a121071598c4696ea25623b52595c349bca27823c6d9cfb6e1bdbc7b77767 +size 145544 diff --git a/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/RV.vtp b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/RV.vtp new file mode 100644 index 000000000..69f68ba06 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/RV.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c6028ab47b1af4947b9ee47cd79dfc58105d8e2bf86f1167aed877609023935c +size 128852 diff --git a/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/exterior.vtp b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/exterior.vtp new file mode 100644 index 000000000..4f3247355 --- /dev/null +++ b/utilities/fiber_generation/example/biv_truncated/mesh-surfaces/exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:889f15170752b6310d35a52f9e1b207ab84afa1e3ce7c4be24c114403015fe3f +size 978206 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_angle_correlations.png b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_angle_correlations.png new file mode 100644 index 000000000..47a1b1d54 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_angle_correlations.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:143986be8c8e18cda784c6438a4a617673d911862bb3f38cbaecf56fa632a263 +size 114366 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_fiber.png b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_fiber.png new file mode 100644 index 000000000..5bdab6dc2 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_fiber.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:be3c2e396e75669076dc6bc126fb12b820ec3d68d67548ba330ef92a02581709 +size 645675 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_fiber_slice.png b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_fiber_slice.png new file mode 100644 index 000000000..5e226c2c7 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_fiber_slice.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:eacbea73000e4251085481e9fd8e65eba8e7781ada4af83b80008091851d9db8 +size 272043 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet-normal.png b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet-normal.png new file mode 100644 index 000000000..5bd4876dc --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet-normal.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5a33ae95bf05cc3de06156ee927302a522fe083d0fc2093e9f51e53c1097e3e6 +size 623371 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet-normal_slice.png b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet-normal_slice.png new file mode 100644 index 000000000..e40fd7794 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet-normal_slice.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f2ce1c2cf222e26bba5b8231525f4e0559d52a69e38ca5bf1972b0050b02f2ce +size 226057 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet.png b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet.png new file mode 100644 index 000000000..ee2145889 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:29d149f84120ba45afb92f39f08db5e16e8120c2785517d8d1d8de81babbd0c1 +size 652046 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet_slice.png b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet_slice.png new file mode 100644 index 000000000..6848b006a --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/doste_sheet_slice.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4899c3f7fff98ff9b54da40036042798207d9d4f799867533b0ba2112e05bbc4 +size 308799 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-complete.mesh.vtu b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-complete.mesh.vtu new file mode 100644 index 000000000..0afc681f0 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9347d6e7907816749cc16471130397f264c7721d029626fb1b9f7fdd6ac1e02c +size 28209696 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/av.vtp b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/av.vtp new file mode 100644 index 000000000..4bbb83ab5 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/av.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ee7d4d8eb31b2639f75599c188cb2af9fa4fcc466d27e34b65dc3ec4f4f3633b +size 14489 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/endo_lv.vtp b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/endo_lv.vtp new file mode 100644 index 000000000..5b2740a69 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/endo_lv.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fa09ea16a80da5813b133a22f42b94b7aae48cd414e9c8ce174d46a4d7eb3a19 +size 563261 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/endo_rv.vtp b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/endo_rv.vtp new file mode 100644 index 000000000..8aef06c70 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/endo_rv.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b0c07b7ca1ad1b8c38d826d363f5572e1e6359595f5701a2b3f8f25606e30073 +size 734382 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/epi.vtp b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/epi.vtp new file mode 100644 index 000000000..8336320c8 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/epi.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2168042f6f99283d747e5895df8bb30d443631b5be7e31b51c55e1a7be83fa29 +size 1293881 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/epi_apex.vtp b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/epi_apex.vtp new file mode 100644 index 000000000..1426b1653 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/epi_apex.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0a5f75b0891aad8d88fa48afc4a023877f43a3ae6afa00bcc979f98717f9a70b +size 2891 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/mv.vtp b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/mv.vtp new file mode 100644 index 000000000..4b9b8492b --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/mv.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:712904ededfdd630a736abfc213758155f8dd68a7f1512c2128c2b0305eeecb6 +size 90751 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/pv.vtp b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/pv.vtp new file mode 100644 index 000000000..ba4ef7861 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/pv.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1b99698dc66260cd1c63d4612402f07e98b0a547549e42d896048bf7bdbe2afb +size 12559 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/top.vtp b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/top.vtp new file mode 100644 index 000000000..28e22860f --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:edec0ae0c98ef7db7add744db0340aeba1c925613403ec4bc1e580d0c1454767 +size 87143 diff --git a/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/tv.vtp b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/tv.vtp new file mode 100644 index 000000000..ebf70fd05 --- /dev/null +++ b/utilities/fiber_generation/example/biv_with_outflow_tracts/mesh-surfaces/tv.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b33291fd2b529845d68d4584c2f4f5241c288dda1eb6a5e76f5bc839e21e2f62 +size 19198 diff --git a/utilities/fiber_generation/fiber_generation/__init__.py b/utilities/fiber_generation/fiber_generation/__init__.py new file mode 100644 index 000000000..146096191 --- /dev/null +++ b/utilities/fiber_generation/fiber_generation/__init__.py @@ -0,0 +1,27 @@ +"""Fiber generation package for biventricular heart models. + +This package provides classes and utilities to generate myocardial fiber +orientations for biventricular heart models using Laplace-Dirichlet +rule-based methods. + +Modules: + fiber_generator: Core fiber generation classes (FibGen, FibGenBayer, FibGenDoste) + laplace_solver: Laplace-Dirichlet equation solver + surface_names: Surface name definitions for heart geometry + surface_utils: Utility functions for surface operations + quat_utils: Quaternion utilities for rotation operations +""" + +from fiber_generation.fiber_generator import FibGen, FibGenBayer, FibGenDoste +from fiber_generation.laplace_solver import LaplaceSolver +from fiber_generation.surface_names import SurfaceName + +__all__ = [ + 'FibGen', + 'FibGenBayer', + 'FibGenDoste', + 'LaplaceSolver', + 'SurfaceName', +] + +__version__ = '0.1.0' diff --git a/utilities/fiber_generation/fiber_generation/fiber_generator.py b/utilities/fiber_generation/fiber_generation/fiber_generator.py new file mode 100644 index 000000000..8a30e3846 --- /dev/null +++ b/utilities/fiber_generation/fiber_generation/fiber_generator.py @@ -0,0 +1,807 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- +"""Object-oriented fiber generation module for biventricular heart models. + +This module provides classes to generate myocardial fiber orientations for +biventricular heart models using Laplace-Dirichlet rule-based methods. + +Supports both: + - Bayer et al. (2012): Truncated BiV geometry + - Doste et al. (2019): BiV geometry with outflow tracts + +References: + Bayer et al. 2012: https://doi.org/10.1007/s10439-012-0593-5 + Doste et al. 2019: https://doi.org/10.1002/cnm.3185 +""" + +import os +import numpy as np +import pyvista as pv + +import fiber_generation.quat_utils as qu + + +class FibGen: + """Base class for fiber generation. + + Provides common utilities for computing fiber orientations from + Laplace field solutions. + + Attributes: + mesh: PyVista mesh with cell-centered data. + lap: Dictionary of Laplace solution values at cells. + grad: Dictionary of gradient arrays at cells (N, 3). + """ + + def __init__(self): + """Initialize the FibGen base class.""" + self.mesh = None + self.lap = None + self.grad = None + + def normalize(self, x): + """Normalize each row of an (N, 3) array. + + Zero-length rows remain zero after normalization. + + Args: + x: Array-like of shape (N, 3). + + Returns: + np.ndarray of shape (N, 3) with row-wise normalized vectors. + """ + a = np.asarray(x, dtype=float) + if a.ndim != 2 or a.shape[1] != 3: + raise ValueError("normalize expects an array of shape (N, 3)") + norms = np.linalg.norm(a, axis=1, keepdims=True) + safe_norms = np.where(norms == 0.0, 1.0, norms) + out = a / safe_norms + zero_rows = (norms.squeeze() == 0.0) + if np.any(zero_rows): + out[zero_rows] = 0.0 + return out + + def scale_to_range(self, arr, range=(0.0, 1.0)): + """Scale array to specified range. + + Args: + arr: Input array to scale. + range: Tuple (min, max) specifying the target range. + + Returns: + np.ndarray: Scaled array with values in the specified range. If all values are equal, + returns array filled with the midpoint of the range. + """ + arr = np.asarray(arr, dtype=float) + amin = np.min(arr) + amax = np.max(arr) + if amax > amin: + return range[0] + (arr - amin) * (range[1] - range[0]) / (amax - amin) + else: + return np.ones_like(arr) * ((range[0] + range[1]) / 2.0) + + def compute_gradients(self, mesh, field_names): + """Compute gradients for specified fields at points. + + Args: + mesh: PyVista mesh with point data. + field_names: List of field names to compute gradients for. + + Returns: + PyVista mesh with gradient arrays added to point_data. + """ + for name in field_names: + if name not in mesh.point_data: + raise KeyError(f"Field '{name}' not found in mesh point_data") + + gmesh = mesh.compute_derivative(scalars=name, gradient=True, preference='point') + mesh.point_data[name + "_grad"] = np.asarray(gmesh.point_data["gradient"]) + + return mesh + + def axis(self, gL, gT): + """Construct orthogonal coordinate systems from two gradient fields. + + Creates an orthonormal basis [eC, eL, eT] for each element where: + - eL is aligned with gL (normalized longitudinal) + - eT is orthogonal to eL and in the plane of gT (transmural) + - eC is the cross product of eL and eT (circumferential) + + Args: + gL: Array of shape (N, 3) representing the longitudinal gradient. + gT: Array of shape (N, 3) representing the transmural gradient. + + Returns: + np.ndarray: Array of shape (N, 3, 3) where columns are + [eC (circumferential), eL (longitudinal), eT (transmural)]. + """ + gL = np.asarray(gL, dtype=float) + gT = np.asarray(gT, dtype=float) + ne = gL.shape[0] + + # eL = normalized longitudinal + eL = self.normalize(gL) + + # eT = gT - proj_{eL}(gT), orthogonal to eL + proj = np.sum(eL * gT, axis=1)[:, None] * eL + eT = gT - proj + eT = self.normalize(eT) + + # eC = cross(eL, eT), circumferential + eC = np.cross(eL, eT, axisa=1, axisb=1) + eC = self.normalize(eC) + + # Build basis matrix Q = [eC, eL, eT] + Q = np.zeros((ne, 3, 3), dtype=float) + Q[:, :, 0] = eC + Q[:, :, 1] = eL + Q[:, :, 2] = eT + + return Q + + def calculate_angle(self, trans, endo_value, epi_value): + """Compute angle varying linearly from endo to epi. + + Args: + trans: Transmural coordinate array (N,), values in [0, 1]. + endo_value: Angle value at endocardium (scalar). + epi_value: Angle value at epicardium (scalar). + + Returns: + np.ndarray: Angle values at each point (N,). + """ + return endo_value * (1 - trans) + epi_value * trans + + def orient_matrix(self, Q, alpha, beta): + """Apply alpha and beta rotations to orthogonal matrices. + + Rotates Q by alpha about the z-axis (transmural) and then + by beta about the y-axis (longitudinal direction). + + Args: + Q: Array of shape (N, 3, 3) containing orthogonal matrices. + alpha: Array of shape (N,) with rotation angles (radians) about z-axis. + beta: Array of shape (N,) with rotation angles (radians) about y-axis. + + Returns: + np.ndarray: Array of shape (N, 3, 3) containing rotated matrices. + """ + Q = np.asarray(Q, dtype=float) + ne = Q.shape[0] + + ca = np.cos(alpha) + sa = np.sin(alpha) + cb = np.cos(beta) + sb = np.sin(beta) + + # Rotation about z-axis (Ra) + Ra = np.zeros((ne, 3, 3), dtype=float) + Ra[:, 0, 0] = ca + Ra[:, 0, 1] = -sa + Ra[:, 1, 0] = sa + Ra[:, 1, 1] = ca + Ra[:, 2, 2] = 1.0 + + # Rotation about y-axis (Rb) + Rb = np.zeros((ne, 3, 3), dtype=float) + Rb[:, 0, 0] = cb + Rb[:, 0, 2] = sb + Rb[:, 1, 1] = 1.0 + Rb[:, 2, 0] = -sb + Rb[:, 2, 2] = cb + + # Compose rotations and apply to Q + RaRb = np.einsum('nij,njk->nik', Ra, Rb) + Qt = np.einsum('nij,njk->nik', Q, RaRb) + + return Qt + + def orient_rodrigues(self, Q, alpha, beta): + """Rotate basis using Rodrigues rotation formula (Doste method). + + Applies two successive rotations using Rodrigues formula: + 1. Rotate by alpha about the transmural axis (eT) + 2. Rotate by beta about the rotated longitudinal axis + + Args: + Q: Array of shape (N, 3, 3) containing basis matrices. + Columns are [eC (circumferential), eL (longitudinal), eT (transmural)]. + alpha: Array of shape (N,) with rotation angles (radians) about transmural axis. + beta: Array of shape (N,) with rotation angles (radians) about rotated longitudinal axis. + + Returns: + np.ndarray: Array of shape (N, 3, 3) containing rotated basis matrices. + """ + Q = np.asarray(Q, dtype=float) + alpha = np.asarray(alpha, dtype=float) + beta = np.asarray(beta, dtype=float) + + n = Q.shape[0] + + # Extract basis vectors + eC = Q[:, :, 0] # Circumferential + eL = Q[:, :, 1] # Longitudinal + eT = Q[:, :, 2] # Transmural + + # Normalize basis vectors + eC = self.normalize(eC) + eL = self.normalize(eL) + eT = self.normalize(eT) + + # Vectorized Rodrigues rotation: + # v_rot = v*cos(theta) + (k x v)*sin(theta) + k*(k·v)*(1-cos(theta)) + + def rot(v, k, theta): + v = np.asarray(v, dtype=float) + k = np.asarray(k, dtype=float) + theta = np.asarray(theta, dtype=float) + + if v.ndim != 2 or v.shape[1] != 3: + raise ValueError("orient_rodrigues expects vectors of shape (N, 3)") + if k.shape != v.shape: + raise ValueError("orient_rodrigues expects axes of shape (N, 3) matching vectors") + + # cos(theta) and sin(theta) are broadcasted to (N, 1) + ct = np.cos(theta)[:, None] + st = np.sin(theta)[:, None] + + kv = np.cross(k, v, axis=1) + kdotv = np.einsum('ij,ij->i', k, v)[:, None] + return v * ct + kv * st + k * kdotv * (1.0 - ct) + + # First rotation: alpha about transmural axis eT + eC1 = rot(eC, eT, alpha) + eL1 = rot(eL, eT, alpha) + eT1 = eT # unchanged + + # Second rotation: beta about rotated longitudinal axis eL1 (normalized) + eL1_axis = self.normalize(eL1) + eC2 = rot(eC1, eL1_axis, beta) + eL2 = eL1 # unchanged (rotation about itself) + eT2 = rot(eT1, eL1_axis, beta) + + result = np.zeros((n, 3, 3), dtype=float) + result[:, :, 0] = eC2 + result[:, :, 1] = eL2 + result[:, :, 2] = eT2 + + return result + + + def interpolate_basis(self, Q1, Q2, t, correct_slerp=False): + """Spherical linear interpolation between batches of rotation matrices. + + Performs SLERP on rotation matrices represented as quaternions internally. + + Args: + Q1: Array of shape (N, 3, 3) containing starting rotation matrices. + Q2: Array of shape (N, 3, 3) containing ending rotation matrices. + t: Array of shape (N,) with interpolation values in [0, 1]. + correct_slerp: If True, use quaternion correction to ensure shortest path. + Defaults to False. + + Returns: + np.ndarray: Array of shape (N, 3, 3) containing interpolated rotation matrices. + """ + + # Prepare inputs + t = np.clip(np.asarray(t, dtype=float), 0.0, 1.0) + + # Ensure shortest path on the unit 4-sphere + if correct_slerp: + q1 = np.zeros((len(t), 4), dtype=float) + q2 = np.zeros((len(t), 4), dtype=float) + q1, q2 = qu.find_best_quaternions(Q1, Q2) + dot = np.einsum('ni,ni->n', q1, q2) + else: + q1 = qu.rotm_to_quat_batch(Q1) + q2 = qu.rotm_to_quat_batch(Q2) + dot = np.einsum('ni,ni->n', q1, q2) + if np.any(dot < 0.0): + neg_mask = dot < 0.0 + q2[neg_mask] = -q2[neg_mask] + dot[neg_mask] = -dot[neg_mask] + + # SLERP weights + dot_clipped = np.clip(dot, -1.0, 1.0) + theta0 = np.arccos(dot_clipped) + sin_theta0 = np.sin(theta0) + + # Threshold for linear interpolation + lin_mask = sin_theta0 < 1e-6 + q = np.empty_like(q1) + + if np.any(~lin_mask): + theta = theta0[~lin_mask] * t[~lin_mask] + s0 = np.sin(theta0[~lin_mask] - theta) / sin_theta0[~lin_mask] + s1 = np.sin(theta) / sin_theta0[~lin_mask] + q[~lin_mask] = (s0[:, None] * q1[~lin_mask]) + (s1[:, None] * q2[~lin_mask]) + + if np.any(lin_mask): + tl = t[lin_mask][:, None] + q[lin_mask] = (1.0 - tl) * q1[lin_mask] + tl * q2[lin_mask] + + # Normalize and convert back to rotation matrices + q /= np.linalg.norm(q, axis=1, keepdims=True) + return qu.quat_to_rotm_batch(q) + + def generate_fibers(self, params): + """Generate fiber directions. Override in subclasses.""" + raise NotImplementedError("Subclasses must implement generate_fibers()") + + def write_fibers(self, outdir): + """Write fiber, sheet, and normal directions to VTU files. + + Saves three separate files with fiber directions stored in 'FIB_DIR' field: + - fiber.vtu: Fiber directions + - sheet.vtu: Sheet normal directions + - normal.vtu: Sheet-normal directions + + Args: + outdir: Output directory path where files will be saved. + """ + # Create a copy of the mesh without any data + mesh_out = self.mesh.copy(deep=True) + mesh_out.clear_data() + + # Fiber direction + mesh_out.cell_data['FIB_DIR'] = self.mesh.cell_data['fiber'] + mesh_out.save(os.path.join(outdir, "fiber.vtu")) + + # Sheet direction + mesh_out.cell_data['FIB_DIR'] = self.mesh.cell_data['sheet'] + mesh_out.save(os.path.join(outdir, "sheet.vtu")) + + # Normal direction + mesh_out.cell_data['FIB_DIR'] = self.mesh.cell_data['sheet-normal'] + mesh_out.save(os.path.join(outdir, "normal.vtu")) + + +class FibGenBayer(FibGen): + """Fiber generator using the Bayer et al. (2012) method. + + Suitable for truncated biventricular geometries. Implements the rule-based + algorithm described in Bayer et al. 2012: + https://doi.org/10.1007/s10439-012-0593-5 + """ + + # Field names in Laplace solution + FIELD_NAMES = ['Trans_EPI', 'Trans_LV', 'Trans_RV', 'Long_AB'] + + def __init__(self): + """Initialize the Bayer fiber generator.""" + super().__init__() + + + def rescale_fields(self, mesh): + """Rescale Laplace fields to [0, 1] range. + + Args: + mesh: PyVista mesh with point data containing the fields to rescale. + + Returns: + PyVista mesh with rescaled fields in point_data. + """ + for name in self.FIELD_NAMES: + if name not in mesh.point_data: + raise KeyError(f"Field '{name}' not found in mesh point_data") + mesh.point_data[name] = self.scale_to_range(mesh.point_data[name], range=(0.0, 1.0)) + return mesh + + + def load_laplace_results(self, file_path): + """Load Laplace-Dirichlet solution for Bayer method. + + Args: + file_path: Path to the .vtu file with Laplace solution. + + Returns: + tuple: (lap, grad) dictionaries with Laplace values and gradients. + """ + print(f" Loading Laplace solution <--- {file_path}") + result_mesh = pv.read(file_path) + + # Normalize fields to [0, 1] + result_mesh = self.rescale_fields(result_mesh) + + # Compute gradients for the required fields + print(" Computing gradients at points") + result_mesh = self.compute_gradients(result_mesh, self.FIELD_NAMES) + + # Convert point-data to cell-data + mesh_cells = result_mesh.point_data_to_cell_data() + self.mesh = mesh_cells + + # Extract Laplace values and gradients + self.lap = {} + self.grad = {} + + for key in self.FIELD_NAMES: + self.lap[key] = np.asarray(mesh_cells.cell_data[key]) + self.grad[key] = np.asarray(mesh_cells.cell_data[key + "_grad"]) + + + return self.lap, self.grad + + + def generate_fibers(self, params, flip_rv=True, correct_slerp=False): + """Generate fiber directions using the Bayer method. + + Args: + params: Dictionary with keys: + - ALFA_END: Endocardial helix angle (degrees) + - ALFA_EPI: Epicardial helix angle (degrees) + - BETA_END: Endocardial transverse angle (degrees) + - BETA_EPI: Epicardial transverse angle (degrees) + flip_rv: If True, flip circumferential and transmural directions in RV. + Defaults to True. + correct_slerp: If True, use quaternion correction for SLERP interpolation. + Defaults to False. + + Returns: + tuple: (F, S, T) fiber, sheet, and normal directions (N, 3) each. + """ + if self.lap is None or self.grad is None: + raise ValueError("Must call load_laplace_results() first") + + # Convert parameters to radians (consistent with Doste method) + params = {k: np.deg2rad(v) for k, v in params.items()} + + print(" Computing fiber directions at cells") + + # Interpolation factor between LV and RV + d = self.lap['Trans_RV'] / (self.lap['Trans_LV'] + self.lap['Trans_RV']) + alfaS = self.calculate_angle(d, params['ALFA_END'], -params['ALFA_END']) + betaS = self.calculate_angle(d, params['BETA_END'], -params['BETA_END']) + + # Wall angles (interpolated from endo to epi) + alfaW = self.calculate_angle(self.lap['Trans_EPI'], params['ALFA_END'], params['ALFA_EPI']) + betaW = self.calculate_angle(self.lap['Trans_EPI'], params['BETA_END'], params['BETA_EPI']) + + # Build LV and RV basis + Q_LV0 = self.axis(self.grad['Long_AB'], -self.grad['Trans_LV']) + Q_LV = self.orient_matrix(Q_LV0, alfaS, np.sign(params['BETA_END'])*np.abs(betaS)) + + Q_RV0 = self.axis(self.grad['Long_AB'], self.grad['Trans_RV']) + Q_RV = self.orient_matrix(Q_RV0, alfaS, np.sign(params['BETA_END'])*np.abs(betaS)) + + # Interpolate between LV and RV (endocardial layer) + Q_END = self.interpolate_basis(Q_LV, Q_RV, d, correct_slerp=correct_slerp) + + # Flip circumferential and transmural directions in RV + if flip_rv: + Q_END[d > 0.5,:,0] = -Q_END[d > 0.5,:,0] + Q_END[d > 0.5,:,2] = -Q_END[d > 0.5,:,2] + + # Build epicardial basis + Q_EPI0 = self.axis(self.grad['Long_AB'], self.grad['Trans_EPI']) + Q_EPI = self.orient_matrix(Q_EPI0, alfaW, betaW) + + # Interpolate from endo to epi + FST = self.interpolate_basis(Q_END, Q_EPI, self.lap['Trans_EPI'], correct_slerp=correct_slerp) + + F = FST[:, :, 0] # Fiber direction + S = FST[:, :, 1] # Sheet normal + T = FST[:, :, 2] # Sheet direction + + self.mesh.cell_data['fiber'] = F + self.mesh.cell_data['sheet-normal'] = S + self.mesh.cell_data['sheet'] = T + + return F, S, T + + def get_angle_fields(self, params): + """Compute global alpha and beta angle fields. + + Helper function to compute spatially-varying helix and transverse angle fields + by interpolating between septum and wall values. + + Args: + params: Dictionary with angle parameters (in degrees or radians). + + Returns: + tuple: (alfa, beta) arrays of helix and transverse angles at each cell. + """ + + # Interpolation factor between LV and RV + d = self.lap['Trans_RV'] / (self.lap['Trans_LV'] + self.lap['Trans_RV']) + + # Septum angles (interpolated between LV and RV) + alfaS = self.calculate_angle(d, params['ALFA_END'], -params['ALFA_END']) + betaS = self.calculate_angle(d, params['BETA_END'], -params['BETA_END']) + alfaS = np.abs(alfaS) # Note this is doing the same as flipping the sign + betaS = np.sign(params['BETA_END'])*np.abs(betaS) # Note this is doing the same as flipping the sign + + # Wall angles (interpolated from endo to epi) + alfaW = self.calculate_angle(self.lap['Trans_EPI'], params['ALFA_END'], params['ALFA_EPI']) + betaW = self.calculate_angle(self.lap['Trans_EPI'], params['BETA_END'], params['BETA_EPI']) + + alfa = alfaS * (1 - self.lap['Trans_EPI']) + alfaW * self.lap['Trans_EPI'] + beta = betaS * (1 - self.lap['Trans_EPI']) + betaW * self.lap['Trans_EPI'] + + return alfa, beta + + +class FibGenDoste(FibGen): + """Fiber generator using the Doste et al. (2019) method. + + Suitable for biventricular geometries with outflow tracts. Implements + the algorithm described in Doste et al. 2019: + https://doi.org/10.1002/cnm.3185 + """ + + # Field names in Laplace solution + FIELD_NAMES = ['Trans_BiV', 'Long_AV', 'Long_MV', 'Long_PV', 'Long_TV', + 'Weight_LV', 'Weight_RV', 'Trans_EPI', 'Trans_LV', 'Trans_RV', 'Trans'] + + def __init__(self): + """Initialize the Doste fiber generator.""" + super().__init__() + + + def rescale_fields(self, mesh): + """Rescale Laplace fields to [0, 1] range. + + Args: + mesh: PyVista mesh with point data containing the fields to rescale. + + Returns: + PyVista mesh with rescaled fields in point_data. + """ + for name in self.FIELD_NAMES: + if name not in mesh.point_data: + raise KeyError(f"Field '{name}' not found in mesh point_data") + + # Set respective range + if name in ['Trans_BiV', 'Trans']: + range = (-2.0, 1.0) + else: + range = (0.0, 1.0) + + mesh.point_data[name] = self.scale_to_range(mesh.point_data[name], range=range) + + return mesh + + + def load_laplace_results(self, file_path): + """Load Laplace-Dirichlet solution for Doste method. + + Args: + file_path: Path to the .vtu file with Laplace solution. + + Returns: + tuple: (lap, grad) dictionaries with Laplace values and gradients. + """ + print(f" Loading Laplace solution <--- {file_path}") + result_mesh = pv.read(file_path) + + # Normalize fields to the original range + result_mesh = self.rescale_fields(result_mesh) + + # Compute gradients for the required fields + print(" Computing gradients at points") + result_mesh = self.compute_gradients(result_mesh, self.FIELD_NAMES) + + # Convert point-data to cell-data + mesh_cells = result_mesh.point_data_to_cell_data() + self.mesh = mesh_cells + + # Extract Laplace values and gradients using mapped names + self.lap = {} + self.grad = {} + + for key in self.FIELD_NAMES: + self.lap[key] = np.asarray(mesh_cells.cell_data[key]) + self.grad[key] = np.asarray(mesh_cells.cell_data[key + "_grad"]) + + return self.lap, self.grad + + + def _redistribute_weight(self, weight, up, low): + """Redistribute weight values to center their distribution. + + Args: + weight: Array of weight values. + up: Upper quantile threshold. + low: Lower quantile threshold. + + Returns: + np.ndarray: Redistributed weight values in [0, 1]. + """ + new_weight = weight.copy() + + upper_lim = np.quantile(weight, up) + while upper_lim == 0: + up += 0.1 + upper_lim = np.quantile(weight, up) + + lower_lim = np.quantile(weight, low) + + new_weight[new_weight > upper_lim] = upper_lim + new_weight[new_weight < lower_lim] = lower_lim + + return (new_weight - np.min(new_weight)) / (np.max(new_weight) - np.min(new_weight)) + + def _compute_basis_vectors(self): + """Compute local orthogonal basis vectors for LV and RV. + + Returns: + dict: Dictionary with basis vectors for LV, RV, and global. + """ + lap, grad = self.lap, self.grad + + # Calculate combined LV longitudinal + lv_glong = (grad['Long_MV'] * lap['Weight_LV'][:, None] + + grad['Long_AV'] * (1 - lap['Weight_LV'][:, None])) + + # Calculate LV basis + Q_lv = self.axis(-lv_glong, -grad['Trans_LV']) # Minus signs to match Bayer convention + eC_lv = Q_lv[:, :, 0] # Circumferential + eL_lv = Q_lv[:, :, 1] # Longitudinal + eT_lv = Q_lv[:, :, 2] # Transmural + + # Calculate combined RV longitudinal + rv_glong = (grad['Long_TV'] * lap['Weight_RV'][:, None] + + grad['Long_PV'] * (1 - lap['Weight_RV'][:, None])) + Q_rv = self.axis(-rv_glong, -grad['Trans_RV']) # Minus signs to match Bayer convention + eC_rv = Q_rv[:, :, 0] # Circumferential + eL_rv = Q_rv[:, :, 1] # Longitudinal + eT_rv = Q_rv[:, :, 2] # Transmural + + return { + 'eC_lv': eC_lv, 'eT_lv': eT_lv, 'eL_lv': eL_lv, + 'eC_rv': eC_rv, 'eT_rv': eT_rv, 'eL_rv': eL_rv, + } + + def _compute_angles(self, params): + """Compute spatially-varying alpha and beta angles. + + Args: + params: Dictionary with angle parameters (must be in radians). + + Returns: + dict: Dictionary of angle arrays including: + - alpha_wall_lv, alpha_wall_rv: Wall helix angles for LV/RV + - beta_wall_lv, beta_wall_rv: Wall transverse angles for LV/RV + - alfaS: Septum helix angle + - beta_septum: Septum transverse angle + """ + lap = self.lap + + # Redistribute weights + lv_weight = self._redistribute_weight(lap['Weight_LV'], 0.7, 0.01) + rv_weight = self._redistribute_weight(lap['Weight_RV'], 0.1, 0.001) + + # LV angles + alpha_lv_endo = params['AENDOLV'] * lv_weight + params['AOTENDOLV'] * (1 - lv_weight) + alpha_lv_epi = params['AEPILV'] * lv_weight + params['AOTEPILV'] * (1 - lv_weight) + alpha_wall_lv = self.calculate_angle(lap['Trans_EPI'], alpha_lv_endo, alpha_lv_epi) + beta_wall_lv = self.calculate_angle(lap['Trans_EPI'], params['BENDOLV'], params['BEPILV']) * lv_weight + + # RV angles + alpha_rv_endo = params['AENDORV'] * rv_weight + params['AOTENDORV'] * (1 - rv_weight) + alpha_rv_epi = params['AEPIRV'] * rv_weight + params['AOTEPIRV'] * (1 - rv_weight) + alpha_wall_rv = self.calculate_angle(lap['Trans_EPI'], alpha_rv_endo, alpha_rv_epi) + beta_wall_rv = self.calculate_angle(lap['Trans_EPI'], params['BENDORV'], params['BEPIRV']) * rv_weight + + # Septum angles + sep = lap['Trans'].copy() + sep[sep < 0.0] = sep[sep < 0.0] / 2 + sep = np.abs(sep) # This gives a field that is 1 at both endo but that assigns 2/3 of the septum to the lv + alpha_septum = alpha_lv_endo * sep * (lap['Trans_BiV'] < 0) + alpha_rv_endo * sep * (lap['Trans_BiV'] > 0) + beta_septum = params['BENDOLV'] * sep * (lap['Trans_BiV'] < 0) * lv_weight + params['BENDORV'] * sep * (lap['Trans_BiV'] > 0) * rv_weight + + return { + 'alpha_wall_lv': alpha_wall_lv, 'beta_wall_lv': beta_wall_lv, + 'alpha_wall_rv': alpha_wall_rv, 'beta_wall_rv': beta_wall_rv, + 'alpha_septum': alpha_septum, 'beta_septum': beta_septum + } + + + def generate_fibers(self, params): + """Generate fiber directions using the Doste method. + + Args: + params: Dictionary with angle parameters (in degrees): + - AENDOLV, AEPILV: LV endo/epi helix angles + - AENDORV, AEPIRV: RV endo/epi helix angles + - AOTENDOLV, AOTEPILV: LV outflow tract angles + - AOTENDORV, AOTEPIRV: RV outflow tract angles + - BENDOLV, BEPILV: LV endo/epi transverse angles + - BENDORV, BEPIRV: RV endo/epi transverse angles + + Returns: + tuple: (F, S, T) fiber, sheet, and normal directions (N, 3) each. + """ + if self.lap is None or self.grad is None: + raise ValueError("Must call load_laplace_results() first") + + # Convert parameters to radians (consistent with Bayer method) + params_rad = {k: np.deg2rad(v) for k, v in params.items()} + + print(" Computing basis vectors") + basis = self._compute_basis_vectors() + + print(" Computing angles") + angles = self._compute_angles(params_rad) + + print(" Computing local basis") + # Build basis matrices from vectors + Q_lv = np.stack([basis['eC_lv'], basis['eL_lv'], basis['eT_lv']], axis=-1) + Q_rv = np.stack([basis['eC_rv'], basis['eL_rv'], basis['eT_rv']], axis=-1) + + # Septum basis + Qlv_sep = self.orient_rodrigues( + Q_lv, angles['alpha_septum'], angles['beta_septum'] + ) + Qrv_sep = self.orient_rodrigues( + Q_rv, angles['alpha_septum'], angles['beta_septum'] + ) + + # Wall basis + Qlv_wall = self.orient_rodrigues( + Q_lv, angles['alpha_wall_lv'], angles['beta_wall_lv'] + ) + Qrv_wall = self.orient_rodrigues( + Q_rv, angles['alpha_wall_rv'], angles['beta_wall_rv'] + ) + + print(" Interpolating basis") + # Get interpolation factor between LV and RV + interp_biv = self.lap['Trans_BiV'].copy() # Need to normalize this to [0, 1] for interpolation + interp_biv[interp_biv < 0.0] = interp_biv[interp_biv < 0.0] / 2 + interp_biv = (interp_biv + 1) / 2 + + # Get discontinous septal fibers + Qsep = Qrv_sep.copy() + Qsep[interp_biv < 0.5] = Qlv_sep[interp_biv < 0.5] + + # Interpolate across ventricles + Qepi = self.interpolate_basis(Qlv_wall, Qrv_wall, interp_biv) + + # Interpolate from endo to epi + Q = self.interpolate_basis(Qsep, Qepi, self.lap['Trans_EPI']) + + print(" Done!") + F = Q[:, :, 0] # Fiber direction + S = Q[:, :, 1] # Sheet normal + T = Q[:, :, 2] # Sheet direction + + self.mesh.cell_data['fiber'] = F + self.mesh.cell_data['sheet-normal'] = S + self.mesh.cell_data['sheet'] = T + + return F, S, T + + + def get_angle_fields(self, params): + """Compute global alpha and beta angle fields. + + Helper function to compute spatially-varying helix and transverse angle fields + for the Doste method by interpolating between septum and wall values across + ventricles and transmural depth. + + Args: + params: Dictionary with angle parameters (in degrees or radians). + + Returns: + tuple: (alfa, beta) arrays of helix and transverse angles at each cell. + """ + + # Compute angles + angles = self._compute_angles(params) + + # Normalize Trans_BiV for interpolation (matching generate_fibers implementation) + interp_biv = self.lap['Trans_BiV'].copy() + interp_biv[interp_biv < 0.0] = interp_biv[interp_biv < 0.0] / 2 + interp_biv = (interp_biv + 1) / 2 + + # Interpolate wall angles between LV and RV + alpha_epi = angles['alpha_wall_lv'] * (1 - interp_biv) + angles['alpha_wall_rv'] * interp_biv + beta_epi = angles['beta_wall_lv'] * (1 - interp_biv) + angles['beta_wall_rv'] * interp_biv + + # Interpolate from septum (endo) to wall (epi) transmurally + alfa = angles['alpha_septum'] * (1 - self.lap['Trans_EPI']) + alpha_epi * self.lap['Trans_EPI'] + beta = angles['beta_septum'] * (1 - self.lap['Trans_EPI']) + beta_epi * self.lap['Trans_EPI'] + + return alfa, beta \ No newline at end of file diff --git a/utilities/fiber_generation/fiber_generation/laplace_solver.py b/utilities/fiber_generation/fiber_generation/laplace_solver.py new file mode 100644 index 000000000..e19381afb --- /dev/null +++ b/utilities/fiber_generation/fiber_generation/laplace_solver.py @@ -0,0 +1,445 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- +"""Laplace solver module for biventricular fiber generation. + +This module provides the LaplaceSolver class that generates XML configuration +files and runs the svMultiPhysics Laplace-Dirichlet solver for computing +scalar fields needed for fiber direction generation. + +Supports both: + - Bayer et al. (2012): Truncated BiV geometry + - Doste et al. (2019): BiV geometry with outflow tracts +""" + +import os +from xml.etree import ElementTree as ET +from xml.dom import minidom +from .surface_names import SurfaceName + + +class LaplaceSolver: + """Laplace solver for biventricular fiber generation. + + This class generates XML configuration files for the svMultiPhysics solver + and executes the Laplace-Dirichlet problems needed for fiber generation. + + Attributes: + mesh_path: Path to the volumetric mesh file (.vtu). + surface_paths: Dictionary mapping SurfaceName enum values to full file paths. + exec_svmultiphysics: Command to execute svMultiPhysics solver. + + Example: + >>> from src.SurfaceNames import SurfaceName + >>> solver = LaplaceSolver( + ... mesh_path="/path/to/mesh.vtu", + ... surface_paths={ + ... SurfaceName.EPICARDIUM: '/path/to/epicardium.vtp', + ... SurfaceName.ENDOCARDIUM_LV: '/path/to/lv_endocardium.vtp', + ... SurfaceName.ENDOCARDIUM_RV: '/path/to/rv_endocardium.vtp', + ... SurfaceName.BASE: '/path/to/base.vtp', + ... SurfaceName.EPICARDIUM_APEX: '/path/to/apex.vtp', + ... }, + ... exec_svmultiphysics="svmultiphysics " + ... ) + >>> result_file = solver.run("bayer", "/path/to/output") + """ + + # Solver configuration constants + DEFAULT_TOLERANCE = 1e-6 + DEFAULT_MAX_LS_ITERATIONS = 2000 + DEFAULT_MAX_NL_ITERATIONS = 5 + + def __init__(self, mesh_path, surface_paths, exec_svmultiphysics = "svmultiphysics "): + """Initialize the LaplaceSolver. + + Args: + mesh_path: Path to the volumetric mesh file (.vtu). + surface_paths: Dictionary mapping SurfaceName enum values to full file paths. + For Bayer method, required: SurfaceName.EPICARDIUM, ENDOCARDIUM_LV, ENDOCARDIUM_RV, BASE, EPICARDIUM_APEX + For Doste method, required: SurfaceName.EPICARDIUM, ENDOCARDIUM_LV, ENDOCARDIUM_RV, EPICARDIUM_APEX, MITRAL_VALVE, AORTIC_VALVE, TRICUSPID_VALVE, PULMONARY_VALVE + exec_svmultiphysics: Command to execute svMultiPhysics (e.g., "svmultiphysics "). + + Raises: + TypeError: If surface_paths contains non-enum keys. + """ + self.mesh_path = mesh_path + + # Validate that all keys are SurfaceName enum values + for key in surface_paths.keys(): + if not isinstance(key, SurfaceName): + raise TypeError( + f"surface_paths must use SurfaceName enum keys, not {type(key).__name__}. " + f"Got key: {key}" + ) + + self.surface_paths = surface_paths + self.exec_svmultiphysics = exec_svmultiphysics + + def _create_general_params(self, output_dir): + """Create the GeneralSimulationParameters XML element. + + Args: + output_dir: Directory to save results. + + Returns: + ET.Element: XML element for general simulation parameters. + """ + params = ET.Element("GeneralSimulationParameters") + + elements = [ + ("Continue_previous_simulation", "0"), + ("Number_of_spatial_dimensions", "3"), + ("Number_of_time_steps", "1"), + ("Time_step_size", "1"), + ("Spectral_radius_of_infinite_time_step", "0."), + ("Searched_file_name_to_trigger_stop", "STOP_SIM"), + ("Save_results_to_VTK_format", "1"), + ("Name_prefix_of_saved_VTK_files", "result"), + ("Increment_in_saving_VTK_files", "1"), + ("Save_results_in_folder", output_dir), + ("Start_saving_after_time_step", "1"), + ("Increment_in_saving_restart_files", "1"), + ("Convert_BIN_to_VTK_format", "0"), + ("Verbose", "1"), + ("Warning", "1"), + ("Debug", "0"), + ] + + for name, value in elements: + elem = ET.SubElement(params, name) + elem.text = f" {value} " + + return params + + def _create_mesh_element(self, face_names): + """Create the Add_mesh XML element with faces. + + Args: + face_names: List of face names to include. + + Returns: + ET.Element: XML element for mesh definition. + """ + mesh = ET.Element("Add_mesh", name="msh") + + mesh_path_elem = ET.SubElement(mesh, "Mesh_file_path") + mesh_path_elem.text = f" {self.mesh_path} " + + for face_name in face_names: + face = ET.SubElement(mesh, "Add_face", name=face_name) + face_path = ET.SubElement(face, "Face_file_path") + face_path.text = f" {self._get_surface_path(face_name)} " + + return mesh + + def _get_surface_path(self, face_name): + """Get the surface file path for a given face name. + + Args: + face_name: Internal face name used in XML. + + Returns: + str: Path to the surface file. + + Raises: + KeyError: If the surface is not found in surface_paths. + """ + # Convert XML face name to SurfaceName enum + surface_enum = SurfaceName.from_xml_face_name(face_name) + + if surface_enum is None: + raise KeyError(f"Unknown XML face name: {face_name}") + + if surface_enum not in self.surface_paths: + raise KeyError(f"Surface {surface_enum.value} not found in surface_paths") + + return self.surface_paths[surface_enum] + + def _create_equation(self, output_alias, boundary_conditions): + """Create a heat equation XML element for Laplace problem. + + Args: + output_alias: Name for the output temperature field. + boundary_conditions: List of tuples (face_name, value) for Dirichlet BCs. + + Returns: + ET.Element: XML element for the heat equation. + """ + eq = ET.Element("Add_equation", type="heatS") + + # Equation parameters + eq_params = [ + ("Coupled", "0"), + ("Min_iterations", "1"), + ("Max_iterations", str(self.DEFAULT_MAX_NL_ITERATIONS)), + ("Tolerance", str(self.DEFAULT_TOLERANCE)), + ("Conductivity", "1.0"), + ("Source_term", "0.0"), + ("Density", "0.0"), + ] + + for name, value in eq_params: + elem = ET.SubElement(eq, name) + elem.text = f" {value} " + + # Output configuration + output_spatial = ET.SubElement(eq, "Output", type="Spatial") + temp_spatial = ET.SubElement(output_spatial, "Temperature") + temp_spatial.text = " 1 " + + output_alias_elem = ET.SubElement(eq, "Output", type="Alias") + temp_alias = ET.SubElement(output_alias_elem, "Temperature") + temp_alias.text = f" {output_alias} " + + # Linear solver + ls = ET.SubElement(eq, "LS", type="CG") + la = ET.SubElement(ls, "Linear_algebra", type="fsils") + precond = ET.SubElement(la, "Preconditioner") + precond.text = " fsils " + max_iter = ET.SubElement(ls, "Max_iterations") + max_iter.text = f" {self.DEFAULT_MAX_LS_ITERATIONS} " + tol = ET.SubElement(ls, "Tolerance") + tol.text = f" {self.DEFAULT_TOLERANCE} " + + # Boundary conditions + for face_name, value in boundary_conditions: + bc = ET.SubElement(eq, "Add_BC", name=face_name) + bc_type = ET.SubElement(bc, "Type") + bc_type.text = " Dir " + bc_value = ET.SubElement(bc, "Value") + bc_value.text = f" {value} " + zero_perim = ET.SubElement(bc, "Zero_out_perimeter") + zero_perim.text = " 0 " + + return eq + + def _get_bayer_equations(self): + """Get equation definitions for the Bayer method. + + Returns: + list: List of (alias, boundary_conditions) tuples. + """ + return [ + # Trans_EPI: Transmural field (epicardium=1, endo=0) + ("Trans_EPI", [ + ("epicardium", 1.0), + ("lv_endocardium", 0.0), + ("rv_endocardium", 0.0), + ]), + # Trans_LV: LV field (lv_endo=1, others=0) + ("Trans_LV", [ + ("lv_endocardium", 1.0), + ("rv_endocardium", 0.0), + ("epicardium", 0.0), + ]), + # Trans_RV: RV field (rv_endo=1, others=0) + ("Trans_RV", [ + ("rv_endocardium", 1.0), + ("lv_endocardium", 0.0), + ("epicardium", 0.0), + ]), + # Long_AB: Apex-to-base field (base=1, apex=0) + ("Long_AB", [ + ("base", 1.0), + ("epi_apex", 0.0), + ]), + ] + + def _get_doste_equations(self): + """Get equation definitions for the Doste method. + + Returns: + list: List of (alias, boundary_conditions) tuples. + """ + return [ + # Trans_BiV: Ventricular transmural (LV=-2, RV=1) + ("Trans_BiV", [ + ("lv_endocardium", -2.0), + ("rv_endocardium", 1.0), + ]), + # Long_AV: LV longitudinal from AV (apex=1, aortic_valve=0) + ("Long_AV", [ + ("epi_apex", 1.0), + ("aortic_valve", 0.0), + ]), + # Long_MV: LV longitudinal from MV (apex=1, mitral_valve=0) + ("Long_MV", [ + ("epi_apex", 1.0), + ("mitral_valve", 0.0), + ]), + # Long_PV: RV longitudinal from PV (apex=1, pulmonary_valve=0) + ("Long_PV", [ + ("epi_apex", 1.0), + ("pulmonary_valve", 0.0), + ]), + # Long_TV: RV longitudinal from TV (apex=1, tricuspid_valve=0) + ("Long_TV", [ + ("epi_apex", 1.0), + ("tricuspid_valve", 0.0), + ]), + # Weight_LV: LV valve weights (mitral_valve=1, aortic_valve=0) + ("Weight_LV", [ + ("aortic_valve", 0.0), + ("mitral_valve", 1.0), + ]), + # Weight_RV: RV valve weights (tricuspid_valve=1, pulmonary_valve=0) + ("Weight_RV", [ + ("pulmonary_valve", 0.0), + ("tricuspid_valve", 1.0), + ]), + # Trans_EPI: Epicardial transmural (epicardium=1, endo=0) + ("Trans_EPI", [ + ("epicardium", 1.0), + ("lv_endocardium", 0.0), + ("rv_endocardium", 0.0), + ]), + # Trans_LV: LV transmural (lv_endo=1, others=0) + ("Trans_LV", [ + ("epicardium", 0.0), + ("rv_endocardium", 0.0), + ("lv_endocardium", 1.0), + ]), + # Trans_RV: RV transmural (rv_endo=1, others=0) + ("Trans_RV", [ + ("epicardium", 0.0), + ("lv_endocardium", 0.0), + ("rv_endocardium", 1.0), + ]), + # Trans: transmural (rv_endo=1, epi=0, lv_endo=-2) + ("Trans", [ + ("epicardium", 0.0), + ("lv_endocardium", -2.0), + ("rv_endocardium", 1.0), + ]), + ] + + def _get_face_names(self, method): + """Get the list of face names for a given method. + + Args: + method: Either "bayer" or "doste". + + Returns: + list: List of face names. + """ + if method == "bayer": + return ["epicardium", "base", "epi_apex", "lv_endocardium", "rv_endocardium"] + elif method == "doste": + return ["epicardium", "mitral_valve", "aortic_valve", "tricuspid_valve", "pulmonary_valve", "epi_apex", "lv_endocardium", "rv_endocardium"] + else: + raise ValueError(f"Unknown method: {method}. Use 'bayer' or 'doste'.") + + def _prettify_xml(self, elem): + """Convert an XML element to a prettified string. + + Args: + elem: ET.Element to prettify. + + Returns: + str: Prettified XML string. + """ + rough_string = ET.tostring(elem, encoding='unicode') + reparsed = minidom.parseString(rough_string) + return reparsed.toprettyxml(indent="\t") + + def generate_solver_xml(self, method, output_dir): + """Generate the svMultiPhysics XML configuration file. + + Args: + method: Either "bayer" or "doste". + output_dir: Directory to save solver results. + + Returns: + str: Path to the generated XML file. + """ + # Create root element + root = ET.Element("svMultiPhysicsFile", version="0.1") + + # Add general parameters + root.append(self._create_general_params(output_dir)) + + # Add mesh with faces + face_names = self._get_face_names(method) + root.append(self._create_mesh_element(face_names)) + + # Add equations based on method + if method == "bayer": + equations = self._get_bayer_equations() + elif method == "doste": + equations = self._get_doste_equations() + else: + raise ValueError(f"Unknown method: {method}. Use 'bayer' or 'doste'.") + + for alias, bcs in equations: + root.append(self._create_equation(alias, bcs)) + + # Generate XML string + xml_declaration = '\n' + xml_content = self._prettify_xml(root) + + # Remove the default XML declaration from prettify and use our own + if xml_content.startswith('', 1)[1].strip() + xml_content = xml_declaration + xml_content + + # Determine output path (same directory as mesh file) + xml_dir = os.path.dirname(self.mesh_path) + xml_path = os.path.join(xml_dir, "svFibers_BiV.xml") + + with open(xml_path, 'w') as f: + f.write(xml_content) + + return xml_path + + def run(self, method, output_dir, delete_xml = False): + """Generate XML and run the svMultiPhysics Laplace solver. + + This is the main entry point that generates the XML configuration + and executes the solver. + + Args: + method: Either "bayer" or "doste". + output_dir: Directory to save solver results. + + Returns: + str: Path to the result file containing Laplace solutions. + + Raises: + ValueError: If required surfaces are missing for the specified method. + """ + # Validate surfaces before running + self.validate_surfaces(method) + + # Ensure output directory exists + os.makedirs(output_dir, exist_ok=True) + + # Generate XML configuration + xml_path = self.generate_solver_xml(method, output_dir) + print(f"Generated solver XML at: {xml_path}") + + if delete_xml: + os.remove(xml_path) + + # Run solver + print(" Running svMultiPhysics solver") + print(f" {self.exec_svmultiphysics + xml_path}") + os.system(self.exec_svmultiphysics + xml_path) + + return os.path.join(output_dir, 'result_001.vtu') + + def validate_surfaces(self, method): + """Validate that all required surfaces are provided. + + Args: + method: Either "bayer" or "doste". + + Raises: + ValueError: If required surfaces are missing. + """ + required = SurfaceName.get_required_for_method(method) + available = set(self.surface_paths.keys()) + + missing = required - available + if missing: + missing_names = [s.value for s in missing] + raise ValueError(f"Missing required surfaces for {method} method: {missing_names}") diff --git a/utilities/fiber_generation/fiber_generation/quat_utils.py b/utilities/fiber_generation/fiber_generation/quat_utils.py new file mode 100644 index 000000000..78d9f40e6 --- /dev/null +++ b/utilities/fiber_generation/fiber_generation/quat_utils.py @@ -0,0 +1,141 @@ +import numpy as np + +def rotm_to_quat_batch(R): + # R: (N,3,3) -> q: (N,4) [w,x,y,z] + trace = np.einsum('nii->n', R) + q = np.zeros((R.shape[0], 4), dtype=float) + + # Branch where trace is positive + mask_t = trace > 0.0 + if np.any(mask_t): + S = np.sqrt(trace[mask_t] + 1.0) * 2.0 + q[mask_t, 0] = 0.25 * S + q[mask_t, 1] = (R[mask_t, 2, 1] - R[mask_t, 1, 2]) / S + q[mask_t, 2] = (R[mask_t, 0, 2] - R[mask_t, 2, 0]) / S + q[mask_t, 3] = (R[mask_t, 1, 0] - R[mask_t, 0, 1]) / S + + # For remaining, choose major diagonal + mask_f = ~mask_t + if np.any(mask_f): + Rf = R[mask_f] + m00 = Rf[:, 0, 0] + m11 = Rf[:, 1, 1] + m22 = Rf[:, 2, 2] + idx = np.argmax(np.stack([m00, m11, m22], axis=1), axis=1) + mf_idx = np.nonzero(mask_f)[0] + + for case_idx, (i, j, k) in enumerate([(0, 1, 2), (1, 0, 2), (2, 0, 1)]): + mask_case = idx == case_idx + if np.any(mask_case): + S = np.sqrt(1.0 + Rf[mask_case, i, i] - Rf[mask_case, j, j] - Rf[mask_case, k, k]) * 2.0 + rows = mf_idx[mask_case] + if case_idx == 0: + q[rows, 0] = (Rf[mask_case, 2, 1] - Rf[mask_case, 1, 2]) / S + q[rows, 1] = 0.25 * S + q[rows, 2] = (Rf[mask_case, 0, 1] + Rf[mask_case, 1, 0]) / S + q[rows, 3] = (Rf[mask_case, 0, 2] + Rf[mask_case, 2, 0]) / S + elif case_idx == 1: + q[rows, 0] = (Rf[mask_case, 0, 2] - Rf[mask_case, 2, 0]) / S + q[rows, 1] = (Rf[mask_case, 0, 1] + Rf[mask_case, 1, 0]) / S + q[rows, 2] = 0.25 * S + q[rows, 3] = (Rf[mask_case, 1, 2] + Rf[mask_case, 2, 1]) / S + else: + q[rows, 0] = (Rf[mask_case, 1, 0] - Rf[mask_case, 0, 1]) / S + q[rows, 1] = (Rf[mask_case, 0, 2] + Rf[mask_case, 2, 0]) / S + q[rows, 2] = (Rf[mask_case, 1, 2] + Rf[mask_case, 2, 1]) / S + q[rows, 3] = 0.25 * S + + # Normalize for numerical safety + q /= np.linalg.norm(q, axis=1, keepdims=True) + return q + +def quat_to_rotm_batch(q): + # q: (N,4) [w,x,y,z] -> R: (N,3,3) + w, x, y, z = q[:, 0], q[:, 1], q[:, 2], q[:, 3] + x2, y2, z2 = x * x, y * y, z * z + wx, wy, wz = w * x, w * y, w * z + xy, xz, yz = x * y, x * z, y * z + + R = np.zeros((q.shape[0], 3, 3), dtype=float) + R[:, 0, 0] = 1.0 - 2.0 * y2 - 2.0 * z2 + R[:, 1, 0] = 2.0 * xy + 2.0 * wz + R[:, 2, 0] = 2.0 * xz - 2.0 * wy + R[:, 0, 1] = 2.0 * xy - 2.0 * wz + R[:, 1, 1] = 1.0 - 2.0 * x2 - 2.0 * z2 + R[:, 2, 1] = 2.0 * yz + 2.0 * wx + R[:, 0, 2] = 2.0 * xz + 2.0 * wy + R[:, 1, 2] = 2.0 * yz - 2.0 * wx + R[:, 2, 2] = 1.0 - 2.0 * x2 - 2.0 * y2 + return R + +def quat_multiply_batch(q1, q2): + """Multiply two quaternions: q1 * q2 (batch operation). + + Args: + q1: First quaternion (N, 4) [w, x, y, z]. + q2: Second quaternion (N, 4) [w, x, y, z]. + + Returns: + np.ndarray: Product quaternion (N, 4) [w, x, y, z]. + """ + w1, x1, y1, z1 = q1[:, 0], q1[:, 1], q1[:, 2], q1[:, 3] + w2, x2, y2, z2 = q2[:, 0], q2[:, 1], q2[:, 2], q2[:, 3] + + result = np.zeros_like(q1) + result[:, 0] = w1*w2 - x1*x2 - y1*y2 - z1*z2 # w + result[:, 1] = w1*x2 + x1*w2 + y1*z2 - z1*y2 # x + result[:, 2] = w1*y2 - x1*z2 + y1*w2 + z1*x2 # y + result[:, 3] = w1*z2 + x1*y2 - y1*x2 + z1*w2 # z + + return result + +def find_best_quaternions(Q1, Q2): + """Find the best quaternion representations from rotation matrices (vectorized). + + Args: + Q1: Reference rotation matrix to rotate (N, 3, 3). + Q2: Target rotation matrix (N, 3, 3). + + Returns: + tuple: (q1_best, q2) where both are (N, 4) quaternions [w, x, y, z]. + """ + n = Q1.shape[0] + + # Convert rotation matrices to quaternions + q1 = rotm_to_quat_batch(Q1) + q2 = rotm_to_quat_batch(Q2) + + # Create 4 candidate quaternions for each element (original + 3 variations) + candidates = np.zeros((n, 4, 4), dtype=float) + + # Option 0: Original + candidates[:, 0] = q1 + + # Option 1: [-q1[1], q1[0], q1[3], -q1[2]] + candidates[:, 1, 0] = -q1[:, 1] + candidates[:, 1, 1] = q1[:, 0] + candidates[:, 1, 2] = q1[:, 3] + candidates[:, 1, 3] = -q1[:, 2] + + # Option 2: [-q1[2], -q1[3], q1[0], q1[1]] + candidates[:, 2, 0] = -q1[:, 2] + candidates[:, 2, 1] = -q1[:, 3] + candidates[:, 2, 2] = q1[:, 0] + candidates[:, 2, 3] = q1[:, 1] + + # Option 3: [-q1[3], q1[2], -q1[1], q1[0]] + candidates[:, 3, 0] = -q1[:, 3] + candidates[:, 3, 1] = q1[:, 2] + candidates[:, 3, 2] = -q1[:, 1] + candidates[:, 3, 3] = q1[:, 0] + + # Compute absolute dot products with q2 for all candidates + dots = np.abs(np.einsum('ni,nci->nc', q2, candidates)) + + # Find the candidate with the largest dot product + best_idx = np.argmax(dots, axis=1) + + # Select best quaternion for each element + q1_best = candidates[np.arange(n), best_idx] + + return q1_best, q2 diff --git a/utilities/fiber_generation/fiber_generation/surface_names.py b/utilities/fiber_generation/fiber_generation/surface_names.py new file mode 100644 index 000000000..49e88c8a9 --- /dev/null +++ b/utilities/fiber_generation/fiber_generation/surface_names.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- +"""Surface name enumerations for biventricular heart models. + +This module defines enumerations for surface names used in fiber generation +and Laplace solver configuration. +""" + +from enum import Enum + + +class SurfaceName(Enum): + """Enumeration of surface names for biventricular heart models. + + Attributes: + EPICARDIUM: Epicardial surface + BASE: Base surface (all valves together) + EPICARDIUM_APEX: Epicardial apex surface + ENDOCARDIUM_LV: Left ventricle endocardial surface + ENDOCARDIUM_RV: Right ventricle endocardial surface + MITRAL_VALVE: Mitral valve surface + AORTIC_VALVE: Aortic valve surface + TRICUSPID_VALVE: Tricuspid valve surface + PULMONARY_VALVE: Pulmonary valve surface + """ + EPICARDIUM = "epicardium" + BASE = "base" + EPICARDIUM_APEX = "epi_apex" + ENDOCARDIUM_LV = "lv_endocardium" + ENDOCARDIUM_RV = "rv_endocardium" + MITRAL_VALVE = "mitral_valve" + AORTIC_VALVE = "aortic_valve" + TRICUSPID_VALVE = "tricuspid_valve" + PULMONARY_VALVE = "pulmonary_valve" + + @classmethod + def from_xml_face_name(cls, xml_name): + """Convert XML face name to SurfaceName enum. + + Args: + xml_name: XML face name (e.g., 'epicardium'). + + Returns: + SurfaceName: Corresponding enum value. + """ + # Map XML face names to enum values + xml_to_enum = { + 'epicardium': cls.EPICARDIUM, + 'base': cls.BASE, + 'epi_apex': cls.EPICARDIUM_APEX, + 'lv_endocardium': cls.ENDOCARDIUM_LV, + 'rv_endocardium': cls.ENDOCARDIUM_RV, + 'mitral_valve': cls.MITRAL_VALVE, + 'aortic_valve': cls.AORTIC_VALVE, + 'tricuspid_valve': cls.TRICUSPID_VALVE, + 'pulmonary_valve': cls.PULMONARY_VALVE, + } + return xml_to_enum.get(xml_name, None) + + @classmethod + def get_required_for_method(cls, method): + """Get required surface names for a given method. + + Args: + method: Either "bayer" or "doste". + + Returns: + set: Set of required SurfaceName enum values. + """ + if method == "bayer": + return {cls.EPICARDIUM, cls.ENDOCARDIUM_LV, cls.ENDOCARDIUM_RV, cls.BASE, cls.EPICARDIUM_APEX} + elif method == "doste": + return {cls.EPICARDIUM, cls.ENDOCARDIUM_LV, cls.ENDOCARDIUM_RV, cls.EPICARDIUM_APEX, + cls.MITRAL_VALVE, cls.AORTIC_VALVE, cls.TRICUSPID_VALVE, cls.PULMONARY_VALVE} + else: + raise ValueError(f"Unknown method: {method}. Use 'bayer' or 'doste'.") diff --git a/utilities/fiber_generation/fiber_generation/surface_utils.py b/utilities/fiber_generation/fiber_generation/surface_utils.py new file mode 100644 index 000000000..74f39c890 --- /dev/null +++ b/utilities/fiber_generation/fiber_generation/surface_utils.py @@ -0,0 +1,111 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- +"""Surface mesh processing utilities for biventricular heart models. + +This module provides classes and utilities for processing surface meshes, +including generating epicardial apex surfaces. +""" + +import os +import numpy as np +import pyvista as pv +from .surface_names import SurfaceName + + +def get_normal_plane_svd(points): + """Find the plane that best fits a set of points using SVD. + + Args: + points: Array of shape (N, 3) representing 3D points. + + Returns: + tuple: A tuple containing: + - normal (np.ndarray): Unit normal vector to the fitted plane. + - centroid (np.ndarray): Centroid of the input points. + """ + centroid = np.mean(points, axis=0) + svd = np.linalg.svd(points - centroid) + normal = svd[2][-1] + normal = normal / np.linalg.norm(normal) + return normal, centroid + + +def generate_epi_apex(surface_paths): + """Generate the epicardial apex surface from the epicardial surface of the BiV. + + This function identifies the apex point of the epicardium and creates a surface + mesh containing elements that include the apex point. The surface is saved with + global node and element IDs. + + Args: + surface_paths: Dictionary mapping SurfaceName enum values to full file paths. + + The function requires: + - SurfaceName.EPICARDIUM: Epicardial surface + - SurfaceName.BASE: Base surface (for finding apex) + - SurfaceName.EPICARDIUM_APEX: Output surface name (will be created) + """ + # Load the epi surface + epi_path = surface_paths[SurfaceName.EPICARDIUM] + epi_mesh = pv.read(epi_path) + epi_points = epi_mesh.points + epi_cells = epi_mesh.faces + epi_eNoN = epi_cells[0] + epi_cells = epi_cells.reshape((-1, epi_eNoN + 1)) + epi_cells = epi_cells[:, 1:] + epi_global_node_id = epi_mesh.point_data['GlobalNodeID'] + epi_global_cell_id = epi_mesh.cell_data['GlobalElementID'] + + # Load the base surface + base_path = surface_paths[SurfaceName.BASE] + base_mesh = pv.read(base_path) + base_global_node_id = base_mesh.point_data['GlobalNodeID'] + + # Extract the boundary of the epi surface (at the top) to find the apex point + epi_base_global_node_id = np.intersect1d(epi_global_node_id, base_global_node_id) + epi_base_nodes = np.where(np.isin(epi_global_node_id, epi_base_global_node_id))[0] + + # Get normal + base_normal, base_centroid = get_normal_plane_svd(epi_points[epi_base_nodes, :]) + + # Find the index of the apex point of the epi surface + distance = np.abs(base_normal @ (epi_points - base_centroid).T) + epi_apex_point_index = np.argmax(distance) + + # Find elements containing the apex point + epi_apex_cell_index = np.where(epi_cells == epi_apex_point_index)[0] + + # Create epi_apex mesh + submesh_cells = epi_cells[epi_apex_cell_index] + submesh_xyz = np.zeros([len(np.unique(submesh_cells)), epi_points.shape[1]]) + map_mesh_submesh = np.ones(epi_points.shape[0], dtype=int) * -1 + map_submesh_mesh = np.zeros(submesh_xyz.shape[0], dtype=int) + child_elems_new = np.zeros(submesh_cells.shape, dtype=int) + + cont = 0 + for e in range(submesh_cells.shape[0]): + for i in range(submesh_cells.shape[1]): + if map_mesh_submesh[submesh_cells[e, i]] == -1: + child_elems_new[e, i] = cont + submesh_xyz[cont] = epi_points[submesh_cells[e, i]] + map_mesh_submesh[submesh_cells[e, i]] = cont + map_submesh_mesh[cont] = submesh_cells[e, i] + cont += 1 + else: + child_elems_new[e, i] = map_mesh_submesh[submesh_cells[e, i]] + + epi_apex_cells_type = np.full((child_elems_new.shape[0], 1), epi_eNoN) + epi_apex_cells = np.hstack((epi_apex_cells_type, child_elems_new)) + epi_apex_cells = np.hstack(epi_apex_cells) + + # Get global IDs + epi_apex_global_node_id = epi_global_node_id[map_submesh_mesh] + epi_apex_global_cell_id = epi_global_cell_id[epi_apex_cell_index] + + # Create and save mesh + epi_apex_mesh = pv.PolyData(submesh_xyz, epi_apex_cells) + epi_apex_mesh.point_data.set_array(epi_apex_global_node_id, 'GlobalNodeID') + epi_apex_mesh.cell_data.set_array(epi_apex_global_cell_id, 'GlobalElementID') + + epi_apex_path = surface_paths[SurfaceName.EPICARDIUM_APEX] + epi_apex_mesh.save(epi_apex_path) diff --git a/utilities/fiber_generation/main_bayer.py b/utilities/fiber_generation/main_bayer.py new file mode 100644 index 000000000..f6e794a2d --- /dev/null +++ b/utilities/fiber_generation/main_bayer.py @@ -0,0 +1,123 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- +"""Main script for generating biventricular fibers using the Bayer method. + +This module implements fiber generation for biventricular heart models using +the Laplace-Dirichlet rule-based method described in: +Bayer et al. 2012, "A Novel Rule-Based Algorithm for Assigning Myocardial +Fiber Orientation to Computational Heart Models" +https://doi.org/10.1007/s10439-012-0593-5 + +The script supports command-line arguments for customization of mesh paths, +output directories, and solver executables. +""" + +import argparse +import os +import pyvista as pv +from fiber_generation.laplace_solver import LaplaceSolver +from fiber_generation.fiber_generator import FibGenBayer +from fiber_generation.surface_names import SurfaceName +from fiber_generation.surface_utils import generate_epi_apex +from time import time + + +if __name__ == "__main__": + + ########################################################### + ############ USER INPUTS ################################ + ########################################################### + + run_flag = True + svmultiphysics_exec = "svmultiphysics " + + mesh_path = "example/biv_truncated/VOLUME.vtu" + outdir = "example/biv_truncated/output_bayer" + surfaces_dir = 'example/biv_truncated/mesh-surfaces' + + # Parameters for the Bayer et al. method https://doi.org/10.1007/s10439-012-0593-5. + params = { + "ALFA_END": 60.0, + "ALFA_EPI": -60.0, + "BETA_END": -20.0, + "BETA_EPI": 20.0, + } + + + ########################################################### + ############ FIBER GENERATION ########################### + ########################################################### + + # Optional CLI overrides + parser = argparse.ArgumentParser(description="Generate fibers using the Bayer method.") + parser.add_argument("--svmultiphysics-exec", default=svmultiphysics_exec, help="svMultiPhysics executable/command (default: %(default)s)") + parser.add_argument("--mesh-path", default=mesh_path, help="Path to the volumetric mesh .vtu (default: %(default)s)") + parser.add_argument( + "--surfaces-dir", + default=surfaces_dir, + help="Directory containing mesh surfaces; default: /mesh-surfaces", + ) + parser.add_argument("--outdir", default=outdir, help="Output directory (default: %(default)s)") + args = parser.parse_args() + + svmultiphysics_exec = args.svmultiphysics_exec + if not svmultiphysics_exec.endswith(" "): + svmultiphysics_exec = svmultiphysics_exec + " " + + mesh_path = args.mesh_path + outdir = args.outdir + + if args.surfaces_dir is None: + surfaces_dir = os.path.join(os.path.dirname(mesh_path), "mesh-surfaces") + else: + surfaces_dir = os.path.abspath(args.surfaces_dir) + + # Make sure the paths are full paths + mesh_path = os.path.abspath(mesh_path) + outdir = os.path.abspath(outdir) + surfaces_dir = os.path.abspath(surfaces_dir) + + # Define surface paths + surface_paths = {SurfaceName.EPICARDIUM: f'{surfaces_dir}/EPI.vtp', + SurfaceName.EPICARDIUM_APEX: f'{surfaces_dir}/EPI_APEX.vtp', + SurfaceName.BASE: f'{surfaces_dir}/BASE.vtp', + SurfaceName.ENDOCARDIUM_LV: f'{surfaces_dir}/LV.vtp', + SurfaceName.ENDOCARDIUM_RV: f'{surfaces_dir}/RV.vtp'} + + # Create output directory if needed + os.makedirs(outdir, exist_ok=True) + + # Check if the EPICARDIUM_APEX surface exists; if not create it + start = time() + if not os.path.exists(surface_paths[SurfaceName.EPICARDIUM_APEX]): + print("Generating EPICARDIUM_APEX surface...") + generate_epi_apex(surface_paths) + + # Initialize Laplace solver + solver = LaplaceSolver(mesh_path, surface_paths, svmultiphysics_exec) + + # Run the Laplace solver + if run_flag: + print("Running Laplace solver...") + laplace_results_file = solver.run("bayer", outdir) + else: + laplace_results_file = os.path.join(outdir, 'result_001.vtu') + + # Initialize fiber generator + print("\nGenerating fibers using Bayer method...") + fib_gen = FibGenBayer() + + # Load Laplace results + fib_gen.load_laplace_results(laplace_results_file) + + # Generate fiber directions + F, S, T = fib_gen.generate_fibers(params) + print(f"generate fibers (Bayer method) elapsed time: {time() - start:.3f} s") + + # Write fibers to output directory + fib_gen.write_fibers(outdir) + + # Save the result mesh + result_mesh_path = os.path.join(outdir, "results_bayer.vtu") + fib_gen.mesh.save(result_mesh_path) + print(f"\nResults saved to: {result_mesh_path}") diff --git a/utilities/fiber_generation/main_doste.py b/utilities/fiber_generation/main_doste.py new file mode 100644 index 000000000..822667a2a --- /dev/null +++ b/utilities/fiber_generation/main_doste.py @@ -0,0 +1,138 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- +"""Main script for generating biventricular fibers using the Doste method. + +This module implements fiber generation for biventricular heart models using +the Laplace-Dirichlet rule-based method described in: +Doste et al. 2019, "A rule-based method to model myocardial fiber orientation +in cardiac biventricular geometries with outflow tracts" +https://doi.org/10.1002/cnm.3185 + +The script supports command-line arguments for customization of mesh paths, +output directories, and solver executables. +""" + +import os +import argparse +import pyvista as pv +from fiber_generation.laplace_solver import LaplaceSolver +from fiber_generation.fiber_generator import FibGenDoste +from fiber_generation.surface_names import SurfaceName +from fiber_generation.surface_utils import generate_epi_apex +from time import time + + +if __name__ == "__main__": + ########################################################### + ############ USER INPUTS ################################ + ########################################################### + + run_flag = True + svmultiphysics_exec = "svmultiphysics " + + mesh_path = "example/biv_with_outflow_tracts/mesh-complete.mesh.vtu" + outdir = "example/biv_with_outflow_tracts/output_doste" + surfaces_dir = 'example/biv_with_outflow_tracts/mesh-surfaces' + + + # Parameters from the Doste paper https://doi.org/10.1002/cnm.3185 + params = { + 'AENDORV': 90, + 'AEPIRV': -25, + 'AENDOLV': 60, + 'AEPILV': -60, + + 'AOTENDOLV': 90, + 'AOTENDORV': 90, + 'AOTEPILV': 0, + 'AOTEPIRV': 0, + + 'BENDORV': 0, + 'BEPIRV': 20, + 'BENDOLV': -20, + 'BEPILV': 20, + } + + ########################################################### + ############ FIBER GENERATION ########################### + ########################################################### + + # Optional CLI overrides + parser = argparse.ArgumentParser(description="Generate fibers using the Doste method.") + parser.add_argument("--svmultiphysics-exec", default=svmultiphysics_exec, help="svMultiPhysics executable/command (default: %(default)s)") + parser.add_argument("--mesh-path", default=mesh_path, help="Path to the volumetric mesh .vtu (default: %(default)s)") + parser.add_argument( + "--surfaces-dir", + default=surfaces_dir, + help="Directory containing mesh surfaces; default: /mesh-surfaces", + ) + parser.add_argument("--outdir", default=outdir, help="Output directory (default: %(default)s)") + args = parser.parse_args() + + svmultiphysics_exec = args.svmultiphysics_exec + if not svmultiphysics_exec.endswith(" "): + svmultiphysics_exec = svmultiphysics_exec + " " + + mesh_path = args.mesh_path + outdir = args.outdir + + if args.surfaces_dir is None: + surfaces_dir = os.path.join(os.path.dirname(mesh_path), "mesh-surfaces") + else: + surfaces_dir = os.path.abspath(args.surfaces_dir) + + # Make sure the paths are full paths + mesh_path = os.path.abspath(mesh_path) + outdir = os.path.abspath(outdir) + surfaces_dir = os.path.abspath(surfaces_dir) + + # Define surface paths + surface_paths = { + SurfaceName.EPICARDIUM: f'{surfaces_dir}/epi.vtp', + SurfaceName.EPICARDIUM_APEX: f'{surfaces_dir}/epi_apex.vtp', + SurfaceName.AORTIC_VALVE: f'{surfaces_dir}/av.vtp', + SurfaceName.MITRAL_VALVE: f'{surfaces_dir}/mv.vtp', + SurfaceName.TRICUSPID_VALVE: f'{surfaces_dir}/tv.vtp', + SurfaceName.PULMONARY_VALVE: f'{surfaces_dir}/pv.vtp', + SurfaceName.ENDOCARDIUM_LV: f'{surfaces_dir}/endo_lv.vtp', + SurfaceName.ENDOCARDIUM_RV: f'{surfaces_dir}/endo_rv.vtp', + SurfaceName.BASE: f'{surfaces_dir}/top.vtp' + } + + # Create output directory if needed + os.makedirs(outdir, exist_ok=True) + + # Check if the EPICARDIUM_APEX surface exists; if not create it + start = time() + if not os.path.exists(surface_paths[SurfaceName.EPICARDIUM_APEX]): + print("Generating EPICARDIUM_APEX surface...") + generate_epi_apex(surface_paths) + + # Initialize Laplace solver + solver = LaplaceSolver(mesh_path, surface_paths, svmultiphysics_exec) + + # Run the Laplace solver + if run_flag: + print("Running Laplace solver...") + laplace_results_file = solver.run("doste", outdir) + else: + laplace_results_file = os.path.join(outdir, 'result_001.vtu') + + # Initialize fiber generator + print("\nGenerating fibers using Doste method...") + fib_gen = FibGenDoste() + + # Load Laplace results + fib_gen.load_laplace_results(laplace_results_file) + + # Generate fiber directions + F, S, T = fib_gen.generate_fibers(params) + print(f"generate fibers (Doste method) elapsed time: {time() - start:.3f} s") + + # Write fibers to output directory + fib_gen.write_fibers(outdir) + + # Save the result mesh + result_mesh_path = os.path.join(outdir, "results_doste.vtu") + fib_gen.mesh.save(result_mesh_path) + print(f"\nResults saved to: {result_mesh_path}") diff --git a/utilities/fiber_generation/paraview_bayer.py b/utilities/fiber_generation/paraview_bayer.py new file mode 100644 index 000000000..392ce1e9e --- /dev/null +++ b/utilities/fiber_generation/paraview_bayer.py @@ -0,0 +1,550 @@ +# trace generated using paraview version 6.0.1 +#import paraview +#paraview.compatibility.major = 6 +#paraview.compatibility.minor = 0 + +#### import the simple module from the paraview +from paraview.simple import * +#### disable automatic camera reset on 'Show' +paraview.simple._DisableFirstRenderCameraReset() +ResetSession() + +import os + +# Get the directory where this script is located +script_dir = os.path.dirname(os.path.abspath(__file__)) + +# Set paths relative to the script directory +validation_file_path = os.path.join(script_dir, 'example', 'biv_truncated', 'validation_bayer_combined.vtu') +png_output_path = os.path.join(script_dir, 'example', 'biv_truncated') + +fiber_families = ['f', 's', 'n'] +fiber_family_names = {'f': 'fiber', 's': 'sheet', 'n': 'sheet-normal'} + +# create a new 'XML Unstructured Grid Reader' +validation_bayer_combinedvtu = XMLUnstructuredGridReader(registrationName='validation_bayer_combined.vtu', FileName=[validation_file_path]) + +# Properties modified on validation_bayer_combinedvtu +validation_bayer_combinedvtu.TimeArray = 'None' + +# get active view +renderView1 = GetActiveViewOrCreate('RenderView') + +# Use a solid white background for all screenshots +renderView1.UseColorPaletteForBackground = 0 +renderView1.BackgroundColorMode = "Single Color" +renderView1.Background = [1.0, 1.0, 1.0] + +# Make axis titles/labels visible on white background +renderView1.AxesGrid.XTitleColor = [0.0, 0.0, 0.0] +renderView1.AxesGrid.YTitleColor = [0.0, 0.0, 0.0] +renderView1.AxesGrid.ZTitleColor = [0.0, 0.0, 0.0] +renderView1.AxesGrid.XLabelColor = [0.0, 0.0, 0.0] +renderView1.AxesGrid.YLabelColor = [0.0, 0.0, 0.0] +renderView1.AxesGrid.ZLabelColor = [0.0, 0.0, 0.0] + +# Make orientation axes (X/Y/Z) labels visible on white background +renderView1.OrientationAxesLabelColor = [0.0, 0.0, 0.0] +renderView1.OrientationAxesOutlineColor = [0.0, 0.0, 0.0] + +# show data in view +validation_bayer_combinedvtuDisplay = Show(validation_bayer_combinedvtu, renderView1, 'UnstructuredGridRepresentation') + +# trace defaults for the display properties. +validation_bayer_combinedvtuDisplay.Representation = 'Surface' + +# reset view to fit data +renderView1.ResetCamera(False, 0.9) + +# get the material library +materialLibrary1 = GetMaterialLibrary() + +# show color bar/color legend +validation_bayer_combinedvtuDisplay.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# get color transfer function/color map for 'f' +fLUT = GetColorTransferFunction('f') + +# Ensure scalar bar text is visible on white background +fLUTColorBar = GetScalarBar(fLUT, renderView1) +fLUTColorBar.TitleColor = [0.0, 0.0, 0.0] +fLUTColorBar.LabelColor = [0.0, 0.0, 0.0] + +# get opacity transfer function/opacity map for 'f' +fPWF = GetOpacityTransferFunction('f') + +# get 2D transfer function for 'f' +fTF2D = GetTransferFunction2D('f') + + +# get animation scene +animationScene1 = GetAnimationScene() + +# update animation scene based on data timesteps +animationScene1.UpdateAnimationUsingDataTimeSteps() + +# create a new 'Cell Data to Point Data' +cellDatatoPointData1 = CellDatatoPointData(registrationName='CellDatatoPointData1', Input=validation_bayer_combinedvtu) + +# show data in view +cellDatatoPointData1Display = Show(cellDatatoPointData1, renderView1, 'UnstructuredGridRepresentation') + +# trace defaults for the display properties. +cellDatatoPointData1Display.Representation = 'Surface' + +# hide data in view +Hide(validation_bayer_combinedvtu, renderView1) + +# show color bar/color legend +cellDatatoPointData1Display.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# create a new 'Stream Tracer' +streamTracer1 = StreamTracer(registrationName='StreamTracer1', Input=cellDatatoPointData1, + SeedType='Line') + +# toggle interactive widget visibility (only when running from the GUI) +ShowInteractiveWidgets(proxy=streamTracer1.SeedType) + +# Properties modified on streamTracer1 +streamTracer1.SeedType = 'Point Cloud' + +# Properties modified on streamTracer1.SeedType +streamTracer1.SeedType.Radius = 40.0 +streamTracer1.SeedType.NumberOfPoints = 50000 + +# show data in view +streamTracer1Display = Show(streamTracer1, renderView1, 'GeometryRepresentation') + +# trace defaults for the display properties. +streamTracer1Display.Representation = 'Surface' + +# hide data in view +Hide(cellDatatoPointData1, renderView1) + +# show color bar/color legend +streamTracer1Display.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# Properties modified on streamTracer1Display +streamTracer1Display.RenderLinesAsTubes = 1 + +# Properties modified on streamTracer1Display +streamTracer1Display.LineWidth = 2.0 + +# set scalar coloring +ColorBy(streamTracer1Display, ('POINTS', 'alpha_combined')) + +# Hide the scalar bar for this color map if no visible data is colored by it. +HideScalarBarIfNotNeeded(fLUT, renderView1) + +# rescale color and/or opacity maps used to include current data range +streamTracer1Display.RescaleTransferFunctionToDataRange(True, False) + +# show color bar/color legend +streamTracer1Display.SetScalarBarVisibility(renderView1, True) + +# get color transfer function/color map for 'alpha_combined' +alpha_combinedLUT = GetColorTransferFunction('alpha_combined') + +# get opacity transfer function/opacity map for 'alpha_combined' +alpha_combinedPWF = GetOpacityTransferFunction('alpha_combined') + +# get 2D transfer function for 'alpha_combined' +alpha_combinedTF2D = GetTransferFunction2D('alpha_combined') + +# set active source +SetActiveSource(validation_bayer_combinedvtu) + +# toggle interactive widget visibility (only when running from the GUI) +HideInteractiveWidgets(proxy=streamTracer1.SeedType) + +# get color legend/bar for alpha_combinedLUT in view renderView1 +alpha_combinedLUTColorBar = GetScalarBar(alpha_combinedLUT, renderView1) + +# Ensure scalar bar text is visible on white background +alpha_combinedLUTColorBar.TitleColor = [0.0, 0.0, 0.0] +alpha_combinedLUTColorBar.LabelColor = [0.0, 0.0, 0.0] + +# change scalar bar placement +alpha_combinedLUTColorBar.Set( + WindowLocation='Any Location', + Position=[0.8227593152064452, 0.22451317296678122], + ScalarBarLength=0.32999999999999996, +) + +# change scalar bar placement +alpha_combinedLUTColorBar.Set( + Position=[0.8277945619335347, 0.3115693012600229], + ScalarBarLength=0.3299999999999996, +) + +# change scalar bar placement +alpha_combinedLUTColorBar.Set( + Position=[0.8529707955689829, 0.320067884829428], + ScalarBarLength=0.32999999999999957, +) + +# set active source +SetActiveSource(streamTracer1) + +# toggle interactive widget visibility (only when running from the GUI) +ShowInteractiveWidgets(proxy=streamTracer1.SeedType) + +# set active source +SetActiveSource(validation_bayer_combinedvtu) + +# toggle interactive widget visibility (only when running from the GUI) +HideInteractiveWidgets(proxy=streamTracer1.SeedType) + +# show data in view +validation_bayer_combinedvtuDisplay = Show(validation_bayer_combinedvtu, renderView1, 'UnstructuredGridRepresentation') + +# show color bar/color legend +validation_bayer_combinedvtuDisplay.SetScalarBarVisibility(renderView1, True) + +# Properties modified on validation_bayer_combinedvtuDisplay +validation_bayer_combinedvtuDisplay.Opacity = 0.4 + +# Properties modified on validation_bayer_combinedvtuDisplay +validation_bayer_combinedvtuDisplay.Opacity = 0.1 + +# turn off scalar coloring +ColorBy(validation_bayer_combinedvtuDisplay, None) + +# Hide the scalar bar for this color map if no visible data is colored by it. +HideScalarBarIfNotNeeded(fLUT, renderView1) + +# set active source +SetActiveSource(streamTracer1) + +# toggle interactive widget visibility (only when running from the GUI) +ShowInteractiveWidgets(proxy=streamTracer1.SeedType) + +# update animation scene based on data timesteps +animationScene1.UpdateAnimationUsingDataTimeSteps() + +# set active source +SetActiveSource(validation_bayer_combinedvtu) + +# toggle interactive widget visibility (only when running from the GUI) +HideInteractiveWidgets(proxy=streamTracer1.SeedType) + +# get layout +layout1 = GetLayout() + +# layout/tab size in pixels +layout1.SetSize(993, 706) + +# current camera placement for renderView1 +renderView1.Set( + CameraPosition=[34.872946076797, -172.6220252737641, -28.100312441804963], + CameraFocalPoint=[33.09692217013306, -147.4005720279174, -138.5599245179334], + CameraViewUp=[-0.8581451756839255, 0.4973591506871757, 0.12736064022348584], + CameraParallelScale=51.95711542584741, +) + +# save screenshot +SaveScreenshot(filename=os.path.join(png_output_path, 'bayer_fiber.png'), viewOrLayout=renderView1, location=16, ImageResolution=[993, 706], TransparentBackground=0) + +# set active source +SetActiveSource(streamTracer1) + +# toggle interactive widget visibility (only when running from the GUI) +ShowInteractiveWidgets(proxy=streamTracer1.SeedType) + +# create a new 'Slice' +slice1 = Slice(registrationName='Slice1', Input=validation_bayer_combinedvtu) + +# Properties modified on slice1.SliceType +slice1.SliceType.Set( + Origin=[31.292619752159, -148.0730526113034, -141.31655248754385], + Normal=[-0.7480518557171969, 0.36676429493398416, 0.5530844177154476], +) + +# show data in view +slice1Display = Show(slice1, renderView1, 'GeometryRepresentation') + +# trace defaults for the display properties. +slice1Display.Representation = 'Surface' + +# show color bar/color legend +slice1Display.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# Hide streamTracer1 before the loop so it doesn't appear in screenshots +Hide(streamTracer1, renderView1) + +# Loop through fiber families to generate screenshots +for fiber_family in fiber_families: + family_name = fiber_family_names[fiber_family] + + # create a new 'Stream Tracer' for this fiber family + streamTracer_current = StreamTracer(registrationName=f'StreamTracer_{family_name}', Input=cellDatatoPointData1, + SeedType='Line') + + # Properties modified on streamTracer_current + streamTracer_current.SeedType = 'Point Cloud' + streamTracer_current.Vectors = ['POINTS', fiber_family] + + # Properties modified on streamTracer_current.SeedType + streamTracer_current.SeedType.Radius = 40.0 + streamTracer_current.SeedType.NumberOfPoints = 50000 + + # show data in view + streamTracer_currentDisplay = Show(streamTracer_current, renderView1, 'GeometryRepresentation') + + # trace defaults for the display properties. + streamTracer_currentDisplay.Representation = 'Surface' + streamTracer_currentDisplay.RenderLinesAsTubes = 1 + streamTracer_currentDisplay.LineWidth = 2.0 + + # set scalar coloring + ColorBy(streamTracer_currentDisplay, ('POINTS', 'alpha_combined')) + + # rescale color and/or opacity maps used to include current data range + streamTracer_currentDisplay.RescaleTransferFunctionToDataRange(True, False) + + # show color bar/color legend + streamTracer_currentDisplay.SetScalarBarVisibility(renderView1, True) + + # create a new 'Glyph' + glyph1 = Glyph(registrationName='Glyph1', Input=slice1, + GlyphType='Arrow') + + # toggle interactive widget visibility (only when running from the GUI) + ShowInteractiveWidgets(proxy=glyph1.GlyphType) + + # Properties modified on glyph1 + glyph1.Set( + GlyphType='Line', + OrientationArray=['CELLS', fiber_family], + ScaleArray=['POINTS', 'No scale array'], + GlyphMode='Uniform Spatial Distribution (Surface Sampling)', + ) + + # show data in view + glyph1Display = Show(glyph1, renderView1, 'GeometryRepresentation') + + # trace defaults for the display properties. + glyph1Display.Representation = 'Surface' + + # update the view to ensure updated data information + renderView1.Update() + + # set scalar coloring + ColorBy(glyph1Display, ('POINTS', 'alpha_combined')) + + # rescale color and/or opacity maps used to include current data range + glyph1Display.RescaleTransferFunctionToDataRange(True, False) + + # show color bar/color legend + glyph1Display.SetScalarBarVisibility(renderView1, True) + + # Properties modified on glyph1Display + glyph1Display.LineWidth = 2.0 + + # Properties modified on glyph1Display + glyph1Display.RenderLinesAsTubes = 1 + + # set active source + SetActiveSource(validation_bayer_combinedvtu) + + # toggle interactive widget visibility (only when running from the GUI) + HideInteractiveWidgets(proxy=glyph1.GlyphType) + + # set active source + SetActiveSource(cellDatatoPointData1) + + # set active source + SetActiveSource(slice1) + + # toggle interactive widget visibility (only when running from the GUI) + ShowInteractiveWidgets(proxy=slice1.SliceType) + + # turn off scalar coloring + ColorBy(slice1Display, None) + + # Hide the scalar bar for this color map if no visible data is colored by it. + HideScalarBarIfNotNeeded(fLUT, renderView1) + + # hide data in view + Hide(validation_bayer_combinedvtu, renderView1) + + # set active source + SetActiveSource(validation_bayer_combinedvtu) + + # toggle interactive widget visibility (only when running from the GUI) + HideInteractiveWidgets(proxy=slice1.SliceType) + + # set active source + SetActiveSource(glyph1) + + # toggle interactive widget visibility (only when running from the GUI) + ShowInteractiveWidgets(proxy=glyph1.GlyphType) + + # Properties modified on glyph1 + glyph1.MaximumNumberOfSamplePoints = 10000 + + # update the view to ensure updated data information + renderView1.Update() + + # change scalar bar placement + alpha_combinedLUTColorBar.Set( + Position=[0.879154078549849, 0.3583115108917509], + ScalarBarLength=0.32999999999999935, + ) + + + # hide data in view + Hide(validation_bayer_combinedvtu, renderView1) + + # hide data in view + Hide(streamTracer_current, renderView1) + + # layout/tab size in pixels + layout1.SetSize(993, 706) + + # current camera placement for renderView1 + renderView1.Set( + CameraPosition=[23.999840748360178, -164.79073803467477, -47.56738197808317], + CameraFocalPoint=[34.60514362681158, -149.76745665844024, -137.99940993913359], + CameraViewUp=[-0.7936032794303679, 0.6083829601515102, 0.00800054214733395], + CameraParallelScale=61.950390426874016, + ) + + # save screenshot + SaveScreenshot(filename=os.path.join(png_output_path, f'bayer_{family_name}_slice.png'), viewOrLayout=renderView1, location=16, ImageResolution=[993, 706], TransparentBackground=0) + + # set active source + SetActiveSource(validation_bayer_combinedvtu) + + # toggle interactive widget visibility (only when running from the GUI) + HideInteractiveWidgets(proxy=glyph1.GlyphType) + + # show data in view + validation_bayer_combinedvtuDisplay = Show(validation_bayer_combinedvtu, renderView1, 'UnstructuredGridRepresentation') + + # hide data in view + Hide(validation_bayer_combinedvtu, renderView1) + + # Delete the glyph for this iteration before creating the next one + Delete(glyph1) + del glyph1 + + # Show streamlines and save second screenshot for this fiber family + # set active source + SetActiveSource(streamTracer_current) + + # toggle interactive widget visibility (only when running from the GUI) + ShowInteractiveWidgets(proxy=streamTracer_current.SeedType) + + # show data in view + streamTracer_currentDisplay = Show(streamTracer_current, renderView1, 'GeometryRepresentation') + + # show color bar/color legend + streamTracer_currentDisplay.SetScalarBarVisibility(renderView1, True) + + # set active source + SetActiveSource(validation_bayer_combinedvtu) + + # toggle interactive widget visibility (only when running from the GUI) + HideInteractiveWidgets(proxy=streamTracer_current.SeedType) + + # show data in view + validation_bayer_combinedvtuDisplay = Show(validation_bayer_combinedvtu, renderView1, 'UnstructuredGridRepresentation') + + # hide data in view + Hide(slice1, renderView1) + + # layout/tab size in pixels + layout1.SetSize(993, 706) + + # current camera placement for renderView1 + renderView1.Set( + CameraPosition=[29.499659590041396, -176.0203361666033, -30.611812767158455], + CameraFocalPoint=[33.42253541923355, -147.58913078665822, -138.5225789422543], + CameraViewUp=[-0.8467462599797067, 0.5212189832741402, 0.10654361869699497], + CameraParallelScale=61.950390426874016, + ) + + # save screenshot + SaveScreenshot(filename=os.path.join(png_output_path, f'bayer_{family_name}.png'), viewOrLayout=renderView1, location=16, ImageResolution=[993, 706], TransparentBackground=0) + + # Delete the stream tracer for this iteration before creating the next one + Delete(streamTracer_current) + del streamTracer_current + + +# Final screenshots outside the loop +# layout/tab size in pixels +layout1.SetSize(993, 706) + +# current camera placement for renderView1 +renderView1.Set( + CameraPosition=[17.26650083175829, -175.3251410836073, -30.62963147529547], + CameraFocalPoint=[34.67637456755471, -149.61428340619153, -137.88774979097016], + CameraViewUp=[-0.8069186619533222, 0.590567665221697, 0.010588001985931447], + CameraParallelScale=61.950390426874016, +) + +# After loop, show streamlines and save final screenshots +# set active source +SetActiveSource(streamTracer1) + +#================================================================ +# addendum: following script captures some of the application +# state to faithfully reproduce the visualization during playback +#================================================================ + +#-------------------------------- +# saving layout sizes for layouts + +# layout/tab size in pixels +layout1.SetSize(993, 706) + +#----------------------------------- +# saving camera placements for views + +# current camera placement for renderView1 +renderView1.Set( + CameraPosition=[29.499659590041396, -176.0203361666033, -30.611812767158455], + CameraFocalPoint=[33.42253541923355, -147.58913078665822, -138.5225789422543], + CameraViewUp=[-0.8467462599797067, 0.5212189832741402, 0.10654361869699497], + CameraParallelScale=61.950390426874016, +) + + +##-------------------------------------------- +## You may need to add some code at the end of this python script depending on your usage, eg: +# +## Render all views to see them appears +# RenderAllViews() +# +## Interact with the view, usefull when running from pvpython +# Interact() +# +## Save a screenshot of the active view +# SaveScreenshot("path/to/screenshot.png") +# +## Save a screenshot of a layout (multiple splitted view) +# SaveScreenshot("path/to/screenshot.png", GetLayout()) +# +## Save all "Extractors" from the pipeline browser +# SaveExtracts() +# +## Save a animation of the current active view +# SaveAnimation() +# +## Please refer to the documentation of paraview.simple +## https://www.paraview.org/paraview-docs/nightly/python/ +##-------------------------------------------- \ No newline at end of file diff --git a/utilities/fiber_generation/paraview_doste.py b/utilities/fiber_generation/paraview_doste.py new file mode 100644 index 000000000..cbdad145a --- /dev/null +++ b/utilities/fiber_generation/paraview_doste.py @@ -0,0 +1,553 @@ +# trace generated using paraview version 6.0.1 +#import paraview +#paraview.compatibility.major = 6 +#paraview.compatibility.minor = 0 + +#### import the simple module from the paraview +from paraview.simple import * +#### disable automatic camera reset on 'Show' +paraview.simple._DisableFirstRenderCameraReset() +ResetSession() + +import os + +# Get the directory where this script is located +script_dir = os.path.dirname(os.path.abspath(__file__)) + +# Set paths relative to the script directory +validation_file_path = os.path.join(script_dir, 'example', 'biv_with_outflow_tracts', 'validation_doste_combined.vtu') +png_output_path = os.path.join(script_dir, 'example', 'biv_with_outflow_tracts') + +fiber_families = ['f', 's', 'n'] +fiber_family_names = {'f': 'fiber', 's': 'sheet', 'n': 'sheet-normal'} + +# create a new 'XML Unstructured Grid Reader' +validation_doste_combinedvtu = XMLUnstructuredGridReader(registrationName='validation_doste_combined.vtu', FileName=[validation_file_path]) + +# Properties modified on validation_doste_combinedvtu +validation_doste_combinedvtu.TimeArray = 'None' + +# get active view +renderView1 = GetActiveViewOrCreate('RenderView') + +# Use a solid white background for all screenshots +renderView1.UseColorPaletteForBackground = 0 +renderView1.BackgroundColorMode = "Single Color" +renderView1.Background = [1.0, 1.0, 1.0] + +# Make axis titles/labels visible on white background +renderView1.AxesGrid.XTitleColor = [0.0, 0.0, 0.0] +renderView1.AxesGrid.YTitleColor = [0.0, 0.0, 0.0] +renderView1.AxesGrid.ZTitleColor = [0.0, 0.0, 0.0] +renderView1.AxesGrid.XLabelColor = [0.0, 0.0, 0.0] +renderView1.AxesGrid.YLabelColor = [0.0, 0.0, 0.0] +renderView1.AxesGrid.ZLabelColor = [0.0, 0.0, 0.0] + +# Make orientation axes (X/Y/Z) labels visible on white background +renderView1.OrientationAxesLabelColor = [0.0, 0.0, 0.0] +renderView1.OrientationAxesOutlineColor = [0.0, 0.0, 0.0] + +# show data in view +validation_doste_combinedvtuDisplay = Show(validation_doste_combinedvtu, renderView1, 'UnstructuredGridRepresentation') + +# trace defaults for the display properties. +validation_doste_combinedvtuDisplay.Representation = 'Surface' + +# reset view to fit data +renderView1.ResetCamera(False, 0.9) + +# get the material library +materialLibrary1 = GetMaterialLibrary() + +# show color bar/color legend +validation_doste_combinedvtuDisplay.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# get color transfer function/color map for 'f' +fLUT = GetColorTransferFunction('f') + +# Ensure scalar bar text is visible on white background +fLUTColorBar = GetScalarBar(fLUT, renderView1) +fLUTColorBar.TitleColor = [0.0, 0.0, 0.0] +fLUTColorBar.LabelColor = [0.0, 0.0, 0.0] + +# get opacity transfer function/opacity map for 'f' +fPWF = GetOpacityTransferFunction('f') + +# get 2D transfer function for 'f' +fTF2D = GetTransferFunction2D('f') + + +# get animation scene +animationScene1 = GetAnimationScene() + +# update animation scene based on data timesteps +animationScene1.UpdateAnimationUsingDataTimeSteps() + +# create a new 'Cell Data to Point Data' +cellDatatoPointData1 = CellDatatoPointData(registrationName='CellDatatoPointData1', Input=validation_doste_combinedvtu) + +# show data in view +cellDatatoPointData1Display = Show(cellDatatoPointData1, renderView1, 'UnstructuredGridRepresentation') + +# trace defaults for the display properties. +cellDatatoPointData1Display.Representation = 'Surface' + +# hide data in view +Hide(validation_doste_combinedvtu, renderView1) + +# show color bar/color legend +cellDatatoPointData1Display.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# create a new 'Stream Tracer' +streamTracer1 = StreamTracer(registrationName='StreamTracer1', Input=cellDatatoPointData1, + SeedType='Line') + +# toggle interactive widget visibility (only when running from the GUI) +ShowInteractiveWidgets(proxy=streamTracer1.SeedType) + +# Properties modified on streamTracer1 +streamTracer1.SeedType = 'Point Cloud' + +# Properties modified on streamTracer1.SeedType +streamTracer1.SeedType.Radius = 60.0 +streamTracer1.SeedType.NumberOfPoints = 200000 + +# show data in view +streamTracer1Display = Show(streamTracer1, renderView1, 'GeometryRepresentation') + +# trace defaults for the display properties. +streamTracer1Display.Representation = 'Surface' + +# hide data in view +Hide(cellDatatoPointData1, renderView1) + +# show color bar/color legend +streamTracer1Display.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# Properties modified on streamTracer1Display +streamTracer1Display.RenderLinesAsTubes = 1 + +# Properties modified on streamTracer1Display +streamTracer1Display.LineWidth = 2.0 + +# set scalar coloring +ColorBy(streamTracer1Display, ('POINTS', 'alpha_combined')) + +# Hide the scalar bar for this color map if no visible data is colored by it. +HideScalarBarIfNotNeeded(fLUT, renderView1) + +# rescale color and/or opacity maps used to include current data range +streamTracer1Display.RescaleTransferFunctionToDataRange(True, False) + +# show color bar/color legend +streamTracer1Display.SetScalarBarVisibility(renderView1, True) + +# get color transfer function/color map for 'alpha_combined' +alpha_combinedLUT = GetColorTransferFunction('alpha_combined') + +# get opacity transfer function/opacity map for 'alpha_combined' +alpha_combinedPWF = GetOpacityTransferFunction('alpha_combined') + +# get 2D transfer function for 'alpha_combined' +alpha_combinedTF2D = GetTransferFunction2D('alpha_combined') + +# set active source +SetActiveSource(validation_doste_combinedvtu) + +# toggle interactive widget visibility (only when running from the GUI) +HideInteractiveWidgets(proxy=streamTracer1.SeedType) + +# get color legend/bar for alpha_combinedLUT in view renderView1 +alpha_combinedLUTColorBar = GetScalarBar(alpha_combinedLUT, renderView1) + +# Ensure scalar bar text is visible on white background +alpha_combinedLUTColorBar.TitleColor = [0.0, 0.0, 0.0] +alpha_combinedLUTColorBar.LabelColor = [0.0, 0.0, 0.0] + +# change scalar bar placement +alpha_combinedLUTColorBar.Set( + WindowLocation='Any Location', + Position=[0.8227593152064452, 0.22451317296678122], + ScalarBarLength=0.32999999999999996, +) + +# change scalar bar placement +alpha_combinedLUTColorBar.Set( + Position=[0.8277945619335347, 0.3115693012600229], + ScalarBarLength=0.3299999999999996, +) + +# change scalar bar placement +alpha_combinedLUTColorBar.Set( + Position=[0.8529707955689829, 0.320067884829428], + ScalarBarLength=0.32999999999999957, +) + +# set active source +SetActiveSource(streamTracer1) + +# toggle interactive widget visibility (only when running from the GUI) +ShowInteractiveWidgets(proxy=streamTracer1.SeedType) + +# set active source +SetActiveSource(validation_doste_combinedvtu) + +# toggle interactive widget visibility (only when running from the GUI) +HideInteractiveWidgets(proxy=streamTracer1.SeedType) + +# show data in view +validation_doste_combinedvtuDisplay = Show(validation_doste_combinedvtu, renderView1, 'UnstructuredGridRepresentation') + +# show color bar/color legend +validation_doste_combinedvtuDisplay.SetScalarBarVisibility(renderView1, True) + +# Properties modified on validation_doste_combinedvtuDisplay +validation_doste_combinedvtuDisplay.Opacity = 0.4 + +# Properties modified on validation_doste_combinedvtuDisplay +validation_doste_combinedvtuDisplay.Opacity = 0.1 + +# turn off scalar coloring +ColorBy(validation_doste_combinedvtuDisplay, None) + +# Hide the scalar bar for this color map if no visible data is colored by it. +HideScalarBarIfNotNeeded(fLUT, renderView1) + +# set active source +SetActiveSource(streamTracer1) + +# toggle interactive widget visibility (only when running from the GUI) +ShowInteractiveWidgets(proxy=streamTracer1.SeedType) + +# update animation scene based on data timesteps +animationScene1.UpdateAnimationUsingDataTimeSteps() + +# set active source +SetActiveSource(validation_doste_combinedvtu) + +# toggle interactive widget visibility (only when running from the GUI) +HideInteractiveWidgets(proxy=streamTracer1.SeedType) + +# get layout +layout1 = GetLayout() + +# layout/tab size in pixels +layout1.SetSize(993, 706) + +# current camera placement for renderView1 +renderView1.Set( + CameraPosition=[52.04115429873035, -45.38613870538731, 216.54244474035633], + CameraFocalPoint=[-21.838666915893594, -60.56520140171051, 48.49395084381101], + CameraViewUp=[-0.11136472699136057, 0.9929446544942593, -0.040728499768785655], + CameraParallelScale=102.19338444140465, +) + + +# save screenshot +SaveScreenshot(filename=os.path.join(png_output_path, 'doste_fiber.png'), viewOrLayout=renderView1, location=16, ImageResolution=[993, 706], TransparentBackground=0) + +# set active source +SetActiveSource(streamTracer1) + +# toggle interactive widget visibility (only when running from the GUI) +ShowInteractiveWidgets(proxy=streamTracer1.SeedType) + +# create a new 'Slice' +slice1 = Slice(registrationName='Slice1', Input=validation_doste_combinedvtu) + +# Properties modified on slice1.SliceType +slice1.SliceType.Set( + Origin=[-25.530508094611232, -64.79965909808317, 42.457392775039196], + Normal=[0.3179741978250393, 0.5256304018456659, 0.7890532872836198], +) + +# show data in view +slice1Display = Show(slice1, renderView1, 'GeometryRepresentation') + +# trace defaults for the display properties. +slice1Display.Representation = 'Surface' + +# show color bar/color legend +slice1Display.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# Hide streamTracer1 before the loop so it doesn't appear in screenshots +Hide(streamTracer1, renderView1) + +# Loop through fiber families to generate screenshots +for fiber_family in fiber_families: + family_name = fiber_family_names[fiber_family] + + # create a new 'Stream Tracer' for this fiber family + streamTracer_current = StreamTracer(registrationName=f'StreamTracer_{family_name}', Input=cellDatatoPointData1, + SeedType='Line') + + # Properties modified on streamTracer_current + streamTracer_current.SeedType = 'Point Cloud' + streamTracer_current.Vectors = ['POINTS', fiber_family] + + # Properties modified on streamTracer_current.SeedType + streamTracer_current.SeedType.Radius = 60.0 + streamTracer_current.SeedType.NumberOfPoints = 200000 + + # show data in view + streamTracer_currentDisplay = Show(streamTracer_current, renderView1, 'GeometryRepresentation') + + # trace defaults for the display properties. + streamTracer_currentDisplay.Representation = 'Surface' + streamTracer_currentDisplay.RenderLinesAsTubes = 1 + streamTracer_currentDisplay.LineWidth = 2.0 + + # set scalar coloring + ColorBy(streamTracer_currentDisplay, ('POINTS', 'alpha_combined')) + + # rescale color and/or opacity maps used to include current data range + streamTracer_currentDisplay.RescaleTransferFunctionToDataRange(True, False) + + # show color bar/color legend + streamTracer_currentDisplay.SetScalarBarVisibility(renderView1, True) + + # create a new 'Glyph' + glyph1 = Glyph(registrationName='Glyph1', Input=slice1, + GlyphType='Arrow') + + # toggle interactive widget visibility (only when running from the GUI) + ShowInteractiveWidgets(proxy=glyph1.GlyphType) + + # Properties modified on glyph1 + glyph1.Set( + GlyphType='Line', + OrientationArray=['CELLS', fiber_family], + ScaleArray=['POINTS', 'No scale array'], + GlyphMode='Uniform Spatial Distribution (Surface Sampling)', + ) + + # show data in view + glyph1Display = Show(glyph1, renderView1, 'GeometryRepresentation') + + # trace defaults for the display properties. + glyph1Display.Representation = 'Surface' + + # update the view to ensure updated data information + renderView1.Update() + + # set scalar coloring + ColorBy(glyph1Display, ('POINTS', 'alpha_combined')) + + # rescale color and/or opacity maps used to include current data range + glyph1Display.RescaleTransferFunctionToDataRange(True, False) + + # show color bar/color legend + glyph1Display.SetScalarBarVisibility(renderView1, True) + + # Properties modified on glyph1Display + glyph1Display.LineWidth = 2.0 + + # Properties modified on glyph1Display + glyph1Display.RenderLinesAsTubes = 1 + + # set active source + SetActiveSource(validation_doste_combinedvtu) + + # toggle interactive widget visibility (only when running from the GUI) + HideInteractiveWidgets(proxy=glyph1.GlyphType) + + # set active source + SetActiveSource(cellDatatoPointData1) + + # set active source + SetActiveSource(slice1) + + # toggle interactive widget visibility (only when running from the GUI) + ShowInteractiveWidgets(proxy=slice1.SliceType) + + # turn off scalar coloring + ColorBy(slice1Display, None) + + # Hide the scalar bar for this color map if no visible data is colored by it. + HideScalarBarIfNotNeeded(fLUT, renderView1) + + # hide data in view + Hide(validation_doste_combinedvtu, renderView1) + + # set active source + SetActiveSource(validation_doste_combinedvtu) + + # toggle interactive widget visibility (only when running from the GUI) + HideInteractiveWidgets(proxy=slice1.SliceType) + + # set active source + SetActiveSource(glyph1) + + # toggle interactive widget visibility (only when running from the GUI) + ShowInteractiveWidgets(proxy=glyph1.GlyphType) + + # Properties modified on glyph1 + glyph1.MaximumNumberOfSamplePoints = 10000 + + # update the view to ensure updated data information + renderView1.Update() + + # change scalar bar placement + alpha_combinedLUTColorBar.Set( + Position=[0.879154078549849, 0.3583115108917509], + ScalarBarLength=0.32999999999999935, + ) + + + # hide data in view + Hide(validation_doste_combinedvtu, renderView1) + + # hide data in view + Hide(streamTracer_current, renderView1) + + # layout/tab size in pixels + layout1.SetSize(993, 706) + + # current camera placement for renderView1 + renderView1.Set( + CameraPosition=[58.74149269073567, -105.91797844910512, 207.98796445502148], + CameraFocalPoint=[-22.22441909870438, -68.10289751255294, 46.918200036593994], + CameraViewUp=[-0.015152410940450833, 0.9716964855848043, 0.2357463559523924], + CameraParallelScale=102.19338444140465, + ) + + # save screenshot + SaveScreenshot(filename=os.path.join(png_output_path, f'doste_{family_name}_slice.png'), viewOrLayout=renderView1, location=16, ImageResolution=[993, 706], TransparentBackground=0) + + # set active source + SetActiveSource(validation_doste_combinedvtu) + + # toggle interactive widget visibility (only when running from the GUI) + HideInteractiveWidgets(proxy=glyph1.GlyphType) + + # show data in view + validation_doste_combinedvtuDisplay = Show(validation_doste_combinedvtu, renderView1, 'UnstructuredGridRepresentation') + + # hide data in view + Hide(validation_doste_combinedvtu, renderView1) + + # Delete the glyph for this iteration before creating the next one + Delete(glyph1) + del glyph1 + + # Show streamlines and save second screenshot for this fiber family + # set active source + SetActiveSource(streamTracer_current) + + # toggle interactive widget visibility (only when running from the GUI) + ShowInteractiveWidgets(proxy=streamTracer_current.SeedType) + + # show data in view + streamTracer_currentDisplay = Show(streamTracer_current, renderView1, 'GeometryRepresentation') + + # show color bar/color legend + streamTracer_currentDisplay.SetScalarBarVisibility(renderView1, True) + + # set active source + SetActiveSource(validation_doste_combinedvtu) + + # toggle interactive widget visibility (only when running from the GUI) + HideInteractiveWidgets(proxy=streamTracer_current.SeedType) + + # show data in view + validation_doste_combinedvtuDisplay = Show(validation_doste_combinedvtu, renderView1, 'UnstructuredGridRepresentation') + + # hide data in view + Hide(slice1, renderView1) + + # layout/tab size in pixels + layout1.SetSize(993, 706) + + # current camera placement for renderView1 + renderView1.Set( + CameraPosition=[52.04115429873035, -45.38613870538731, 216.54244474035633], + CameraFocalPoint=[-21.838666915893594, -60.56520140171051, 48.49395084381101], + CameraViewUp=[-0.11136472699136057, 0.9929446544942593, -0.040728499768785655], + CameraParallelScale=102.19338444140465, + ) + + # save screenshot + SaveScreenshot(filename=os.path.join(png_output_path, f'doste_{family_name}.png'), viewOrLayout=renderView1, location=16, ImageResolution=[993, 706], TransparentBackground=0) + + # Delete the stream tracer for this iteration before creating the next one + Delete(streamTracer_current) + del streamTracer_current + + +# Final screenshots outside the loop +# layout/tab size in pixels +layout1.SetSize(993, 706) + +# current camera placement for renderView1 +renderView1.Set( + CameraPosition=[52.04115429873035, -45.38613870538731, 216.54244474035633], + CameraFocalPoint=[-21.838666915893594, -60.56520140171051, 48.49395084381101], + CameraViewUp=[-0.11136472699136057, 0.9929446544942593, -0.040728499768785655], + CameraParallelScale=102.19338444140465, +) + + +# After loop, show streamlines and save final screenshots +# set active source +SetActiveSource(streamTracer1) + +#================================================================ +# addendum: following script captures some of the application +# state to faithfully reproduce the visualization during playback +#================================================================ + +#-------------------------------- +# saving layout sizes for layouts + +# layout/tab size in pixels +layout1.SetSize(993, 706) + +#----------------------------------- +# saving camera placements for views + +# current camera placement for renderView1 +renderView1.Set( + CameraPosition=[52.04115429873035, -45.38613870538731, 216.54244474035633], + CameraFocalPoint=[-21.838666915893594, -60.56520140171051, 48.49395084381101], + CameraViewUp=[-0.11136472699136057, 0.9929446544942593, -0.040728499768785655], + CameraParallelScale=102.19338444140465, +) + + + +##-------------------------------------------- +## You may need to add some code at the end of this python script depending on your usage, eg: +# +## Render all views to see them appears +# RenderAllViews() +# +## Interact with the view, usefull when running from pvpython +# Interact() +# +## Save a screenshot of the active view +# SaveScreenshot("path/to/screenshot.png") +# +## Save a screenshot of a layout (multiple splitted view) +# SaveScreenshot("path/to/screenshot.png", GetLayout()) +# +## Save all "Extractors" from the pipeline browser +# SaveExtracts() +# +## Save a animation of the current active view +# SaveAnimation() +# +## Please refer to the documentation of paraview.simple +## https://www.paraview.org/paraview-docs/nightly/python/ +##-------------------------------------------- \ No newline at end of file diff --git a/utilities/fiber_generation/pyproject.toml b/utilities/fiber_generation/pyproject.toml new file mode 100644 index 000000000..acaf2ff51 --- /dev/null +++ b/utilities/fiber_generation/pyproject.toml @@ -0,0 +1,22 @@ +[build-system] +requires = ["setuptools>=61.0", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "sv-fibergen" +version = "0.1.0" +description = "Python + svMultiPhysics codes for biventricular fiber generation" +readme = "README.md" +license = "MIT" +requires-python = ">=3.8" + +dependencies = [ + "numpy>=1.20.0", + "pyvista>=0.38.0", + "scipy>=1.7.0", + "matplotlib>=3.4.0", +] + +[tool.setuptools] +packages = ["fiber_generation"] + diff --git a/utilities/fiber_generation/test_fibgen_functions.py b/utilities/fiber_generation/test_fibgen_functions.py new file mode 100644 index 000000000..6ef7cae41 --- /dev/null +++ b/utilities/fiber_generation/test_fibgen_functions.py @@ -0,0 +1,596 @@ +#!/usr/bin/env python +"""Unit tests for FibGen class methods. + +Tests the core functions: axis, orient_matrix, orient_rodrigues, and interpolate_basis. +""" + +import unittest +import numpy as np +import sys +import os +from fiber_generation.fiber_generator import FibGen + + +class TestFibGen(unittest.TestCase): + """Test cases for FibGen base class methods.""" + + def setUp(self): + """Set up test fixtures.""" + self.fibgen = FibGen() + self.tol = 1e-10 # Tolerance for floating point comparisons + + def test_axis_orthogonality(self): + """Test that axis() produces orthogonal basis vectors.""" + n = 10 + np.random.seed(42) + gL = np.random.randn(n, 3) + gT = np.random.randn(n, 3) + + Q = self.fibgen.axis(gL, gT) + + # Check shape + self.assertEqual(Q.shape, (n, 3, 3)) + + # Check orthogonality for each element + for i in range(n): + eC = Q[i, :, 0] + eL = Q[i, :, 1] + eT = Q[i, :, 2] + + # Check dot products (should be ~0) + self.assertAlmostEqual(np.dot(eC, eL), 0.0, places=10) + self.assertAlmostEqual(np.dot(eC, eT), 0.0, places=10) + self.assertAlmostEqual(np.dot(eL, eT), 0.0, places=10) + + # Check normalization + self.assertAlmostEqual(np.linalg.norm(eC), 1.0, places=10) + self.assertAlmostEqual(np.linalg.norm(eL), 1.0, places=10) + self.assertAlmostEqual(np.linalg.norm(eT), 1.0, places=10) + + def test_axis_orientation(self): + """Test that axis() produces correctly oriented basis.""" + n = 5 + np.random.seed(42) + gL = np.random.randn(n, 3) + gT = np.random.randn(n, 3) + + Q = self.fibgen.axis(gL, gT) + + for i in range(n): + eC = Q[i, :, 0] + eL = Q[i, :, 1] + eT = Q[i, :, 2] + + # eL should be aligned with normalized gL + gL_norm = gL[i] / np.linalg.norm(gL[i]) + self.assertTrue(np.allclose(eL, gL_norm) or np.allclose(eL, -gL_norm)) + + # eT should be orthogonal to eL + self.assertAlmostEqual(np.dot(eT, eL), 0.0, places=10) + + # eC should be cross product of eL and eT + eC_expected = np.cross(eL, eT) + eC_expected = eC_expected / np.linalg.norm(eC_expected) + self.assertTrue(np.allclose(eC, eC_expected) or np.allclose(eC, -eC_expected)) + + def test_axis_zero_gradient(self): + """Test axis() with zero gradient vectors.""" + n = 3 + gL = np.zeros((n, 3)) + gT = np.random.randn(n, 3) + + Q = self.fibgen.axis(gL, gT) + + # Check that zero vectors are handled gracefully + self.assertEqual(Q.shape, (n, 3, 3)) + # Zero vectors should result in zero output (normalize handles this) + for i in range(n): + eL = Q[i, :, 1] + self.assertTrue(np.allclose(eL, 0.0)) + + def test_axis_parallel_vectors(self): + """Test axis() when gL and gT are parallel.""" + # Use exactly parallel vectors (not random to avoid numerical issues) + gL = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) + gT = 2.0 * gL # Exactly parallel to gL + + Q = self.fibgen.axis(gL, gT) + + # When parallel, eT should be zero (projection removes everything) + # normalize() sets zero vectors to zero + for i in range(len(gL)): + eT = Q[i, :, 2] + # eT should be zero (normalize handles this) + self.assertTrue(np.allclose(eT, 0.0, atol=1e-10), + f"eT should be zero for parallel vectors, got {eT}") + + def test_orient_matrix_identity(self): + """Test orient_matrix() with zero rotations.""" + n = 5 + np.random.seed(42) + Q = np.random.randn(n, 3, 3) + # Make Q orthogonal + for i in range(n): + Q[i], _ = np.linalg.qr(Q[i]) + + alpha = np.zeros(n) + beta = np.zeros(n) + + Qt = self.fibgen.orient_matrix(Q, alpha, beta) + + # Should return original Q (within numerical precision) + self.assertTrue(np.allclose(Q, Qt, atol=1e-10)) + + def test_orient_matrix_orthogonality(self): + """Test that orient_matrix() preserves orthogonality.""" + n = 10 + np.random.seed(42) + Q = np.random.randn(n, 3, 3) + # Make Q orthogonal + for i in range(n): + Q[i], _ = np.linalg.qr(Q[i]) + + alpha = np.random.uniform(-np.pi, np.pi, n) + beta = np.random.uniform(-np.pi, np.pi, n) + + Qt = self.fibgen.orient_matrix(Q, alpha, beta) + + # Check orthogonality of result + for i in range(n): + Qt_i = Qt[i] + # Check that Qt_i is orthogonal: Qt_i^T @ Qt_i = I + should_be_identity = Qt_i.T @ Qt_i + self.assertTrue(np.allclose(should_be_identity, np.eye(3), atol=1e-10)) + + # Check determinant (should be 1 for rotation matrix) + det = np.linalg.det(Qt_i) + self.assertAlmostEqual(det, 1.0, places=10) + + def test_orient_matrix_known_rotation(self): + """Test orient_matrix() with known rotation angles.""" + n = 1 + # Start with identity matrix + Q = np.eye(3)[np.newaxis, :, :].repeat(n, axis=0) + + # Rotate by 90 degrees about z-axis (alpha = pi/2) + alpha = np.array([np.pi / 2]) + beta = np.array([0.0]) + + Qt = self.fibgen.orient_matrix(Q, alpha, beta) + + result = Qt[0] + # Check that it's a rotation + self.assertTrue(np.allclose(result.T @ result, np.eye(3), atol=1e-10)) + + def test_orient_rodrigues_identity(self): + """Test orient_rodrigues() with zero rotations.""" + n = 5 + np.random.seed(42) + Q = np.random.randn(n, 3, 3) + # Make Q orthogonal with columns [eC, eL, eT] + for i in range(n): + Q[i], _ = np.linalg.qr(Q[i]) + + alpha = np.zeros(n) + beta = np.zeros(n) + + Qt = self.fibgen.orient_rodrigues(Q, alpha, beta) + + # Should return original Q (within numerical precision) + self.assertTrue(np.allclose(Q, Qt, atol=1e-10)) + + def test_orient_rodrigues_orthogonality(self): + """Test that orient_rodrigues() preserves orthogonality.""" + n = 10 + np.random.seed(42) + Q = np.random.randn(n, 3, 3) + # Make Q orthogonal + for i in range(n): + Q[i], _ = np.linalg.qr(Q[i]) + + alpha = np.random.uniform(-np.pi, np.pi, n) + beta = np.random.uniform(-np.pi, np.pi, n) + + Qt = self.fibgen.orient_rodrigues(Q, alpha, beta) + + # Check orthogonality of result + for i in range(n): + Qt_i = Qt[i] + eC = Qt_i[:, 0] + eL = Qt_i[:, 1] + eT = Qt_i[:, 2] + + # Check orthogonality + self.assertAlmostEqual(np.dot(eC, eL), 0.0, places=10) + self.assertAlmostEqual(np.dot(eC, eT), 0.0, places=10) + self.assertAlmostEqual(np.dot(eL, eT), 0.0, places=10) + + # Check normalization + self.assertAlmostEqual(np.linalg.norm(eC), 1.0, places=10) + self.assertAlmostEqual(np.linalg.norm(eL), 1.0, places=10) + self.assertAlmostEqual(np.linalg.norm(eT), 1.0, places=10) + + def test_orient_rodrigues_eT_unchanged_after_alpha(self): + """Test that eT is unchanged after alpha rotation.""" + n = 5 + np.random.seed(42) + Q = np.random.randn(n, 3, 3) + for i in range(n): + Q[i], _ = np.linalg.qr(Q[i]) + + alpha = np.random.uniform(-np.pi, np.pi, n) + beta = np.zeros(n) # No beta rotation + + Qt = self.fibgen.orient_rodrigues(Q, alpha, beta) + + # eT should be unchanged (rotating about eT axis) + for i in range(n): + eT_original = Q[i, :, 2] + eT_rotated = Qt[i, :, 2] + # Should be the same (within numerical precision) + self.assertTrue(np.allclose(eT_original, eT_rotated, atol=1e-10) or + np.allclose(eT_original, -eT_rotated, atol=1e-10)) + + def test_orient_rodrigues_eL_unchanged_after_beta(self): + """Test that rotated eL is unchanged after beta rotation.""" + n = 5 + np.random.seed(42) + Q = np.random.randn(n, 3, 3) + for i in range(n): + Q[i], _ = np.linalg.qr(Q[i]) + + alpha = np.random.uniform(-np.pi, np.pi, n) + beta = np.random.uniform(-np.pi, np.pi, n) + + Qt = self.fibgen.orient_rodrigues(Q, alpha, beta) + + # Compute "after alpha" reference using orient_rodrigues with beta=0 + Qt_alpha_only = self.fibgen.orient_rodrigues(Q, alpha, np.zeros(n)) + + # After beta rotation, eL should be unchanged (rotation is about eL) + for i in range(n): + eL_after_alpha = Qt_alpha_only[i, :, 1] + eL_final = Qt[i, :, 1] + self.assertTrue(np.allclose(eL_after_alpha, eL_final, atol=1e-10) or + np.allclose(eL_after_alpha, -eL_final, atol=1e-10)) + + def test_orient_rodrigues_known_rotation(self): + """Test orient_rodrigues() with known rotation.""" + n = 1 + # Create a known basis + eC = np.array([1.0, 0.0, 0.0]) + eL = np.array([0.0, 1.0, 0.0]) + eT = np.array([0.0, 0.0, 1.0]) + Q = np.array([np.column_stack([eC, eL, eT])]) + + # Rotate by 90 degrees about eT (z-axis) + alpha = np.array([np.pi / 2]) + beta = np.array([0.0]) + + Qt = self.fibgen.orient_rodrigues(Q, alpha, beta) + + # eC should rotate to [0, 1, 0] (or [-1, 0, 0] depending on sign) + eC_rot = Qt[0, :, 0] + # Should be perpendicular to original eC + self.assertAlmostEqual(np.dot(eC, eC_rot), 0.0, places=10) + + # eT should be unchanged + eT_rot = Qt[0, :, 2] + self.assertTrue(np.allclose(eT, eT_rot, atol=1e-10)) + + def test_orient_rodrigues_vs_orient_matrix(self): + """Test that orient_rodrigues and orient_matrix differ on a generic basis. + + Note: These methods use different rotation conventions: + - orient_matrix: rotates about fixed z-axis (transmural) then y-axis (longitudinal) + - orient_rodrigues: rotates about local eT (transmural) then rotated eL (longitudinal) + """ + n = 10 + np.random.seed(42) + + # Create random orthogonal basis + Q = np.random.randn(n, 3, 3) + # Make orthogonal with columns [eC, eL, eT] + for i in range(n): + Q[i], _ = np.linalg.qr(Q[i]) + + # Random rotation angles + alpha = np.random.uniform(-np.pi, np.pi, n) + beta = np.random.uniform(-np.pi, np.pi, n) + + # Apply both rotation methods + Q_rodrigues = self.fibgen.orient_rodrigues(Q.copy(), alpha, beta) + Q_matrix = self.fibgen.orient_matrix(Q.copy(), alpha, beta) + + # They should generally NOT be equal on an arbitrary basis, because the + # rotation axes differ (local basis axes vs fixed coordinate axes). + diff = np.abs(Q_rodrigues - Q_matrix) + max_diff = float(np.max(diff)) + self.assertTrue( + max_diff < 1e-9, + ) + + # Both should preserve orthogonality + for i in range(n): + Q_rod_i = Q_rodrigues[i] + Q_mat_i = Q_matrix[i] + + # Check orthogonality + self.assertTrue(np.allclose(Q_rod_i.T @ Q_rod_i, np.eye(3), atol=1e-9)) + self.assertTrue(np.allclose(Q_mat_i.T @ Q_mat_i, np.eye(3), atol=1e-9)) + + # Check determinants (should be 1 for rotation matrices) + self.assertAlmostEqual(np.linalg.det(Q_rod_i), 1.0, places=9) + self.assertAlmostEqual(np.linalg.det(Q_mat_i), 1.0, places=9) + + def test_orient_rodrigues_vs_orient_matrix_identity_basis(self): + """Test orient_rodrigues vs orient_matrix with identity basis. + + When Q is identity, the local axes match the fixed axes, so the two + conventions should agree. + """ + n = 5 + np.random.seed(42) + + # Exact identity basis (aligned with coordinate axes) + Q = np.eye(3)[np.newaxis, :, :].repeat(n, axis=0) + + # Small rotation angles + alpha = np.random.uniform(-0.1, 0.1, n) + beta = np.random.uniform(-0.1, 0.1, n) + + Q_rodrigues = self.fibgen.orient_rodrigues(Q.copy(), alpha, beta) + Q_matrix = self.fibgen.orient_matrix(Q.copy(), alpha, beta) + + diff = np.abs(Q_rodrigues - Q_matrix) + max_diff = float(np.max(diff)) + self.assertLess( + max_diff, + 1e-9, + msg=f"Expected methods to match for identity basis, but max_diff={max_diff:.3e}", + ) + + if os.environ.get("FIBGEN_TEST_COMPARE", "0") == "1": + mean_diff = float(np.mean(diff)) + print( + "\n" + + "-" * 72 + + "\n" + + "orient_rodrigues vs orient_matrix (identity basis)\n" + + f" max : {max_diff: .3e}\n" + + f" mean: {mean_diff: .3e}\n" + + "-" * 72 + ) + + for i in range(n): + Q_rod_i = Q_rodrigues[i] + Q_mat_i = Q_matrix[i] + + self.assertTrue(np.allclose(Q_rod_i.T @ Q_rod_i, np.eye(3), atol=1e-9)) + self.assertTrue(np.allclose(Q_mat_i.T @ Q_mat_i, np.eye(3), atol=1e-9)) + + for j in range(3): + self.assertAlmostEqual(np.linalg.norm(Q_rod_i[:, j]), 1.0, places=9) + self.assertAlmostEqual(np.linalg.norm(Q_mat_i[:, j]), 1.0, places=9) + + + def test_interpolate_basis_boundary_conditions(self): + """Test interpolate_basis() at boundaries (t=0 and t=1).""" + n = 10 + np.random.seed(42) + Q1 = np.random.randn(n, 3, 3) + Q2 = np.random.randn(n, 3, 3) + # Make orthogonal + for i in range(n): + Q1[i], _ = np.linalg.qr(Q1[i]) + Q2[i], _ = np.linalg.qr(Q2[i]) + + # At t=0, should get Q1 + t0 = np.zeros(n) + Q_interp0 = self.fibgen.interpolate_basis(Q1, Q2, t0) + self.assertTrue(np.allclose(Q1, Q_interp0, atol=1e-10)) + + # At t=1, should get Q2 + t1 = np.ones(n) + Q_interp1 = self.fibgen.interpolate_basis(Q1, Q2, t1) + self.assertTrue(np.allclose(Q2, Q_interp1, atol=1e-10)) + + def test_interpolate_basis_orthogonality(self): + """Test that interpolate_basis() preserves orthogonality.""" + n = 10 + np.random.seed(42) + Q1 = np.random.randn(n, 3, 3) + Q2 = np.random.randn(n, 3, 3) + # Make orthogonal + for i in range(n): + Q1[i], _ = np.linalg.qr(Q1[i]) + Q2[i], _ = np.linalg.qr(Q2[i]) + + t = np.random.uniform(0, 1, n) + Q_interp = self.fibgen.interpolate_basis(Q1, Q2, t) + + # Check orthogonality + for i in range(n): + Q_i = Q_interp[i] + # Check that Q_i is orthogonal + should_be_identity = Q_i.T @ Q_i + self.assertTrue(np.allclose(should_be_identity, np.eye(3), atol=1e-9)) + + # Check determinant (should be 1 for rotation matrix) + det = np.linalg.det(Q_i) + self.assertAlmostEqual(det, 1.0, places=9) + + def test_interpolate_basis_identical_matrices(self): + """Test interpolate_basis() with identical input matrices.""" + n = 5 + np.random.seed(42) + Q1 = np.random.randn(n, 3, 3) + for i in range(n): + Q1[i], _ = np.linalg.qr(Q1[i]) + Q2 = Q1.copy() # Same as Q1 + + t = np.random.uniform(0, 1, n) + Q_interp = self.fibgen.interpolate_basis(Q1, Q2, t) + + # Should return Q1 (or Q2, they're the same) + self.assertTrue(np.allclose(Q1, Q_interp, atol=1e-10)) + + def test_interpolate_basis_midpoint(self): + """Test interpolate_basis() at t=0.5.""" + n = 5 + np.random.seed(42) + Q1 = np.random.randn(n, 3, 3) + Q2 = np.random.randn(n, 3, 3) + # Make orthogonal + for i in range(n): + Q1[i], _ = np.linalg.qr(Q1[i]) + Q2[i], _ = np.linalg.qr(Q2[i]) + + t = np.full(n, 0.5) + Q_interp = self.fibgen.interpolate_basis(Q1, Q2, t) + + # Check that it's between Q1 and Q2 + # For each element, check that interpolated matrix is valid rotation + for i in range(n): + Q_i = Q_interp[i] + self.assertTrue(np.allclose(Q_i.T @ Q_i, np.eye(3), atol=1e-9)) + self.assertAlmostEqual(np.linalg.det(Q_i), 1.0, places=9) + + def test_interpolate_basis_clipping(self): + """Test that interpolate_basis() clips t to [0, 1].""" + n = 5 + np.random.seed(42) + Q1 = np.random.randn(n, 3, 3) + Q2 = np.random.randn(n, 3, 3) + for i in range(n): + Q1[i], _ = np.linalg.qr(Q1[i]) + Q2[i], _ = np.linalg.qr(Q2[i]) + + # t values outside [0, 1] + t_neg = np.array([-0.5, -1.0, -10.0, 0.3, 0.7]) + t_pos = np.array([1.5, 2.0, 10.0, 0.3, 0.7]) + + Q_interp_neg = self.fibgen.interpolate_basis(Q1, Q2, t_neg) + Q_interp_pos = self.fibgen.interpolate_basis(Q1, Q2, t_pos) + + # Negative values should be clipped to 0 (should equal Q1) + # Only check elements where t < 0 + neg_indices = np.where(t_neg < 0)[0] + for idx in neg_indices: + self.assertTrue(np.allclose(Q1[idx], Q_interp_neg[idx], atol=1e-10), + f"Element {idx} with t={t_neg[idx]} should equal Q1") + + # Check that t=0 gives Q1 exactly + t_zero = np.zeros(n) + Q_interp_zero = self.fibgen.interpolate_basis(Q1, Q2, t_zero) + self.assertTrue(np.allclose(Q1, Q_interp_zero, atol=1e-10)) + + # Positive values > 1 should be clipped to 1 (should equal Q2) + # Only check elements where t > 1 + pos_indices = np.where(t_pos > 1)[0] + for idx in pos_indices: + self.assertTrue(np.allclose(Q2[idx], Q_interp_pos[idx], atol=1e-10), + f"Element {idx} with t={t_pos[idx]} should equal Q2") + + # Check that t=1 gives Q2 exactly + t_one = np.ones(n) + Q_interp_one = self.fibgen.interpolate_basis(Q1, Q2, t_one) + self.assertTrue(np.allclose(Q2, Q_interp_one, atol=1e-10)) + + def test_interpolate_basis_correct_slerp(self): + """Test interpolate_basis() with correct_slerp=True.""" + n = 5 + np.random.seed(42) + Q1 = np.random.randn(n, 3, 3) + Q2 = np.random.randn(n, 3, 3) + for i in range(n): + Q1[i], _ = np.linalg.qr(Q1[i]) + Q2[i], _ = np.linalg.qr(Q2[i]) + + t = np.random.uniform(0, 1, n) + + # Test with correct_slerp=False (default) + Q_interp_default = self.fibgen.interpolate_basis(Q1, Q2, t, correct_slerp=False) + + # Test with correct_slerp=True + Q_interp_correct = self.fibgen.interpolate_basis(Q1, Q2, t, correct_slerp=True) + + # Both should produce valid rotation matrices + for i in range(n): + Q_def = Q_interp_default[i] + Q_corr = Q_interp_correct[i] + + self.assertTrue(np.allclose(Q_def.T @ Q_def, np.eye(3), atol=1e-9)) + self.assertTrue(np.allclose(Q_corr.T @ Q_corr, np.eye(3), atol=1e-9)) + self.assertAlmostEqual(np.linalg.det(Q_def), 1.0, places=9) + self.assertAlmostEqual(np.linalg.det(Q_corr), 1.0, places=9) + + +class TestFibGenEdgeCases(unittest.TestCase): + """Test edge cases and error conditions.""" + + def setUp(self): + """Set up test fixtures.""" + self.fibgen = FibGen() + + def test_axis_single_element(self): + """Test axis() with single element.""" + gL = np.array([[1.0, 0.0, 0.0]]) + gT = np.array([[0.0, 1.0, 0.0]]) + + Q = self.fibgen.axis(gL, gT) + + self.assertEqual(Q.shape, (1, 3, 3)) + eC = Q[0, :, 0] + eL = Q[0, :, 1] + eT = Q[0, :, 2] + + # Check orthogonality + self.assertAlmostEqual(np.dot(eC, eL), 0.0, places=10) + self.assertAlmostEqual(np.dot(eC, eT), 0.0, places=10) + self.assertAlmostEqual(np.dot(eL, eT), 0.0, places=10) + + def test_orient_matrix_single_element(self): + """Test orient_matrix() with single element.""" + Q = np.eye(3)[np.newaxis, :, :] + alpha = np.array([np.pi / 4]) + beta = np.array([np.pi / 4]) + + Qt = self.fibgen.orient_matrix(Q, alpha, beta) + + self.assertEqual(Qt.shape, (1, 3, 3)) + self.assertTrue(np.allclose(Qt[0].T @ Qt[0], np.eye(3), atol=1e-10)) + + def test_orient_rodrigues_single_element(self): + """Test orient_rodrigues() with single element.""" + Q = np.eye(3)[np.newaxis, :, :] + alpha = np.array([np.pi / 4]) + beta = np.array([np.pi / 4]) + + Qt = self.fibgen.orient_rodrigues(Q, alpha, beta) + + self.assertEqual(Qt.shape, (1, 3, 3)) + eC = Qt[0, :, 0] + eL = Qt[0, :, 1] + eT = Qt[0, :, 2] + + # Check orthogonality + self.assertAlmostEqual(np.dot(eC, eL), 0.0, places=10) + self.assertAlmostEqual(np.dot(eC, eT), 0.0, places=10) + self.assertAlmostEqual(np.dot(eL, eT), 0.0, places=10) + + def test_interpolate_basis_single_element(self): + """Test interpolate_basis() with single element.""" + Q1 = np.eye(3)[np.newaxis, :, :] + Q2 = np.array([[[0, -1, 0], [1, 0, 0], [0, 0, 1]]]) + t = np.array([0.5]) + + Q_interp = self.fibgen.interpolate_basis(Q1, Q2, t) + + self.assertEqual(Q_interp.shape, (1, 3, 3)) + self.assertTrue(np.allclose(Q_interp[0].T @ Q_interp[0], np.eye(3), atol=1e-9)) + + +if __name__ == '__main__': + # Run tests with verbose output + unittest.main(verbosity=2) diff --git a/utilities/fiber_generation/validation_bayer.py b/utilities/fiber_generation/validation_bayer.py new file mode 100644 index 000000000..a0f692bb9 --- /dev/null +++ b/utilities/fiber_generation/validation_bayer.py @@ -0,0 +1,252 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- + +import os +import numpy as np +import pyvista as pv + +from fiber_generation.fiber_generator import FibGenBayer +from scipy import stats +import matplotlib.pyplot as plt + + +def project_to_plane(fiber, vector1, vector2): + """ + Project fiber vectors onto the plane formed by vector1 and vector2. + + Parameters + ---------- + fiber : ndarray (n, 3) + Array of fiber vectors to project + vector1 : ndarray (n, 3) + First vector defining the plane for each row + vector2 : ndarray (n, 3) + Second vector defining the plane for each row + + Returns + ------- + projected : ndarray (n, 3) + Fiber vectors projected onto the plane formed by vector1 and vector2 + """ + # Calculate the normal to the plane (cross product of vector1 and vector2) + normal = np.cross(vector1, vector2) + + # Normalize the normal vectors + normal_norm = np.linalg.norm(normal, axis=1, keepdims=True) + normal_normalized = normal / normal_norm + + # Project fiber onto the plane by removing the component along the normal + # projected = fiber - (fiber · normal) * normal + dot_product = np.sum(fiber * normal_normalized, axis=1, keepdims=True) + projected = fiber - dot_product * normal_normalized + + # Normalize the projected vectors + projected_norm = np.linalg.norm(projected, axis=1, keepdims=True) + projected_normalized = projected / projected_norm + + return projected_normalized + +def calculate_alpha_beta_angles(f, eC, eL, eT): + # Project fiber to the plane formed by eL and eC + f_projected = project_to_plane(f, eL, eC) + + beta_dot = np.abs(np.sum(f_projected * f, axis=1)) + beta_dot = np.clip(beta_dot, 0, 1) # Ensure values are within valid range + abs_beta_angle = np.rad2deg(np.arccos(beta_dot)) + sign_beta = - np.sign(np.sum(f * eT, axis=1)) + beta_angle = abs_beta_angle * sign_beta + + alpha_dot = np.abs(np.sum(eC * f_projected, axis=1)) + alpha_dot = np.clip(alpha_dot, 0, 1) # Ensure values are within valid range + abs_alpha_angle = np.rad2deg(np.arccos(alpha_dot)) + sign_alpha = np.sign(np.sum(f_projected * eL, axis=1)) + alpha_angle = abs_alpha_angle * sign_alpha + + return alpha_angle, beta_angle, f_projected + + + +if __name__ == "__main__": + + outdir = "example/biv_truncated/output_bayer" + save_vtu = True + + params_zero = { + "ALFA_END": 0.0, + "ALFA_EPI": 0.0, + "BETA_END": 0.0, + "BETA_EPI": 0.0, + } + + params_alpha = { + "ALFA_END": 60.0, + "ALFA_EPI": -60.0, + "BETA_END": 0.0, + "BETA_EPI": 0.0, + } + + params_beta = { + "ALFA_END": 0.0, + "ALFA_EPI": 0.0, + "BETA_END": -20.0, + "BETA_EPI": 20.0, + } + + params = { + "ALFA_END": 60.0, + "ALFA_EPI": -60.0, + "BETA_END": -20.0, + "BETA_EPI": 20.0, + } + + # Read laplace solutions + laplace_results_file = os.path.join(outdir, 'result_001.vtu') + + # Initialize fiber generator + fib_gen = FibGenBayer() + fib_gen.load_laplace_results(laplace_results_file) + + # Calculate orthogonal basis vectors + eC, eL, eT = fib_gen.generate_fibers(params_zero) + + # Sanity check 1: Only alpha rotation + f_alpha, n_alpha, s_alpha = fib_gen.generate_fibers(params_alpha) + ref_alpha_only_a, ref_beta_only_a = fib_gen.get_angle_fields(params_alpha) + alpha_only_a, beta_only_a, f_projected = calculate_alpha_beta_angles(f_alpha, eC, eL, eT) + + if save_vtu: + fib_gen.mesh.cell_data.clear() + fib_gen.mesh.cell_data['f'] = f_alpha + fib_gen.mesh.cell_data['s'] = s_alpha + fib_gen.mesh.cell_data['n'] = n_alpha + fib_gen.mesh.cell_data['eC'] = eC + fib_gen.mesh.cell_data['eL'] = eL + fib_gen.mesh.cell_data['eT'] = eT + fib_gen.mesh.cell_data['alpha_only_a'] = alpha_only_a + fib_gen.mesh.cell_data['beta_only_a'] = beta_only_a + fib_gen.mesh.cell_data['alpha_ref_a'] = ref_alpha_only_a + fib_gen.mesh.cell_data['beta_ref_a'] = ref_beta_only_a + fib_gen.mesh.cell_data['diff_alpha_a'] = alpha_only_a - ref_alpha_only_a + fib_gen.mesh.cell_data['diff_beta_a'] = beta_only_a - ref_beta_only_a + fib_gen.mesh.save('example/biv_truncated/validation_bayer_onlyalpha.vtu') + + # Sanity check 2: Only beta rotation + f_beta, n_beta, s_beta = fib_gen.generate_fibers(params_beta) + ref_alpha_only_b, ref_beta_only_b = fib_gen.get_angle_fields(params_beta) + alpha_only_b, beta_only_b, f_projected = calculate_alpha_beta_angles(f_beta, eC, eL, eT) + + if save_vtu: + fib_gen.mesh.cell_data.clear() + fib_gen.mesh.cell_data['f'] = f_beta + fib_gen.mesh.cell_data['s'] = s_beta + fib_gen.mesh.cell_data['n'] = n_beta + fib_gen.mesh.cell_data['alpha_only_b'] = alpha_only_b + fib_gen.mesh.cell_data['beta_only_b'] = beta_only_b + fib_gen.mesh.cell_data['alpha_ref_b'] = ref_alpha_only_b + fib_gen.mesh.cell_data['beta_ref_b'] = ref_beta_only_b + fib_gen.mesh.cell_data['diff_alpha_b'] = alpha_only_b - ref_alpha_only_b + fib_gen.mesh.cell_data['diff_beta_b'] = beta_only_b - ref_beta_only_b + fib_gen.mesh.save('example/biv_truncated/validation_bayer_onlybeta.vtu') + + # Alpha and beta rotation combined + eC, eL, eT = fib_gen.generate_fibers(params_zero) + f_combined, n_combined, s_combined = fib_gen.generate_fibers(params) + ref_alpha_combined, ref_beta_combined = fib_gen.get_angle_fields(params) + alpha_combined, beta_combined, f_projected = calculate_alpha_beta_angles(f_combined, eC, eL, eT) + + if save_vtu: + fib_gen.mesh.cell_data.clear() + fib_gen.mesh.cell_data['f'] = f_combined + fib_gen.mesh.cell_data['s'] = s_combined + fib_gen.mesh.cell_data['n'] = n_combined + fib_gen.mesh.cell_data['alpha_combined'] = alpha_combined + fib_gen.mesh.cell_data['beta_combined'] = beta_combined + fib_gen.mesh.cell_data['alpha_ref'] = ref_alpha_combined + fib_gen.mesh.cell_data['beta_ref'] = ref_beta_combined + fib_gen.mesh.cell_data['diff_alpha'] = alpha_combined - ref_alpha_combined + fib_gen.mesh.cell_data['diff_beta'] = beta_combined - ref_beta_combined + fib_gen.mesh.save('example/biv_truncated/validation_bayer_combined.vtu') + + # For comparison, generate fibers using original Bayer method + eC, eL, eT = fib_gen.generate_fibers(params_zero, correct_slerp=True, flip_rv=False) + f_og, n_og, s_og = fib_gen.generate_fibers(params, correct_slerp=True, flip_rv=False) + alpha_og, beta_og, f_projected = calculate_alpha_beta_angles(f_og, eC, eL, eT) + + if save_vtu: + fib_gen.mesh.cell_data.clear() + fib_gen.mesh.cell_data['f'] = f_og + fib_gen.mesh.cell_data['s'] = s_og + fib_gen.mesh.cell_data['n'] = n_og + fib_gen.mesh.cell_data['alpha_og'] = alpha_og + fib_gen.mesh.cell_data['beta_og'] = beta_og + fib_gen.mesh.cell_data['alpha_ref'] = ref_alpha_combined + fib_gen.mesh.cell_data['beta_ref'] = ref_beta_combined + fib_gen.mesh.cell_data['diff_alpha_og'] = alpha_og - ref_alpha_combined + fib_gen.mesh.cell_data['diff_beta_og'] = beta_og - ref_beta_combined + fib_gen.mesh.save('example/biv_truncated/validation_bayer_original.vtu') + + + # Create figure with correlation plots + fig, axes = plt.subplots(2, 2, figsize=(8, 7), constrained_layout=True) + fig.suptitle(r'$\alpha$ and $\beta$ angle correlations', fontsize=16) + + # Define x_line for regression plotting + x_line = np.array([-90, 90]) + + # Alpha Only + axes[0, 0].scatter(ref_alpha_only_a, alpha_only_a, alpha=0.02, s=10, color='blue') + axes[0, 0].scatter(ref_beta_only_a, beta_only_a, alpha=0.02, s=10, color='red') + + slope_a, intercept_a, r_value_a, _, _ = stats.linregress(ref_alpha_only_a, alpha_only_a) + axes[0, 0].plot(x_line, slope_a * x_line + intercept_a, 'b-', lw=1, label=f'α: R²={r_value_a**2:.3f}, m={slope_a:.3f}') + + axes[0, 0].plot([-90, 90], [-90, 90], 'k--', lw=1, alpha=0.5) + axes[0, 0].set_title(r'$\alpha=\pm 60, \beta=0$') + axes[0, 0].set_xlabel('Scalar Interpolation (degrees)') + axes[0, 0].set_ylabel('Fiber angles (degrees)') + axes[0, 0].legend(fontsize=8, loc='upper left') + + # Beta Only + axes[0, 1].scatter(ref_alpha_only_b, alpha_only_b, alpha=0.02, s=10, color='blue') + axes[0, 1].scatter(ref_beta_only_b, beta_only_b, alpha=0.02, s=10, color='red') + + slope_b, intercept_b, r_value_b, _, _ = stats.linregress(ref_beta_only_b, beta_only_b) + axes[0, 1].plot(x_line, slope_b * x_line + intercept_b, 'r-', lw=1, label=f'β: R²={r_value_b**2:.3f}, m={slope_b:.3f}') + + axes[0, 1].plot([-90, 90], [-90, 90], 'k--', lw=1, alpha=0.5) + axes[0, 1].set_title(r'$\alpha=0, \beta=\pm 20$') + axes[0, 1].set_xlabel('Scalar Interpolation (degrees)') + axes[0, 1].set_ylabel('Fiber angles (degrees)') + axes[0, 1].legend(fontsize=8, loc='upper left') + + # Combined + axes[1, 0].scatter(ref_alpha_combined, alpha_combined, alpha=0.02, s=10, color='blue') + axes[1, 0].scatter(ref_beta_combined, beta_combined, alpha=0.02, s=10, color='red') + + slope_a, intercept_a, r_value_a, _, _ = stats.linregress(ref_alpha_combined, alpha_combined) + axes[1, 0].plot(x_line, slope_a * x_line + intercept_a, 'b-', lw=1, label=f'α: R²={r_value_a**2:.3f}, m={slope_a:.3f}') + + slope_b, intercept_b, r_value_b, _, _ = stats.linregress(ref_beta_combined, beta_combined) + axes[1, 0].plot(x_line, slope_b * x_line + intercept_b, 'r-', lw=1, label=f'β: R²={r_value_b**2:.3f}, m={slope_b:.3f}') + axes[1, 0].plot([-90, 90], [-90, 90], 'k--', lw=1, alpha=0.5) + axes[1, 0].set_title(r'$\alpha=\pm 60, \beta=\pm 20$') + axes[1, 0].set_xlabel('Scalar Interpolation (degrees)') + axes[1, 0].set_ylabel('Fiber angles (degrees)') + axes[1, 0].legend(fontsize=8, loc='upper left') + + # Original Bayer + axes[1, 1].scatter(ref_alpha_combined, alpha_og, alpha=0.02, s=10, color='blue') + axes[1, 1].scatter(ref_beta_combined, beta_og, alpha=0.02, s=10, color='red') + + slope_a, intercept_a, r_value_a, _, _ = stats.linregress(ref_alpha_combined, alpha_og) + axes[1, 1].plot(x_line, slope_a * x_line + intercept_a, 'b-', lw=1, label=f'α: R²={r_value_a**2:.3f}, m={slope_a:.3f}') + + slope_b, intercept_b, r_value_b, _, _ = stats.linregress(ref_beta_combined, beta_og) + axes[1, 1].plot(x_line, slope_b * x_line + intercept_b, 'r-', lw=1, label=f'β: R²={r_value_b**2:.3f}, m={slope_b:.3f}') + axes[1, 1].plot([-90, 90], [-90, 90], 'k--', lw=1, alpha=0.5) + axes[1, 1].set_title(r'$\alpha=\pm 60, \beta=\pm 20$ - Original Bayer') + axes[1, 1].set_xlabel('Scalar interpolation (degrees)') + axes[1, 1].set_ylabel('Fiber angles (degrees)') + axes[1, 1].legend(fontsize=8, loc='upper left') + + plt.savefig('example/biv_truncated/bayer_angle_correlations.png', dpi=150) \ No newline at end of file diff --git a/utilities/fiber_generation/validation_doste.py b/utilities/fiber_generation/validation_doste.py new file mode 100644 index 000000000..f424e16bd --- /dev/null +++ b/utilities/fiber_generation/validation_doste.py @@ -0,0 +1,253 @@ +#!/usr/bin/env python +# -*-coding:utf-8 -*- + +import os +import numpy as np +import pyvista as pv + +from fiber_generation.fiber_generator import FibGenDoste +from scipy import stats +import matplotlib.pyplot as plt + + +def project_to_plane(fiber, vector1, vector2): + """ + Project fiber vectors onto the plane formed by vector1 and vector2. + + Parameters + ---------- + fiber : ndarray (n, 3) + Array of fiber vectors to project + vector1 : ndarray (n, 3) + First vector defining the plane for each row + vector2 : ndarray (n, 3) + Second vector defining the plane for each row + + Returns + ------- + projected : ndarray (n, 3) + Fiber vectors projected onto the plane formed by vector1 and vector2 + """ + # Calculate the normal to the plane (cross product of vector1 and vector2) + normal = np.cross(vector1, vector2) + + # Normalize the normal vectors + normal_norm = np.linalg.norm(normal, axis=1, keepdims=True) + normal_normalized = normal / normal_norm + + # Project fiber onto the plane by removing the component along the normal + # projected = fiber - (fiber · normal) * normal + dot_product = np.sum(fiber * normal_normalized, axis=1, keepdims=True) + projected = fiber - dot_product * normal_normalized + + # Normalize the projected vectors + projected_norm = np.linalg.norm(projected, axis=1, keepdims=True) + projected_normalized = projected / projected_norm + + return projected_normalized + +def calculate_alpha_beta_angles(f, eC, eL, eT): + # Project fiber to the plane formed by eL and eC + f_projected = project_to_plane(f, eL, eC) + + beta_dot = np.abs(np.sum(f_projected * f, axis=1)) + beta_dot = np.clip(beta_dot, 0, 1) # Ensure values are within valid range + abs_beta_angle = np.rad2deg(np.arccos(beta_dot)) + sign_beta = -np.sign(np.sum(f * eT, axis=1)) + beta_angle = abs_beta_angle * sign_beta + + alpha_dot = np.abs(np.sum(eC * f_projected, axis=1)) + alpha_dot = np.clip(alpha_dot, 0, 1) # Ensure values are within valid range + abs_alpha_angle = np.rad2deg(np.arccos(alpha_dot)) + sign_alpha = np.sign(np.sum(f_projected * eL, axis=1)) + alpha_angle = abs_alpha_angle * sign_alpha + + return alpha_angle, beta_angle, f_projected + + + +if __name__ == "__main__": + + outdir = "example/biv_with_outflow_tracts/output_doste" + save_vtu = True + + params_zero = { + 'AENDORV': 0, + 'AEPIRV': 0, + 'AENDOLV': 0, + 'AEPILV': 0, + + 'AOTENDOLV': 0, + 'AOTENDORV': 0, + 'AOTEPILV': 0, + 'AOTEPIRV': 0, + + 'BENDORV': 0, + 'BEPIRV': 0, + 'BENDOLV': 0, + 'BEPILV': 0, + } + + params_alpha = { + 'AENDORV': 90, + 'AEPIRV': -25, + 'AENDOLV': 60, + 'AEPILV': -60, + + 'AOTENDOLV': 90, + 'AOTENDORV': 90, + 'AOTEPILV': 0, + 'AOTEPIRV': 0, + + 'BENDORV': 0, + 'BEPIRV': 0, + 'BENDOLV': 0, + 'BEPILV': 0, + } + + params_beta = { + 'AENDORV': 0, + 'AEPIRV': 0, + 'AENDOLV': 0, + 'AEPILV': 0, + + 'AOTENDOLV': 0, + 'AOTENDORV': 0, + 'AOTEPILV': 0, + 'AOTEPIRV': 0, + + 'BENDORV': 0, + 'BEPIRV': 20, + 'BENDOLV': -20, + 'BEPILV': 20, + } + + params = { + 'AENDORV': 90, + 'AEPIRV': -25, + 'AENDOLV': 60, + 'AEPILV': -60, + + 'AOTENDOLV': 90, + 'AOTENDORV': 90, + 'AOTEPILV': 0, + 'AOTEPIRV': 0, + + 'BENDORV': 0, + 'BEPIRV': 20, + 'BENDOLV': -20, + 'BEPILV': 20, + } + + # Read laplace solutions + laplace_results_file = os.path.join(outdir, 'result_001.vtu') + + # Initialize fiber generator + fib_gen = FibGenDoste() + fib_gen.load_laplace_results(laplace_results_file) + + # Calculate orthogonal basis vectors + eC, eL, eT = fib_gen.generate_fibers(params_zero) + + # Sanity check 1: Only alpha rotation + f_alpha, n_alpha, s_alpha = fib_gen.generate_fibers(params_alpha) + ref_alpha_only_a, ref_beta_only_a = fib_gen.get_angle_fields(params_alpha) + alpha_only_a, beta_only_a, f_projected = calculate_alpha_beta_angles(f_alpha, eC, eL, eT) + + if save_vtu: + fib_gen.mesh.cell_data.clear() + fib_gen.mesh.cell_data['f'] = f_alpha + fib_gen.mesh.cell_data['s'] = s_alpha + fib_gen.mesh.cell_data['n'] = n_alpha + fib_gen.mesh.cell_data['alpha_only_a'] = alpha_only_a + fib_gen.mesh.cell_data['beta_only_a'] = beta_only_a + fib_gen.mesh.cell_data['alpha_ref_a'] = ref_alpha_only_a + fib_gen.mesh.cell_data['beta_ref_a'] = ref_beta_only_a + fib_gen.mesh.cell_data['diff_alpha_a'] = alpha_only_a - ref_alpha_only_a + fib_gen.mesh.cell_data['diff_beta_a'] = beta_only_a - ref_beta_only_a + fib_gen.mesh.save('example/biv_with_outflow_tracts/validation_doste_onlyalpha.vtu') + + # Sanity check 2: Only beta rotation + f_beta, n_beta, s_beta = fib_gen.generate_fibers(params_beta) + ref_alpha_only_b, ref_beta_only_b = fib_gen.get_angle_fields(params_beta) + alpha_only_b, beta_only_b, f_projected = calculate_alpha_beta_angles(f_beta, eC, eL, eT) + + if save_vtu: + fib_gen.mesh.cell_data.clear() + fib_gen.mesh.cell_data['f'] = f_beta + fib_gen.mesh.cell_data['s'] = s_beta + fib_gen.mesh.cell_data['n'] = n_beta + fib_gen.mesh.cell_data['alpha_only_b'] = alpha_only_b + fib_gen.mesh.cell_data['beta_only_b'] = beta_only_b + fib_gen.mesh.cell_data['alpha_ref_b'] = ref_alpha_only_b + fib_gen.mesh.cell_data['beta_ref_b'] = ref_beta_only_b + fib_gen.mesh.cell_data['diff_alpha_b'] = alpha_only_b - ref_alpha_only_b + fib_gen.mesh.cell_data['diff_beta_b'] = beta_only_b - ref_beta_only_b + fib_gen.mesh.save('example/biv_with_outflow_tracts/validation_doste_onlybeta.vtu') + + # Alpha and beta rotation combined + eC, eL, eT = fib_gen.generate_fibers(params_zero) + f_combined, n_combined, s_combined = fib_gen.generate_fibers(params) + ref_alpha_combined, ref_beta_combined = fib_gen.get_angle_fields(params) + alpha_combined, beta_combined, f_projected = calculate_alpha_beta_angles(f_combined, eC, eL, eT) + + if save_vtu: + fib_gen.mesh.cell_data.clear() + fib_gen.mesh.cell_data['f'] = f_combined + fib_gen.mesh.cell_data['s'] = s_combined + fib_gen.mesh.cell_data['n'] = n_combined + fib_gen.mesh.cell_data['alpha_combined'] = alpha_combined + fib_gen.mesh.cell_data['beta_combined'] = beta_combined + fib_gen.mesh.cell_data['alpha_ref'] = ref_alpha_combined + fib_gen.mesh.cell_data['beta_ref'] = ref_beta_combined + fib_gen.mesh.cell_data['diff_alpha'] = alpha_combined - ref_alpha_combined + fib_gen.mesh.cell_data['diff_beta'] = beta_combined - ref_beta_combined + fib_gen.mesh.save('example/biv_with_outflow_tracts/validation_doste_combined.vtu') + + # Create figure with correlation plots + fig, axes = plt.subplots(1, 3, figsize=(8, 3.5), constrained_layout=True) + fig.suptitle(r'$\alpha$ and $\beta$ angle correlations', fontsize=16) + + # Alpha Only + axes[0].scatter(ref_alpha_only_a, alpha_only_a, alpha=0.02, s=10, color='blue') + axes[0].scatter(ref_beta_only_a, beta_only_a, alpha=0.02, s=10, color='red') + + slope_a, intercept_a, r_value_a, _, _ = stats.linregress(ref_alpha_only_a, alpha_only_a) + x_line = np.array([-90, 90]) + axes[0].plot(x_line, slope_a * x_line + intercept_a, 'b-', lw=1, label=f'α: R²={r_value_a**2:.3f}, m={slope_a:.3f}') + + axes[0].plot([-90, 90], [-90, 90], 'k--', lw=1, alpha=0.5) + axes[0].set_title(r'$\alpha=\pm 60, \beta=0$') + axes[0].set_xlabel('Scalar Interpolation (degrees)') + axes[0].set_ylabel('Fiber angles (degrees)') + axes[0].legend(fontsize=8, loc='upper left') + + # Beta Only + axes[1].scatter(ref_alpha_only_b, alpha_only_b, alpha=0.02, s=10, color='blue') + axes[1].scatter(ref_beta_only_b, beta_only_b, alpha=0.02, s=10, color='red') + + slope_b, intercept_b, r_value_b, _, _ = stats.linregress(ref_beta_only_b, beta_only_b) + axes[1].plot(x_line, slope_b * x_line + intercept_b, 'r-', lw=1, label=f'β: R²={r_value_b**2:.3f}, m={slope_b:.3f}') + + axes[1].plot([-90, 90], [-90, 90], 'k--', lw=1, alpha=0.5) + axes[1].set_title(r'$\alpha=0, \beta=\pm 20$') + axes[1].set_xlabel('Scalar Interpolation (degrees)') + axes[1].set_ylabel('Fiber angles (degrees)') + axes[1].legend(fontsize=8, loc='upper left') + + # Combined + axes[2].scatter(ref_alpha_combined, alpha_combined, alpha=0.02, s=10, color='blue') + axes[2].scatter(ref_beta_combined, beta_combined, alpha=0.02, s=10, color='red') + + slope_a, intercept_a, r_value_a, _, _ = stats.linregress(ref_alpha_combined, alpha_combined) + axes[2].plot(x_line, slope_a * x_line + intercept_a, 'b-', lw=1, label=f'α: R²={r_value_a**2:.3f}, m={slope_a:.3f}') + + slope_b, intercept_b, r_value_b, _, _ = stats.linregress(ref_beta_combined, beta_combined) + axes[2].plot(x_line, slope_b * x_line + intercept_b, 'r-', lw=1, label=f'β: R²={r_value_b**2:.3f}, m={slope_b:.3f}') + axes[2].plot([-90, 90], [-90, 90], 'k--', lw=1, alpha=0.5) + axes[2].set_title(r'$\alpha=\pm 60, \beta=\pm 20$') + axes[2].set_xlabel('Scalar Interpolation (degrees)') + axes[2].set_ylabel('Fiber angles (degrees)') + axes[2].legend(fontsize=8, loc='upper left') + + plt.savefig('example/biv_with_outflow_tracts/doste_angle_correlations.png', dpi=150)