API Reference

Calibration

FilmCalibration

Source code in app/calibration/calibration.py

class FilmCalibration:

    def __init__(self, groundtruth_image: np.ndarray, bits_per_channel=8, 
                 calibration_type: str = 'single-channel', fitting_function_name: str = 'polynomial',
                 filter_type='median'):
        """
        Initializes the film calibration process by setting the ground truth image 
        and the calibration type.

        Parameters
        ----------
        groundtruth_image : np.ndarray
            The original ground truth image where the irradiated films are located.
            It must be a NumPy array with dimensions [height, width, channels].
        bits_per_channel : int, optional
            Number of bits per channel (default is 8).
        calibration_type : str, optional
            Calibration type ('single-channel' or 'multi-channel'). Default is 'single-channel'.
        fitting_function_name : str, optional
            Name of the fitting function to use. Default is 'polynomial'.
        filter_type : str, optional
            Type of filter to apply to images. Default is 'median'.
        """
        if not isinstance(groundtruth_image, np.ndarray):
            if isinstance(groundtruth_image, str):
                # If a string is provided, attempt to read the image from the file
                groundtruth_image = read_image(groundtruth_image)
            else:
                raise ValueError("groundtruth_image must be a NumPy array.")
        # Store the ground truth image and basic calibration parameters
        self.groundtruth_image = groundtruth_image
        self.bits_per_channel = bits_per_channel
        self.calibration_type = calibration_type
        self.filter_type = filter_type
        # Dictionary to store CalibrationDose instances mapped by dose value
        self.doses = {}
        # List to store the pre-irradiation pixel values for each channel
        self.pixel_values_before = [None, None, None]
        # List to store the mapping between dose and independent variable per channel
        self.dose_to_independent_by_channel = []  
        self.parameters = None
        self.uncertainties = None  # Will store the standard deviations of the fitting parameters
        self.fitting_func_name = fitting_function_name
        # Retrieve the fitting function instance based on the provided function name
        self.fitting_func_instance = get_fitting_function(fitting_function_name)

    def add_roi(self, dose: float, x: int, y: int, width: int, height: int):
        """
        Registers a region of interest (ROI) associated with a given dose. If the dose 
        does not exist, a new CalibrationDose instance is created.

        Parameters
        ----------
        dose : float
            Dose value (in Gy).
        x : int
            X-coordinate (top-left corner) of the ROI.
        y : int
            Y-coordinate (top-left corner) of the ROI.
        size : int
            Length of the square defining the ROI.
        """
        # If the dose is not already registered, create a new CalibrationDose instance
        if dose not in self.doses:
            self.doses[dose] = CalibrationDose(dose, calibration=self)
        # Add the ROI to the corresponding dose
        self.doses[dose].add_roi(x, y, width, height)

    def get_total_roi_count(self) -> int:
        """
        Returns the total number of ROIs registered across all doses.

        Returns
        -------
        int
            Total number of ROIs.
        """
        # Sum the number of ROIs for each dose
        return sum(calib_dose.get_roi_count() for calib_dose in self.doses.values())

    def get_total_dose_count(self) -> int:
        """
        Returns the total number of doses registered in the calibration.

        Returns
        -------
        int
            Total number of doses.
        """
        # The count of doses is equal to the number of keys in the doses dictionary
        return len(self.doses)

    def get_rois_by_dose(self):
        """
        Prints out the ROIs for each registered dose.
        """
        # Iterate over each dose and its corresponding CalibrationDose object to print the ROIs
        for dose, calib_dose in self.doses.items():
            print(f"Dose: {calib_dose.value} Gy")
            print(f"ROIs: {calib_dose.rois} \n")


    def compute_channel_independent_values(self, channel: int = 0) -> dict:
        """
        Computes the independent variable (e.g., netOD or netT) for each dose using 
        the global pre-irradiation pixel value (PV_before).

        Parameters
        ----------
        channel : int, optional
            The channel index for which to compute the independent variable (default is 0).

        Returns
        -------
        dict
            A dictionary mapping dose values to the computed independent variable.
        """
        # Ensure that dose 0 is defined to set the pre-irradiation pixel value (PV_before)
        if 0 not in self.doses:
            raise ValueError("Dose 0 must be defined in order to obtain pixel_values_before.")

        # Compute the average pixel value for all doses for the specified channel
        for dose, calib_dose in self.doses.items():
            calib_dose.compute_average_pv(channel)

        # Set the global PV_before using dose 0
        self.pixel_values_before[channel] = self.doses[0].pixel_values_after.get(channel, None)
        if self.pixel_values_before[channel] is None:
            raise ValueError(f"PV_before (from dose 0) for channel {channel} is not defined.")

        dose_to_value = {}
        # Compute the independent variable (e.g., netOD or netT) for each dose
        for dose, calib_dose in self.doses.items():
            dose_to_value[dose] = calib_dose.compute_independent_value(self.pixel_values_before, channel)

        # Return the mapping sorted by dose value
        return dict(sorted(dose_to_value.items()))


    def calibrate(self, calibration_type: str = 'single-channel'):
        """
        Calibrates the film based on the specified calibration type.
        Currently, only single-channel calibration is implemented, which includes:

          1. Computing channel averages.
          2. Computing netOD (or another independent variable) for each dose.
          3. Fitting the calibration function using curve_fit.

        The fitting uses the independent variable as the independent value and the dose as the dependent value.

        Parameters
        ----------
        calibration_type : str, optional
            Calibration type ('single-channel' or 'multi-channel'). Default is 'single-channel'.

        Returns
        -------
        list
            A list of optimal parameters for each channel obtained from curve_fit.

        Raises
        ------
        ValueError
            If the pre-irradiation pixel value (PV_before) is not defined.
        """
        # Only single-channel calibration is currently implemented
        if calibration_type != 'single-channel':
            raise NotImplementedError("Only single-channel calibration is implemented.")

        # Ensure that dose 0 exists for setting PV_before
        if 0 not in self.doses:
            raise ValueError("Dose 0 must be defined to set pixel_values_before.")

        dose_to_independent_by_channel = []
        parameters = []
        uncertainties = []

        # Iterate over the three channels of the image
        for channel in range(0, 3):
            # Calculate the independent variable mapping for the given channel
            dose_to_independent = self.compute_channel_independent_values(channel)
            dose_to_independent_by_channel.append(dose_to_independent)

            # Prepare lists for curve fitting: dose values and corresponding independent variable values
            dose_list = np.array(list(dose_to_independent.keys()))
            independent_list = np.array(list(dose_to_independent.values()))

            # Get initial parameter guess for the fitting function based on expected number of parameters
            num_params = len(self.fitting_func_instance.param_names)
            p0 = self.fitting_func_instance.initial_param_guess

            # Perform the curve fitting using the fitting function and the data
            popt, pcov = curve_fit(self.fitting_func_instance.func, independent_list, dose_list,
                                   p0=p0, maxfev=10000)
            parameters.append(popt)
            # Calculate uncertainties as the square root of the diagonal elements of the covariance matrix
            uncertainties.append(np.sqrt(np.diag(pcov)))

        # Store the computed mappings, parameters, and uncertainties for later use
        self.dose_to_independent_by_channel = dose_to_independent_by_channel
        self.parameters = parameters
        self.uncertainties = uncertainties

        return parameters

    def get_metric(self, metric_name: str, channel: int = None):
        """
        Computes a specified metric for the calibration curve(s) and returns a tuple 
        containing the metric value and its formatted name.

        Parameters
        ----------
        metric_name : str
            The metric to compute. Options are:
            - 'r2' or 'r^2' for the coefficient of determination,
            - 'rmse' for the root mean square error,
            - 'mse' for the mean square error,
            - 'chi2', 'chi-squared', or 'chi^2' for the chi-squared value.
        channel : int, optional
            The channel index for which to compute the metric. If None, returns a
            dictionary mapping each channel index to a (value, formatted name) tuple.

        Returns
        -------
        tuple or dict
            A tuple (metric_value, formatted_metric_name) if a channel is specified,
            otherwise a dictionary mapping channel indices to such tuples.
        """
        # Ensure that calibration has been performed and the necessary parameters exist
        if self.parameters is None or not self.dose_to_independent_by_channel:
            raise ValueError("Calibration parameters have not been computed. Please run calibrate() first.")

        def compute_for_channel(ch):
            # Get dose values and the corresponding independent variable values for the channel
            dose_to_x = self.dose_to_independent_by_channel[ch]
            dose_list = np.array(list(dose_to_x.keys()))
            x_list = np.array(list(dose_to_x.values()))
            popt = self.parameters[ch]
            # Compute the predicted dose using the fitting function with the optimized parameters
            dose_predicted = self.fitting_func_instance.func(x_list, *popt)

            metric = metric_name.lower()
            # Calculate the requested metric
            if metric in ['r2', 'r^2']:
                value = r2_score(dose_list, dose_predicted)
                formatted = "R²"
            elif metric in ['rmse', 'RMSE']:
                value = root_mean_squared_error(dose_list, dose_predicted)
                formatted = "RMSE"
            elif metric in ['mse', 'MSE']:
                value = mean_squared_error(dose_list, dose_predicted)
                formatted = "MSE"
            elif metric in ['chi2', 'chi-squared', 'chi squared', 'chi^2']:
                value = np.sum((dose_list - dose_predicted) ** 2)
                formatted = "χ²"
            else:
                raise ValueError("Metric not recognized. Available metrics: 'r2', 'rmse', 'mse', and 'chi2'.")
            return value, formatted

        if channel is not None:
            # Return the metric for the specified channel
            return compute_for_channel(channel)
        else:
            # Otherwise, compute the metric for each channel and return in a dictionary
            metrics = {}
            for ch in range(len(self.parameters)):
                metrics[ch] = compute_for_channel(ch)
            return metrics

    def graph_calibration_curve(self, metric_name='r2'):
        """
        Graphs the calibration curves for each channel.
        The x-axis represents the independent variable (e.g., netOD or netT) and
        the y-axis represents the dose.

        Parameters
        ----------
        metric_name : str, optional
            The metric to display on the curve legend (default is 'r2').
        """
        colors = ['r', 'g', 'b']
        plt.figure(figsize=(8, 6))

        # Iterate over each channel to plot experimental data points and the fitted curve
        for i, dose_to_x in enumerate(self.dose_to_independent_by_channel):

            dose_list = np.array(list(dose_to_x.keys()))
            x_list = np.array(list(dose_to_x.values()))

            popt = self.parameters[i]
            uncert = self.uncertainties[i]

            # Define a smooth range for the independent variable for curve plotting
            x_fit = np.linspace(min(x_list), max(x_list) + 0.012, 300)
            dose_fit = self.fitting_func_instance.func(x_fit, *popt)

            # Predict the dose values at the measured independent variable values
            dose_predicted = self.fitting_func_instance.func(x_list, *popt)

            # Compute the specified metric for the current channel
            metric_value, metric_text = self.get_metric(metric_name, channel=i)

            # Create a label displaying optimized parameters with their uncertainties and the metric
            label_text = "\n".join(
                f"{name}={p:.3f}±{u:.3f}" 
                for name, p, u in zip(self.fitting_func_instance.param_names, popt, uncert)
            )
            label_text += f"\n{metric_text}={metric_value:.3f}"

            # Plot experimental data points and the fitted curve
            plt.scatter(x_list, dose_list, color=colors[i])
            plt.plot(x_fit, dose_fit, color=colors[i], linestyle='--', label=label_text)

        plt.title(f"Calibration Curves {self.fitting_func_instance.description}", fontsize=14)
        plt.ylabel("Dose (Gy)", fontsize=12)
        plt.xlabel(self.fitting_func_instance.independent_variable, fontsize=12)
        # Place the legend outside the plot area for clarity
        plt.legend(bbox_to_anchor=(1.04, 0.5), loc='center left')
        plt.grid(True)
        #plt.show()
        return plt.gcf()  

    def graph_response_curve(self, metric_name='r2'):
        """
        Graphs the response curves for each channel.
        The x-axis corresponds to the dose (Gy) and the y-axis corresponds to the film response 
        (e.g., netOD or netT), which is the independent variable.

        For each channel, this method numerically inverts the calibration function (which 
        returns the dose for a given response) to obtain the independent variable for a range of doses.

        Parameters
        ----------
        metric_name : str, optional
            The metric to display on the curve legend (default is 'r2').
        """
        colors = ['r', 'g', 'b']
        plt.figure(figsize=(8, 6))

        for i, dose_to_x in enumerate(self.dose_to_independent_by_channel):

            dose_list = np.array(list(dose_to_x.keys()))
            x_list = np.array(list(dose_to_x.values()))

            popt = self.parameters[i]
            uncert = self.uncertainties[i]

            # Define a smooth range for the independent variable inversion
            x_fit = np.linspace(min(x_list), max(x_list) + 0.015, 300)
            dose_fit = self.fitting_func_instance.func(x_fit, *popt)

            # Predict the dose values at the measured independent variable values
            dose_predicted = self.fitting_func_instance.func(x_list, *popt)

            # Obtain the specified metric for the current channel to include in the legend
            metric_value, metric_text = self.get_metric(metric_name, channel=i)

            # Create a label displaying optimized parameters with uncertainties and the metric
            label_text = "\n".join(
                f"{name}={p:.3f}±{u:.3f}"
                for name, p, u in zip(self.fitting_func_instance.param_names, popt, uncert)
            )
            label_text += f"\n{metric_text}={metric_value:.3f}"

            plt.scatter(dose_list, x_list, color=colors[i])
            plt.plot(dose_fit, x_fit, color=colors[i], linestyle='--', label=label_text)

        plt.title(f"Response Curves {self.fitting_func_instance.description}", fontsize=14)
        plt.xlabel("Dose (Gy)", fontsize=12)
        plt.ylabel(self.fitting_func_instance.independent_variable, fontsize=12)
        # Place the legend outside the plot area for clarity
        plt.legend(bbox_to_anchor=(1.04, 0.5), loc='center left')
        plt.grid(True)
        plt.show()

    def to_json(self, filename: str):
        """
        Exports the FilmCalibration instance to a JSON file.
        Only the following attributes are saved:
            - groundtruth_image (converted to list)
            - bits_per_channel
            - calibration_type
            - filter_type
            - pixel_values_before
            - dose_to_independent_by_channel
            - parameters (converted to lists)
            - uncertainties (converted to lists)
            - fitting_func_name
        Note: The doses, pre-irradiation pixel values (pixel_values_before), and the fitting function
              instance are not saved.

        Parameters
        ----------
        filename : str
            The path to the JSON file where the instance will be saved.
        """
        # Prepare a dictionary with the data to export, converting NumPy arrays to lists
        data = {
            "groundtruth_image": self.groundtruth_image.tolist(),
            "bits_per_channel": self.bits_per_channel,
            "calibration_type": self.calibration_type,
            "filter_type": self.filter_type,
            "pixel_values_before": [None if x is None else float(x) for x in self.pixel_values_before],
            "dose_to_independent_by_channel": self.dose_to_independent_by_channel,
            "parameters": [p.tolist() if isinstance(p, (np.ndarray, list)) else p for p in self.parameters] if self.parameters is not None else None,
            "uncertainties": [u.tolist() if isinstance(u, (np.ndarray, list)) else u for u in self.uncertainties] if self.uncertainties is not None else None,
            "fitting_func_name": self.fitting_func_name
        }
        # Write the dictionary to the specified JSON file
        with open(filename, 'w') as f:
            json.dump(data, f)

    @classmethod
    def from_json(cls, filename: str):
        """
        Loads a FilmCalibration instance from a JSON file.
        The JSON must contain the keys:
            - groundtruth_image
            - bits_per_channel
            - calibration_type
            - filter_type
            - dose_to_independent_by_channel
            - parameters
            - uncertainties
            - fitting_func_name
        Note: The doses, pre-irradiation pixel values (pixel_values_before), and the fitting function
              instance are not stored; the fitting function instance is re-initialized using the fitting_func_name.

        Parameters
        ----------
        filename : str
            The path to the JSON file to load.

        Returns
        -------
        FilmCalibration
            The reconstructed FilmCalibration instance.
        """
        # Open and read the JSON file to load the data
        with open(filename, 'r') as f:
            data = json.load(f)

        # Reconstruct the ground truth image from the stored list
        groundtruth_image = np.array(data["groundtruth_image"])
        # Create an instance using the basic stored parameters
        instance = cls(
            groundtruth_image=groundtruth_image,
            bits_per_channel=data.get("bits_per_channel"),
            calibration_type=data.get("calibration_type"),
            fitting_function_name=data.get("fitting_func_name"),
            filter_type=data.get("filter_type")
        )

        # Load the global pre-irradiation pixel values, converting each to np.float64 if necessary
        instance.pixel_values_before = [
            None if x is None else np.float64(x) for x in data.get("pixel_values_before", [None, None, None])
        ]

        # Reconstruct the mapping of dose to independent variable for each channel
        dose_to_independent_by_channel = data.get("dose_to_independent_by_channel", [])
        instance.dose_to_independent_by_channel = [
            {float(k): np.float64(v) for k, v in d.items()} for d in dose_to_independent_by_channel
        ]

        # Reconstruct the fitting parameters if available
        parameters = data.get("parameters")
        if parameters is not None:
            instance.parameters = [np.array(p) for p in parameters]
        else:
            instance.parameters = None

        # Reconstruct the uncertainties if available
        uncertainties = data.get("uncertainties")
        if uncertainties is not None:
            instance.uncertainties = [np.array(u) for u in uncertainties]
        else:
            instance.uncertainties = None

        return instance

    def compute_dose_map(self, film_file: str, channel: int = 0, new_size=512) -> np.ndarray:
        """
        Loads an irradiated film from a .tif file and computes the dose map using the calibrated model
        and the single-channel method.

        Parameters
        ----------
        film_file : str
            Path to the .tif file that contains the irradiated film.
        channel : int, optional
            The channel of the image to be used for calculating the dose map (default is 0).
        new_size : int, optional
            A new size for resizing the image (not used in the current implementation).

        Returns
        -------
        np.ndarray
            A 2D array representing the calculated dose map for the film.

        Raises
        ------
        ValueError
            If the global pre-irradiation pixel value (PV_before) for the specified channel is not defined
            or if the independent variable used is not supported.
        """
        # Load the film image from the file
        film_image = read_image(film_file)

        # Optionally resize the image using OpenCV (commented out to preserve original logic)
        # film_image = cv2.resize(film_image, (new_size, new_size), interpolation=cv2.INTER_NEAREST)

        # If the image has multiple channels, extract the specified channel; otherwise use the image itself
        if film_image.ndim == 3:
            film_channel = film_image[:, :, channel]
        else:
            film_channel = film_image

        # Apply the specified filter to reduce noise if a filter type is provided
        if self.filter_type is not None:
            film_channel = filter_image(film_channel, self.filter_type)

        # Retrieve the global pre-irradiation pixel value for the specified channel (set during calibration)
        PV_before = self.pixel_values_before[channel]
        if PV_before is None:
            raise ValueError(f"The global PV_before for channel {channel} is not defined. Please ensure dose 0 is calibrated.")

        # Calculate the independent variable (netOD or netT) for each pixel
        independent_variable = self.fitting_func_instance.independent_variable
        if independent_variable == "netOD":
            # Replace zeros to avoid division by zero
            film_channel_safe = np.where(film_channel == 0, 1e-6, film_channel)
            x_map = np.log10(PV_before / film_channel_safe)
        elif independent_variable == "netT":
            x_map = (film_channel - PV_before)
        else:
            raise ValueError(f"Unsupported independent variable type: {independent_variable}")

        # Retrieve the calibrated parameters for the specified channel
        popt = self.parameters[channel]

        # Apply the calibration function using the optimized parameters to obtain the dose map
        dose_map = self.fitting_func_instance.func(x_map, *popt)

        # Replace any NaN values with 0
        dose_map = np.nan_to_num(dose_map)

        return dose_map

