Coverage for src/precon3d/resampler.py: 26%
356 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-23 03:50 +0000
1"""Resample slices"""
3# Default python modules
4from enum import Enum
5from numpy.typing import NDArray
6from typing import NamedTuple, Final, Tuple, List
7from pathlib import Path
8import matplotlib.pyplot as plt
9from mpl_toolkits.mplot3d import Axes3D
11# from tqdm import tqdm
12from rich.progress import Progress, BarColumn, TextColumn, TimeElapsedColumn
13from datetime import datetime
14import typer
16import csv
17import os
18import re
19import sys
20import shutil
21import time
23# 3rd party modules
24import numpy as np
25import pandas as pd
26import matplotlib.pyplot as plt
27from scipy import optimize
28from aicspylibczi import CziFile
29import seaborn as sns
31# Local scripts
32import precon3d.czi_info as ci
33import precon3d.utility as ut
34import precon3d.factory as gen_types
36# pylint: disable=wildcard-import
37from precon3d.custom_types import *
38from precon3d._my_typer_cli import CustomCLIGroup, CustomCLICommand
40# CLI
41app = typer.Typer(
42 cls=CustomCLIGroup,
43 no_args_is_help=True,
44 short_help="Resample slices",
45)
47UM_TO_MM: Final[float] = 1.0e-3
50class ResamplePlotOptions(NamedTuple):
51 """Displays the resample plot with the specified plot options."""
53 display: bool
54 raw: bool
55 skip_removed: bool
56 resampled: bool
57 fit: bool
58 legend: bool
59 # segments_legend: bool = False
62class ResampleMethod(Enum):
63 """Defines the available resample methods as an enumeration."""
65 SIMPLE = 1
66 ROBUST = 2
69# Example using the ResampleMethod enum:
70# Person = namedtuple('Person', ['name', 'age', 'method'])
71# Andrew = Person("Andrew", 33, ResampleMethod.ROBUST)
74def tileregion_support_points(czi_file_path: Path) -> List[SupportPoints]:
75 """Extract scene support points data from CZI metadata.
77 Args:
78 czi_file_path (Path): Path to the `.czi` file.
80 Returns:
81 List[SupportPoints]: List of scene support points data.
82 """
83 parent_parent_dir_name = czi_file_path.parent.parent.name
84 parent_dir_name = czi_file_path.parent.name
85 file_basename = czi_file_path.stem
87 czi_fname = f"{parent_parent_dir_name}_{parent_dir_name}_{file_basename}"
89 # Remove all whitespace from the string
90 czi_fname = re.sub(r"\s+", "", czi_fname)
92 # Regular expression to match numerical values
93 def zero_pad(match):
94 return match.group().zfill(3) # Pad the matched number with zeros
96 # Replace all numerical values with zero-padded versions
97 czi_fname_zeropadded = re.sub(r"\d+", zero_pad, czi_fname)
99 metadata = CziFile(czi_file_path).meta
101 all_scene_sp = []
102 for each_scene in metadata.findall(".//TileRegion"):
103 if each_scene.find(".//IsUsedForAcquisition").text == "true":
104 scene_name = each_scene.attrib["Name"]
105 scene_z = each_scene.find(".//Z").text
107 sp_list = []
108 for each_sp in each_scene.findall(".//SupportPoints//SupportPoint"):
109 sp_list.append(
110 SupportPoint(
111 x=float(each_sp.find("X").text),
112 y=float(each_sp.find("Y").text),
113 z=float(each_sp.find("Z").text),
114 )
115 )
117 all_scene_sp.append(
118 SupportPoints(
119 z_height=scene_z,
120 positions=sp_list,
121 czi_fname=czi_fname_zeropadded,
122 scene_name=scene_name,
123 )
124 )
126 return all_scene_sp
129def fit_plane(points: list[SupportPoints]) -> tuple[np.ndarray, float]:
130 """
131 Fits a plane to the given support points and calculates the normal vector and tilt angle.
133 Args:
134 points (list[tuple[float, float, float]]): List of support points as (x, y, z) tuples.
136 Returns:
137 tuple[np.ndarray, float]: A tuple containing the unit normal vector (as a NumPy array)
138 """
139 # Extract coordinates from the support points
140 X = np.array([p.x for p in points])
141 Y = np.array([p.y for p in points])
142 Z = np.array([p.z for p in points])
144 # Fit a plane using least squares
145 A = np.c_[X, Y, np.ones(X.shape)]
146 C, _, _, _ = np.linalg.lstsq(A, Z, rcond=None)
148 # Plane equation: Z = C[0]*X + C[1]*Y + C[2]
149 normal_vector = np.array([C[0], C[1], -1])
150 normal_vector_magnitude = np.linalg.norm(normal_vector)
151 normal_vector_unit = normal_vector / normal_vector_magnitude
153 # Round the vector to remove extra whitespace and ensure clean formatting
154 normal_vector_unit = np.round(normal_vector_unit, decimals=6)
156 # Calculate the tilt angle from the z-axis
157 z_axis = np.array([0, 0, -1])
158 cos_theta = np.dot(normal_vector_unit, z_axis)
159 tilt_angle = np.arccos(cos_theta) * (180 / np.pi) # Convert to degrees
161 # Return the unit normal vector and the tilt angle
162 return normal_vector_unit, float(tilt_angle)
165def fit_plane_and_plot(points: list[SupportPoints]) -> None:
166 """
167 Fits a plane to the given support points, plots the plane, normal vector, and calculates the tilt angle from the z-axis.
169 Args:
170 points (list[SupportPoint]): List of support points with x, y, z coordinates.
172 Returns:
173 None
174 """
175 # Extract coordinates from the support points
176 X = np.array([p.x for p in points])
177 Y = np.array([p.y for p in points])
178 Z = np.array([p.z for p in points])
180 # Fit a plane using least squares
181 A = np.c_[X, Y, np.ones(X.shape)]
182 C, _, _, _ = np.linalg.lstsq(A, Z, rcond=None)
184 # Plane equation: Z = C[0]*X + C[1]*Y + C[2]
185 normal_vector = np.array([C[0], C[1], -1])
186 normal_vector_magnitude = np.linalg.norm(normal_vector)
187 normal_vector_unit = normal_vector / normal_vector_magnitude
189 # Calculate the tilt angle from the z-axis
190 z_axis = np.array([0, 0, -1])
191 cos_theta = np.dot(normal_vector_unit, z_axis)
192 tilt_angle = np.arccos(cos_theta) * (180 / np.pi) # Convert to degrees
194 # Create a grid to plot the plane
195 x_range = np.linspace(X.min(), X.max(), 10)
196 y_range = np.linspace(Y.min(), Y.max(), 10)
197 X_grid, Y_grid = np.meshgrid(x_range, y_range)
198 Z_grid = C[0] * X_grid + C[1] * Y_grid + C[2]
200 # Plotting
201 fig = plt.figure(figsize=(10, 8))
202 ax = fig.add_subplot(111, projection="3d")
204 # Plot the original points
205 ax.scatter(X, Y, Z, color="r", label="Support Points")
207 # Plot the fitted plane
208 ax.plot_surface(
209 X_grid, Y_grid, Z_grid, alpha=0.25, color="b", rstride=100, cstride=100
210 )
212 # Plot the normal vector
213 origin = np.mean(X), np.mean(Y), np.mean(Z)
214 ax.quiver(
215 *origin,
216 *normal_vector_unit,
217 length=20,
218 color="k",
219 label="Normal Vector",
220 )
222 ax.quiver(*origin, 1, 0, 0, length=100, color="r", label="X")
223 ax.quiver(*origin, 0, 1, 0, length=100, color="g", label="Y")
224 ax.quiver(*origin, 0, 0, 1, length=100, color="b", label="Z")
226 # Set labels
227 ax.set_xlabel("X")
228 ax.set_ylabel("Y")
229 ax.set_zlabel("Z")
230 ax.set_title("Fitted Plane and Normal Vector")
231 ax.legend()
233 ax.set_box_aspect((1, 1, 0.1))
234 # ax.auto_scale_xyz([0, 500], [0, 500], [0, 0.15])
236 # Display the equation of the plane and tilt angle
237 print(f"Plane equation: Z = {C[0]:.4f}*X + {C[1]:.4f}*Y + {C[2]:.4f}")
238 print(f"Magnitude of tilt from Z-axis: {tilt_angle:.4f} degrees")
240 plt.show()
243@app.command(
244 cls=CustomCLICommand,
245 name="main",
246 short_help="save support point information",
247)
248def save_czi_tileregion_support_points(configfile: Path):
249 """
250 Extracts support points from CZI files in the input directory and saves them to a CSV file.
252 Parameters
253 ----------
254 params : Path to config file with general attrs
255 """
257 # cast to Path object if not already a Path object
258 if not isinstance(configfile, Path):
259 configfile = Path(configfile)
261 config_dict = ut.read_config(configfile)
262 params = gen_types.create_resampler_parameters(config_dict)
264 # Validate file extension
265 if params.file_extension != ".czi":
266 raise ValueError("At this time, only .czi file extension is supported")
268 # Input and output directories
269 czi_input_directory = params.input_directory
270 output_directory = params.output_directory
271 output_directory.mkdir(parents=True, exist_ok=True)
273 # Sort the CZI files in the specified input directory
274 sorted_czi_files = ut.sorted_files(
275 directory=czi_input_directory, file_extension=params.file_extension
276 )
278 # Prepare the output file path
279 output_file_name = re.sub(
280 r"\s+",
281 "",
282 f"{czi_input_directory.parent.parent.name}_{czi_input_directory.parent.name}_{czi_input_directory.name}_support_points.csv",
283 )
284 output_file_path = output_directory / output_file_name
286 print(f"Found {len(sorted_czi_files)} CZI files in {czi_input_directory}\n")
287 print(f"Support point information will be saved in {output_file_path}\n")
288 print(f'{" Reading the files ":*^28}')
290 # Progress bar setup
291 total_files = len(sorted_czi_files)
292 with Progress(
293 TextColumn("[progress.description]{task.description}"),
294 BarColumn(),
295 TextColumn("[progress.percentage]{task.percentage:>3.0f}%"),
296 TextColumn("[bold blue]{task.completed}/{task.total}"),
297 ) as progress:
298 task = progress.add_task("Reading...", total=total_files)
300 # Open the CSV file for writing
301 with open(output_file_path, mode="w", newline="") as csv_file:
302 csv_writer = csv.writer(csv_file)
304 # Write the header row to the CSV file
305 csv_writer.writerow(
306 [
307 "File Name",
308 "Scene",
309 "z height from metadata",
310 "Support Point #",
311 "Support Point X",
312 "Support Point Y",
313 "Support Point Z",
314 "normal vector",
315 "angle from z",
316 ]
317 )
319 # Process each CZI file
320 for each_czi_path in sorted_czi_files:
321 try:
322 # Extract support points from the current CZI file
323 tr_sp = tileregion_support_points(each_czi_path)
325 # Write support points to the CSV file
326 for scene_sp in tr_sp:
327 normal_vec, angle_from_z = fit_plane(scene_sp.positions)
328 for idx, each_sp in enumerate(scene_sp.positions, start=1):
329 csv_writer.writerow(
330 [
331 scene_sp.czi_fname,
332 scene_sp.scene_name,
333 scene_sp.z_height,
334 idx,
335 each_sp.x,
336 each_sp.y,
337 each_sp.z,
338 normal_vec,
339 angle_from_z,
340 ]
341 )
343 except RuntimeError as error:
344 # Print an error message if the file cannot be processed
345 print(f"{each_czi_path} has an error: {error}")
346 continue
347 finally:
348 progress.update(task, advance=1)
350 print(f"Support points successfully saved to {output_file_path}")
353def plot_support_points(csv_file: Path):
354 """
355 Reads support point data from a CSV file and plots the z-positions of the support points
356 for each slice using Seaborn.
358 Parameters
359 ----------
360 csv_file : str
361 Path to the CSV file containing support point data.
363 Returns
364 -------
365 None
366 """
367 # Read the CSV file into a Pandas DataFrame
368 df = pd.read_csv(csv_file)
370 # Ensure the data is sorted by slice and support point number
371 df = df.sort_values(by=["File Name", "Support Point #"])
373 # Plot the z-positions for each slice
374 plt.figure(figsize=(12, 6))
375 sns.lineplot(
376 data=df,
377 x="File Name", # Slice (x-axis)
378 y="Support Point Z", # Z position (y-axis)
379 hue="Support Point #", # Support Point # (overlaid lines)
380 marker="o",
381 )
383 # Customize the plot
384 plt.title("Z Positions of Support Points Across Slices", fontsize=16)
385 plt.xlabel("Slice (File Name)", fontsize=14)
386 plt.ylabel("Z Position", fontsize=14)
387 plt.xticks(rotation=45, fontsize=10)
388 plt.legend(title="Support Point #", fontsize=10)
389 plt.grid(True, linestyle="--", alpha=0.5)
391 # Show the plot
392 plt.tight_layout()
393 plt.show()
396def plot_support_points_and_metadata(csv_file: Path):
397 """
398 Reads support point data from a CSV file and creates multiple plots:
399 1. Z positions of support points across slices.
400 2. Z height from metadata across slices.
401 3. Angle from z-axis across slices.
402 4. Polar plot visualizing normal vectors looking down the z-axis.
404 Parameters
405 ----------
406 csv_file : str
407 Path to the CSV file containing support point data.
409 Returns
410 -------
411 None
412 """
413 # Read the CSV file into a Pandas DataFrame
414 df = pd.read_csv(csv_file)
416 # Ensure the data is sorted by slice and support point number
417 df = df.sort_values(by=["File Name", "Support Point #"])
419 # Preprocess the normal vector column to fix missing commas
420 def parse_normal_vector(vec: str):
421 """
422 Parses a string representation of a normal vector into a NumPy array,
423 ensuring no leading commas and replacing internal whitespace with commas.
424 """
426 # Remove leading and trailing whitespace
427 vec = vec.strip()
429 # Use regular expressions to replace internal whitespace with commas
430 vec = re.sub(r"\s+", ",", vec) # Replace all consecutive whitespace with commas
432 if vec.startswith("[,"):
433 vec = "[" + vec[2:]
435 return np.array(eval(vec)) # Convert string to NumPy array
437 df["normal vector"] = df["normal vector"].apply(parse_normal_vector)
439 # Create subplots
440 fig, axes = plt.subplots(
441 2, 2, figsize=(14, 10), gridspec_kw={"width_ratios": [2, 1]}
442 )
443 fig.suptitle(f"{csv_file.stem}\nSupport Points Analysis", fontsize=16)
445 def reduce_filename_length_vectorized(filenames: pd.Series) -> pd.Series:
446 """
447 Reduces filename length by removing repeated substrings between underscores
448 using vectorized operations.
450 Parameters
451 ----------
452 filenames : pd.Series
453 Pandas Series containing filenames.
455 Returns
456 -------
457 pd.Series
458 Pandas Series with reduced filenames.
459 """
460 # Extract all substrings between underscores
461 substrings = filenames.str.extractall(r"_(.*?)_")[0]
463 # Find common substrings across all filenames
464 common_substrings = substrings.value_counts()[
465 substrings.value_counts() == len(filenames)
466 ].index
468 # Remove common substrings from filenames using vectorized string replacement
469 reduced_filenames = filenames
470 for common in common_substrings:
471 reduced_filenames = reduced_filenames.str.replace(
472 f"_{common}_", "_", regex=False
473 )
475 return reduced_filenames
477 # Reduce filename length
478 df["Reduced File Name"] = reduce_filename_length_vectorized(df["File Name"])
480 # Adjust x-ticks spacing
481 x_ticks = df["Reduced File Name"].unique()
483 # Plot 1: Z positions of support points across slices
484 sns.lineplot(
485 data=df,
486 x="Reduced File Name",
487 y="Support Point Z",
488 hue="Support Point #",
489 marker=".",
490 ax=axes[0, 0],
491 )
492 axes[0, 0].set_title("Z Positions of Support Points Across Slices", fontsize=14)
493 axes[0, 0].set_xlabel("Slice (File Name)", fontsize=12)
494 axes[0, 0].set_ylabel("Z Position", fontsize=12)
495 axes[0, 0].tick_params(axis="x", rotation=45)
496 axes[0, 0].grid(True, linestyle="--", alpha=0.5)
498 # Adjust x-ticks spacing
499 axes[0, 0].set_xticks(x_ticks[::50]) # Show every second slice for better spacing
501 # Plot 2: Z height from metadata across slices
502 sns.lineplot(
503 data=df,
504 x="Reduced File Name",
505 y="z height from metadata",
506 markers=False,
507 ax=axes[0, 1],
508 )
509 axes[0, 1].scatter(
510 df["Reduced File Name"][::10], # Select every 2nd point for x-axis
511 df["z height from metadata"][::10], # Select every 2nd point for y-axis
512 color="blue",
513 alpha=0.7,
514 s=10,
515 )
516 axes[0, 1].set_title("Z Height from Metadata Across Slices", fontsize=14)
517 axes[0, 1].set_xlabel("Slice (File Name)", fontsize=12)
518 axes[0, 1].set_ylabel("Z Height", fontsize=12)
519 axes[0, 1].tick_params(axis="x", rotation=45)
520 axes[0, 1].grid(True, linestyle="--", alpha=0.5)
522 # Adjust x-ticks spacing
523 axes[0, 1].set_xticks(x_ticks[::50]) # Show every second slice for better spacing
525 # Plot 3: Angle from z-axis across slices
526 sns.lineplot(
527 data=df,
528 x="Reduced File Name",
529 y="angle from z",
530 marker=".",
531 ax=axes[1, 0],
532 )
533 axes[1, 0].set_title("Angle from Z-Axis Across Slices", fontsize=14)
534 axes[1, 0].set_xlabel("Slice (File Name)", fontsize=12)
535 axes[1, 0].set_ylabel("Angle (Degrees)", fontsize=12)
536 axes[1, 0].tick_params(axis="x", rotation=45)
537 axes[1, 0].grid(True, linestyle="--", alpha=0.5)
539 # Adjust x-ticks spacing
540 axes[1, 0].set_xticks(x_ticks[::50]) # Show every second slice for better spacing
542 # Plot 4: Polar plot of normal vectors
543 axes[1, 1] = plt.subplot(2, 2, 4, polar=True)
544 # for _, row in df.iterrows():
545 # normal_vec = row["normal vector"]
546 # angle = np.arctan2(normal_vec[1], normal_vec[0]) # Angle in radians
547 # magnitude = np.linalg.norm(normal_vec[:2]) # Magnitude in the XY plane
548 # axes[1, 1].scatter(angle, magnitude, color='r', s=5)
549 # Extract normal vectors as a NumPy
551 # Randomly sample 100 rows (or fewer if the DataFrame has less than 100 rows)
552 sampled_df = df.sample(n=min(100, len(df)), random_state=42)
554 # Extract normal vectors as a NumPy array
555 normal_vectors = np.array(sampled_df["normal vector"].tolist())
557 # Calculate angles and magnitudes using vectorized operations
558 angles = np.arctan2(normal_vectors[:, 1], normal_vectors[:, 0]) # Angles in radians
559 magnitudes = np.linalg.norm(
560 normal_vectors[:, :2], axis=1
561 ) # Magnitudes in the XY plane
563 axes[1, 1].scatter(angles, magnitudes, color="r", s=5)
564 axes[1, 1].spines["polar"].set_visible(False) # This is key for polar plots
565 axes[1, 1].set_title(
566 "Polar Plot of Normal Vectors (Looking Down Z-Axis)", fontsize=14
567 )
568 axes[1, 1].set_theta_zero_location("N") # North points up
569 axes[1, 1].set_theta_direction(-1) # Clockwise direction
570 # axes[1, 1].axis('off')
572 # axes[1, 1].legend(loc="upper right", fontsize=8)
574 # Adjust layout and show the plots
575 plt.tight_layout(rect=[0, 0, 1, 0.95])
576 plt.show()
579def save_support_point_info(params: GeneralAttrs) -> Path:
580 """
581 Save support point information extracted from CZI files.
583 This function reads CZI files from the specified input directory, extracts support points,
584 and saves the information into a CSV file in the output directory. Currently, only files
585 with the '.czi' extension are supported.
587 Parameters
588 ----------
589 params : ut.GeneralAttrs
590 An object containing the following attributes:
591 - file_extension (str): The file extension to filter files (must be '.czi').
592 - input_directory (Path): The directory containing the input CZI files.
593 - output_directory (Path): The directory where the output CSV file will be saved.
595 Returns
596 -------
597 Path
598 The path to the CSV file containing the support point information.
600 Raises
601 ------
602 ValueError
603 If the file extension is not '.czi'.
605 Notes
606 -----
607 - The output CSV file will be named 'support_points.csv'.
608 - If an error occurs while processing a CZI file, a message will be printed, and the
609 function will continue with the next file.
610 - The output directory will be created if it does not exist.
611 """
613 ## Tile reading, CZI only
614 file_extension = params.file_extension
615 if file_extension != ".czi":
616 raise ValueError("At this time, only .czi file extension is supported")
618 czi_input_directory = params.input_directory
619 # Sort the CZI files in the specified input directory
620 sorted_czi_files = ut.sorted_files(
621 directory=czi_input_directory, file_extension=file_extension
622 )
624 output_directory = params.output_directory
625 # Create the output directory if it doesn't exist
626 output_directory.mkdir(parents=True, exist_ok=True)
628 output_file_name = "support_points.csv"
629 file_path_czi_sp_out = output_directory.joinpath(output_file_name)
630 file_path_czi_info = output_directory.joinpath("czi_info.csv")
632 print(f"Found {len(sorted_czi_files)} czi files in {czi_input_directory}\n")
633 print(f"Support point information will be saved in {file_path_czi_sp_out}\n")
635 print(f'{" Reading the files ":*^28}')
637 total_files = len(sorted_czi_files)
639 with Progress(
640 TextColumn("[progress.description]{task.description}"),
641 BarColumn(),
642 TextColumn("[progress.percentage]{task.percentage:>3.0f}%"),
643 TextColumn("[bold blue]{task.completed}/{task.total}"),
644 ) as progress:
645 task = progress.add_task("Reading...", total=total_files)
647 for each_czi_path in sorted_czi_files:
648 try:
649 # Extract support points from the current CZI file
650 all_scene_sp = ci.scene_support_points(each_czi_path)
652 # check if empty
653 if all_scene_sp: # case one, there are tile regions and tile arrays
654 if all_scene_sp[0].z_heights.size == 0:
655 all_scene_sp = ci.scene_support_points_global(each_czi_path)
656 else: # case two, empty, only tile arrays
657 all_scene_sp = ci.scene_support_points_global(each_czi_path)
658 except RuntimeError as error:
659 # Print an error message if the file cannot be processed
660 # This scenario usually occurs if the czi is corrupted
661 print(f"{each_czi_path} has an error: {error}")
662 continue
663 finally:
664 progress.update(task, advance=1)
665 # Dump the support point data to a CSV file
666 dump_czi_support_point_info(all_scene_sp, file_path_czi_sp_out)
667 dump_file_info_to_csv(each_czi_path, file_path_czi_info)
668 # Uncomment the following lines to dump scene-specific support point data
669 # for each_scene_sp in all_scene_sp:
670 # # Determine the file name based on the scene name
671 # file_path = output_directory / f"{each_scene_sp.scene}_raw_sp_info.csv"
672 # dump_scene_support_point_info(each_scene_sp, file_path)
674 return file_path_czi_sp_out
677def dump_czi_support_point_info(
678 all_sp_info: list[SceneSupportPoints], csv_out_file: Path
679):
680 """to come"""
681 czi_fname = all_sp_info[0].czi_fname
683 all_z_heights = np.array([])
684 for each_scene_sp in all_sp_info:
685 # Ensure each_scene_sp.z_heights is a NumPy array
686 each_scene_sp_z_height = np.array(each_scene_sp.z_heights)
687 # Concatenate the current z_heights with the accumulated ones
688 all_z_heights = np.concatenate((all_z_heights, each_scene_sp_z_height))
690 # compute slice average
691 average_z_height = np.mean(all_z_heights)
693 # Check if the file exists and is empty to decide on writing the header
694 file_exists = csv_out_file.exists()
696 with open(csv_out_file, "a", newline="", encoding="utf-8") as csvfile:
697 sp_writer = csv.writer(csvfile)
699 # If the file didn't exist or was empty, write the header
700 if not file_exists or csv_out_file.stat().st_size == 0:
701 sp_writer.writerow(
702 [
703 "File name",
704 "Average height (um)",
705 "All heights (um)",
706 ]
707 )
709 # Write the data row
710 sp_writer.writerow([czi_fname, average_z_height, all_z_heights])
713def dump_file_info_to_csv(file_path: Path, csv_file: Path):
714 # Ensure the file_path is a Path object
715 # file_path = Path(file_path)
717 # Check if the file exists
718 if not file_path.exists():
719 print(f"The file {file_path} does not exist.")
720 return
722 # Extracting information
723 parent1 = file_path.parent.name
724 parent2 = file_path.parent.parent.name if file_path.parent.parent else ""
725 file_name = file_path.name
726 file_size = file_path.stat().st_size
727 # creation_time = datetime.fromtimestamp(file_path.stat().st_ctime).strftime('%Y-%m-%d %H:%M:%S')
728 creation_time = ci.acquisition_date_and_time(file_path)
730 # Data to write
731 data = [parent2, parent1, file_name, file_size, creation_time]
733 # Check if the file exists and is empty to decide on writing the header
734 file_exists = csv_file.exists()
735 # Writing to CSV
736 with open(csv_file, mode="a", newline="", encoding="utf-8") as file:
737 writer = csv.writer(file)
739 # If the file didn't exist or was empty, write the header
740 if not file_exists or csv_file.stat().st_size == 0:
741 writer.writerow(
742 [
743 "Parent Dir2",
744 "Parent Dir1",
745 "Filename",
746 "File size",
747 "Creation time",
748 ]
749 )
750 writer.writerow(data)
752 # print(f"Data for {file_name} has been written to {csv_file}")
755def dump_scene_support_point_info(sp_info: SceneSupportPoints, csv_out_file: Path):
756 """Dump scene support point information to a CSV file.
758 Parameters
759 ----------
760 sp_info : SceneSupportPoints
761 An object containing information about the scene, including the filename,
762 scene name, and z heights of support points.
763 csv_out_file : Path
764 The path to the CSV file where the information will be written.
766 Notes
767 -----
768 If the specified CSV file does not exist or is empty, a header row will be
769 written. Each subsequent call will append a new row with the scene information.
770 """
771 czi_fname = sp_info.czi_fname
772 scene_name = sp_info.scene
773 if len(sp_info.z_heights) == 0: # empty
774 average_z_height = sp_info.z
775 else:
776 average_z_height = np.mean(sp_info.z_heights)
777 # breakpoint()
779 # Check if the file exists and is empty to decide on writing the header
780 file_exists = csv_out_file.exists()
782 with open(csv_out_file, "a", newline="", encoding="utf-8") as csvfile:
783 sp_writer = csv.writer(csvfile)
785 # If the file didn't exist or was empty, write the header
786 if not file_exists or csv_out_file.stat().st_size == 0:
787 sp_writer.writerow(
788 [
789 "File name",
790 "Scene name",
791 "Average height (um)",
792 "All heights (um)",
793 # "All Support Points (x,y,z)",
794 ]
795 )
797 # Write the data row
798 sp_writer.writerow([czi_fname, scene_name, average_z_height, sp_info.z_heights])
801def raw_z_heights(csv_file_path: Path, pol_images: bool) -> np.array:
802 """Returns the z height position from the csv file.
804 Args:
805 csv_file_path: The fully pathed csv file.
807 Returns:
808 The z-height positions.
809 """
811 df = pd.read_csv(csv_file_path)
813 # Check if the "Average height (um)" column exists
814 if "Average height (um)" in df.columns:
815 height_data = df.get("Average height (um)")
816 else:
817 # If not found, get unique z heights for each unique filename
818 unique_heights = df[["File Name", "z height from metadata"]].drop_duplicates()
819 height_data = unique_heights["z height from metadata"]
821 height_arr = height_data.to_numpy()
823 if pol_images:
824 return height_arr[1::13]
826 return height_arr
829def process_scene_resample(
830 raw_z_csv: Path,
831 target_thickness: float, # units of micron
832 plot_opts: ResamplePlotOptions,
833 pol_images: bool,
834 skip_slice_idx: np.array = np.array(0),
835) -> bool:
836 """Resamples an individual scene for a dataset to provide 3D slices with
837 nearly uniform slice thickness.
839 Assumes we are working with directories of images, not HDF files to store the origin data
841 Resample the slice images to be near uniform thickness
843 Args:
844 raw_z_csv: The csv file containing the raw z values.
846 Returns:
847 True if the function succeeded, False otherwise.
848 """
850 # 1
852 processed = False # by default, the fuction has not yet been processed
854 all_slice_z = raw_z_heights(raw_z_csv, pol_images)
856 # all_slice_czi_paths = ut.sorted_files(Path(czi_input_dir), "czi")
857 # all_slice_tif_paths = io.sorted_files(scene_input_dir, "tif")
859 # check the file indices match properly, same sizes
860 # assert len(all_slice_czi_paths) == len(
861 # all_slice_tif_paths
862 # ), "number of czi files used to calculate z heights do not match the number of image slices attempting to resample."
864 # IGNORE FOR NOW
865 # assert all_slice_z.shape[0] == len(
866 # all_slice_czi_paths
867 # ), "number of z heights in the csv do not match the number of czi files used to calculate z heights."
869 # assert all_slice_z.shape[0] == len(
870 # all_slice_tif_paths
871 # ), "number of z heights in the csv do not match the number of image slices attempting to resample."
873 # 2 return new slice indices
874 resampled_idx = scene_resample_simple(all_slice_z, target_thickness, skip_slice_idx)
876 # plot and verify depending on plot options
878 # for each channel in the scene, place files into new directory
879 # (could be with mv or cp depending on how much we trust ourselves)
881 if plot_opts.display:
882 fig, ax = plt.subplots(figsize=(6, 5), dpi=150)
883 ax.set_ylabel("Z Position (mm)", fontweight="bold")
884 ax.set_xlabel("Slice Number", labelpad=5, fontweight="bold")
886 if plot_opts.raw:
887 plt.scatter(
888 np.arange(0, len(all_slice_z)),
889 all_slice_z * UM_TO_MM,
890 s=3,
891 label="All Slices",
892 color="k",
893 )
894 if plot_opts.skip_removed:
895 skip_remove_z = np.delete(all_slice_z, skip_slice_idx)
896 plt.scatter(
897 np.arange(0, len(skip_remove_z)),
898 skip_remove_z * UM_TO_MM,
899 s=3,
900 label="Skipped Slice Removed",
901 color="orange",
902 alpha=0.5,
903 )
904 if plot_opts.resampled:
905 plt.scatter(
906 np.arange(0, len(resampled_idx)),
907 all_slice_z[resampled_idx] * UM_TO_MM,
908 s=3,
909 label="Resampled Height",
910 color="gray",
911 )
912 if plot_opts.fit:
913 q = fit_linear_segments(
914 x_values=np.arange(0, len(resampled_idx)),
915 y_values=all_slice_z[resampled_idx] * UM_TO_MM,
916 segment_count=1,
917 )
918 qx = np.array([q[0][1], q[0][0]])
919 qy = np.array([q[1][1], q[1][0]])
920 resampled_slope = (q[1][1] - q[1][0]) / (q[0][1] - q[0][0])
922 plt.plot(
923 qx,
924 qy,
925 label=f"Overall {np.around(resampled_slope, 4)*(1/(UM_TO_MM))} um/slice",
926 color="k",
927 linestyle="--",
928 )
929 title = f"Resampling for target thickness of\n{target_thickness} micron"
930 if plot_opts.legend:
931 plt.legend()
932 plt.grid()
933 ax.set_axisbelow(True)
934 plt.title(title, fontweight="bold")
935 plt.tight_layout()
936 plt.show()
938 if not ut.yes_no("\tContinue with resampling?"):
939 print("Exiting data processor...")
940 sys.exit(0)
942 # new_slice_idx = 0
943 # for slice_idx in resampled_idx:
944 # file_dest = os.path.join(scene_output_dir, f"Slice_{new_slice_idx:>04}")
946 # # COPY FILES
947 # shutil.copyfile(all_slice_tif_paths[slice_idx], file_dest)
949 # new_slice_idx += 1
951 processed = True # overwrite, if we get to this point, the functiion succeede
953 return processed
956def scene_resample_simple(
957 z_heights: NDArray, # Array of original z heights
958 target_thickness: float, # Desired thickness of the slices in microns
959 skip_slice_idx: NDArray = np.array([], dtype=int), # Indices of slices to skip
960) -> NDArray: # Indices of the resampled heights in the original z_heights array
961 """Resamples an individual scene for a dataset to provide 3D slices with nearly uniform slice thickness.
963 Args:
964 z_heights (np.ndarray): Array of original z heights.
965 target_thickness (float): Desired thickness of the slices in microns.
966 skip_slice_idx (np.ndarray, optional): Array of indices to skip during resampling. Defaults to an empty array.
968 Returns:
969 np.ndarray: Indices of the resampled heights in the original z_heights array.
970 """
972 # Ensure skip_slice_idx is a 1D array
973 skip_slice_idx = np.asarray(skip_slice_idx).flatten()
975 # Remove skipped slice indices from z_heights
976 if skip_slice_idx.size > 0:
977 sel_z_heights = np.delete(z_heights, skip_slice_idx)
978 else:
979 sel_z_heights = z_heights
981 # Generate target heights based on the original z_heights range
982 target_heights = np.arange(
983 sel_z_heights[0], sel_z_heights[-1] + target_thickness, target_thickness
984 )
986 # Resample heights based on the target heights
987 resampled_heights = []
988 for height in target_heights:
989 diff = np.abs(sel_z_heights - height)
990 index = diff.argmin()
991 resampled_heights.append(sel_z_heights[index])
993 # Convert the result to a NumPy array
994 resampled_heights = np.asarray(resampled_heights)
996 # Get the indices of the resampled heights in the original z_heights
997 resampled_heights_indices = np.array(
998 [np.argwhere(z_heights == val)[0][0] for val in resampled_heights]
999 )
1001 return resampled_heights_indices
1004def fit_linear_segments(
1005 x_values: NDArray[np.float64],
1006 y_values: NDArray[np.float64],
1007 segment_count: int,
1008) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
1009 """
1010 Fit linear segments to a set of points.
1011 """
1012 if segment_count >= len(x_values):
1013 raise ValueError(
1014 "The number of segments must be less than the number of data points."
1015 )
1017 # Calculate the range of x values
1018 x_min = x_values.min()
1019 x_max = x_values.max()
1021 # Improved initial guess for segment boundaries
1022 # Distribute segments evenly across the x range
1023 segment_boundaries_x = np.linspace(x_min, x_max, segment_count + 1)
1025 # For y-values at segment boundaries, interpolate based on the existing points
1026 segment_boundaries_y = np.interp(segment_boundaries_x, x_values, y_values)
1028 def error_function(params):
1029 """
1030 Calculate the mean squared error between the model and the actual y-values.
1031 """
1032 # Reconstruct the y-values at segment boundaries from params
1033 y_values_at_boundaries = params
1034 # Use fixed x segment boundaries and interpolated y-values for error calculation
1035 y_predicted = np.interp(x_values, segment_boundaries_x, y_values_at_boundaries)
1036 return np.mean((y_values - y_predicted) ** 2)
1038 # Optimize only the y-values at segment boundaries
1039 initial_guess = segment_boundaries_y
1040 result = optimize.minimize(error_function, x0=initial_guess, method="Nelder-Mead")
1042 # Use the optimized y-values with the original x segment boundaries
1043 optimized_y = result.x
1045 return segment_boundaries_x, optimized_y
1048@app.callback()
1049def callback():
1050 """
1051 precon3d.resampler tries to evenly space out the slices collected based on the microscope focus height.
1052 """