From ab4c3773dc9ed9a53abbd68a24479de76078ac19 Mon Sep 17 00:00:00 2001 From: josephbrunet Date: Wed, 17 Jun 2026 14:54:03 +0200 Subject: [PATCH 1/6] add video --- pyproject.toml | 2 + .../scripts/visualize_streamlines.py | 113 ++++++++- .../fury_plotting_streamlines.py | 230 ++++++++++++++---- src/cardiotensor/visualization/streamlines.py | 129 +++++++--- 4 files changed, 389 insertions(+), 85 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 4e7b196..2616710 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -36,6 +36,7 @@ dependencies = [ "opencv-python-headless", "pandas", "psutil", + "pyvista", "scikit-image", "structure-tensor>=0.3.0", "tifffile", @@ -77,6 +78,7 @@ git = [ cardio-analysis = "cardiotensor.scripts.gui_analysis_tool:script" cardio-analysis-histo = "cardiotensor.scripts.analysis_histograms:script" cardio-analysis-streamlines = "cardiotensor.scripts.analysis_streamlines:script" +cardio-create-movie = "cardiotensor.scripts.create_movie:script" cardio-generate-streamlines = "cardiotensor.scripts.generate_streamlines:script" cardio-tensor = "cardiotensor.scripts.compute_orientation:script" cardio-tensor-slurm = "cardiotensor.scripts.slurm_launcher:script" diff --git a/src/cardiotensor/scripts/visualize_streamlines.py b/src/cardiotensor/scripts/visualize_streamlines.py index faa0c75..1340149 100644 --- a/src/cardiotensor/scripts/visualize_streamlines.py +++ b/src/cardiotensor/scripts/visualize_streamlines.py @@ -81,7 +81,7 @@ def script() -> None: "--line-width", type=float, default=4.0, - help="Tube width", + help="Tube radius in tube mode, or screen line width in line mode", ) parser.add_argument( "--subsample", @@ -131,7 +131,7 @@ def script() -> None: parser.add_argument( "--no-interactive", action="store_true", - help="Disable interactive window, useful when only saving a screenshot", + help="Disable interactive window, useful when only saving a screenshot or video", ) parser.add_argument( "--screenshot", @@ -139,11 +139,33 @@ def script() -> None: default=None, help="Path to save a PNG screenshot, no file is saved if omitted", ) + parser.add_argument( + "--video", + type=str, + default=None, + metavar="PATH", + help=( + "Path for an orbit video (e.g. orbit.mp4 or orbit.gif). " + "Implies --no-interactive. FURY MP4 uses OpenCV; PyVista MP4 may require imageio-ffmpeg." + ), + ) + parser.add_argument( + "--video-fps", + type=int, + default=30, + help="Frames per second for the orbit video (default: 30)", + ) + parser.add_argument( + "--video-frames", + type=int, + default=120, + help="Total number of frames in the orbit video (default: 120 = 4 s at 30 fps)", + ) parser.add_argument( "--spline-subdiv", type=int, - default=16, - help="Spline subdivisions for tube rendering (higher = smoother, heavier)", + default=2, + help="Spline/display subdivisions for smoother streamlines (higher = smoother, heavier)", ) parser.add_argument("--width", type=int, default=800, help="Window width in pixels") parser.add_argument( @@ -155,6 +177,61 @@ def script() -> None: default=None, help="Colormap name. By default, choose one from --color-by: helix_angle for HA/IA, viridis for EL, hsv for AZ.", ) + parser.add_argument( + "--backend", + choices=("fury", "pyvista"), + default="fury", + help="Rendering backend", + ) + parser.add_argument( + "--mode", + choices=("tube", "line"), + default="tube", + help="Render explicit tube geometry or faster line geometry", + ) + parser.add_argument( + "--background-color", + type=str, + default=None, + help="Renderer background color. Default: white for PyVista, black for FURY", + ) + parser.add_argument( + "--tube-sides", + type=int, + default=9, + help="Number of sides for tube geometry", + ) + parser.add_argument( + "--opacity", + type=float, + default=1.0, + help="Streamline opacity for the PyVista backend", + ) + parser.add_argument( + "--shadows", + action="store_true", + help="Enable PyVista shadow rendering at startup", + ) + parser.add_argument( + "--no-shadows", + action="store_true", + help=argparse.SUPPRESS, + ) + parser.add_argument( + "--hide-axes", + action="store_true", + help="Hide the PyVista orientation axes", + ) + parser.add_argument( + "--show-bounds", + action="store_true", + help="Show the PyVista bounds grid at startup", + ) + parser.add_argument( + "--hide-bounds", + action="store_true", + help=argparse.SUPPRESS, + ) args = parser.parse_args() @@ -163,7 +240,8 @@ def script() -> None: # List available color-by and exit if requested available_dpp = _discover_trk_color_fields(trk_path) - available = ["elevation"] + available_dpp + computed_color_options = ["elevation", "azimuth", "az", "el"] + available = computed_color_options + available_dpp if args.list_color_by: print("Available color-by options:") for k in available: @@ -185,12 +263,14 @@ def script() -> None: print(f"Auto color-by resolved to: {color_by}") # Validate choice - if color_by != "elevation" and color_by not in dpp_lower: + if color_by not in computed_color_options and color_by not in dpp_lower: print( f"Requested color-by '{args.color_by}' not found. Available: {', '.join(available)}" ) sys.exit(2) + selected_color_by = dpp_lower[color_by] if color_by in dpp_lower else color_by + # Choose colormap only when explicitly requested; otherwise the visualizer # picks the default based on color_by. chosen_cmap = None @@ -215,24 +295,35 @@ def script() -> None: tuple(args.crop_z or [-float("inf"), float("inf")]), ) + # --video implies off-screen + is_interactive = not args.no_interactive and not args.video + # Call the visualizer - # Always render as tubes as requested earlier visualize_streamlines( streamlines_file=trk_path, - color_by=( - color_by if color_by == "elevation" else dpp_lower[color_by] - ), # pass original key if stored + color_by=selected_color_by, line_width=args.line_width, subsample_factor=args.subsample, filter_min_len=args.min_length, downsample_factor=args.downsample_factor, max_streamlines=args.max_streamlines, crop_bounds=crop_bounds, - interactive=not args.no_interactive, + interactive=is_interactive, screenshot_path=args.screenshot, + video_path=args.video, + video_fps=args.video_fps, + video_frames=args.video_frames, window_size=(args.width, args.height), colormap=chosen_cmap, spline_subdiv=args.spline_subdiv, + backend=args.backend, + mode=args.mode, + background_color=args.background_color, + tube_sides=args.tube_sides, + pyvista_opacity=args.opacity, + pyvista_show_axes=not args.hide_axes, + pyvista_show_bounds=args.show_bounds and not args.hide_bounds, + pyvista_shadows=args.shadows and not args.no_shadows, ) diff --git a/src/cardiotensor/visualization/fury_plotting_streamlines.py b/src/cardiotensor/visualization/fury_plotting_streamlines.py index 2668981..7ff1d2a 100644 --- a/src/cardiotensor/visualization/fury_plotting_streamlines.py +++ b/src/cardiotensor/visualization/fury_plotting_streamlines.py @@ -1,7 +1,11 @@ #!/usr/bin/env python3 from __future__ import annotations +import datetime +import inspect import random +import tempfile +from pathlib import Path import fury import matplotlib.pyplot as plt @@ -41,6 +45,71 @@ def downsample_streamline(streamline: np.ndarray, factor: int = 2) -> np.ndarray return streamline if len(streamline) < 3 or factor <= 1 else streamline[::factor] +def _supported_actor_kwargs(func, **kwargs): + try: + params = inspect.signature(func).parameters + except (TypeError, ValueError): + return kwargs + + if any(param.kind == inspect.Parameter.VAR_KEYWORD for param in params.values()): + return kwargs + return {name: value for name, value in kwargs.items() if name in params} + + +def _write_frame_sequence_to_video( + frame_paths: list[Path], output_path: Path, fps: int +) -> None: + if not frame_paths: + raise ValueError("No FURY frames were recorded for the video.") + + output_path.parent.mkdir(parents=True, exist_ok=True) + suffix = output_path.suffix.lower() + if suffix == ".gif": + try: + import imageio.v2 as imageio + except ImportError as err: + raise RuntimeError( + "GIF export requires imageio: pip install imageio" + ) from err + + with imageio.get_writer( + str(output_path), mode="I", duration=1.0 / fps + ) as writer: + for frame_path in frame_paths: + writer.append_data(imageio.imread(frame_path)) + return + + try: + import cv2 + except ImportError as err: + raise RuntimeError( + "MP4 export requires OpenCV: pip install opencv-python" + ) from err + + first = cv2.imread(str(frame_paths[0]), cv2.IMREAD_COLOR) + if first is None: + raise RuntimeError(f"Could not read recorded frame: {frame_paths[0]}") + + height, width = first.shape[:2] + fourcc = cv2.VideoWriter_fourcc(*"mp4v") + writer = cv2.VideoWriter( + str(output_path), fourcc, float(fps), (width, height), isColor=True + ) + if not writer.isOpened(): + raise RuntimeError(f"Could not open video writer for {output_path}") + + try: + for frame_path in frame_paths: + frame = cv2.imread(str(frame_path), cv2.IMREAD_COLOR) + if frame is None: + raise RuntimeError(f"Could not read recorded frame: {frame_path}") + if frame.shape[:2] != (height, width): + frame = cv2.resize(frame, (width, height), interpolation=cv2.INTER_AREA) + writer.write(frame) + finally: + writer.release() + + def _compute_streamline_bounds( streamlines_xyz: list[np.ndarray], ) -> tuple[np.ndarray, np.ndarray]: @@ -122,6 +191,7 @@ def __init__( lut, background_color="black", spline_subdiv=16, + tube_sides: int = 9, color_range: tuple[float, float] | None = None, color_label: str = "Angle", ): @@ -140,6 +210,7 @@ def __init__( # thickness state self.linewidth = max(1.0, float(line_width)) # used for both line and tube self.spline_subdiv = spline_subdiv + self.tube_sides = max(3, int(tube_sides)) self.scale_bar = None self.scale_bar_on = False @@ -189,10 +260,8 @@ def _build_pipeline(self): linewidth=self.linewidth, spline_subdiv=self.spline_subdiv, lookup_colormap=self.lut, - tube_sides=20, - # lod=False, # <— - # lod_points=20000, # optional, tune - # lod_points_size=2, # optional, tune + tube_sides=self.tube_sides, + **_supported_actor_kwargs(actor.streamtube, lod=False), ) else: self.actor0 = actor.line( @@ -200,6 +269,7 @@ def _build_pipeline(self): colors=self.flat_vals, # scalars linewidth=self.linewidth, lookup_colormap=self.lut, + **_supported_actor_kwargs(actor.line, lod=False), ) self.scene.add(self.actor0) @@ -211,8 +281,8 @@ def _build_pipeline(self): colors=self.flat_vals, # pass scalars linewidth=1.0, # very light lookup_colormap=self.lut, - # lod=False, # make it deterministic fake_tube=True, # a bit of shading to hint tubes + **_supported_actor_kwargs(actor.line, lod=False), ) self.actor_fast.SetVisibility(False) self.scene.add(self.actor_fast) @@ -238,10 +308,10 @@ def _build_pipeline(self): def _style_streamline_actor(self): prop = self.actor0.GetProperty() prop.SetInterpolationToPhong() - prop.SetAmbient(0.1) - prop.SetDiffuse(0.95) - prop.SetSpecular(0.25) - prop.SetSpecularPower(12) + prop.SetAmbient(0.45) # raised from 0.1 → brighter unlit sides, less dark + prop.SetDiffuse(0.75) + prop.SetSpecular(0.15) + prop.SetSpecularPower(10) def _add_origin_marker(self): bounds_size = np.array( @@ -347,11 +417,9 @@ def _rebuild_unclipped_actor(self): colors=self.flat_vals, linewidth=self.linewidth, spline_subdiv=self.spline_subdiv, - tube_sides=20, + tube_sides=self.tube_sides, lookup_colormap=self.lut, - # lod=False, # <— - # lod_points=20000, # optional, tune - # lod_points_size=2, # optional, tune + **_supported_actor_kwargs(actor.streamtube, lod=False), ) else: self.actor0 = actor.line( @@ -359,6 +427,7 @@ def _rebuild_unclipped_actor(self): colors=self.flat_vals, linewidth=self.linewidth, lookup_colormap=self.lut, + **_supported_actor_kwargs(actor.line, lod=False), ) self.scene.add(self.actor0) @@ -420,16 +489,13 @@ def _on_keypress(self, obj, evt): self._render_now() elif key == "p": - import datetime - import os - ts = datetime.datetime.now().strftime("%Y%m%d_%H%M%S") - out_path = os.path.abspath(f"view_{ts}_hires.png") + out_path = Path.cwd() / f"view_{ts}_hires.png" try: highres_size = (2000, 2000) fury.window.record( scene=self.scene, - out_path=out_path, + out_path=str(out_path), size=highres_size, reset_camera=False, ) @@ -439,6 +505,10 @@ def _on_keypress(self, obj, evt): except Exception as e: print(f"Failed to save screenshot: {e}") + elif key == "v": + ts = datetime.datetime.now().strftime("%Y%m%d_%H%M%S") + self._record_orbit(Path.cwd() / f"fury_orbit_{ts}.mp4", 120, 30) + elif key == "r": self.plane_rep.SetOrigin(*self.center) self.plane_rep.SetNormal(1, 0, 0) @@ -473,35 +543,91 @@ def _lod_on(self, obj=None, evt=None): # Apply clipping only if it was previously enabled if self.clipping_active: - self.actor_fast.GetMapper().AddClippingPlane( - self.plane_fn - ) # Apply clipping to the fast actor + mapper = self.actor_fast.GetMapper() + mapper.RemoveAllClippingPlanes() + mapper.AddClippingPlane(self.plane_fn) except Exception: pass self._render_now() def _lod_off(self, obj=None, evt=None): """Switch back to full-res actor, applying clipping if it was previously enabled.""" - # Rebuild actor to ensure it shows at full resolution with clipping applied - self._rebuild_unclipped_actor() - try: self.actor_fast.SetVisibility(False) # Hide the fast actor self.actor0.SetVisibility(True) # Show the full-res actor # Apply clipping to the full-res actor if it was previously enabled if self.clipping_active: - self.actor0.GetMapper().AddClippingPlane( - self.plane_fn - ) # Reapply clipping + mapper = self.actor0.GetMapper() + mapper.RemoveAllClippingPlanes() + mapper.AddClippingPlane(self.plane_fn) except Exception: pass self._render_now() + def _set_default_camera(self) -> None: + self.scene.reset_camera() + self.scene.azimuth(15) + self.scene.elevation(10) + self.scene.zoom(1.1) + + def _record_orbit( + self, video_path: str | Path, video_frames: int, video_fps: int + ) -> None: + out_path = Path(video_path).resolve() + n_frames = max(1, int(video_frames)) + fps = max(1, int(video_fps)) + step_degrees = 360.0 / float(n_frames) + + try: + self.actor_fast.SetVisibility(False) + self.actor0.SetVisibility(True) + except Exception: + pass + + self._set_default_camera() + print( + f"Recording FURY orbit video: {n_frames} frames at {fps} fps -> {out_path}" + ) + + with tempfile.TemporaryDirectory(prefix="cardiotensor_fury_orbit_") as tmpdir: + tmpdir_path = Path(tmpdir) + frame_paths = [] + try: + from tqdm import tqdm + + frame_iter = tqdm( + range(n_frames), total=n_frames, unit="frame", desc="Recording" + ) + except ImportError: + frame_iter = range(n_frames) + + for frame_idx in frame_iter: + frame_path = tmpdir_path / f"frame_{frame_idx:05d}.png" + fury.window.record( + scene=self.scene, + out_path=str(frame_path), + size=self.window_size, + reset_camera=False, + ) + frame_paths.append(frame_path) + self.scene.azimuth(step_degrees) + + _write_frame_sequence_to_video(frame_paths, out_path, fps) + + print(f"Saved FURY orbit video to {out_path}") + # --------------------------- # main entry # --------------------------- - def run(self, interactive: bool, screenshot_path: str | None): + def run( + self, + interactive: bool, + screenshot_path: str | None, + video_path: str | None = None, + video_fps: int = 30, + video_frames: int = 120, + ): if interactive: self.showm = window.ShowManager( scene=self.scene, size=self.window_size, reset_camera=False @@ -540,23 +666,29 @@ def run(self, interactive: bool, screenshot_path: str | None): self.showm.iren.AddObserver("KeyPressEvent", self._on_keypress) - self.scene.reset_camera() - - self.scene.azimuth(15) - self.scene.elevation(10) - self.scene.zoom(1.1) + self._set_default_camera() print( - "Keys: O toggle plane, H hide gizmo, I flip side, R reset plane, +/- thickness, B background, S scale bar, P save PNG" + "Keys: O toggle plane, H hide gizmo, I flip side, R reset plane, +/- thickness, B background, S scale bar, P save PNG, V record orbit" ) self.showm.start() else: - if not screenshot_path: - raise ValueError("Must specify screenshot_path when interactive=False.") - self.scene.reset_camera() - fury.window.record( - scene=self.scene, out_path=screenshot_path, size=self.window_size - ) + if not screenshot_path and not video_path: + raise ValueError( + "Must specify screenshot_path or video_path when interactive=False." + ) + self._set_default_camera() + if screenshot_path: + screenshot = Path(screenshot_path).resolve() + screenshot.parent.mkdir(parents=True, exist_ok=True) + fury.window.record( + scene=self.scene, + out_path=str(screenshot), + size=self.window_size, + reset_camera=False, + ) + if video_path: + self._record_orbit(video_path, video_frames, video_fps) # =========================== @@ -569,16 +701,21 @@ def show_streamlines( line_width: float = 4, interactive: bool = True, screenshot_path: str | None = None, + video_path: str | None = None, + video_fps: int = 30, + video_frames: int = 120, window_size: tuple[int, int] = (800, 800), downsample_factor: int = 2, max_streamlines: int | None = None, filter_min_len: int | None = None, subsample_factor: int = 1, - crop_bounds: tuple[tuple[float, float], tuple[float, float], tuple[float, float]] - | None = None, + crop_bounds: ( + tuple[tuple[float, float], tuple[float, float], tuple[float, float]] | None + ) = None, colormap=None, background_color: str | tuple[float, float, float] = "black", spline_subdiv: int = 16, + tube_sides: int = 9, color_range: tuple[float, float] | None = None, color_label: str = "Angle (deg)", ): @@ -691,10 +828,17 @@ def show_streamlines( lut, background_color=background_color, spline_subdiv=spline_subdiv, + tube_sides=tube_sides, color_range=lut_range, color_label=color_label, ) - viewer.run(interactive=interactive, screenshot_path=screenshot_path) + viewer.run( + interactive=interactive, + screenshot_path=screenshot_path, + video_path=video_path, + video_fps=video_fps, + video_frames=video_frames, + ) # --------------------------- diff --git a/src/cardiotensor/visualization/streamlines.py b/src/cardiotensor/visualization/streamlines.py index 8950fbd..ca2356a 100644 --- a/src/cardiotensor/visualization/streamlines.py +++ b/src/cardiotensor/visualization/streamlines.py @@ -5,7 +5,6 @@ from cardiotensor.colormaps.helix_angle import helix_angle_cmap from cardiotensor.utils.streamlines_io_utils import load_trk_streamlines -from cardiotensor.visualization.fury_plotting_streamlines import show_streamlines ANGLE_RANGES = { "HA": (-90.0, 90.0), @@ -61,13 +60,31 @@ def visualize_streamlines( crop_bounds: tuple | None = None, # ((xmin,xmax),(ymin,ymax),(zmin,zmax)) interactive: bool = True, screenshot_path: str | None = None, + video_path: str | None = None, + video_fps: int = 30, + video_frames: int = 120, window_size: tuple[int, int] = (800, 800), colormap=None, - spline_subdiv: int = 16, + spline_subdiv: int = 2, + backend: str = "fury", + mode: str = "tube", + background_color: str | tuple[float, float, float] | None = None, + tube_sides: int = 9, + pyvista_opacity: float = 1.0, + pyvista_show_axes: bool = True, + pyvista_show_bounds: bool = False, + pyvista_shadows: bool = False, ): """ Visualize .trk streamlines with per-point angle-based coloring. - Always renders in tube mode. + + Parameters + ---------- + backend + Rendering backend: "fury" for tractography-oriented interaction or + "pyvista" for polished scientific plots and screenshots. + mode + "tube" for explicit tube geometry or "line" for faster line rendering. """ p = Path(streamlines_file) if not p.exists(): @@ -98,13 +115,13 @@ def visualize_streamlines( print("🧭 You can use: color_by = ha, ia, az, el, elevation, azimuth\n") # Decide the color scalar - mode = color_by.lower().strip() + color_mode = color_by.lower().strip() color_values: list[np.ndarray] | None = None color_range: tuple[float, float] | None = None color_label = "Angle (deg)" - if mode in {"ha", "ia", "az", "el"}: - key = mode.upper() + if color_mode in {"ha", "ia", "az", "el"}: + key = color_mode.upper() color_range = ANGLE_RANGES[key] color_label = f"{key} (deg)" if key in attrs: @@ -119,39 +136,89 @@ def visualize_streamlines( ) else: raise ValueError(f"No per-point attribute '{key}' found in .trk") - elif mode in {"elevation", "azimuth"}: + elif color_mode in {"elevation", "azimuth"}: az_list, el_list = _compute_az_el_from_streamlines(streamlines_xyz) - color_values = el_list if mode == "elevation" else az_list - color_range = (-90.0, 90.0) if mode == "elevation" else ANGLE_RANGES["AZ"] - color_label = "Elevation (deg)" if mode == "elevation" else "Azimuth (deg)" + color_values = el_list if color_mode == "elevation" else az_list + color_range = (-90.0, 90.0) if color_mode == "elevation" else ANGLE_RANGES["AZ"] + color_label = ( + "Elevation (deg)" if color_mode == "elevation" else "Azimuth (deg)" + ) else: raise ValueError("color_by must be one of: ha, ia, az, el, elevation, azimuth") # Default colormap selection if colormap is None: - if mode == "el": + if color_mode == "el": colormap = cm.viridis - elif mode in {"ha", "ia", "elevation"}: + elif color_mode in {"ha", "ia", "elevation"}: colormap = helix_angle_cmap else: colormap = cm.hsv - # Always use tube mode - show_streamlines( - streamlines_xyz=streamlines_xyz, - color_values=color_values, - mode="tube", - line_width=line_width, - interactive=interactive, - screenshot_path=screenshot_path, - window_size=window_size, - downsample_factor=downsample_factor, - max_streamlines=max_streamlines, - filter_min_len=filter_min_len, - subsample_factor=subsample_factor, - crop_bounds=crop_bounds, - colormap=colormap, - spline_subdiv=spline_subdiv, - color_range=color_range, - color_label=color_label, - ) + backend = backend.lower().strip() + render_mode = mode.lower().strip() + if render_mode not in {"tube", "line"}: + raise ValueError("mode must be one of: tube, line") + + if backend == "fury": + from cardiotensor.visualization.fury_plotting_streamlines import ( + show_streamlines, + ) + + show_streamlines( + streamlines_xyz=streamlines_xyz, + color_values=color_values, + mode=render_mode, + line_width=line_width, + interactive=interactive, + screenshot_path=screenshot_path, + video_path=video_path, + video_fps=video_fps, + video_frames=video_frames, + window_size=window_size, + downsample_factor=downsample_factor, + max_streamlines=max_streamlines, + filter_min_len=filter_min_len, + subsample_factor=subsample_factor, + crop_bounds=crop_bounds, + colormap=colormap, + background_color=background_color or "black", + spline_subdiv=spline_subdiv, + tube_sides=tube_sides, + color_range=color_range, + color_label=color_label, + ) + elif backend == "pyvista": + from cardiotensor.visualization.pyvista_plotting_streamlines import ( + show_streamlines_pyvista, + ) + + show_streamlines_pyvista( + streamlines_xyz=streamlines_xyz, + color_values=color_values, + mode=render_mode, + line_width=line_width, + interactive=interactive, + screenshot_path=screenshot_path, + video_path=video_path, + video_fps=video_fps, + video_frames=video_frames, + window_size=window_size, + downsample_factor=downsample_factor, + max_streamlines=max_streamlines, + filter_min_len=filter_min_len, + subsample_factor=subsample_factor, + crop_bounds=crop_bounds, + colormap=colormap, + background_color=background_color or "white", + tube_sides=tube_sides, + color_range=color_range, + color_label=color_label, + opacity=pyvista_opacity, + show_axes=pyvista_show_axes, + show_bounds=pyvista_show_bounds, + shadows=pyvista_shadows, + spline_subdiv=spline_subdiv, + ) + else: + raise ValueError("backend must be one of: fury, pyvista") From 3f5f0df61893c51be620bcb4e8d612891a6d29a6 Mon Sep 17 00:00:00 2001 From: josephbrunet Date: Wed, 17 Jun 2026 14:55:11 +0200 Subject: [PATCH 2/6] fix helical name --- docs/advanced/python-api.md | 14 +++---- docs/getting-started/examples.md | 2 +- docs/getting-technical/conventions.md | 4 +- docs/getting-technical/index.md | 2 +- docs/getting-technical/structure_tensor.md | 2 +- docs/getting-technical/tractography.md | 4 +- docs/reference/cli.md | 2 +- examples/README.md | 2 +- src/cardiotensor/__init__.py | 4 +- .../analysis/gui_analysis_tool.py | 2 +- src/cardiotensor/orientation/__init__.py | 4 +- .../visualization/vector_field.py | 2 +- .../test_orientation_computation_functions.py | 42 +++++++++++++++---- tests/test_utils.py | 2 + 14 files changed, 59 insertions(+), 29 deletions(-) diff --git a/docs/advanced/python-api.md b/docs/advanced/python-api.md index 6c92ff1..fd5be12 100644 --- a/docs/advanced/python-api.md +++ b/docs/advanced/python-api.md @@ -55,12 +55,12 @@ print(fa.shape) # (Y, X), values in [0, 1] --- -## Helix and Transverse Angles +## Helical and Intrusion Angles ```python import numpy as np import matplotlib.pyplot as plt -from cardiotensor import compute_helix_and_transverse_angles +from cardiotensor import compute_helical_and_intrusion_angles from cardiotensor.orientation import rotate_vectors_to_new_axis # Use one slice of the eigenvector field computed above @@ -76,17 +76,17 @@ center_point = ( z_index, ) -helix_angle, intrusion_angle = compute_helix_and_transverse_angles( +helical_angle, intrusion_angle = compute_helical_and_intrusion_angles( rotated_vectors, center_point, ) # Both arrays shape: (Y, X), values in degrees -# The second output corresponds to the transverse/intrusion companion angle map +# The second output corresponds to the intrusion angle map fig, axes = plt.subplots(1, 2, figsize=(10, 4), constrained_layout=True) -im0 = axes[0].imshow(helix_angle, cmap="viridis", vmin=-90, vmax=90) -axes[0].set_title("Helix angle") +im0 = axes[0].imshow(helical_angle, cmap="viridis", vmin=-90, vmax=90) +axes[0].set_title("Helical angle") axes[0].axis("off") fig.colorbar(im0, ax=axes[0], label="degrees") @@ -113,7 +113,7 @@ from cardiotensor import ( calculate_structure_tensor, compute_azimuth_and_elevation, compute_fraction_anisotropy, - compute_helix_and_transverse_angles, + compute_helical_and_intrusion_angles, # Tractography generate_streamlines_from_vector_field, generate_streamlines_from_params, diff --git a/docs/getting-started/examples.md b/docs/getting-started/examples.md index 47bd574..37acae7 100644 --- a/docs/getting-started/examples.md +++ b/docs/getting-started/examples.md @@ -73,7 +73,7 @@ The `./examples/` directory contains: 3. Outputs will be saved in the `./output` directory with the following structure: ``` ./output - ├── HA # Helix angle results + ├── HA # Helical angle results ├── IA # Intrusion angle results ├── FA # Fractional anisotropy results └── eigen_vec # 3rd Eigenvectors diff --git a/docs/getting-technical/conventions.md b/docs/getting-technical/conventions.md index a5db575..7abed3f 100644 --- a/docs/getting-technical/conventions.md +++ b/docs/getting-technical/conventions.md @@ -13,13 +13,13 @@ Eigenvector maps are saved in the shape `(3, Z, Y, X)`: - 3 = x, y, z vector components - Each voxel contains the orientation of the 3rd eigenvector (principal myocyte axis) -These vector fields are used for computing helix/intrusion angles, streamlines, and for visualization. +These vector fields are used for computing helical/intrusion angles, streamlines, and for visualization. ## Units of measurement - Lengths are expressed in **pixels** (or voxels for 3D). -- Angles (helix, intrusion) are in **degrees**. +- Angles (helical, intrusion) are in **degrees**. - Fractional anisotropy is **dimensionless**, ranging from 0 (isotropic) to 1 (highly anisotropic). diff --git a/docs/getting-technical/index.md b/docs/getting-technical/index.md index 3c69093..5f37ec4 100644 --- a/docs/getting-technical/index.md +++ b/docs/getting-technical/index.md @@ -8,7 +8,7 @@ This section provides a deeper understanding of how `cardiotensor` works, includ - [Structure Tensor](./structure_tensor.md): Mathematical background on the structure tensor and eigen decomposition -- [Angle Definitions](./angles.md): Detailed explanation of helix and intrusion angles +- [Angle Definitions](./angles.md): Detailed explanation of helical and intrusion angles - [Fractional Anisotropy](./fractional_anisotropy.md) diff --git a/docs/getting-technical/structure_tensor.md b/docs/getting-technical/structure_tensor.md index 88608e0..b60f5e6 100644 --- a/docs/getting-technical/structure_tensor.md +++ b/docs/getting-technical/structure_tensor.md @@ -71,7 +71,7 @@ Each structure tensor $S$ is decomposed into eigenvalues $\lambda_1 \leq \lambda ## Output and Interpretation * The third eigenvector $\vec{v}_1$ is stored as the local myocyte orientation. -* Helix and intrusion angles are computed from $\vec{v}_1$ after transforming it into cylindrical coordinates. +* Helical and intrusion angles are computed from $\vec{v}_1$ after transforming it into cylindrical coordinates. * Fractional Anisotropy (FA) is computed using the three eigenvalues (see [FA section](./fractional_anisotropy.md)). --- diff --git a/docs/getting-technical/tractography.md b/docs/getting-technical/tractography.md index d9a9c79..aa6213b 100644 --- a/docs/getting-technical/tractography.md +++ b/docs/getting-technical/tractography.md @@ -56,7 +56,7 @@ Streamlines shorter than `min_length_pts` (default: 10 points) are discarded. Yo The streamline results are saved in `.trk` format as: * `streamlines`: a list of arrays, each of shape (N, 3), where N is the number of points in that streamline -* `ha_values`: sampled HA (helix angle) values along each point +* `ha_values`: sampled HA (helical angle) values along each point These can be loaded in Python or exported to `.vtk` for 3D visualization in ParaView. @@ -74,7 +74,7 @@ cardio-generate config.conf --seeds 20000 --bin 2 --step 0.5 --fa-threshold 0.15 * Streamlines are computed from the 3rd eigenvector of the structure tensor, corresponding to the myocyte axis. * FA is computed once from the structure tensor and optionally downsampled for speed. -* Helix angle (HA) is sampled at each point along the streamline and saved for further analysis. +* Helical angle (HA) is sampled at each point along the streamline and saved for further analysis. --- diff --git a/docs/reference/cli.md b/docs/reference/cli.md index 343e474..d7f36a3 100644 --- a/docs/reference/cli.md +++ b/docs/reference/cli.md @@ -7,7 +7,7 @@ Cardiotensor provides powerful features for analyzing 3D cardiac imaging data, f Compute myocyte orientation from a 3D volume using a configuration file. - `cardio-tensor` - Computes structure tensor, helix/transverse angle, FA, and eigenvectors. + Computes structure tensor, helical/intrusion angle, FA, and eigenvectors. - `cardio-tensor-slurm` Submits chunked `cardio-tensor` runs as SLURM array jobs. diff --git a/examples/README.md b/examples/README.md index e61e712..19a09b5 100644 --- a/examples/README.md +++ b/examples/README.md @@ -73,7 +73,7 @@ To process the entire volume: 3. The results will be saved in the `./output` directory with the following structure: ``` ./output - ├── HA # Helix angle results + ├── HA # Helical angle results ├── IA # Intrusion angle results └── FA # Fractional anisotropy results ``` diff --git a/src/cardiotensor/__init__.py b/src/cardiotensor/__init__.py index 5ad55b8..bcb7c71 100644 --- a/src/cardiotensor/__init__.py +++ b/src/cardiotensor/__init__.py @@ -45,7 +45,7 @@ calculate_structure_tensor, compute_azimuth_and_elevation, compute_fraction_anisotropy, - compute_helix_and_transverse_angles, + compute_helical_and_intrusion_angles, compute_orientation, ) from cardiotensor.tractography import ( @@ -72,7 +72,7 @@ "calculate_structure_tensor", "compute_azimuth_and_elevation", "compute_fraction_anisotropy", - "compute_helix_and_transverse_angles", + "compute_helical_and_intrusion_angles", # Tractography "generate_streamlines_from_vector_field", "generate_streamlines_from_params", diff --git a/src/cardiotensor/analysis/gui_analysis_tool.py b/src/cardiotensor/analysis/gui_analysis_tool.py index ad47f38..f8fc0fb 100644 --- a/src/cardiotensor/analysis/gui_analysis_tool.py +++ b/src/cardiotensor/analysis/gui_analysis_tool.py @@ -105,7 +105,7 @@ def plot_label_and_limits(mode: str): if m == "FA": return "Fractional Anisotropy", 0.0, 1.0 if m in {"HA", "IA", "EL"}: - name = {"HA": "Helix Angle", "IA": "Intrusion Angle", "EL": "Elevation Angle"}[ + name = {"HA": "Helical Angle", "IA": "Intrusion Angle", "EL": "Elevation Angle"}[ m ] return f"{name} (°)", -90.0, 90.0 diff --git a/src/cardiotensor/orientation/__init__.py b/src/cardiotensor/orientation/__init__.py index ed28aff..91f6b34 100644 --- a/src/cardiotensor/orientation/__init__.py +++ b/src/cardiotensor/orientation/__init__.py @@ -4,7 +4,7 @@ calculate_structure_tensor, compute_azimuth_and_elevation, compute_fraction_anisotropy, - compute_helix_and_transverse_angles, + compute_helical_and_intrusion_angles, orient_vectors_z_positive, rotate_vectors_to_new_axis, ) @@ -16,7 +16,7 @@ "calculate_structure_tensor", "compute_azimuth_and_elevation", "compute_fraction_anisotropy", - "compute_helix_and_transverse_angles", + "compute_helical_and_intrusion_angles", "compute_orientation", "orient_vectors_z_positive", "rotate_vectors_to_new_axis", diff --git a/src/cardiotensor/visualization/vector_field.py b/src/cardiotensor/visualization/vector_field.py index 090f621..d16183b 100644 --- a/src/cardiotensor/visualization/vector_field.py +++ b/src/cardiotensor/visualization/vector_field.py @@ -49,7 +49,7 @@ def visualize_vector_field( Path to the 3D vector field (directory or file). The vector field must be stored as (3, Z, Y, X) or (Z, Y, X, 3) numpy arrays. color_volume_path : str or Path, optional - Path to a scalar volume used to color the vectors (e.g., helix angles). + Path to a scalar volume used to color the vectors (e.g., helical angles). If RGB, the channels will be averaged to a single scalar map. mask_path : str or Path, optional Path to a binary mask volume. Vectors outside the mask are set to NaN diff --git a/tests/test_orientation_computation_functions.py b/tests/test_orientation_computation_functions.py index f1799d4..c36b829 100644 --- a/tests/test_orientation_computation_functions.py +++ b/tests/test_orientation_computation_functions.py @@ -8,7 +8,7 @@ calculate_structure_tensor, compute_azimuth_and_elevation, compute_fraction_anisotropy, - compute_helix_and_transverse_angles, + compute_helical_and_intrusion_angles, interpolate_points, orient_vectors_z_positive, plot_images, @@ -127,28 +127,56 @@ def test_compute_azimuth_and_elevation_uses_z_positive_convention(): assert np.allclose(azimuth, [[0, 0]], atol=1e-6) -def test_compute_helix_and_transverse_angles_orients_circumferential_positive(): +def test_compute_helical_and_intrusion_angles_orients_circumferential_positive(): vector_slice = np.zeros((3, 1, 2), dtype=np.float32) vector_slice[:, 0, 0] = [0, 1, 1] vector_slice[:, 0, 1] = [0, -1, -1] - helix, transverse = compute_helix_and_transverse_angles( + helical, intrusion = compute_helical_and_intrusion_angles( vector_slice, center_point=(0, 0, 0) ) - assert np.allclose(helix, [[45, 45]], atol=1e-6) - assert np.allclose(transverse, [[0, 0]], atol=1e-6) + assert np.allclose(helical, [[45, 45]], atol=1e-6) + assert np.allclose(intrusion, [[0, 0]], atol=1e-6) + + +def test_compute_angles_use_full_vector_without_projection(): + vector_slice = np.zeros((3, 1, 1), dtype=np.float32) + vector_slice[:, 0, 0] = [1, 2, 3] + + helical, intrusion = compute_helical_and_intrusion_angles( + vector_slice, center_point=(0, 0, 0) + ) + + expected_helical = np.rad2deg(np.arctan2(3, np.hypot(1, 2))) + expected_intrusion = np.rad2deg(np.arctan2(1, np.hypot(2, 3))) + assert np.allclose(helical, [[expected_helical]], atol=1e-6) + assert np.allclose(intrusion, [[expected_intrusion]], atol=1e-6) + + +def test_compute_angles_can_return_projected_intrusion_values(): + vector_slice = np.zeros((3, 1, 1), dtype=np.float32) + vector_slice[:, 0, 0] = [1, 2, 3] + + helical, intrusion = compute_helical_and_intrusion_angles( + vector_slice, center_point=(0, 0, 0), projected=True + ) + + expected_helical = np.rad2deg(np.arctan2(3, 2)) + expected_intrusion = np.rad2deg(np.arctan2(1, 2)) + assert np.allclose(helical, [[expected_helical]], atol=1e-6) + assert np.allclose(intrusion, [[expected_intrusion]], atol=1e-6) def test_write_images_and_vectors(tmp_path: Path): # --- Prepare dummy 2D slices --- - img_helix = np.ones((5, 5), dtype=np.float32) + img_helical = np.ones((5, 5), dtype=np.float32) img_intrusion = np.ones((5, 5), dtype=np.float32) img_FA = np.ones((5, 5), dtype=np.float32) # --- Test write_images --- write_images( - img_helix, + img_helical, img_intrusion, img_FA, start_index=0, diff --git a/tests/test_utils.py b/tests/test_utils.py index 157184d..1212c6f 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -36,6 +36,7 @@ def test_read_conf_file(tmp_path): "REVERSE = False\n\n" "[ANGLE CALCULATION]\n" "WRITE_ANGLES = True\n" + "PROJECTED_ANGLES = True\n" "AXIS_POINTS = (0,0,0), (1,1,1)\n\n" "[TEST]\n" "TEST = True\n" @@ -65,6 +66,7 @@ def test_read_conf_file(tmp_path): # Angles assert config["WRITE_ANGLES"] is True + assert config["PROJECTED_ANGLES"] is True assert isinstance(config["AXIS_POINTS"], list) assert config["AXIS_POINTS"] == [(0, 0, 0), (1, 1, 1)] From f39c970e97f86c23a30ace0a40ab477b4b64a72c Mon Sep 17 00:00:00 2001 From: josephbrunet Date: Wed, 17 Jun 2026 14:57:00 +0200 Subject: [PATCH 3/6] remove projection --- .../orientation_computation_functions.py | 52 +++++++++++++------ .../orientation_computation_pipeline.py | 39 ++++++++++---- .../scripts/compute_orientation.py | 9 ++++ 3 files changed, 73 insertions(+), 27 deletions(-) diff --git a/src/cardiotensor/orientation/orientation_computation_functions.py b/src/cardiotensor/orientation/orientation_computation_functions.py index 55d0371..4341d8b 100644 --- a/src/cardiotensor/orientation/orientation_computation_functions.py +++ b/src/cardiotensor/orientation/orientation_computation_functions.py @@ -372,19 +372,31 @@ def orient_vectors_z_positive(vector_field_slice: np.ndarray) -> np.ndarray: return oriented -def compute_helix_and_transverse_angles( +def compute_helical_and_intrusion_angles( vector_field_2d: np.ndarray, center_point: tuple[int, int, int], + projected: bool = False, ) -> tuple[np.ndarray, np.ndarray]: """ - Computes helix and transverse angles from a 2D vector field. + Computes helical and intrusion angles from a 2D vector field. + + By default, both angles use the full local vector. The helical angle retains + radial contribution instead of first projecting into the circumferential/ + longitudinal plane, and the intrusion angle retains longitudinal + contribution instead of first projecting into the radial/circumferential + plane. + + Set projected=True to return the legacy projected helical angle and + projected intrusion angle for comparison with projection-based literature. Args: vector_field_2d (np.ndarray): 2D orientation vector field. center_point (Tuple[int, int, int]): Coordinates of the center point. + projected (bool): If True, compute projected helical/intrusion angles. + Defaults to False for unprojected 3D helical/intrusion angles. Returns: - Tuple[np.ndarray, np.ndarray]: Helix and transverse angle arrays. + Tuple[np.ndarray, np.ndarray]: Helical and intrusion angle arrays. """ center = center_point[0:2] # Replace with actual values rows, cols = vector_field_2d.shape[1:3] @@ -419,19 +431,27 @@ def compute_helix_and_transverse_angles( reshaped_rotated_vector_field ) - # Calculate helix and transverse angles - helix_angle = np.arctan2( - reshaped_rotated_vector_field[2, :, :], - reshaped_rotated_vector_field[1, :, :], - ) - transverse_angle = np.arctan2( - reshaped_rotated_vector_field[0, :, :], - reshaped_rotated_vector_field[1, :, :], - ) - helix_angle = np.rad2deg(helix_angle) - transverse_angle = np.rad2deg(transverse_angle) + radial_component = reshaped_rotated_vector_field[0, :, :] + circumferential_component = reshaped_rotated_vector_field[1, :, :] + longitudinal_component = reshaped_rotated_vector_field[2, :, :] + + # Calculate helical and intrusion angles + if projected: + helical_angle = np.arctan2(longitudinal_component, circumferential_component) + intrusion_angle = np.arctan2(radial_component, circumferential_component) + else: + helical_angle = np.arctan2( + longitudinal_component, + np.hypot(radial_component, circumferential_component), + ) + intrusion_angle = np.arctan2( + radial_component, + np.hypot(circumferential_component, longitudinal_component), + ) + helical_angle = np.rad2deg(helical_angle) + intrusion_angle = np.rad2deg(intrusion_angle) - return helix_angle, transverse_angle + return helical_angle, intrusion_angle def compute_azimuth_and_elevation( @@ -539,7 +559,7 @@ def plot_images( colormap_angle=None, colormap_angle2=None, colormap_FA=None, - angle1_title: str = "Helix Angle", + angle1_title: str = "Helical Angle", angle2_title: str = "Intrusion Angle", angle_ranges: tuple[tuple[float, float], tuple[float, float]] | None = None, save_path: str | None = None, diff --git a/src/cardiotensor/orientation/orientation_computation_pipeline.py b/src/cardiotensor/orientation/orientation_computation_pipeline.py index 3c8a048..65114c8 100644 --- a/src/cardiotensor/orientation/orientation_computation_pipeline.py +++ b/src/cardiotensor/orientation/orientation_computation_pipeline.py @@ -18,7 +18,7 @@ calculate_structure_tensor, compute_azimuth_and_elevation, compute_fraction_anisotropy, - compute_helix_and_transverse_angles, + compute_helical_and_intrusion_angles, interpolate_points, plot_images, remove_padding, @@ -32,7 +32,7 @@ # --- small helpers --- def _resolve_colormap(colormap: str | None): - """Resolve colormap names, including the project-specific helix angle map.""" + """Resolve colormap names, including the project-specific helical angle map.""" if colormap is None: return None cmap_name = colormap.strip() @@ -161,6 +161,7 @@ def compute_orientation( colormap: str | None = None, colormap_angle1: str | None = None, colormap_angle2: str | None = None, + projected: bool = False, ) -> None: """ Compute the orientation for a volume dataset. @@ -178,6 +179,7 @@ def compute_orientation( vertical_padding: Padding slices for tensor computation. write_vectors: Whether to save eigenvectors. Ignored in test mode. write_angles: Whether to save HA/IA/FA maps. + projected: If True in ha_ia mode, write projected HA/IA legacy maps. use_gpu: Use GPU acceleration for tensor computation. is_test: If True, runs in test mode and outputs plots. n_slice_test: Number of slices to process in test mode. @@ -194,7 +196,9 @@ def compute_orientation( raise ValueError("sigma must be <= rho") if angle_mode.lower() == "ha_ia": - angle_names = ("HA", "IA") + angle_names = ( + ("HA_projected", "IA_projected") if projected else ("HA", "IA") + ) elif angle_mode.lower() == "az_el": angle_names = ("AZ", "EL") else: @@ -210,6 +214,10 @@ def compute_orientation( ) write_vectors = False + projected_status = ( + projected if angle_mode.lower().strip() == "ha_ia" else "[n/a]" + ) + print(f""" Parameters: - Volume path: {volume_path} @@ -221,6 +229,7 @@ def compute_orientation( - truncate: {truncate} - Write angles: {write_angles} - Angle mode: {angle_mode} -> {angle_names[0]}, {angle_names[1]} + - Projected HA/IA:{projected_status} - Write vectors: {write_vectors} - Use GPU: {use_gpu} - Test mode: {is_test} @@ -392,6 +401,7 @@ def update_bar(_): colormap, colormap_angle1, colormap_angle2, + projected, ), callback=update_bar, ) @@ -440,6 +450,7 @@ def update_bar(_): colormap, colormap_angle1, colormap_angle2, + projected, ) if is_test: @@ -470,6 +481,7 @@ def compute_slice_angles_and_anisotropy( colormap: str | None = None, colormap_angle1: str | None = None, colormap_angle2: str | None = None, + projected: bool = False, ) -> tuple[float, float, bool]: """ Compute either HA/IA or Azimuth/Elevation plus FA for a single slice, @@ -478,7 +490,9 @@ def compute_slice_angles_and_anisotropy( # Decide angle labels and ranges based on mode mode = angle_mode.lower().strip() if mode == "ha_ia": - angle_names = ("HA", "IA") + angle_names = ( + ("HA_projected", "IA_projected") if projected else ("HA", "IA") + ) angle_ranges = ((-90.0, 90.0), (-90.0, 90.0)) elif mode == "az_el": angle_names = ("AZ", "EL") @@ -547,8 +561,8 @@ def compute_slice_angles_and_anisotropy( ) if mode == "ha_ia": - img_angle1, img_angle2 = compute_helix_and_transverse_angles( - vector_field_slice_rotated, center_point + img_angle1, img_angle2 = compute_helical_and_intrusion_angles( + vector_field_slice_rotated, center_point, projected=projected ) else: # "az_el" img_angle1, img_angle2 = compute_azimuth_and_elevation( @@ -560,11 +574,14 @@ def compute_slice_angles_and_anisotropy( # Test mode: visualize a 2x2 figure and write to test subfolder if is_test: - titles = ( - ("Helix Angle", "Intrusion Angle") - if mode == "ha_ia" - else ("Azimuth", "Elevation") - ) + if mode == "ha_ia": + titles = ( + ("Projected Helical Angle", "Projected Intrusion Angle") + if projected + else ("Helical Angle", "Intrusion Angle") + ) + else: + titles = ("Azimuth", "Elevation") test_output_dir = os.path.join(output_dir, "test_slice") plot_images( img_slice, diff --git a/src/cardiotensor/scripts/compute_orientation.py b/src/cardiotensor/scripts/compute_orientation.py index 381d2c9..7c730a4 100644 --- a/src/cardiotensor/scripts/compute_orientation.py +++ b/src/cardiotensor/scripts/compute_orientation.py @@ -67,6 +67,11 @@ def script() -> None: default=None, help="Colormap for the second angle output, for example IA or EL.", ) + parser.add_argument( + "--projected", + action="store_true", + help="In ha_ia mode, write legacy projected HA/IA maps instead of 3D HA/IA.", + ) args = parser.parse_args() conf_file_path = args.conf_file_path @@ -77,6 +82,7 @@ def script() -> None: colormap = args.colormap colormap_angle1 = args.colormap_angle1 colormap_angle2 = args.colormap_angle2 + projected = args.projected # --- Read configuration --- try: @@ -106,6 +112,7 @@ def script() -> None: is_test = force_test or params.get("TEST", False) n_slice_test = params.get("N_SLICE_TEST", None) show_quiver = params.get("SHOW_QUIVER", True) + projected = projected or params.get("PROJECTED_ANGLES", False) n_chunk = params.get("N_CHUNK", 100) colormap = colormap or params.get("COLORMAP", None) colormap_angle1 = colormap_angle1 or params.get("COLORMAP_ANGLE1", None) @@ -152,6 +159,7 @@ def script() -> None: colormap=colormap, colormap_angle1=colormap_angle1, colormap_angle2=colormap_angle2, + projected=projected, ) print(f"--- {time.time() - t0:.1f} seconds (TEST mode) ---") return @@ -199,6 +207,7 @@ def script() -> None: colormap=colormap, colormap_angle1=colormap_angle1, colormap_angle2=colormap_angle2, + projected=projected, ) elapsed = time.time() - t0 From 6fe36e674447fcddd58a5f02eced8499d34b9986 Mon Sep 17 00:00:00 2001 From: josephbrunet Date: Wed, 17 Jun 2026 14:57:44 +0200 Subject: [PATCH 4/6] add param --- docs/reference/configuration.md | 6 ++-- examples/parameters_example.conf | 58 ++++++++++++++++---------------- src/cardiotensor/utils/utils.py | 7 ++++ 3 files changed, 40 insertions(+), 31 deletions(-) diff --git a/docs/reference/configuration.md b/docs/reference/configuration.md index a81ad1b..adfeef3 100644 --- a/docs/reference/configuration.md +++ b/docs/reference/configuration.md @@ -61,7 +61,7 @@ REVERSE = False WRITE_ANGLES = True # Choose which angle pair to compute: -# ha_ia → Helix Angle / Intrusion Angle (standard myocardial architecture) +# ha_ia → Helical Angle / Intrusion Angle (standard myocardial architecture) # az_el → Azimuth / Elevation (generic vector orientation in 3D) ANGLE_MODE = ha_ia @@ -145,7 +145,9 @@ OUTPUT_TYPE = 8bit --- ### `[ANGLE CALCULATION]` -- **`WRITE_ANGLES`**: Save helix and intrusion angles and fractional anisotropy values. +- **`WRITE_ANGLES`**: Save helical and intrusion angles and fractional anisotropy values. +- **`ANGLE_MODE`**: Choose `ha_ia` for cardiac HA/IA or `az_el` for generic azimuth/elevation maps. +- **`PROJECTED_ANGLES`**: In `ha_ia` mode, keep `False` for unprojected 3D HA/IA. Set `True` only to write legacy projected `HA_projected`/`IA_projected` comparison maps. - **`AXIS_POINTS`**: List of 3D points `[X, Y, Z]` along the ventricle axis. Typically, the first and last points correspond to the mitral valve and apex. Intermediate points help create a curved axis via interpolation. --- diff --git a/examples/parameters_example.conf b/examples/parameters_example.conf index ea2c98a..811cc20 100644 --- a/examples/parameters_example.conf +++ b/examples/parameters_example.conf @@ -49,7 +49,7 @@ REVERSE = False WRITE_ANGLES = True # Choose which angle pair to compute: -# ha_ia → Helix Angle / Intrusion Angle (standard myocardial architecture) +# ha_ia → Helical Angle / Intrusion Angle (standard myocardial architecture) # az_el → Azimuth / Elevation (generic vector orientation in 3D) ANGLE_MODE = ha_ia @@ -60,22 +60,22 @@ ANGLE_MODE = ha_ia AXIS_POINTS = [104,110,116], [41,87,210], [68,95,162] -[TEST] -# Enable test mode: -# - True: Process and plot only a single slice for testing. -# - False: Perform the full processing on the entire volume. -TEST = True - -# Specify the slice number to process when test mode is enabled. -N_SLICE_TEST = 155 - -# Overlay the in-plane vector field as a quiver plot on the grayscale test slice. -# - True: show the arrows -# - False: keep only the grayscale slice and angle/FA panels -SHOW_QUIVER = True - - -[OUTPUT] +[TEST] +# Enable test mode: +# - True: Process and plot only a single slice for testing. +# - False: Perform the full processing on the entire volume. +TEST = True + +# Specify the slice number to process when test mode is enabled. +N_SLICE_TEST = 155 + +# Overlay the in-plane vector field as a quiver plot on the grayscale test slice. +# - True: show the arrows +# - False: keep only the grayscale slice and angle/FA panels +SHOW_QUIVER = True + + +[OUTPUT] # Path to the folder where the results will be saved OUTPUT_PATH =./output @@ -83,15 +83,15 @@ OUTPUT_PATH =./output # Default format is jp2 OUTPUT_FORMAT = tif -# Type of pixel values in the output file: -# - "8bit" for grayscale 8-bit images -# - "rgb" for 3-channel color images -OUTPUT_TYPE = 8bit - -# Optional colormaps for RGB output and test plots. -# Use any Matplotlib colormap name, or helix_angle for the CardioTensor map. -# COLORMAP applies to both angles unless COLORMAP_ANGLE1 or COLORMAP_ANGLE2 is set. -# For ANGLE_MODE = az_el, defaults are helix_angle for AZ and viridis for EL. -# COLORMAP = -# COLORMAP_ANGLE1 = -# COLORMAP_ANGLE2 = +# Type of pixel values in the output file: +# - "8bit" for grayscale 8-bit images +# - "rgb" for 3-channel color images +OUTPUT_TYPE = 8bit + +# Optional colormaps for RGB output and test plots. +# Use any Matplotlib colormap name, or helix_angle for the CardioTensor map. +# COLORMAP applies to both angles unless COLORMAP_ANGLE1 or COLORMAP_ANGLE2 is set. +# For ANGLE_MODE = az_el, defaults are helix_angle for AZ and viridis for EL. +# COLORMAP = +# COLORMAP_ANGLE1 = +# COLORMAP_ANGLE2 = diff --git a/src/cardiotensor/utils/utils.py b/src/cardiotensor/utils/utils.py index aed5d19..59ff4c1 100644 --- a/src/cardiotensor/utils/utils.py +++ b/src/cardiotensor/utils/utils.py @@ -189,6 +189,13 @@ def parse_coordinates(section: str, option: str, fallback: str = ""): "ANGLE_MODE": config.get( "ANGLE CALCULATION", "ANGLE_MODE", fallback="ha_ia" ).strip(), + "PROJECTED_ANGLES": config.getboolean( + "ANGLE CALCULATION", + "PROJECTED_ANGLES", + fallback=config.getboolean( + "ANGLE CALCULATION", "PROJECTED", fallback=False + ), + ), "AXIS_POINTS": parse_coordinates("ANGLE CALCULATION", "AXIS_POINTS"), # TEST "TEST": config.getboolean("TEST", "TEST", fallback=False), From 0650ae192c8b086b54f07264f427ee20b197ceca Mon Sep 17 00:00:00 2001 From: josephbrunet Date: Wed, 17 Jun 2026 14:57:49 +0200 Subject: [PATCH 5/6] fix doc --- docs/getting-technical/angles.md | 31 ++++++++++++++++++++++++++----- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/docs/getting-technical/angles.md b/docs/getting-technical/angles.md index a95a098..724b3c5 100644 --- a/docs/getting-technical/angles.md +++ b/docs/getting-technical/angles.md @@ -1,6 +1,8 @@ # Angle Definitions -This page explains how helix and intrusion angles are calculated from the 3D eigenvector field derived by `cardiotensor`. +This page explains how helical and intrusion angles are calculated from the 3D eigenvector field derived by `cardiotensor`. + +By default, Cardiotensor reports **unprojected 3D angles**. The primary eigenvector is not first flattened onto a 2D plane before the angle is measured, because projection discards one component of the vector and can bias the resulting angle. Projection-based angles are available only for legacy comparison. ## Coordinate System @@ -10,12 +12,17 @@ A transformation to a cylindrical coordinate system is defined for each voxel ba - **Circumferential (θ)**: tangential around the ventricle - **Longitudinal (z)**: base to apex direction -To compute local fiber angles consistently, all eigenvectors are first rotated into this cylindrical coordinate frame. This alignment is performed using Rodrigues' rotation formula, which computes the minimal-angle rotation that maps the global reference axis (here the z-axis) onto the local longitudinal axis at each point. This allows a robust comparison of orientations across the myocardium. +To compute local fiber angles consistently, all eigenvectors are first rotated into this cylindrical coordinate frame. This alignment is performed using the Rodrigues rotation formula, which computes the minimal-angle rotation that maps the global reference axis (here the z-axis) onto the local longitudinal axis at each point. This allows a robust comparison of orientations across the myocardium. + +## Helical Angle (HA) +The helical angle is defined as the angle between the primary myocyte-axis eigenvector \\( \vec{v}_1 \\) and the local circumferential plane. -## Helix Angle (HA) +In local cylindrical coordinates, with radial component \\(R\\), circumferential component \\(C\\), and longitudinal component \\(L\\), the unprojected helical angle is: -The helix angle is defined as the angle between the third eigenvector \\( \vec{v}_3 \\) (smallest eigenvalue direction) and the **local circumferential plane**. +\\[ +\mathrm{HA} = \arctan2\left(L, \sqrt{R^2 + C^2}\right) +\\] It captures the transmural variation of fiber orientation from epicardium to endocardium. @@ -26,10 +33,24 @@ Typical pattern: ## Intrusion Angle (IA) -The intrusion angle is the angle between \\( \vec{v}_3 \\) and the **tangential plane** (longitudinal + circumferential). +The intrusion angle is the angle between the primary myocyte-axis eigenvector \\( \vec{v}_1 \\) and the **tangential plane** (longitudinal + circumferential). + +Using the same local components, the unprojected intrusion angle is: + +\\[ +\mathrm{IA} = \arctan2\left(R, \sqrt{C^2 + L^2}\right) +\\] It captures radial deviation of fiber aggregates and can help identify wall thickening or microstructural disruptions. +## Projection Bias + +Conventional projected angles are computed after removing one vector component, for example \\(\arctan2(L, C)\\) for projected helical angle or \\(\arctan2(R, C)\\) for projected intrusion angle. These projected quantities can differ from the true 3D orientation when the discarded component is large. + +Set `PROJECTED_ANGLES = True` only when you need legacy projected `HA_projected` and `IA_projected` maps for comparison with literature. + +This convention follows the projection-error discussion in Agger et al., "Anatomically correct assessment of the orientation of the cardiomyocytes using diffusion tensor imaging", *NMR in Biomedicine* (2020), https://doi.org/10.1002/nbm.4205. + ## Angle Ranges Both angles are reported in degrees: From 0c13b0dcb65b4046d59b046cc4a158581c13c58d Mon Sep 17 00:00:00 2001 From: josephbrunet Date: Wed, 17 Jun 2026 17:37:52 +0200 Subject: [PATCH 6/6] fix assert issue --- tests/test_orientation_computation_functions.py | 12 ++++++------ tests/test_utils.py | 4 +++- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/tests/test_orientation_computation_functions.py b/tests/test_orientation_computation_functions.py index c36b829..e1e9bcd 100644 --- a/tests/test_orientation_computation_functions.py +++ b/tests/test_orientation_computation_functions.py @@ -136,8 +136,8 @@ def test_compute_helical_and_intrusion_angles_orients_circumferential_positive() vector_slice, center_point=(0, 0, 0) ) - assert np.allclose(helical, [[45, 45]], atol=1e-6) - assert np.allclose(intrusion, [[0, 0]], atol=1e-6) + np.testing.assert_allclose(helical, [[45, 45]], atol=1e-6) + np.testing.assert_allclose(intrusion, [[0, 0]], atol=1e-6) def test_compute_angles_use_full_vector_without_projection(): @@ -150,8 +150,8 @@ def test_compute_angles_use_full_vector_without_projection(): expected_helical = np.rad2deg(np.arctan2(3, np.hypot(1, 2))) expected_intrusion = np.rad2deg(np.arctan2(1, np.hypot(2, 3))) - assert np.allclose(helical, [[expected_helical]], atol=1e-6) - assert np.allclose(intrusion, [[expected_intrusion]], atol=1e-6) + np.testing.assert_allclose(helical, [[expected_helical]], atol=1e-6) + np.testing.assert_allclose(intrusion, [[expected_intrusion]], atol=1e-6) def test_compute_angles_can_return_projected_intrusion_values(): @@ -164,8 +164,8 @@ def test_compute_angles_can_return_projected_intrusion_values(): expected_helical = np.rad2deg(np.arctan2(3, 2)) expected_intrusion = np.rad2deg(np.arctan2(1, 2)) - assert np.allclose(helical, [[expected_helical]], atol=1e-6) - assert np.allclose(intrusion, [[expected_intrusion]], atol=1e-6) + np.testing.assert_allclose(helical, [[expected_helical]], atol=1e-6) + np.testing.assert_allclose(intrusion, [[expected_intrusion]], atol=1e-6) def test_write_images_and_vectors(tmp_path: Path): diff --git a/tests/test_utils.py b/tests/test_utils.py index 1212c6f..b4940fa 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,3 +1,5 @@ +from unittest import TestCase + import numpy as np import pytest @@ -66,7 +68,7 @@ def test_read_conf_file(tmp_path): # Angles assert config["WRITE_ANGLES"] is True - assert config["PROJECTED_ANGLES"] is True + TestCase().assertIs(config["PROJECTED_ANGLES"], True) assert isinstance(config["AXIS_POINTS"], list) assert config["AXIS_POINTS"] == [(0, 0, 0), (1, 1, 1)]