__init__(groundtruth_image, bits_per_channel=8, calibration_type='single-channel', fitting_function_name='polynomial', filter_type='median')

Initializes the film calibration process by setting the ground truth image and the calibration type.

Parameters:

Name	Type	Description	Default
groundtruth_image	ndarray	The original ground truth image where the irradiated films are located. It must be a NumPy array with dimensions [height, width, channels].	required
bits_per_channel	int	Number of bits per channel (default is 8).	8
calibration_type	str	Calibration type ('single-channel' or 'multi-channel'). Default is 'single-channel'.	'single-channel'
fitting_function_name	str	Name of the fitting function to use. Default is 'polynomial'.	'polynomial'
filter_type	str	Type of filter to apply to images. Default is 'median'.	'median'

Source code in app/calibration/calibration.py

def __init__(self, groundtruth_image: np.ndarray, bits_per_channel=8, 
             calibration_type: str = 'single-channel', fitting_function_name: str = 'polynomial',
             filter_type='median'):
    """
    Initializes the film calibration process by setting the ground truth image 
    and the calibration type.

    Parameters
    ----------
    groundtruth_image : np.ndarray
        The original ground truth image where the irradiated films are located.
        It must be a NumPy array with dimensions [height, width, channels].
    bits_per_channel : int, optional
        Number of bits per channel (default is 8).
    calibration_type : str, optional
        Calibration type ('single-channel' or 'multi-channel'). Default is 'single-channel'.
    fitting_function_name : str, optional
        Name of the fitting function to use. Default is 'polynomial'.
    filter_type : str, optional
        Type of filter to apply to images. Default is 'median'.
    """
    if not isinstance(groundtruth_image, np.ndarray):
        if isinstance(groundtruth_image, str):
            # If a string is provided, attempt to read the image from the file
            groundtruth_image = read_image(groundtruth_image)
        else:
            raise ValueError("groundtruth_image must be a NumPy array.")
    # Store the ground truth image and basic calibration parameters
    self.groundtruth_image = groundtruth_image
    self.bits_per_channel = bits_per_channel
    self.calibration_type = calibration_type
    self.filter_type = filter_type
    # Dictionary to store CalibrationDose instances mapped by dose value
    self.doses = {}
    # List to store the pre-irradiation pixel values for each channel
    self.pixel_values_before = [None, None, None]
    # List to store the mapping between dose and independent variable per channel
    self.dose_to_independent_by_channel = []  
    self.parameters = None
    self.uncertainties = None  # Will store the standard deviations of the fitting parameters
    self.fitting_func_name = fitting_function_name
    # Retrieve the fitting function instance based on the provided function name
    self.fitting_func_instance = get_fitting_function(fitting_function_name)

add_roi(dose, x, y, width, height)

Registers a region of interest (ROI) associated with a given dose. If the dose does not exist, a new CalibrationDose instance is created.

Parameters:

Name	Type	Description	Default
dose	float	Dose value (in Gy).	required
x	int	X-coordinate (top-left corner) of the ROI.	required
y	int	Y-coordinate (top-left corner) of the ROI.	required
size	int	Length of the square defining the ROI.	required

Source code in app/calibration/calibration.py

def add_roi(self, dose: float, x: int, y: int, width: int, height: int):
    """
    Registers a region of interest (ROI) associated with a given dose. If the dose 
    does not exist, a new CalibrationDose instance is created.

    Parameters
    ----------
    dose : float
        Dose value (in Gy).
    x : int
        X-coordinate (top-left corner) of the ROI.
    y : int
        Y-coordinate (top-left corner) of the ROI.
    size : int
        Length of the square defining the ROI.
    """
    # If the dose is not already registered, create a new CalibrationDose instance
    if dose not in self.doses:
        self.doses[dose] = CalibrationDose(dose, calibration=self)
    # Add the ROI to the corresponding dose
    self.doses[dose].add_roi(x, y, width, height)

