# Load all images to create the hyperspectral stack
hyperspectral_stack = []
for tif_file in tqdm(tiff_files, desc="Loading TIFF files", leave=False):
    hyperspectral_stack.append(imread(tif_file))

print(type(hyperspectral_stack))
print(len(hyperspectral_stack))
print(np.array(hyperspectral_stack).shape)
# Stack the images along a new dimension
print(f"")

new_hyperspectal_stack = np.swapaxes(hyperspectral_stack, 0, 2)
print(type(new_hyperspectal_stack))
print(len(new_hyperspectal_stack))
print(np.array(new_hyperspectal_stack).shape)


hyperspectral_stack = np.stack(hyperspectral_stack, axis=-1)
print(type(hyperspectral_stack))
print(len(hyperspectral_stack))
print(np.array(hyperspectral_stack).shape)

<class 'list'>
2782
(2782, 512, 512)

<class 'numpy.ndarray'>
512
(512, 512, 2782)

<class 'numpy.ndarray'>
512
(512, 512, 2782)

import numpy as np
import matplotlib.pyplot as plt
from scipy.datasets import electrocardiogram
from scipy.signal import find_peaks
x = electrocardiogram()[2000:4000]
peaks, _ = find_peaks(x, height=0)
plt.plot(x)
plt.plot(peaks, x[peaks], "x")
plt.plot(np.zeros_like(x), "--", color="gray")
plt.show()

class CylinderGeometry(BaseModel):
    """
    Cylindrical geometry in pixel coordinates (origin: top-left).
    Only store primitives; derive the rest via properties.
    """
    left_edge: int = Field(..., ge=0)
    right_edge: int = Field(..., ge=0)
    top_edge: int = Field(..., ge=0)
    bottom_edge: int = Field(..., ge=0)
    center_x: int = Field(..., ge=0)
    center_y: int = Field(..., ge=0)

    @model_validator(mode="after")
    def _validate_bounds(self):
        if self.right_edge <= self.left_edge:
            raise ValueError("right_edge must be > left_edge")
        if self.bottom_edge <= self.top_edge:
            raise ValueError("bottom_edge must be > top_edge")
        return self

    @property
    def width(self) -> int:
        return self.right_edge - self.left_edge

    @property
    def height(self) -> int:
        return self.bottom_edge - self.top_edge + 1

    @property
    def radius(self) -> int:
        return self.width // 2

    def __str__(self) -> str:
        return (
            "CylinderGeometry(\n"
            f"  Boundaries: left={self.left_edge}, right={self.right_edge}, "
            f"top={self.top_edge}, bottom={self.bottom_edge}\n"
            f"  Center: ({self.center_x}, {self.center_y})\n"
            f"  Dimensions: width={self.width}, height={self.height}, radius={self.radius}\n"
            ")"
        )

import numpy as np
x = np.arange(6).reshape(2,3)
x
array([[0, 1, 2],
       [3, 4, 5]])
np.argwhere(x>1)
array([[0, 2],
       [1, 0],
       [1, 1],
       [1, 2]])

import numpy as np
a = np.array([[10., 7., 4.], [3., 2., 1.]])
a[0][1] = np.nan
a
np.percentile(a, 50)
np.nanpercentile(a, 50)
np.nanpercentile(a, 50, axis=0)
np.nanpercentile(a, 50, axis=1, keepdims=True)
m = np.nanpercentile(a, 50, axis=0)
out = np.zeros_like(m)
np.nanpercentile(a, 50, axis=0, out=out)
m

def visualize_hyperspectral_radiographs(hyperspectral_stack: np.ndarray, 
                                       selected_indices: Optional[list[int]] = None, 
                                       figsize: tuple[int, int] = (15, 7), 
                                       cmap: str = "gray", 
                                       percentile_range: tuple[float, float] = (2, 98)) -> plt.Figure:
    """
    Visualize selected radiographs from a hyperspectral stack.
    
    Parameters:
    -----------
    hyperspectral_stack : numpy.ndarray
        3D array with shape (height, width, n_spectral_bins)
    selected_indices : list of int, optional
        Indices of radiographs to visualize. If None, uses beginning, middle, and end.
    figsize : tuple, optional
        Figure size (width, height). Default is (15, 7).
    cmap : str, optional
        Colormap for displaying images. Default is "gray".
    percentile_range : tuple, optional
        Percentile range for intensity scaling (min, max). Default is (2, 98).
    
    Returns:
    --------
    fig : matplotlib.figure.Figure
        The created figure object
    """
    # Default to beginning, middle, and end if no indices provided
    if selected_indices is None:
        n_images = hyperspectral_stack.shape[-1]
        selected_indices = [0, n_images//2, n_images-1]
    
    # Create figure with gridspec layout
    fig = plt.figure(figsize=figsize)
    gs = gridspec.GridSpec(2, len(selected_indices), height_ratios=[5, 0.3], hspace=0.01)
    
    # Calculate dynamic range from selected subset
    sub_selection = hyperspectral_stack[..., selected_indices]
    vmin = np.nanpercentile(sub_selection, percentile_range[0])
    vmax = np.nanpercentile(sub_selection, percentile_range[1])
    
    # Create subplots for images in the top row
    img = None  # Will store the last image for colorbar
    for i, idx in enumerate(selected_indices):
        ax = fig.add_subplot(gs[0, i])
        ax.set_title(f"Radiograph {idx}")
        img = ax.imshow(hyperspectral_stack[..., idx], cmap=cmap, origin="lower", 
                       vmin=vmin, vmax=vmax)
        ax.axis("off")
    
    # Create colorbar in the bottom row, spanning all columns
    cbar_ax = fig.add_subplot(gs[1, :])
    cbar = plt.colorbar(img, cax=cbar_ax, orientation='horizontal')
    cbar.set_label('Intensity', fontsize=12)

    return fig

_ = visualize_hyperspectral_radiographs(hyperspectral_stack)

file = open("my_file.png", "rb")
image = file.read()
img = widgets.Image(value=image, format='png')
display(img)

display(widgets.HTML("<h1>Instrument Information</h1> <hr style=height:5px; border:none; color:#333;background-color:#333'>"))

$ lsof -ti:<port number.
20584
$ kill -9 20584