add crop during visu

Author: JosephBrunet

src/cardiotensor/orientation/orientation_computation_pipeline.py

@@ -40,7 +40,6 @@ def compute_orientation(
         conf_file_path (str): Path to the configuration file.
         start_index (int, optional): Start index for processing. Default is 0.
         end_index (int | None): End index for processing. Default is None.
-        use_gpu (bool, optional): Whether to use GPU for calculations. Default is False.
 
     Returns:
         None

src/cardiotensor/scripts/generate_streamlines.py

@@ -83,6 +83,11 @@ def script():
         print(f"⚠️ No eigenvector directory at {eigen_dir}")
         sys.exit(1)
 
+
+    #---------------------------------
+    # BINNING
+
+
     if bin_factor > 1:
         downsample_vector_volume(eigen_dir, bin_factor, OUTPUT_DIR)
         vec_load_dir = OUTPUT_DIR / f"bin{bin_factor}" / "eigen_vec"
@@ -93,64 +98,86 @@ def script():
         start_binned = start_idx
         end_binned = end_idx
 
+    FA_dir = OUTPUT_DIR / "FA"
+    if bin_factor > 1:
+        downsample_volume(
+            input_path=FA_dir,
+            bin_factor=bin_factor,
+            output_dir=OUTPUT_DIR,
+            subfolder="FA",
+            out_ext="tif",
+            min_value=0,
+            max_value=255,
+        )
+        fa_load_dir = OUTPUT_DIR / f"bin{bin_factor}" / "FA"
+    else:
+        fa_load_dir = FA_dir
+        
+    HA_dir = OUTPUT_DIR / "HA"
+    if bin_factor > 1:
+        downsample_volume(
+            input_path=HA_dir,
+            bin_factor=bin_factor,
+            output_dir=OUTPUT_DIR,
+            subfolder="HA",
+            out_ext="tif",
+            min_value=0,
+            max_value=255,
+        )
+        ha_load_dir = OUTPUT_DIR / f"bin{bin_factor}" / "HA"
+    else:
+        ha_load_dir = HA_dir
+        
+    #---------------------------------
+
+
+
     print("Loading vector field...")
     vec_reader = DataReader(vec_load_dir)
-    vector_slices = vec_reader.load_volume(
+    vector_field = vec_reader.load_volume(
         start_index=start_binned, end_index=end_binned
     )
 
-    if vector_slices.ndim == 4 and vector_slices.shape[0] == 3:
-        vector_field = vector_slices
-    elif vector_slices.ndim == 4 and vector_slices.shape[-1] == 3:
-        vector_field = np.moveaxis(vector_slices, -1, 0)
+    if vector_field.ndim == 4 and vector_field.shape[0] == 3:
+        vector_field = vector_field
+    elif vector_field.ndim == 4 and vector_field.shape[-1] == 3:
+        vector_field = np.moveaxis(vector_field, -1, 0)
     else:
-        print(f"⚠️ Unexpected vector field shape: {vector_slices.shape}")
+        print(f"⚠️ Unexpected vector field shape: {vector_field.shape}")
         sys.exit(1)
 
     print("Ensuring Z-components are positive...")
     neg_mask = vector_field[2] < 0
     vector_field[:, neg_mask] *= -1
+    del neg_mask
 
     if MASK_PATH:
         print("Applying mask from config...")
         mask_reader = DataReader(MASK_PATH)
 
         # Load the corresponding mask volume, resampled to match vector field shape
-        mask_volume = mask_reader.load_volume(
+        mask = mask_reader.load_volume(
             start_index=start_binned,
             end_index=end_binned,
             unbinned_shape=vec_reader.shape[1:],  # (Z, Y, X)
         )
 
-        mask = (mask_volume > 0).astype(np.uint8)
-
+        mask = (mask > 0).astype(np.uint8)
         vector_field[:, mask == 0] = np.nan
 
     # Seed mask from FA
     print("Generating seed mask using FA > 0.4...")
-    FA_dir = OUTPUT_DIR / "FA"
-    if bin_factor > 1:
-        downsample_volume(
-            input_path=FA_dir,
-            bin_factor=bin_factor,
-            output_dir=OUTPUT_DIR,
-            subfolder="FA",
-            out_ext="tif",
-            min_value=0,
-            max_value=255,
-        )
-        fa_load_dir = OUTPUT_DIR / f"bin{bin_factor}" / "FA"
-    else:
-        fa_load_dir = FA_dir
 
     if not fa_load_dir.exists():
         print(f"⚠️ No FA directory found at {fa_load_dir}")
         sys.exit(1)
 
     print("Loading FA volume...")
     fa_reader = DataReader(fa_load_dir)
-    fa_volume = fa_reader.load_volume(start_index=start_binned, end_index=end_binned)
-    seed_mask = (fa_volume > 0.4).astype(np.uint8)
+    # FA vlume is stored as 8 bits (0-255), so we need to scale it
+    fa_volume = fa_reader.load_volume(start_index=start_binned, end_index=end_binned)    
+    
+    seed_mask = (fa_volume > (0.25 * 255))
 
     print(f"Selecting {num_seeds} random seed points from mask...")
     valid_indices = np.argwhere(seed_mask > 0)
@@ -161,6 +188,7 @@ def script():
     chosen_indices = valid_indices[
         np.random.choice(valid_indices.shape[0], num_seeds, replace=False)
     ]
+    del valid_indices
 
     # Trace streamlines
     streamlines = generate_streamlines_from_vector_field(
@@ -175,25 +203,10 @@ def script():
     )
 
     # Load HA for sampling
-    HA_dir = OUTPUT_DIR / "HA"
     if not HA_dir.exists():
         print(f"⚠️ No HA directory found at {HA_dir}")
         sys.exit(1)
 
-    if bin_factor > 1:
-        downsample_volume(
-            input_path=HA_dir,
-            bin_factor=bin_factor,
-            output_dir=OUTPUT_DIR,
-            subfolder="HA",
-            out_ext="tif",
-            min_value=0,
-            max_value=255,
-        )
-        ha_load_dir = OUTPUT_DIR / f"bin{bin_factor}" / "HA"
-    else:
-        ha_load_dir = HA_dir
-
     ha_reader = DataReader(ha_load_dir)
     HA_volume = ha_reader.load_volume(start_index=start_binned, end_index=end_binned)
 

src/cardiotensor/scripts/visualize_streamlines.py

@@ -34,6 +34,20 @@
 from cardiotensor.utils.utils import read_conf_file
 
 
+def compute_streamline_bounds(streamlines):
+    """
+    Compute min and max bounds in X, Y, Z for a list of 3D streamlines.
+
+    Returns:
+        ((x_min, x_max), (y_min, y_max), (z_min, z_max))
+    """
+    all_points = np.concatenate(streamlines, axis=0)
+    x_min, y_min, z_min = all_points.min(axis=0)
+    x_max, y_max, z_max = all_points.max(axis=0)
+    return (x_min, x_max), (y_min, y_max), (z_min, z_max)
+
+
+
 def script():
     parser = argparse.ArgumentParser(
         description="Visualize streamlines from .npz file."
@@ -90,6 +104,19 @@ def script():
         default=None,
         help="Maximum number of streamlines to render. If not set, no limit is applied.",
     )
+    
+    parser.add_argument(
+        "--crop-z", nargs=2, type=float, metavar=("ZMIN", "ZMAX"),
+        help="Crop streamlines to this Z range (after conversion to (x,y,z))."
+    )
+    parser.add_argument(
+        "--crop-y", nargs=2, type=float, metavar=("YMIN", "YMAX"),
+        help="Crop streamlines to this Y range."
+    )
+    parser.add_argument(
+        "--crop-x", nargs=2, type=float, metavar=("XMIN", "XMAX"),
+        help="Crop streamlines to this X range."
+    )
 
     parser.add_argument(
         "--no-interactive",
@@ -117,6 +144,7 @@ def script():
         help="Window or screenshot height in pixels. Default: 800.",
     )
 
+
     args = parser.parse_args()
     conf_path = Path(args.conf_file)
 
@@ -166,7 +194,23 @@ def script():
         # Convert from [0, 255] to [-90, 90]
         ha_values = ha_values.astype(np.float32)
         color_values = (ha_values / 255.0) * 180.0 - 90.0
-
+        
+    # Build crop bounds with open-ended defaults
+    crop_z = tuple(args.crop_z) if args.crop_z else (-float("inf"), float("inf"))
+    crop_y = tuple(args.crop_y) if args.crop_y else (-float("inf"), float("inf"))
+    crop_x = tuple(args.crop_x) if args.crop_x else (-float("inf"), float("inf"))
+
+    if any([args.crop_z, args.crop_y, args.crop_x]):
+        crop_bounds = (crop_x, crop_y, crop_z)
+    else:
+        crop_bounds = None
+
+    print(f"Original bounds for streamlines:")
+    x_bounds, y_bounds, z_bounds = compute_streamline_bounds(streamlines_xyz)
+    print(f"X bounds: {int(x_bounds[0])} - {int(x_bounds[1])}")
+    print(f"Y bounds: {int(y_bounds[0])} - {int(y_bounds[1])}")
+    print(f"Z bounds: {int(z_bounds[0])} - {int(z_bounds[1])}")
+        
     show_streamlines(
         streamlines_xyz=streamlines_xyz,
         color_values=color_values,
@@ -179,6 +223,7 @@ def script():
         max_streamlines=args.max_streamlines,
         filter_min_len=args.min_length,
         subsample_factor=args.subsample,
+        crop_bounds=crop_bounds,
     )
 
 

src/cardiotensor/tractography/visualize_streamlines.py

@@ -39,6 +39,7 @@ def show_streamlines(
     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,
 ):
     print(f"Initial number of streamlines: {len(streamlines_xyz)}")
     if filter_min_len:
@@ -50,6 +51,50 @@ def show_streamlines(
     if max_streamlines:
         print(f"Limiting to max {max_streamlines} streamlines")
 
+
+
+
+
+
+    # Crop streamlines and color values (point-wise) if bounds are provided
+    print(f"Cropping streamlines within bounds: {crop_bounds}" if crop_bounds else "No cropping applied.")
+    if crop_bounds is not None:
+        z_min, z_max = crop_bounds[2]
+        y_min, y_max = crop_bounds[1]
+        x_min, x_max = crop_bounds[0]
+
+        new_streamlines = []
+        new_color_values = []
+        color_idx = 0
+
+        for sl in streamlines_xyz:
+            n_pts = len(sl)
+            cl = color_values[color_idx : color_idx + n_pts]
+
+            sl = np.asarray(sl)
+            cl = np.asarray(cl)
+
+            within = (
+                (sl[:, 0] >= x_min) & (sl[:, 0] <= x_max) &
+                (sl[:, 1] >= y_min) & (sl[:, 1] <= y_max) &
+                (sl[:, 2] >= z_min) & (sl[:, 2] <= z_max)
+            )
+
+            if np.any(within):  # Keep only remaining points
+                new_sl = sl[within]
+                new_cl = cl[within]
+                if len(new_sl) > 0:
+                    new_streamlines.append(new_sl)
+                    new_color_values.append(new_cl)
+
+            color_idx += n_pts
+
+        streamlines_xyz = new_streamlines
+        color_values = np.concatenate(new_color_values) if new_color_values else np.array([])
+    else:
+        print("No cropping applied.")
+
+
     # Downsample and filter
     downsampled_streamlines = []
     downsampled_colors = []
@@ -88,6 +133,11 @@ def show_streamlines(
 
     print(f"Final number of streamlines to render: {len(streamlines_xyz)}")
 
+
+
+    if not color_values:
+        raise ValueError("❌ No streamlines left after filtering and cropping. Adjust parameters like --crop or --min-length.")
+    
     flat_colors = np.concatenate(color_values)
     print(f"Coloring mode: min={flat_colors.min():.2f}, max={flat_colors.max():.2f}")
     print(f"Rendering mode: {mode}")

src/cardiotensor/utils/downsampling.py

@@ -246,8 +246,6 @@ def process_image_block(
     downsampling in XY, converting to 8-bit, and writing to disk.
 
     Args:
-        z_start (int): Start index of the Z-block.
-        z_end (int): End index of the Z-block (exclusive).
         file_list (list): List of file paths (entire volume stack).
         bin_factor (int): Binning factor for XY downsampling.
         out_file (Path): Output file path for the downsampled image.