calibrate(calibration_type='single-channel')

Calibrates the film based on the specified calibration type. Currently, only single-channel calibration is implemented, which includes:

1. Computing channel averages.
2. Computing netOD (or another independent variable) for each dose.
3. Fitting the calibration function using curve_fit.

The fitting uses the independent variable as the independent value and the dose as the dependent value.

Parameters:

Name	Type	Description	Default
calibration_type	str	Calibration type ('single-channel' or 'multi-channel'). Default is 'single-channel'.	'single-channel'

Returns:

Type	Description
list	A list of optimal parameters for each channel obtained from curve_fit.

Raises:

Type	Description
ValueError	If the pre-irradiation pixel value (PV_before) is not defined.

Source code in app/calibration/calibration.py

def calibrate(self, calibration_type: str = 'single-channel'):
    """
    Calibrates the film based on the specified calibration type.
    Currently, only single-channel calibration is implemented, which includes:

      1. Computing channel averages.
      2. Computing netOD (or another independent variable) for each dose.
      3. Fitting the calibration function using curve_fit.

    The fitting uses the independent variable as the independent value and the dose as the dependent value.

    Parameters
    ----------
    calibration_type : str, optional
        Calibration type ('single-channel' or 'multi-channel'). Default is 'single-channel'.

    Returns
    -------
    list
        A list of optimal parameters for each channel obtained from curve_fit.

    Raises
    ------
    ValueError
        If the pre-irradiation pixel value (PV_before) is not defined.
    """
    # Only single-channel calibration is currently implemented
    if calibration_type != 'single-channel':
        raise NotImplementedError("Only single-channel calibration is implemented.")

    # Ensure that dose 0 exists for setting PV_before
    if 0 not in self.doses:
        raise ValueError("Dose 0 must be defined to set pixel_values_before.")

    dose_to_independent_by_channel = []
    parameters = []
    uncertainties = []

    # Iterate over the three channels of the image
    for channel in range(0, 3):
        # Calculate the independent variable mapping for the given channel
        dose_to_independent = self.compute_channel_independent_values(channel)
        dose_to_independent_by_channel.append(dose_to_independent)

        # Prepare lists for curve fitting: dose values and corresponding independent variable values
        dose_list = np.array(list(dose_to_independent.keys()))
        independent_list = np.array(list(dose_to_independent.values()))

        # Get initial parameter guess for the fitting function based on expected number of parameters
        num_params = len(self.fitting_func_instance.param_names)
        p0 = self.fitting_func_instance.initial_param_guess

        # Perform the curve fitting using the fitting function and the data
        popt, pcov = curve_fit(self.fitting_func_instance.func, independent_list, dose_list,
                               p0=p0, maxfev=10000)
        parameters.append(popt)
        # Calculate uncertainties as the square root of the diagonal elements of the covariance matrix
        uncertainties.append(np.sqrt(np.diag(pcov)))

    # Store the computed mappings, parameters, and uncertainties for later use
    self.dose_to_independent_by_channel = dose_to_independent_by_channel
    self.parameters = parameters
    self.uncertainties = uncertainties

    return parameters

compute_channel_independent_values(channel=0)

Computes the independent variable (e.g., netOD or netT) for each dose using the global pre-irradiation pixel value (PV_before).

Parameters:

Name	Type	Description	Default
channel	int	The channel index for which to compute the independent variable (default is 0).	0

Returns:

Type	Description
dict	A dictionary mapping dose values to the computed independent variable.

Source code in app/calibration/calibration.py

def compute_channel_independent_values(self, channel: int = 0) -> dict:
    """
    Computes the independent variable (e.g., netOD or netT) for each dose using 
    the global pre-irradiation pixel value (PV_before).

    Parameters
    ----------
    channel : int, optional
        The channel index for which to compute the independent variable (default is 0).

    Returns
    -------
    dict
        A dictionary mapping dose values to the computed independent variable.
    """
    # Ensure that dose 0 is defined to set the pre-irradiation pixel value (PV_before)
    if 0 not in self.doses:
        raise ValueError("Dose 0 must be defined in order to obtain pixel_values_before.")

    # Compute the average pixel value for all doses for the specified channel
    for dose, calib_dose in self.doses.items():
        calib_dose.compute_average_pv(channel)

    # Set the global PV_before using dose 0
    self.pixel_values_before[channel] = self.doses[0].pixel_values_after.get(channel, None)
    if self.pixel_values_before[channel] is None:
        raise ValueError(f"PV_before (from dose 0) for channel {channel} is not defined.")

    dose_to_value = {}
    # Compute the independent variable (e.g., netOD or netT) for each dose
    for dose, calib_dose in self.doses.items():
        dose_to_value[dose] = calib_dose.compute_independent_value(self.pixel_values_before, channel)

    # Return the mapping sorted by dose value
    return dict(sorted(dose_to_value.items()))

compute_dose_map(film_file, channel=0, new_size=512)

Loads an irradiated film from a .tif file and computes the dose map using the calibrated model and the single-channel method.

Parameters:

Name	Type	Description	Default
film_file	str	Path to the .tif file that contains the irradiated film.	required
channel	int	The channel of the image to be used for calculating the dose map (default is 0).	0
new_size	int	A new size for resizing the image (not used in the current implementation).	512

Returns:

Type	Description
ndarray	A 2D array representing the calculated dose map for the film.

Raises:

Type	Description
ValueError	If the global pre-irradiation pixel value (PV_before) for the specified channel is not defined or if the independent variable used is not supported.

Source code in app/calibration/calibration.py

def compute_dose_map(self, film_file: str, channel: int = 0, new_size=512) -> np.ndarray:
    """
    Loads an irradiated film from a .tif file and computes the dose map using the calibrated model
    and the single-channel method.

    Parameters
    ----------
    film_file : str
        Path to the .tif file that contains the irradiated film.
    channel : int, optional
        The channel of the image to be used for calculating the dose map (default is 0).
    new_size : int, optional
        A new size for resizing the image (not used in the current implementation).

    Returns
    -------
    np.ndarray
        A 2D array representing the calculated dose map for the film.

    Raises
    ------
    ValueError
        If the global pre-irradiation pixel value (PV_before) for the specified channel is not defined
        or if the independent variable used is not supported.
    """
    # Load the film image from the file
    film_image = read_image(film_file)

    # Optionally resize the image using OpenCV (commented out to preserve original logic)
    # film_image = cv2.resize(film_image, (new_size, new_size), interpolation=cv2.INTER_NEAREST)

    # If the image has multiple channels, extract the specified channel; otherwise use the image itself
    if film_image.ndim == 3:
        film_channel = film_image[:, :, channel]
    else:
        film_channel = film_image

    # Apply the specified filter to reduce noise if a filter type is provided
    if self.filter_type is not None:
        film_channel = filter_image(film_channel, self.filter_type)

    # Retrieve the global pre-irradiation pixel value for the specified channel (set during calibration)
    PV_before = self.pixel_values_before[channel]
    if PV_before is None:
        raise ValueError(f"The global PV_before for channel {channel} is not defined. Please ensure dose 0 is calibrated.")

    # Calculate the independent variable (netOD or netT) for each pixel
    independent_variable = self.fitting_func_instance.independent_variable
    if independent_variable == "netOD":
        # Replace zeros to avoid division by zero
        film_channel_safe = np.where(film_channel == 0, 1e-6, film_channel)
        x_map = np.log10(PV_before / film_channel_safe)
    elif independent_variable == "netT":
        x_map = (film_channel - PV_before)
    else:
        raise ValueError(f"Unsupported independent variable type: {independent_variable}")

    # Retrieve the calibrated parameters for the specified channel
    popt = self.parameters[channel]

    # Apply the calibration function using the optimized parameters to obtain the dose map
    dose_map = self.fitting_func_instance.func(x_map, *popt)

    # Replace any NaN values with 0
    dose_map = np.nan_to_num(dose_map)

    return dose_map

from_json(filename) classmethod

Loads a FilmCalibration instance from a JSON file. The JSON must contain the keys: - groundtruth_image - bits_per_channel - calibration_type - filter_type - dose_to_independent_by_channel - parameters - uncertainties - fitting_func_name Note: The doses, pre-irradiation pixel values (pixel_values_before), and the fitting function instance are not stored; the fitting function instance is re-initialized using the fitting_func_name.

Parameters:

Name	Type	Description	Default
filename	str	The path to the JSON file to load.	required

Returns:

Type	Description
FilmCalibration	The reconstructed FilmCalibration instance.

Source code in app/calibration/calibration.py

@classmethod
def from_json(cls, filename: str):
    """
    Loads a FilmCalibration instance from a JSON file.
    The JSON must contain the keys:
        - groundtruth_image
        - bits_per_channel
        - calibration_type
        - filter_type
        - dose_to_independent_by_channel
        - parameters
        - uncertainties
        - fitting_func_name
    Note: The doses, pre-irradiation pixel values (pixel_values_before), and the fitting function
          instance are not stored; the fitting function instance is re-initialized using the fitting_func_name.

    Parameters
    ----------
    filename : str
        The path to the JSON file to load.

    Returns
    -------
    FilmCalibration
        The reconstructed FilmCalibration instance.
    """
    # Open and read the JSON file to load the data
    with open(filename, 'r') as f:
        data = json.load(f)

    # Reconstruct the ground truth image from the stored list
    groundtruth_image = np.array(data["groundtruth_image"])
    # Create an instance using the basic stored parameters
    instance = cls(
        groundtruth_image=groundtruth_image,
        bits_per_channel=data.get("bits_per_channel"),
        calibration_type=data.get("calibration_type"),
        fitting_function_name=data.get("fitting_func_name"),
        filter_type=data.get("filter_type")
    )

    # Load the global pre-irradiation pixel values, converting each to np.float64 if necessary
    instance.pixel_values_before = [
        None if x is None else np.float64(x) for x in data.get("pixel_values_before", [None, None, None])
    ]

    # Reconstruct the mapping of dose to independent variable for each channel
    dose_to_independent_by_channel = data.get("dose_to_independent_by_channel", [])
    instance.dose_to_independent_by_channel = [
        {float(k): np.float64(v) for k, v in d.items()} for d in dose_to_independent_by_channel
    ]

    # Reconstruct the fitting parameters if available
    parameters = data.get("parameters")
    if parameters is not None:
        instance.parameters = [np.array(p) for p in parameters]
    else:
        instance.parameters = None

    # Reconstruct the uncertainties if available
    uncertainties = data.get("uncertainties")
    if uncertainties is not None:
        instance.uncertainties = [np.array(u) for u in uncertainties]
    else:
        instance.uncertainties = None

    return instance

get_metric(metric_name, channel=None)

Computes a specified metric for the calibration curve(s) and returns a tuple containing the metric value and its formatted name.

Parameters:

Name	Type	Description	Default
metric_name	str	The metric to compute. Options are: - 'r2' or 'r^2' for the coefficient of determination, - 'rmse' for the root mean square error, - 'mse' for the mean square error, - 'chi2', 'chi-squared', or 'chi^2' for the chi-squared value.	required
channel	int	The channel index for which to compute the metric. If None, returns a dictionary mapping each channel index to a (value, formatted name) tuple.	None

Returns:

Type	Description
tuple or dict	A tuple (metric_value, formatted_metric_name) if a channel is specified, otherwise a dictionary mapping channel indices to such tuples.

Source code in app/calibration/calibration.py

def get_metric(self, metric_name: str, channel: int = None):
    """
    Computes a specified metric for the calibration curve(s) and returns a tuple 
    containing the metric value and its formatted name.

    Parameters
    ----------
    metric_name : str
        The metric to compute. Options are:
        - 'r2' or 'r^2' for the coefficient of determination,
        - 'rmse' for the root mean square error,
        - 'mse' for the mean square error,
        - 'chi2', 'chi-squared', or 'chi^2' for the chi-squared value.
    channel : int, optional
        The channel index for which to compute the metric. If None, returns a
        dictionary mapping each channel index to a (value, formatted name) tuple.

    Returns
    -------
    tuple or dict
        A tuple (metric_value, formatted_metric_name) if a channel is specified,
        otherwise a dictionary mapping channel indices to such tuples.
    """
    # Ensure that calibration has been performed and the necessary parameters exist
    if self.parameters is None or not self.dose_to_independent_by_channel:
        raise ValueError("Calibration parameters have not been computed. Please run calibrate() first.")

    def compute_for_channel(ch):
        # Get dose values and the corresponding independent variable values for the channel
        dose_to_x = self.dose_to_independent_by_channel[ch]
        dose_list = np.array(list(dose_to_x.keys()))
        x_list = np.array(list(dose_to_x.values()))
        popt = self.parameters[ch]
        # Compute the predicted dose using the fitting function with the optimized parameters
        dose_predicted = self.fitting_func_instance.func(x_list, *popt)

        metric = metric_name.lower()
        # Calculate the requested metric
        if metric in ['r2', 'r^2']:
            value = r2_score(dose_list, dose_predicted)
            formatted = "R²"
        elif metric in ['rmse', 'RMSE']:
            value = root_mean_squared_error(dose_list, dose_predicted)
            formatted = "RMSE"
        elif metric in ['mse', 'MSE']:
            value = mean_squared_error(dose_list, dose_predicted)
            formatted = "MSE"
        elif metric in ['chi2', 'chi-squared', 'chi squared', 'chi^2']:
            value = np.sum((dose_list - dose_predicted) ** 2)
            formatted = "χ²"
        else:
            raise ValueError("Metric not recognized. Available metrics: 'r2', 'rmse', 'mse', and 'chi2'.")
        return value, formatted

    if channel is not None:
        # Return the metric for the specified channel
        return compute_for_channel(channel)
    else:
        # Otherwise, compute the metric for each channel and return in a dictionary
        metrics = {}
        for ch in range(len(self.parameters)):
            metrics[ch] = compute_for_channel(ch)
        return metrics

get_rois_by_dose()

Prints out the ROIs for each registered dose.

Source code in app/calibration/calibration.py

def get_rois_by_dose(self):
    """
    Prints out the ROIs for each registered dose.
    """
    # Iterate over each dose and its corresponding CalibrationDose object to print the ROIs
    for dose, calib_dose in self.doses.items():
        print(f"Dose: {calib_dose.value} Gy")
        print(f"ROIs: {calib_dose.rois} \n")

get_total_dose_count()

Returns the total number of doses registered in the calibration.

Returns:

Type	Description
int	Total number of doses.

Source code in app/calibration/calibration.py

def get_total_dose_count(self) -> int:
    """
    Returns the total number of doses registered in the calibration.

    Returns
    -------
    int
        Total number of doses.
    """
    # The count of doses is equal to the number of keys in the doses dictionary
    return len(self.doses)

get_total_roi_count()

Returns the total number of ROIs registered across all doses.

Returns:

Type	Description
int	Total number of ROIs.

Source code in app/calibration/calibration.py

def get_total_roi_count(self) -> int:
    """
    Returns the total number of ROIs registered across all doses.

    Returns
    -------
    int
        Total number of ROIs.
    """
    # Sum the number of ROIs for each dose
    return sum(calib_dose.get_roi_count() for calib_dose in self.doses.values())

graph_calibration_curve(metric_name='r2')

Graphs the calibration curves for each channel. The x-axis represents the independent variable (e.g., netOD or netT) and the y-axis represents the dose.

Parameters:

Name	Type	Description	Default
metric_name	str	The metric to display on the curve legend (default is 'r2').	'r2'

Source code in app/calibration/calibration.py

def graph_calibration_curve(self, metric_name='r2'):
    """
    Graphs the calibration curves for each channel.
    The x-axis represents the independent variable (e.g., netOD or netT) and
    the y-axis represents the dose.

    Parameters
    ----------
    metric_name : str, optional
        The metric to display on the curve legend (default is 'r2').
    """
    colors = ['r', 'g', 'b']
    plt.figure(figsize=(8, 6))

    # Iterate over each channel to plot experimental data points and the fitted curve
    for i, dose_to_x in enumerate(self.dose_to_independent_by_channel):

        dose_list = np.array(list(dose_to_x.keys()))
        x_list = np.array(list(dose_to_x.values()))

        popt = self.parameters[i]
        uncert = self.uncertainties[i]

        # Define a smooth range for the independent variable for curve plotting
        x_fit = np.linspace(min(x_list), max(x_list) + 0.012, 300)
        dose_fit = self.fitting_func_instance.func(x_fit, *popt)

        # Predict the dose values at the measured independent variable values
        dose_predicted = self.fitting_func_instance.func(x_list, *popt)

        # Compute the specified metric for the current channel
        metric_value, metric_text = self.get_metric(metric_name, channel=i)

        # Create a label displaying optimized parameters with their uncertainties and the metric
        label_text = "\n".join(
            f"{name}={p:.3f}±{u:.3f}" 
            for name, p, u in zip(self.fitting_func_instance.param_names, popt, uncert)
        )
        label_text += f"\n{metric_text}={metric_value:.3f}"

        # Plot experimental data points and the fitted curve
        plt.scatter(x_list, dose_list, color=colors[i])
        plt.plot(x_fit, dose_fit, color=colors[i], linestyle='--', label=label_text)

    plt.title(f"Calibration Curves {self.fitting_func_instance.description}", fontsize=14)
    plt.ylabel("Dose (Gy)", fontsize=12)
    plt.xlabel(self.fitting_func_instance.independent_variable, fontsize=12)
    # Place the legend outside the plot area for clarity
    plt.legend(bbox_to_anchor=(1.04, 0.5), loc='center left')
    plt.grid(True)
    #plt.show()
    return plt.gcf()  

graph_response_curve(metric_name='r2')

Graphs the response curves for each channel. The x-axis corresponds to the dose (Gy) and the y-axis corresponds to the film response (e.g., netOD or netT), which is the independent variable.

For each channel, this method numerically inverts the calibration function (which returns the dose for a given response) to obtain the independent variable for a range of doses.

Parameters:

Name	Type	Description	Default
metric_name	str	The metric to display on the curve legend (default is 'r2').	'r2'

Source code in app/calibration/calibration.py

def graph_response_curve(self, metric_name='r2'):
    """
    Graphs the response curves for each channel.
    The x-axis corresponds to the dose (Gy) and the y-axis corresponds to the film response 
    (e.g., netOD or netT), which is the independent variable.

    For each channel, this method numerically inverts the calibration function (which 
    returns the dose for a given response) to obtain the independent variable for a range of doses.

    Parameters
    ----------
    metric_name : str, optional
        The metric to display on the curve legend (default is 'r2').
    """
    colors = ['r', 'g', 'b']
    plt.figure(figsize=(8, 6))

    for i, dose_to_x in enumerate(self.dose_to_independent_by_channel):

        dose_list = np.array(list(dose_to_x.keys()))
        x_list = np.array(list(dose_to_x.values()))

        popt = self.parameters[i]
        uncert = self.uncertainties[i]

        # Define a smooth range for the independent variable inversion
        x_fit = np.linspace(min(x_list), max(x_list) + 0.015, 300)
        dose_fit = self.fitting_func_instance.func(x_fit, *popt)

        # Predict the dose values at the measured independent variable values
        dose_predicted = self.fitting_func_instance.func(x_list, *popt)

        # Obtain the specified metric for the current channel to include in the legend
        metric_value, metric_text = self.get_metric(metric_name, channel=i)

        # Create a label displaying optimized parameters with uncertainties and the metric
        label_text = "\n".join(
            f"{name}={p:.3f}±{u:.3f}"
            for name, p, u in zip(self.fitting_func_instance.param_names, popt, uncert)
        )
        label_text += f"\n{metric_text}={metric_value:.3f}"

        plt.scatter(dose_list, x_list, color=colors[i])
        plt.plot(dose_fit, x_fit, color=colors[i], linestyle='--', label=label_text)

    plt.title(f"Response Curves {self.fitting_func_instance.description}", fontsize=14)
    plt.xlabel("Dose (Gy)", fontsize=12)
    plt.ylabel(self.fitting_func_instance.independent_variable, fontsize=12)
    # Place the legend outside the plot area for clarity
    plt.legend(bbox_to_anchor=(1.04, 0.5), loc='center left')
    plt.grid(True)
    plt.show()

to_json(filename)

Exports the FilmCalibration instance to a JSON file. Only the following attributes are saved: - groundtruth_image (converted to list) - bits_per_channel - calibration_type - filter_type - pixel_values_before - dose_to_independent_by_channel - parameters (converted to lists) - uncertainties (converted to lists) - fitting_func_name Note: The doses, pre-irradiation pixel values (pixel_values_before), and the fitting function instance are not saved.

Parameters:

Name	Type	Description	Default
filename	str	The path to the JSON file where the instance will be saved.	required

Source code in app/calibration/calibration.py

def to_json(self, filename: str):
    """
    Exports the FilmCalibration instance to a JSON file.
    Only the following attributes are saved:
        - groundtruth_image (converted to list)
        - bits_per_channel
        - calibration_type
        - filter_type
        - pixel_values_before
        - dose_to_independent_by_channel
        - parameters (converted to lists)
        - uncertainties (converted to lists)
        - fitting_func_name
    Note: The doses, pre-irradiation pixel values (pixel_values_before), and the fitting function
          instance are not saved.

    Parameters
    ----------
    filename : str
        The path to the JSON file where the instance will be saved.
    """
    # Prepare a dictionary with the data to export, converting NumPy arrays to lists
    data = {
        "groundtruth_image": self.groundtruth_image.tolist(),
        "bits_per_channel": self.bits_per_channel,
        "calibration_type": self.calibration_type,
        "filter_type": self.filter_type,
        "pixel_values_before": [None if x is None else float(x) for x in self.pixel_values_before],
        "dose_to_independent_by_channel": self.dose_to_independent_by_channel,
        "parameters": [p.tolist() if isinstance(p, (np.ndarray, list)) else p for p in self.parameters] if self.parameters is not None else None,
        "uncertainties": [u.tolist() if isinstance(u, (np.ndarray, list)) else u for u in self.uncertainties] if self.uncertainties is not None else None,
        "fitting_func_name": self.fitting_func_name
    }
    # Write the dictionary to the specified JSON file
    with open(filename, 'w') as f:
        json.dump(data, f)

CalibrationDose

Source code in app/calibration/dose.py

class CalibrationDose:
    def __init__(self, dose_value: float, calibration=None):
        """
        Initializes an instance for a specific dose.

        Parameters
        ----------
        dose_value : float
            Dose value (in Gy).
        calibration : FilmCalibration, optional
            The FilmCalibration instance that this dose is part of.
        """
        self.calibration = calibration
        self.value = dose_value
        self.rois = []  # List of ROIs associated with this dose; each ROI is defined as a tuple (x, y, size)
        # Dictionary to store the average pixel value after exposure (PV_after) per channel.
        self.pixel_values_after = {}  # {channel: PV_after}
        # Dictionary to store the computed independent variable (e.g., netOD or netT) per channel.
        self.independent_values = {}  # {channel: computed independent variable}

    def add_roi(self, x: int, y: int, width: int, height: int):
        """
        Adds a square region of interest (ROI) associated with this dose.

        Parameters
        ----------
        x : int
            X-coordinate (top-left corner) of the ROI.
        y : int
            Y-coordinate (top-left corner) of the ROI.
        size : int
            The side length of the square ROI.
        """
        self.rois.append((x, y, width, height))

    def compute_average_pv(self, channel: int = 0) -> float:
        """
        Computes the average post-exposure pixel value (PV_after) for all ROIs associated with this dose.
        It applies a median filter on the specified channel of the ground truth image.

        Parameters
        ----------
        channel : int, optional
            The image channel to process (default is 0).

        Returns
        -------
        float
            The computed average pixel value for the specified channel or None if no ROI exists.
        """
        roi_values = []
        for (x, y, w, h) in self.rois:
            # Extract the ROI from the ground truth image for the specified channel.
            # use width and height to extract the ROI
            roi = self.calibration.groundtruth_image[y:y+h, x:x+w, channel]
            # Apply the median filter to the ROI to reduce noise.
            filtered_roi = filter_image(roi, filter_type=self.calibration.filter_type)
            # Compute the average pixel value of the filtered ROI.
            average_value = np.mean(filtered_roi)
            roi_values.append(average_value)
        # If there are computed ROI values, calculate and store the overall average.
        if roi_values:
            self.pixel_values_after[channel] = np.mean(roi_values)
        else:
            # TODO: Consider raising an exception if no ROIs were added for this dose.
            pass
        return self.pixel_values_after.get(channel, None)

    def compute_independent_value(self, pixel_values_before, channel: int) -> float:
        """
        Computes the independent variable for this dose based on the pre-exposure pixel value (PV_before)
        and the average post-exposure pixel value (PV_after). The type of independent variable (e.g., netOD or netT)
        is determined by the configuration of the calibration function.

        Parameters
        ----------
        pixel_values_before : dict
            A dictionary mapping each channel to its pre-exposure pixel value.
        channel : int
            The image channel to use for the calculation.

        Returns
        -------
        float
            The computed independent variable value for this dose on the specified channel.

        Raises
        ------
        ValueError
            If the average post-exposure pixel value for the given channel has not been computed.
        """
        independent_variable = self.calibration.fitting_func_instance.independent_variable

        # Check that the average pixel value after exposure has been computed for the channel.
        if channel not in self.pixel_values_after or self.pixel_values_after[channel] is None:
            raise ValueError(f"The average pixel value after exposure for channel {channel} has not been computed yet.")

        if independent_variable == 'netOD':
            # Compute net optical density: netOD = log10(PV_before / PV_after)
            value = np.log10(pixel_values_before[channel] / self.pixel_values_after[channel])
        elif independent_variable == 'netT':
            bits_per_channel = self.calibration.bits_per_channel
            # Compute net transmission (netT) as the absolute difference between PV_after and PV_before.
            # Optionally, this difference may be normalized by 2^(bits_per_channel) if needed.
            value = np.abs(self.pixel_values_after[channel] - pixel_values_before[channel])
        else:
            raise ValueError(f"Unsupported independent variable type: {independent_variable}")

        self.independent_values[channel] = value
        return value

    def get_roi_count(self) -> int:
        """
        Returns the number of regions of interest (ROIs) registered for this dose.

        Returns
        -------
        int
            The number of ROIs.
        """
        return len(self.rois)

    def __repr__(self):
        return f"CalibrationDose(value={self.value}, num_rois={self.get_roi_count()})"

__init__(dose_value, calibration=None)

Initializes an instance for a specific dose.

Parameters:

Name	Type	Description	Default
dose_value	float	Dose value (in Gy).	required
calibration	FilmCalibration	The FilmCalibration instance that this dose is part of.	None

Source code in app/calibration/dose.py

def __init__(self, dose_value: float, calibration=None):
    """
    Initializes an instance for a specific dose.

    Parameters
    ----------
    dose_value : float
        Dose value (in Gy).
    calibration : FilmCalibration, optional
        The FilmCalibration instance that this dose is part of.
    """
    self.calibration = calibration
    self.value = dose_value
    self.rois = []  # List of ROIs associated with this dose; each ROI is defined as a tuple (x, y, size)
    # Dictionary to store the average pixel value after exposure (PV_after) per channel.
    self.pixel_values_after = {}  # {channel: PV_after}
    # Dictionary to store the computed independent variable (e.g., netOD or netT) per channel.
    self.independent_values = {}  # {channel: computed independent variable}

add_roi(x, y, width, height)

Adds a square region of interest (ROI) associated with this dose.

Parameters:

Name	Type	Description	Default
x	int	X-coordinate (top-left corner) of the ROI.	required
y	int	Y-coordinate (top-left corner) of the ROI.	required
size	int	The side length of the square ROI.	required

Source code in app/calibration/dose.py

def add_roi(self, x: int, y: int, width: int, height: int):
    """
    Adds a square region of interest (ROI) associated with this dose.

    Parameters
    ----------
    x : int
        X-coordinate (top-left corner) of the ROI.
    y : int
        Y-coordinate (top-left corner) of the ROI.
    size : int
        The side length of the square ROI.
    """
    self.rois.append((x, y, width, height))

compute_average_pv(channel=0)

Computes the average post-exposure pixel value (PV_after) for all ROIs associated with this dose. It applies a median filter on the specified channel of the ground truth image.

Parameters:

Name	Type	Description	Default
channel	int	The image channel to process (default is 0).	0

Returns:

Type	Description
float	The computed average pixel value for the specified channel or None if no ROI exists.

Source code in app/calibration/dose.py

def compute_average_pv(self, channel: int = 0) -> float:
    """
    Computes the average post-exposure pixel value (PV_after) for all ROIs associated with this dose.
    It applies a median filter on the specified channel of the ground truth image.

    Parameters
    ----------
    channel : int, optional
        The image channel to process (default is 0).

    Returns
    -------
    float
        The computed average pixel value for the specified channel or None if no ROI exists.
    """
    roi_values = []
    for (x, y, w, h) in self.rois:
        # Extract the ROI from the ground truth image for the specified channel.
        # use width and height to extract the ROI
        roi = self.calibration.groundtruth_image[y:y+h, x:x+w, channel]
        # Apply the median filter to the ROI to reduce noise.
        filtered_roi = filter_image(roi, filter_type=self.calibration.filter_type)
        # Compute the average pixel value of the filtered ROI.
        average_value = np.mean(filtered_roi)
        roi_values.append(average_value)
    # If there are computed ROI values, calculate and store the overall average.
    if roi_values:
        self.pixel_values_after[channel] = np.mean(roi_values)
    else:
        # TODO: Consider raising an exception if no ROIs were added for this dose.
        pass
    return self.pixel_values_after.get(channel, None)

compute_independent_value(pixel_values_before, channel)

Computes the independent variable for this dose based on the pre-exposure pixel value (PV_before) and the average post-exposure pixel value (PV_after). The type of independent variable (e.g., netOD or netT) is determined by the configuration of the calibration function.

Parameters:

Name	Type	Description	Default
pixel_values_before	dict	A dictionary mapping each channel to its pre-exposure pixel value.	required
channel	int	The image channel to use for the calculation.	required

Returns:

Type	Description
float	The computed independent variable value for this dose on the specified channel.

Raises:

Type	Description
ValueError	If the average post-exposure pixel value for the given channel has not been computed.

Source code in app/calibration/dose.py

def compute_independent_value(self, pixel_values_before, channel: int) -> float:
    """
    Computes the independent variable for this dose based on the pre-exposure pixel value (PV_before)
    and the average post-exposure pixel value (PV_after). The type of independent variable (e.g., netOD or netT)
    is determined by the configuration of the calibration function.

    Parameters
    ----------
    pixel_values_before : dict
        A dictionary mapping each channel to its pre-exposure pixel value.
    channel : int
        The image channel to use for the calculation.

    Returns
    -------
    float
        The computed independent variable value for this dose on the specified channel.

    Raises
    ------
    ValueError
        If the average post-exposure pixel value for the given channel has not been computed.
    """
    independent_variable = self.calibration.fitting_func_instance.independent_variable

    # Check that the average pixel value after exposure has been computed for the channel.
    if channel not in self.pixel_values_after or self.pixel_values_after[channel] is None:
        raise ValueError(f"The average pixel value after exposure for channel {channel} has not been computed yet.")

    if independent_variable == 'netOD':
        # Compute net optical density: netOD = log10(PV_before / PV_after)
        value = np.log10(pixel_values_before[channel] / self.pixel_values_after[channel])
    elif independent_variable == 'netT':
        bits_per_channel = self.calibration.bits_per_channel
        # Compute net transmission (netT) as the absolute difference between PV_after and PV_before.
        # Optionally, this difference may be normalized by 2^(bits_per_channel) if needed.
        value = np.abs(self.pixel_values_after[channel] - pixel_values_before[channel])
    else:
        raise ValueError(f"Unsupported independent variable type: {independent_variable}")

    self.independent_values[channel] = value
    return value

get_roi_count()

Returns the number of regions of interest (ROIs) registered for this dose.

Returns:

Type	Description
int	The number of ROIs.

Source code in app/calibration/dose.py

def get_roi_count(self) -> int:
    """
    Returns the number of regions of interest (ROIs) registered for this dose.

    Returns
    -------
    int
        The number of ROIs.
    """
    return len(self.rois)

functions.py

This module defines a FittingFunction class that encapsulates a fitting function, its name, descriptive text, the names of its fitting parameters, and the possible types of independent variables ("x") that can be used (e.g., netOD, netT, reflectance).

FittingFunction

Source code in app/calibration/functions.py

class FittingFunction:
    def __init__(self, name: str, func, initial_param_guess: list, description: str, param_names: list, independent_variable: str):
        """
        Initializes a FittingFunction instance.

        Parameters
        ----------
        name : str
            The name of the fitting function.
        func : callable
            The fitting function. The independent variable (x) should be its first argument.
        initial_param_guess : list
            A list of initial parameter guesses for the function.
        description : str
            A descriptive text (e.g., a LaTeX string) representing the fitting function.
        param_names : list
            A list of fitting parameter names (e.g., ["a", "b", "n"]).
        independent_variable : str
            The independent variable of the function (e.g., "netOD", "netT").
        """
        self.name = name
        self.func = func
        self.initial_param_guess = initial_param_guess
        self.description = description
        self.param_names = param_names
        self.independent_variable = independent_variable

    def __repr__(self):
        return (f"FittingFunction(name={self.name}, "
                f"param_names={self.param_names}, independent_variable={self.independent_variable})")

__init__(name, func, initial_param_guess, description, param_names, independent_variable)

Initializes a FittingFunction instance.

Parameters:

Name	Type	Description	Default
name	str	The name of the fitting function.	required
func	callable	The fitting function. The independent variable (x) should be its first argument.	required
initial_param_guess	list	A list of initial parameter guesses for the function.	required
description	str	A descriptive text (e.g., a LaTeX string) representing the fitting function.	required
param_names	list	A list of fitting parameter names (e.g., ["a", "b", "n"]).	required
independent_variable	str	The independent variable of the function (e.g., "netOD", "netT").	required

Source code in app/calibration/functions.py

def __init__(self, name: str, func, initial_param_guess: list, description: str, param_names: list, independent_variable: str):
    """
    Initializes a FittingFunction instance.

    Parameters
    ----------
    name : str
        The name of the fitting function.
    func : callable
        The fitting function. The independent variable (x) should be its first argument.
    initial_param_guess : list
        A list of initial parameter guesses for the function.
    description : str
        A descriptive text (e.g., a LaTeX string) representing the fitting function.
    param_names : list
        A list of fitting parameter names (e.g., ["a", "b", "n"]).
    independent_variable : str
        The independent variable of the function (e.g., "netOD", "netT").
    """
    self.name = name
    self.func = func
    self.initial_param_guess = initial_param_guess
    self.description = description
    self.param_names = param_names
    self.independent_variable = independent_variable

cuadratic(x, a, b, c)

Quadratic fitting function.

Parameters:

Name	Type	Description	Default
x	float or array - like	The independent variable (e.g., netOD).	required
a	float	Fitting parameters.	required
b	float	Fitting parameters.	required
c	float	Fitting parameters.	required

Returns:

Type	Description
float or array - like	Calculated dose.

Source code in app/calibration/functions.py

def cuadratic(x, a, b, c):
    """
    Quadratic fitting function.

    Parameters
    ----------
    x : float or array-like
        The independent variable (e.g., netOD).
    a, b, c : float
        Fitting parameters.

    Returns
    -------
    float or array-like
        Calculated dose.
    """
    return a * x**2 + b * x + c 

get_fitting_function(name)

Retrieves the FittingFunction instance by name.

Parameters:

Name	Type	Description	Default
name	str	The name of the fitting function (e.g., "polynomial").	required

Returns:

Type	Description
FittingFunction	The corresponding FittingFunction instance if found; otherwise, None.

Source code in app/calibration/functions.py

def get_fitting_function(name: str) -> FittingFunction:
    """
    Retrieves the FittingFunction instance by name.

    Parameters
    ----------
    name : str
        The name of the fitting function (e.g., "polynomial").

    Returns
    -------
    FittingFunction
        The corresponding FittingFunction instance if found; otherwise, None.
    """
    return fitting_functions.get(name)

polynomial(x, a, b, n)

Polynomial fitting function.

Parameters:

Name	Type	Description	Default
x	float or array - like	The independent variable (e.g., netOD).	required
a	float	Fitting parameters.	required
b	float	Fitting parameters.	required
n	float	Fitting parameters.	required

Returns:

Type	Description
float or array - like	Calculated dose.

Source code in app/calibration/functions.py

def polynomial(x, a, b, n):
    """
    Polynomial fitting function.

    Parameters
    ----------
    x : float or array-like
        The independent variable (e.g., netOD).
    a, b, n : float
        Fitting parameters.

    Returns
    -------
    float or array-like
        Calculated dose.
    """
    return a * x + b * (x ** n)

rational(x, a, b)

Rational fitting function.

Parameters:

Name	Type	Description	Default
x	float or array - like	The independent variable (e.g., netOD).	required
a	float	Fitting parameters.	required
b	float	Fitting parameters.	required

Returns:

Type	Description
float or array - like	Calculated dose.

Source code in app/calibration/functions.py

def rational(x, a, b):
    """
    Rational fitting function.

    Parameters
    ----------
    x : float or array-like
        The independent variable (e.g., netOD).
    a, b : float
        Fitting parameters.

    Returns
    -------
    float or array-like
        Calculated dose.
    """
    return (a * x) / (1 - b * x)

image_processing.py

This module provides functions for processing images, including reading TIFF and DICOM images, applying filters, displaying images, cropping regions of interest (ROIs), and performing template matching.

crop_square_roi(image, x, y, side_length)

Crops a square region of interest (ROI) from the image.

Parameters:

Name	Type	Description	Default
image	ndarray	Input image as a NumPy array.	required
x	int	X-coordinate (column) of the top-left corner of the ROI.	required
y	int	Y-coordinate (row) of the top-left corner of the ROI.	required
side_length	int	Length of one side of the square ROI.	required

Returns:

Type	Description
ndarray	The cropped ROI as a NumPy array.

Raises:

Type	Description
ValueError	If the ROI exceeds the image boundaries.

Source code in app/calibration/image_processing.py

def crop_square_roi(image, x, y, side_length):
    """
    Crops a square region of interest (ROI) from the image.

    Parameters
    ----------
    image : ndarray
        Input image as a NumPy array.
    x : int
        X-coordinate (column) of the top-left corner of the ROI.
    y : int
        Y-coordinate (row) of the top-left corner of the ROI.
    side_length : int
        Length of one side of the square ROI.

    Returns
    -------
    ndarray
        The cropped ROI as a NumPy array.

    Raises
    ------
    ValueError
        If the ROI exceeds the image boundaries.
    """
    # Get image dimensions (height and width)
    max_y, max_x = image.shape[:2]
    if y + side_length > max_y or x + side_length > max_x:
        raise ValueError("The ROI exceeds the image boundaries. Adjust (x, y) or 'side_length'.")
    roi = image[y:y+side_length, x:x+side_length]
    return roi

filter_image(image, filter_type=None, kernel_size=3)

Applies a filtering or preprocessing operation to the provided image.

Parameters:

Name	Type	Description	Default
image	ndarray	Input image data (grayscale).	required
filter_type	str	Type of filter to apply. Possible values include: "none", "gaussian", "median", "sobel", etc. This function currently demonstrates placeholder options.	None
kernel_size	int	The size of the filter kernel.	3

Returns:

Type	Description
ndarray	The filtered image as a NumPy array.

Source code in app/calibration/image_processing.py

def filter_image(image: np.ndarray, filter_type: str = None, kernel_size: int = 3) -> np.ndarray:
    """
    Applies a filtering or preprocessing operation to the provided image.

    Parameters
    ----------
    image : np.ndarray
        Input image data (grayscale).
    filter_type : str, optional
        Type of filter to apply. Possible values include:
          "none", "gaussian", "median", "sobel", etc.
        This function currently demonstrates placeholder options.
    kernel_size : int, optional
        The size of the filter kernel.

    Returns
    -------
    np.ndarray
        The filtered image as a NumPy array.
    """

    def wiener_filter_manual(img, K=30):
        kernel_size = 3
        h = gaussian(kernel_size, kernel_size / 3).reshape(kernel_size, 1)
        h = np.dot(h, h.transpose())
        h /= np.sum(h)
        kernel = h
        kernel /= np.sum(kernel)
        transformed = fft2(np.copy(img))
        kernel = fft2(kernel, s = img.shape)
        kernel = np.conj(kernel) / (np.abs(kernel) ** 2 + K)
        transformed = transformed * kernel
        wiener = np.abs(ifft2(transformed))
        return wiener

    # Ensure the image has only one channel (grayscale)
    if image.ndim > 2 and image.shape[-1] > 1:
        raise ValueError("The input image must be single-channel (grayscale).")

    if filter_type is None:
        filtered = image

    elif filter_type == "median":
        filtered = median_filter(image, size=kernel_size)

    elif filter_type == "wiener-scipy":
        filtered = wiener_scipy(image, mysize=(kernel_size, kernel_size))

    elif filter_type == "wiener-manual":
        filtered = wiener_filter_manual(image)

    elif filter_type == "wiener-skimage-1":
        k = kernel_size
        psf = np.ones((k, k)) / (k * k)
        filtered, _ = unsupervised_wiener(image, psf)

    elif filter_type == "wiener-skimage-2" or filter_type == "wiener":
        k = kernel_size
        psf = np.ones((k, k)) / (k * k)
        filtered = wiener(image, psf, balance=0.35)

    else:
        raise ValueError(f"Unknown filter type: {filter_type}.")

    return filtered

get_real_dimensions(image_path)

Calculates the physical dimensions (in centimeters) of an image from its metadata. Supports TIFF and DICOM image formats.

Parameters:

Name	Type	Description	Default
image_path	str	Path to the image file.	required

Returns:

Type	Description
tuple	(width_cm, height_cm) - the width and height in centimeters.

Raises:

Type	Description
ValueError	If required metadata is missing or the format is unsupported.

Source code in app/calibration/image_processing.py

def get_real_dimensions(image_path):
    """
    Calculates the physical dimensions (in centimeters) of an image from its metadata.
    Supports TIFF and DICOM image formats.

    Parameters
    ----------
    image_path : str
        Path to the image file.

    Returns
    -------
    tuple
        (width_cm, height_cm) - the width and height in centimeters.

    Raises
    ------
    ValueError
        If required metadata is missing or the format is unsupported.
    """
    ext = os.path.splitext(image_path)[-1].lower()

    if ext in [".tif", ".tiff"]:
        # Load the TIFF image using PIL
        img = Image.open(image_path)
        # Get the resolution in DPI (dots per inch)
        dpi = img.info.get("dpi", (300, 300))  # Assume 300 dpi if not specified
        # Get image dimensions in pixels
        width_px, height_px = img.size
        # Convert dimensions to centimeters (1 inch = 2.54 cm)
        width_cm = (width_px / dpi[0]) * 2.54
        height_cm = (height_px / dpi[1]) * 2.54
        return width_cm, height_cm

    elif ext in [".dcm"]:
        # Load the DICOM image
        ds = pydicom.dcmread(image_path)
        rows = ds.Rows
        cols = ds.Columns
        # Check for PixelSpacing information
        if hasattr(ds, "PixelSpacing"):
            pixel_spacing = ds.PixelSpacing  # [row spacing, column spacing] in mm
            width_cm = (cols * float(pixel_spacing[0])) / 10
            height_cm = (rows * float(pixel_spacing[1])) / 10
            return width_cm, height_cm
        else:
            raise ValueError("DICOM - PixelSpacing information not found.")
    else:
        raise ValueError("Unsupported format. Use TIFF or DICOM.")

read_image(image_path)

Reads an image from the specified file path using skimage.io. For TIFF images, a custom reader is used.

Parameters:

Name	Type	Description	Default
image_path	str	Path to the image file.	required

Returns:

Name	Type	Description
image	ndarray	The loaded image as a NumPy array.

Source code in app/calibration/image_processing.py

def read_image(image_path): 
    """
    Reads an image from the specified file path using skimage.io.
    For TIFF images, a custom reader is used.

    Parameters
    ----------
    image_path : str
        Path to the image file.

    Returns
    -------
    image : ndarray
        The loaded image as a NumPy array.
    """
    if image_path.lower().endswith('.tif') or image_path.lower().endswith('.tiff'):
        image = read_image_tif(image_path)
    else:
        image = io.imread(image_path)

    return image

read_image_tif(image_path)

Reads a TIFF image from the specified file path using tifffile. Normalizes the image to the range [0, 1] based on its bit depth.

Parameters:

Name	Type	Description	Default
image_path	str	Path to the TIFF image file.	required

Returns:

Name	Type	Description
image	ndarray	The loaded and normalized image as a NumPy array.

Source code in app/calibration/image_processing.py

def read_image_tif(image_path):
    """
    Reads a TIFF image from the specified file path using tifffile.
    Normalizes the image to the range [0, 1] based on its bit depth.

    Parameters
    ----------
    image_path : str
        Path to the TIFF image file.

    Returns
    -------
    image : ndarray
        The loaded and normalized image as a NumPy array.
    """
    # Read the TIFF image using tifffile
    image = tifffile.imread(image_path)
    # Print the number of bits per channel (e.g., 8, 16, or 32)
    print(f"Bits per channel: {tif_bits_per_channel(image_path)}")
    image = image / (2 ** tif_bits_per_channel(image_path) - 1)
    return image

reflect_image(image, mode='horizontal')

Reflects (flips) an image along the specified axis.

Parameters:

Name	Type	Description	Default
image	ndarray	The input image as a NumPy array.	required
mode	str	'horizontal' to flip left to right, 'vertical' to flip top to bottom.	'horizontal'

Returns:

Type	Description
ndarray	The reflected image.

Source code in app/calibration/image_processing.py

def reflect_image(image, mode='horizontal'):
    """
    Reflects (flips) an image along the specified axis.

    Parameters
    ----------
    image : ndarray
        The input image as a NumPy array.
    mode : str, optional
        'horizontal' to flip left to right, 'vertical' to flip top to bottom.

    Returns
    -------
    ndarray
        The reflected image.
    """
    if mode.lower() == 'horizontal':
        # Horizontal flip: reflect left <-> right
        reflected_image = np.fliplr(image)
    elif mode.lower() == 'vertical':
        # Vertical flip: reflect top <-> bottom
        reflected_image = np.flipud(image)
    else:
        raise ValueError("Mode must be 'horizontal' or 'vertical'.")

    return reflected_image

rotate_image(image, times=1, direction='clockwise')

Rotates an image by 90-degree increments.

Parameters:

Name	Type	Description	Default
image	ndarray	NumPy array representing the image.	required
times	int	Number of 90° rotations (e.g., 1 for 90°, 2 for 180°, etc.).	1
direction	str	'clockwise' for clockwise rotation, or 'counterclockwise' for counterclockwise rotation.	'clockwise'

Returns:

Type	Description
ndarray	The rotated image.

Source code in app/calibration/image_processing.py

def rotate_image(image, times=1, direction='clockwise'):
    """
    Rotates an image by 90-degree increments.

    Parameters
    ----------
    image : ndarray
        NumPy array representing the image.
    times : int, optional
        Number of 90° rotations (e.g., 1 for 90°, 2 for 180°, etc.).
    direction : str, optional
        'clockwise' for clockwise rotation, or 'counterclockwise' for counterclockwise rotation.

    Returns
    -------
    ndarray
        The rotated image.
    """
    # Limit the number of rotations to the range [0, 3]
    times = times % 4

    if direction.lower() == 'clockwise':
        # np.rot90 rotates counterclockwise by default;
        # rotate by a negative number of times to rotate clockwise.
        k = -times
    elif direction.lower() == 'counterclockwise':
        k = times
    else:
        raise ValueError("The 'direction' parameter must be 'clockwise' or 'counterclockwise'.")

    rotated_image = np.rot90(image, k=k)
    return rotated_image

save_image(image, output_path)

Saves an image (NumPy array) to the specified file path.

Parameters:

Name	Type	Description	Default
image	ndarray	NumPy array representing the image to be saved.	required
output_path	str	Full path of the output file (including the file format, e.g., .tiff or .png).	required

Source code in app/calibration/image_processing.py

def save_image(image, output_path):
    """
    Saves an image (NumPy array) to the specified file path.

    Parameters
    ----------
    image : ndarray
        NumPy array representing the image to be saved.
    output_path : str
        Full path of the output file (including the file format, e.g., .tiff or .png).
    """
    io.imsave(output_path, image)

show_image(image, title=None, show_labels=True, show_axis=True)

Displays the provided image along with an optional coordinate system to facilitate region of interest (ROI) selection.

Parameters:

Name	Type	Description	Default
image	ndarray	The image to display.	required
title	str	Title for the displayed image.	None
show_labels	bool	If True, display axis labels.	True
show_axis	bool	If True, display the axis.	True

Returns:

Name	Type	Description
image	ndarray	The same image that was displayed.

Source code in app/calibration/image_processing.py

def show_image(image, title=None, show_labels=True, show_axis=True):
    """
    Displays the provided image along with an optional coordinate system
    to facilitate region of interest (ROI) selection.

    Parameters
    ----------
    image : ndarray
        The image to display.
    title : str, optional
        Title for the displayed image.
    show_labels : bool, optional
        If True, display axis labels.
    show_axis : bool, optional
        If True, display the axis.

    Returns
    -------
    image : ndarray
        The same image that was displayed.
    """
    # Create a figure and display the image (x corresponds to columns, y to rows)
    plt.figure(figsize=(6, 6))
    plt.imshow(image)
    if title is not None:
        plt.title(title)
    if show_labels:
        plt.xlabel("X Coordinate (column)")
        plt.ylabel("Y Coordinate (row)")
    if show_axis:
        plt.axis('on')
    plt.show()

    return image

template_matching(TPS_map_path, film_tif_path, output_path, verbose=False)

Performs template matching between a dose map (TPS map) from a DICOM file and an original film TIFF image. The function aligns the images using template matching, crops the matching region, applies transformations to standardize orientation, and saves the result.

Parameters:

Name	Type	Description	Default
TPS_map_path	str	File path to the TPS map DICOM file.	required
film_tif_path	str	File path to the original film TIFF image.	required
output_path	str	File path where the processed (cropped and transformed) image will be saved.	required

Source code in app/calibration/image_processing.py

def template_matching(TPS_map_path, film_tif_path, output_path, verbose=False):
    """
    Performs template matching between a dose map (TPS map) from a DICOM file and an original film TIFF image.
    The function aligns the images using template matching, crops the matching region,
    applies transformations to standardize orientation, and saves the result.

    Parameters
    ----------
    TPS_map_path : str
        File path to the TPS map DICOM file.
    film_tif_path : str
        File path to the original film TIFF image.
    output_path : str
        File path where the processed (cropped and transformed) image will be saved.
    """
    # Define file paths
    imageA_path = TPS_map_path
    film_tif_full_path = film_tif_path

    # Load the dose map (Image A) from the DICOM file
    imageA = pydicom.dcmread(imageA_path).pixel_array
    imageA = imageA.astype(np.float32)

    # Load the original film TIFF image (preserving its properties)
    film_orig = tifffile.imread(film_tif_full_path)

    # For processing, use the green channel if the image is multichannel; otherwise, use a copy of the original image
    if film_orig.ndim == 3:
        imageB = film_orig[:, :, 1]
    else:
        imageB = film_orig.copy()
    imageB = imageB.astype(np.float32)

    # Verify that the images have been loaded correctly
    if imageA is None or imageA.size == 0:
        print(f'Error loading Image A from {imageA_path}')
    if imageB is None or imageB.size == 0:
        print(f'Error loading Image B from {film_tif_full_path}')

    # Assume the images are 2D arrays
    heightA, widthA = imageA.shape
    heightB, widthB = imageB.shape

    # Obtain real dimensions in centimeters for both images
    widthA_cm, heightA_cm = get_real_dimensions(imageA_path)
    widthB_cm, heightB_cm = get_real_dimensions(film_tif_full_path)
    if verbose:
        print("Image A dimensions (cm):", widthA_cm, heightA_cm)
        print("Image B dimensions (cm):", widthB_cm, heightB_cm)

    # Calculate the resolution (pixels per cm) for each image
    resolutionA_x = widthA / widthA_cm
    resolutionA_y = heightA / heightA_cm
    resolutionB_x = widthB / widthB_cm
    resolutionB_y = heightB / heightB_cm

    if verbose:
        print(f'Image A resolution: {resolutionA_x:.2f} px/cm x {resolutionA_y:.2f} px/cm')
        print(f'Image B resolution: {resolutionB_x:.2f} px/cm x {resolutionB_y:.2f} px/cm')

    # If the resolutions differ, rescale Image A to match the physical scale of Image B
    if abs(resolutionA_x - resolutionB_x) > 1e-2 or abs(resolutionA_y - resolutionB_y) > 1e-2:
        new_widthA = int(widthA_cm * resolutionB_x)
        new_heightA = int(heightA_cm * resolutionB_y)
        imageA = cv2.resize(imageA, (new_widthA, new_heightA), interpolation=cv2.INTER_LINEAR)
        if verbose:
            print(f'Rescaled Image A to: {new_widthA}x{new_heightA} pixels')

    # Normalize both images to the range [0, 1]
    imageA = cv2.normalize(imageA, None, 0, 1, cv2.NORM_MINMAX, dtype=cv2.CV_32F)
    imageB = cv2.normalize(imageB, None, 0, 1, cv2.NORM_MINMAX, dtype=cv2.CV_32F)

    # Transform Image B for template matching:
    # Apply horizontal flip and 180° rotation (equivalent to a vertical flip)
    imageB_trans = cv2.flip(imageB, 1)
    imageB_trans = cv2.rotate(imageB_trans, cv2.ROTATE_180)
    # Invert intensities (dark becomes white and vice versa)
    imageB_trans = 1.0 - imageB_trans

    # Display the transformed images for template matching
    fig, axs = plt.subplots(1, 2, figsize=(12, 6))
    im1 = axs[0].imshow(imageA, cmap='gray')
    axs[0].set_title("Image A (Template)")
    axs[0].axis("off")
    fig.colorbar(im1, ax=axs[0])
    im2 = axs[1].imshow(imageB_trans, cmap='gray')
    axs[1].set_title("Transformed and Inverted Image B")
    axs[1].axis("off")
    fig.colorbar(im2, ax=axs[1])
    #plt.show()

    # Apply template matching using the TM_CCOEFF_NORMED method
    result = cv2.matchTemplate(imageB_trans, imageA, cv2.TM_CCOEFF_NORMED)
    min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(result)
    if verbose:
        print(f'Maximum correlation value: {max_val:.3f}')
        print(f'Top-left location in transformed image: {max_loc}')

    # Get the template size (dimensions of Image A)
    template_height, template_width = imageA.shape

    # Given that the transformation applied to imageB is equivalent to a vertical flip,
    # the coordinate mapping from imageB_trans to the original film image is:
    #   (x, y)_orig = (x, H - y - template_height)
    # where H is the height of the original Image B.
    H = heightB  
    col_t, row_t = max_loc
    orig_top = H - row_t - template_height
    orig_left = col_t
    orig_bottom = H - row_t
    orig_right = col_t + template_width

    if verbose:
        print('Coordinates in the original image for the crop:')
        print(f'  Top-left: ({orig_left}, {orig_top})')
        print(f'  Bottom-right: ({orig_right}, {orig_bottom})')

    # To visualize the bounding box on the original image:
    # Convert film_orig to 8-bit for display (preserving channels if present)
    if film_orig.ndim == 3:
        if film_orig.dtype != np.uint8:
            film_disp = cv2.normalize(film_orig, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
        else:
            film_disp = film_orig.copy()
    else:
        film_disp = cv2.normalize(film_orig, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
        film_disp = cv2.cvtColor(film_disp, cv2.COLOR_GRAY2BGR)

    # Draw the bounding box on the original image
    cv2.rectangle(film_disp, (orig_left, orig_top), (orig_right, orig_bottom), (0, 0, 255), 2)
    plt.figure(figsize=(6,6))
    # Convert BGR to RGB for matplotlib if the image is in color
    if film_disp.ndim == 3:
        film_disp_rgb = cv2.cvtColor(film_disp, cv2.COLOR_BGR2RGB)
        plt.imshow(film_disp_rgb)
    else:
        plt.imshow(film_disp, cmap='gray')
    plt.title("Original Image with Bounding Box")
    plt.axis("off")
    #plt.show()

    # Crop the identified region from the original image (without transformation)
    cropped = film_orig[orig_top:orig_bottom, orig_left:orig_right].copy()

    # To match the orientation of the transformed image, apply the same transformation:
    # horizontal flip and 180° rotation.
    cropped_trans = cv2.flip(cropped, 1)
    cropped_trans = cv2.rotate(cropped_trans, cv2.ROTATE_180)

    # Save the transformed crop as a TIFF file (preserving image properties as much as possible)
    tifffile.imwrite(output_path, cropped_trans)
    print(f'Cropped (flipped and rotated) image saved as {output_path}')

tif_bits_per_channel(image_path)

Retrieves the number of bits per channel for a TIFF image.

Parameters:

Name	Type	Description	Default
image_path	str	Path to the TIFF image.	required

Returns:

Type	Description
int	The number of bits per channel.

Raises:

Type	Description
ValueError	If the BitsPerSample tag is not found or if its values are not identical.

Source code in app/calibration/image_processing.py

def tif_bits_per_channel(image_path: str) -> int:
    """
    Retrieves the number of bits per channel for a TIFF image.

    Parameters
    ----------
    image_path : str
        Path to the TIFF image.

    Returns
    -------
    int
        The number of bits per channel.

    Raises
    ------
    ValueError
        If the BitsPerSample tag is not found or if its values are not identical.
    """
    with Image.open(image_path) as img:
        meta_dict = {TAGS[key]: img.tag[key] for key in img.tag.keys()}
        bits_per_sample = meta_dict.get('BitsPerSample', None)

        if bits_per_sample is None:
            raise ValueError('BitsPerSample tag not found in the image.')
        # If there is more than one value, ensure they are all identical.
        elif len(set(bits_per_sample)) != 1:
            raise ValueError('BitsPerSample values are not identical.')
        return bits_per_sample[0]

Dose Analysis

display_dose_profile(path, x=None, y=None, save_fig=False)

Loads a dose map from a file (.npy or .dcm) and displays interactive profiles: - Dose map image with reference lines. - Horizontal (row) and vertical (column) profiles with interactive sliders.

Parameters:

Name	Type	Description	Default
path	str	File path to the dose map (.npy or .dcm).	required
x	array - like	Horizontal coordinates; if None, pixel indices are used.	None
y	array - like	Vertical coordinates; if None, pixel indices are used.	None
save_fig	bool	If True, saves the figure as a PNG file.	False

Source code in app/dose_analysis/dose_profile.py

def display_dose_profile(path, x=None, y=None, save_fig=False):
    """
    Loads a dose map from a file (.npy or .dcm) and displays interactive profiles:
      - Dose map image with reference lines.
      - Horizontal (row) and vertical (column) profiles with interactive sliders.

    Parameters
    ----------
    path : str
        File path to the dose map (.npy or .dcm).
    x : array-like, optional
        Horizontal coordinates; if None, pixel indices are used.
    y : array-like, optional
        Vertical coordinates; if None, pixel indices are used.
    save_fig : bool, optional
        If True, saves the figure as a PNG file.
    """

    # Load the file based on its extension
    if path.endswith('.dcm'):
        ds = pydicom.dcmread(path)
        dose_map = ds.pixel_array * ds.DoseGridScaling
    elif path.endswith('.npy'):
        dose_map = np.load(path)
    else:
        raise ValueError("Unsupported file format. Use .npy or .dcm")

    # Get the base name of the file for output naming
    output_name = os.path.splitext(os.path.basename(path))[0]

    # Dimensions of the dose map
    nx, ny = dose_map.shape
    if x is None:
        x = np.arange(ny)
    if y is None:
        y = np.arange(nx)

    max_dose = dose_map.max()
    default_row = nx // 2
    default_col = ny // 2

    # For the inverted slider, adjust the initial value:
    # Since the slider value is inverted, the real index is (nx - 1) - slider value.
    default_slider_value = (nx - 1) - default_row

    # Create the main figure with a custom layout
    fig = plt.figure(figsize=(12, 8))
    gs_main = gridspec.GridSpec(1, 2, width_ratios=[3, 2],
                                 left=0.1, right=0.9, top=0.9, bottom=0.3, wspace=0.1)

    # Dose map panel: display the dose map with reference lines
    ax_map = fig.add_subplot(gs_main[0, 0])
    extent = [x[0], x[-1], y[0], y[-1]]
    im = ax_map.imshow(dose_map, cmap='inferno')  # You may add extent=extent if required
    ax_map.set_title("Dose Map")
    ax_map.set_xlabel("X")
    ax_map.set_ylabel("Y")
    # Initial reference lines using default_row and default_col
    line_h = ax_map.axhline(y=y[default_row], color='cyan', lw=2)
    line_v = ax_map.axvline(x=x[default_col], color='lime', lw=2)

    # Profiles panel: subplots for horizontal and vertical profiles
    gs_profiles = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=gs_main[0, 1], hspace=0.4)

    # Horizontal profile (row): using default_row index
    ax_hprofile = fig.add_subplot(gs_profiles[0, 0])
    h_profile_line, = ax_hprofile.plot(x, dose_map[default_row, :], lw=2)
    ax_hprofile.set_title(f"Horizontal Profile (row {default_row})")
    ax_hprofile.set_xlabel("X")
    ax_hprofile.set_ylabel("Dose (Gy)")
    ax_hprofile.grid(True)
    ax_hprofile.set_ylim(0, max_dose)

    # Vertical profile (column): using default_col index
    ax_vprofile = fig.add_subplot(gs_profiles[1, 0])
    v_profile_line, = ax_vprofile.plot(y, dose_map[:, default_col], lw=2, color='green')
    ax_vprofile.set_title(f"Vertical Profile (column {default_col})")
    ax_vprofile.set_xlabel("Y")
    ax_vprofile.set_ylabel("Dose (Gy)")
    ax_vprofile.grid(True)
    ax_vprofile.set_ylim(0, max_dose)

    # Create slider for the horizontal reference (row)
    ax_slider_row = fig.add_axes([0.02, 0.3, 0.03, 0.6])
    # Note: instead of using invert_yaxis(), we map the value in the callback
    slider_row = Slider(ax_slider_row, 'Row (Y)', 0, nx - 1,
                        valinit=default_slider_value, valstep=1, 
                        orientation='vertical', color='dimgray')

    # Create slider for the vertical reference (column)
    ax_slider_col = fig.add_axes([0.13, 0.20, 0.4, 0.03])
    slider_col = Slider(ax_slider_col, 'Column (X)', 0, ny - 1,
                        valinit=default_col, valstep=1, color='dimgray')

    # Update function for the row slider:
    def update_row(val):
        # Apply the inverse transformation: the actual row index is (nx - 1) - slider value
        row = int((nx - 1) - slider_row.val)
        line_h.set_ydata([y[row], y[row]])
        h_profile_line.set_ydata(dose_map[row, :])
        ax_hprofile.set_title(f"Horizontal Profile (row {row})")
        # Manually update the slider text to show the actual row index
        slider_row.valtext.set_text(str(row))
        fig.canvas.draw_idle()

    # Update function for the column slider:
    def update_col(val):
        col = int(slider_col.val)
        line_v.set_xdata([x[col], x[col]])
        v_profile_line.set_ydata(dose_map[:, col])
        ax_vprofile.set_title(f"Vertical Profile (column {col})")
        fig.canvas.draw_idle()

    # Connect the sliders to the update functions
    slider_row.on_changed(update_row)
    slider_col.on_changed(update_col)

    # Save the figure if requested
    if save_fig:
        fig.savefig(f"{output_name}_dose_profiles.png", bbox_inches='tight', dpi=300)

    plt.show()

plot_isodose_map(path, save_fig=False)

Plots an isodose map from a given file and provides interactive sliders to adjust the contour levels and image opacity.

The input file can be either a DICOM (.dcm) file or a NumPy binary (.npy) file. In the case of a DICOM file, the function multiplies the pixel array by the DoseGridScaling attribute. The isodose levels are initially set at 30%, 60%, and 90% of the maximum dose.

Parameters:

Name	Type	Description	Default
path	str	Path to the dose map file (supported formats: .dcm or .npy).	required
save_fig	bool	If True, saves the resulting figure as a PNG file with a name based on the input file.	False

Raises:

Type	Description
ValueError	If the file format is not supported.

Source code in app/dose_analysis/isodoses.py

def plot_isodose_map(path, save_fig=False):
    """
    Plots an isodose map from a given file and provides interactive sliders to adjust the contour levels and image opacity.

    The input file can be either a DICOM (.dcm) file or a NumPy binary (.npy) file.
    In the case of a DICOM file, the function multiplies the pixel array by the DoseGridScaling attribute.
    The isodose levels are initially set at 30%, 60%, and 90% of the maximum dose.

    Parameters
    ----------
    path : str
        Path to the dose map file (supported formats: .dcm or .npy).
    save_fig : bool, optional
        If True, saves the resulting figure as a PNG file with a name based on the input file.

    Raises
    ------
    ValueError
        If the file format is not supported.
    """
    # Define colors for the isodose contours
    curves_colors = ['darkorange', 'turquoise', 'blueviolet']
    # Get the base name (without extension) of the file to use for saving the figure
    output_name = os.path.splitext(os.path.basename(path))[0]

    # Load the file based on its extension
    if path.endswith('.dcm'):
        ds = pydicom.dcmread(path)
        dose_map = ds.pixel_array * ds.DoseGridScaling
    elif path.endswith('.npy'):
        dose_map = np.load(path)
    else:
        raise ValueError("Unsupported file format. Use .dcm or .npy")

    max_dose = dose_map.max()

    # Create a matplotlib figure and axis for displaying the dose map
    fig, ax = plt.subplots()
    plt.subplots_adjust(bottom=0.25, right=0.85)
    im = ax.imshow(dose_map, alpha=0.2, cmap='viridis_r')
    ax.set_title("Isodoses")

    # Set initial contour levels as percentages of the maximum dose
    initial_values = [30, 60, 90]
    initial_levels = [max_dose * p / 100 for p in initial_values]

    # Draw the initial contours using the preset levels and colors
    cs1 = ax.contour(dose_map, levels=[initial_levels[0]], colors=curves_colors[0])
    cs2 = ax.contour(dose_map, levels=[initial_levels[1]], colors=curves_colors[1])
    cs3 = ax.contour(dose_map, levels=[initial_levels[2]], colors=curves_colors[2])
    contours = [cs1, cs2, cs3]

    # Define the color for the slider background
    axcolor = 'lightgoldenrodyellow'
    # Create axes for the three sliders that control the contour levels
    slider_ax1 = plt.axes([0.15, 0.1, 0.65, 0.03], facecolor=axcolor)
    slider_ax2 = plt.axes([0.15, 0.06, 0.65, 0.03], facecolor=axcolor)
    slider_ax3 = plt.axes([0.15, 0.02, 0.65, 0.03], facecolor=axcolor)

    # Create sliders for each contour level (with no label)
    slider1 = Slider(slider_ax1, '', 0, 100, valinit=30, valstep=1,
                     valfmt='%d%%', facecolor=curves_colors[0])
    slider2 = Slider(slider_ax2, '', 0, 100, valinit=60, valstep=1,
                     valfmt='%d%%', facecolor=curves_colors[1])
    slider3 = Slider(slider_ax3, '', 0, 100, valinit=90, valstep=1,
                     valfmt='%d%%', facecolor=curves_colors[2])

    # Hide the vertical lines on the sliders, if they exist
    for slider in [slider1, slider2, slider3]:
        if hasattr(slider, 'vline'):
            slider.vline.set_visible(False)

    # Create an axis for the opacity slider (vertical slider)
    slider_ax_alpha = plt.axes([0.9, 0.25, 0.03, 0.65], facecolor=axcolor)
    alpha_slider = Slider(slider_ax_alpha, 'Opacity', 0, 1, valinit=0.2,
                          orientation='vertical', valfmt='%1.2f', valstep=0.01)

    # Hide the horizontal line of the opacity slider, if it exists
    if hasattr(alpha_slider, 'hline'):
        alpha_slider.hline.set_visible(False)

    def update(val):
        """
        Update callback for the contour level sliders.
        Removes old contour lines and draws new ones based on slider values.
        """
        nonlocal contours
        # Remove existing contours
        for cs in contours:
            cs.remove()
        contours.clear()

        # Compute new levels based on current slider values
        levels = [max_dose * slider.val / 100 for slider in [slider1, slider2, slider3]]
        # Draw new contours with updated levels and colors
        cs1 = ax.contour(dose_map, levels=[levels[0]], colors=curves_colors[0])
        cs2 = ax.contour(dose_map, levels=[levels[1]], colors=curves_colors[1])
        cs3 = ax.contour(dose_map, levels=[levels[2]], colors=curves_colors[2])
        contours = [cs1, cs2, cs3]
        fig.canvas.draw_idle()

    def update_alpha(val):
        """
        Update callback for the opacity slider.
        Adjusts the alpha transparency of the displayed dose map.
        """
        im.set_alpha(val)
        fig.canvas.draw_idle()

    # Register the update callbacks with the sliders
    slider1.on_changed(update)
    slider2.on_changed(update)
    slider3.on_changed(update)
    alpha_slider.on_changed(update_alpha)

    # Save the figure if requested
    if save_fig:
        fig.savefig(f"{output_name}_isodoses.png", dpi=300)
    plt.show()