optimiser.stokes

Stokes-based pixel-wise stress inversion.

This module implements stress field recovery using normalized Stokes components computed from polarimetric measurements. This is the primary and recommended method for most photoelastic analysis tasks.

Key Functions

• compute_stokes_components() - Calculate Stokes parameters from intensity measurements
• compute_normalized_stokes() - Normalize Stokes components for analysis
• recover_stress_map_stokes() - Main function for stress field recovery
• predict_stokes() - Forward model for validation

stokes

Local stress measurement using polarimetric imaging.

This module implements the local stress measurement algorithm using Mueller matrix calculus and multi-wavelength polarimetry to recover the full 2D stress tensor at each pixel from polarimetric images.

Functions

predict_stokes(sigma_xx, sigma_yy, sigma_xy, C, nu, L, wavelength, S_i_hat)

Predict normalized Stokes vector components from stress tensor.

Parameters:

Name	Type	Description	Default
sigma_xx	float	Normal stress component in x direction (Pa).	required
sigma_yy	float	Normal stress component in y direction (Pa).	required
sigma_xy	float	Shear stress component (Pa).	required
C	float	Stress-optic coefficient (1/Pa).	required
nu	float	Solid fraction (use 1.0 for solid samples).	required
L	float	Sample thickness (m).	required
wavelength	float	Wavelength of light (m).	required
S_i_hat	array - like	Incoming normalized Stokes vector [S1_hat, S2_hat] or [S1_hat, S2_hat, S3_hat]. If 2 elements, S3_hat is assumed to be 0 (no circular polarization). If 3 elements, S3_hat represents circular polarization component.	required

Returns:

Name	Type	Description
S_p_hat	ndarray	Predicted normalized Stokes components [S1_hat, S2_hat].

Source code in photoelastimetry/optimiser/stokes.py

def predict_stokes(sigma_xx, sigma_yy, sigma_xy, C, nu, L, wavelength, S_i_hat):
    """
    Predict normalized Stokes vector components from stress tensor.

    Parameters
    ----------
    sigma_xx : float
        Normal stress component in x direction (Pa).
    sigma_yy : float
        Normal stress component in y direction (Pa).
    sigma_xy : float
        Shear stress component (Pa).
    C : float
        Stress-optic coefficient (1/Pa).
    nu : float
        Solid fraction (use 1.0 for solid samples).
    L : float
        Sample thickness (m).
    wavelength : float
        Wavelength of light (m).
    S_i_hat : array-like
        Incoming normalized Stokes vector [S1_hat, S2_hat] or [S1_hat, S2_hat, S3_hat].
        If 2 elements, S3_hat is assumed to be 0 (no circular polarization).
        If 3 elements, S3_hat represents circular polarization component.

    Returns
    -------
    S_p_hat : ndarray
        Predicted normalized Stokes components [S1_hat, S2_hat].
    """
    theta = compute_principal_angle(sigma_xx, sigma_yy, sigma_xy)
    delta = compute_retardance(sigma_xx, sigma_yy, sigma_xy, C, nu, L, wavelength)

    M = mueller_matrix(theta, delta)

    # Extend S_i_hat to full Stokes vector
    S_i_hat = np.asarray(S_i_hat)
    if len(S_i_hat) == 2:
        # Backward compatibility: assume S3 = 0 (no circular polarization)
        S_i_full = np.array([1.0, S_i_hat[0], S_i_hat[1], 0.0])
    elif len(S_i_hat) == 3:
        # Use provided circular component
        S_i_full = np.array([1.0, S_i_hat[0], S_i_hat[1], S_i_hat[2]])
    else:
        raise ValueError(f"S_i_hat must have 2 or 3 elements, got {len(S_i_hat)}")

    # Apply Mueller matrix
    S_m = M @ S_i_full

    # Return normalized components (excluding S0)
    S_p_hat = S_m[1:3]

    return S_p_hat

compute_residual(stress_params, S_m_hat, wavelengths, C_values, nu, L, S_i_hat)

Compute residual between measured and predicted Stokes components.

Parameters:

Name	Type	Description	Default
stress_params	array - like	Stress tensor components [sigma_xx, sigma_yy, sigma_xy].	required
S_m_hat	ndarray	Measured normalized Stokes components, shape (3, 2) for RGB channels.	required
wavelengths	array - like	Wavelengths for R, G, B channels (m).	required
C_values	array - like	Stress-optic coefficients for R, G, B channels (1/Pa).	required
nu	float	Solid fraction (use 1.0 for solid samples).	required
L	float	Sample thickness (m).	required
S_i_hat	array - like	Incoming normalized Stokes vector [S1_hat, S2_hat] or [S1_hat, S2_hat, S3_hat].	required

Returns:

Name	Type	Description
residual	float	Sum of squared residuals across all colour channels.

Source code in photoelastimetry/optimiser/stokes.py

def compute_residual(stress_params, S_m_hat, wavelengths, C_values, nu, L, S_i_hat):
    """
    Compute residual between measured and predicted Stokes components.

    Parameters
    ----------
    stress_params : array-like
        Stress tensor components [sigma_xx, sigma_yy, sigma_xy].
    S_m_hat : ndarray
        Measured normalized Stokes components, shape (3, 2) for RGB channels.
    wavelengths : array-like
        Wavelengths for R, G, B channels (m).
    C_values : array-like
        Stress-optic coefficients for R, G, B channels (1/Pa).
    nu : float
        Solid fraction (use 1.0 for solid samples).
    L : float
        Sample thickness (m).
    S_i_hat : array-like
        Incoming normalized Stokes vector [S1_hat, S2_hat] or [S1_hat, S2_hat, S3_hat].

    Returns
    -------
    residual : float
        Sum of squared residuals across all colour channels.
    """
    sigma_xx, sigma_yy, sigma_xy = stress_params

    residual = 0.0
    for c in range(3):  # R, G, B
        S_p_hat = predict_stokes(
            sigma_xx,
            sigma_yy,
            sigma_xy,
            C_values[c],
            nu,
            L,
            wavelengths[c],
            S_i_hat,
        )
        diff = S_m_hat[c] - S_p_hat
        residual += np.sum(diff**2)

    return residual

recover_stress_tensor(S_m_hat, wavelengths, C_values, nu, L, S_i_hat, initial_guess=None, track_history=False, max_fringes=6)

Recover stress tensor components by minimizing residual.

Parameters:

Name	Type	Description	Default
S_m_hat	ndarray	Measured normalized Stokes components, shape (3, 2) for RGB channels. Each row is [S1_hat, S2_hat] for a colour channel.	required
wavelengths	array - like	Wavelengths for R, G, B channels (m).	required
C_values	array - like	Stress-optic coefficients for R, G, B channels (1/Pa).	required
nu	float	Solid fraction (use 1.0 for solid samples).	required
L	float	Sample thickness (m).	required
S_i_hat	array - like	Incoming normalized Stokes vector [S1_hat, S2_hat] or [S1_hat, S2_hat, S3_hat].	required
initial_guess	array - like	Initial guess for stress tensor [sigma_xx, sigma_yy, sigma_xy]. Default is [1, 1, 1].	None
track_history	bool	If True, track optimization history for debugging plots. Default is False.	False
max_fringes	float	Maximum expected fringe order. Used to set bounds on stress components. Default is 6 fringes, which corresponds to ~1.7 MPa for typical materials.	6

Returns:

Name	Type	Description
stress_tensor	ndarray	Recovered stress tensor components [sigma_xx, sigma_yy, sigma_xy].
success	bool	Whether optimization was successful.
history	(dict, optional)	Only returned if track_history=True. Contains: - 'all_paths': list of dicts, each containing the optimization path from a start point - 'best_path_index': index of the path that led to the best solution

Source code in photoelastimetry/optimiser/stokes.py

def recover_stress_tensor(
    S_m_hat, wavelengths, C_values, nu, L, S_i_hat, initial_guess=None, track_history=False, max_fringes=6
):
    """
    Recover stress tensor components by minimizing residual.

    Parameters
    ----------
    S_m_hat : ndarray
        Measured normalized Stokes components, shape (3, 2) for RGB channels.
        Each row is [S1_hat, S2_hat] for a colour channel.
    wavelengths : array-like
        Wavelengths for R, G, B channels (m).
    C_values : array-like
        Stress-optic coefficients for R, G, B channels (1/Pa).
    nu : float
        Solid fraction (use 1.0 for solid samples).
    L : float
        Sample thickness (m).
    S_i_hat : array-like
        Incoming normalized Stokes vector [S1_hat, S2_hat] or [S1_hat, S2_hat, S3_hat].
    initial_guess : array-like, optional
        Initial guess for stress tensor [sigma_xx, sigma_yy, sigma_xy].
        Default is [1, 1, 1].
    track_history : bool, optional
        If True, track optimization history for debugging plots. Default is False.
    max_fringes : float, optional
        Maximum expected fringe order. Used to set bounds on stress components.
        Default is 6 fringes, which corresponds to ~1.7 MPa for typical materials.

    Returns
    -------
    stress_tensor : ndarray
        Recovered stress tensor components [sigma_xx, sigma_yy, sigma_xy].
    success : bool
        Whether optimization was successful.
    history : dict, optional
        Only returned if track_history=True. Contains:
        - 'all_paths': list of dicts, each containing the optimization path from a start point
        - 'best_path_index': index of the path that led to the best solution
    """
    if initial_guess is None:
        initial_guess = np.array([1.0, 1.0, 0.0])

    # Run core optimization
    if track_history:
        result, all_paths = _optimize_stress_tensor(
            S_m_hat,
            wavelengths,
            C_values,
            nu,
            L,
            S_i_hat,
            initial_guess,
            max_fringes,
            callback=None,
            track_all_paths=True,
        )

        # Find which path was the best
        best_path_index = None
        for i, path in enumerate(all_paths):
            if path["is_best"]:
                best_path_index = i
                break

        history = {"all_paths": all_paths, "best_path_index": best_path_index}
        return result.x, result.success, history
    else:
        result = _optimize_stress_tensor(
            S_m_hat,
            wavelengths,
            C_values,
            nu,
            L,
            S_i_hat,
            initial_guess,
            max_fringes,
            callback=None,
            track_all_paths=False,
        )
        return result.x, result.success

recover_stress_tensor_live(S_m_hat, wavelengths, C_values, nu, L, S_i_hat, initial_guess=None, update_interval=5, max_fringes=6)

Recover stress tensor with live plotting of optimization progress.

This function is useful for debugging and understanding the optimization process for a single pixel. It creates a live-updating plot that shows how the stress components and predicted Stokes parameters evolve during optimization.

Parameters:

Name	Type	Description	Default
S_m_hat	ndarray	Measured normalized Stokes components, shape (3, 2) for RGB channels.	required
wavelengths	array - like	Wavelengths for R, G, B channels (m).	required
C_values	array - like	Stress-optic coefficients for R, G, B channels (1/Pa).	required
nu	float	Solid fraction (use 1.0 for solid samples).	required
L	float	Sample thickness (m).	required
S_i_hat	array - like	Incoming normalized Stokes vector [S1_hat, S2_hat] or [S1_hat, S2_hat, S3_hat].	required
initial_guess	array - like	Initial guess for stress tensor [sigma_xx, sigma_yy, sigma_xy].	None
update_interval	int	Update plot every N iterations. Default is 5.	5
max_fringes	float	Maximum expected fringe order for setting bounds. Default is 6.	6

Returns:

Name	Type	Description
stress_tensor	ndarray	Recovered stress tensor components [sigma_xx, sigma_yy, sigma_xy].
success	bool	Whether optimization was successful.
history	dict	Optimization history for further analysis.
fig	Figure	The figure object (will be kept open).

Source code in photoelastimetry/optimiser/stokes.py

def recover_stress_tensor_live(
    S_m_hat, wavelengths, C_values, nu, L, S_i_hat, initial_guess=None, update_interval=5, max_fringes=6
):
    """
    Recover stress tensor with live plotting of optimization progress.

    This function is useful for debugging and understanding the optimization process
    for a single pixel. It creates a live-updating plot that shows how the stress
    components and predicted Stokes parameters evolve during optimization.

    Parameters
    ----------
    S_m_hat : ndarray
        Measured normalized Stokes components, shape (3, 2) for RGB channels.
    wavelengths : array-like
        Wavelengths for R, G, B channels (m).
    C_values : array-like
        Stress-optic coefficients for R, G, B channels (1/Pa).
    nu : float
        Solid fraction (use 1.0 for solid samples).
    L : float
        Sample thickness (m).
    S_i_hat : array-like
        Incoming normalized Stokes vector [S1_hat, S2_hat] or [S1_hat, S2_hat, S3_hat].
    initial_guess : array-like, optional
        Initial guess for stress tensor [sigma_xx, sigma_yy, sigma_xy].
    update_interval : int, optional
        Update plot every N iterations. Default is 5.
    max_fringes : float, optional
        Maximum expected fringe order for setting bounds. Default is 6.

    Returns
    -------
    stress_tensor : ndarray
        Recovered stress tensor components [sigma_xx, sigma_yy, sigma_xy].
    success : bool
        Whether optimization was successful.
    history : dict
        Optimization history for further analysis.
    fig : matplotlib.figure.Figure
        The figure object (will be kept open).
    """
    import matplotlib.pyplot as plt

    from photoelastimetry.plotting import plot_optimization_history_live

    if initial_guess is None:
        initial_guess = np.array([1.0, 1.0, 0.0])

    # Storage for tracking
    history = {"stress_params": [], "S_predicted": [], "residuals": []}
    fig, axes = None, None
    iteration_count = [0]  # Use list to allow modification in nested function

    def callback_func(xk):
        """Callback with live plotting."""
        # Store current state
        history["stress_params"].append(xk.copy())

        S_pred_all = np.zeros((3, 2))
        for c in range(3):
            S_pred_all[c] = predict_stokes(xk[0], xk[1], xk[2], C_values[c], nu, L, wavelengths[c], S_i_hat)
        history["S_predicted"].append(S_pred_all.copy())

        residual = compute_residual(xk, S_m_hat, wavelengths, C_values, nu, L, S_i_hat)
        history["residuals"].append(residual)

        # Update plot periodically
        iteration_count[0] += 1
        if iteration_count[0] % update_interval == 0 or iteration_count[0] == 1:
            nonlocal fig, axes
            # Convert to arrays for plotting
            hist_for_plot = {
                "stress_params": np.array(history["stress_params"]),
                "S_predicted": np.array(history["S_predicted"]),
                "residuals": np.array(history["residuals"]),
            }
            fig, axes = plot_optimization_history_live(hist_for_plot, S_m_hat, fig, axes)
            plt.pause(0.01)

    # Run core optimization with live plotting callback
    result = _optimize_stress_tensor(
        S_m_hat, wavelengths, C_values, nu, L, S_i_hat, initial_guess, max_fringes, callback=callback_func
    )

    # Final update
    history["stress_params"] = np.array(history["stress_params"])
    history["S_predicted"] = np.array(history["S_predicted"])
    history["residuals"] = np.array(history["residuals"])

    fig, axes = plot_optimization_history_live(history, S_m_hat, fig, axes)
    plt.ioff()  # Turn off interactive mode

    return result.x, result.success, history, fig

compute_solid_fraction(S0, S_ref, mu, L)

Compute solid fraction from intensity using Beer-Lambert law.

Parameters:

Name	Type	Description	Default
S0	array - like	Measured intensity (from colour channel with absorptive dye).	required
S_ref	float	Reference light intensity before passing through sample.	required
mu	float	Absorption coefficient for the colour channel (calibrated parameter).	required
L	float	Sample thickness (m).	required

Returns:

Name	Type	Description
nu	array - like	Solid fraction values.

Source code in photoelastimetry/optimiser/stokes.py

def compute_solid_fraction(S0, S_ref, mu, L):
    """
    Compute solid fraction from intensity using Beer-Lambert law.

    Parameters
    ----------
    S0 : array-like
        Measured intensity (from colour channel with absorptive dye).
    S_ref : float
        Reference light intensity before passing through sample.
    mu : float
        Absorption coefficient for the colour channel (calibrated parameter).
    L : float
        Sample thickness (m).

    Returns
    -------
    nu : array-like
        Solid fraction values.
    """
    # Beer-Lambert: S0 = S_ref * exp(-mu * nu * L)
    # Solving for nu: nu = -ln(S0 / S_ref) / (mu * L)
    S0_safe = np.maximum(S0, 1e-10)
    nu = -np.log(S0_safe / S_ref) / (mu * L)
    return nu

recover_stress_map_stokes(image_stack, wavelengths, C_values, nu, L, S_i_hat, initial_guess_map=None, n_jobs=-1)

Recover full 2D stress tensor map from polarimetric image stack using Stokes method.

Parameters:

Name	Type	Description	Default
image_stack	ndarray	Image stack of shape [H, W, 3, 4] where: - H, W are image dimensions - 3 colour channels (R, G, B) - 4 polarisation angles (0, 45, 90, 135 degrees)	required
wavelengths	array - like	Wavelengths for R, G, B channels (m).	required
C_values	array - like	Stress-optic coefficients for R, G, B channels (1/Pa).	required
nu	float or ndarray	Solid fraction. Use 1.0 for solid samples. Can be scalar or array matching image dimensions.	required
L	float	Sample thickness (m).	required
S_i_hat	array - like	Incoming normalized Stokes vector [S1_hat, S2_hat] or [S1_hat, S2_hat, S3_hat].	required
initial_guess_map	ndarray	Initial guess stress map [H, W, 3].	None
n_jobs	int	Number of parallel jobs. -1 uses all available cores (default: -1).	-1

Returns:

Name	Type	Description
stress_map	ndarray	Array of shape [H, W, 3] containing [sigma_xx, sigma_yy, sigma_xy] in Pa.

Source code in photoelastimetry/optimiser/stokes.py

def recover_stress_map_stokes(
    image_stack,
    wavelengths,
    C_values,
    nu,
    L,
    S_i_hat,
    initial_guess_map=None,
    n_jobs=-1,
):
    """
    Recover full 2D stress tensor map from polarimetric image stack using Stokes method.

    Parameters
    ----------
    image_stack : ndarray
        Image stack of shape [H, W, 3, 4] where:
        - H, W are image dimensions
        - 3 colour channels (R, G, B)
        - 4 polarisation angles (0, 45, 90, 135 degrees)
    wavelengths : array-like
        Wavelengths for R, G, B channels (m).
    C_values : array-like
        Stress-optic coefficients for R, G, B channels (1/Pa).
    nu : float or ndarray
        Solid fraction. Use 1.0 for solid samples.
        Can be scalar or array matching image dimensions.
    L : float
        Sample thickness (m).
    S_i_hat : array-like
        Incoming normalized Stokes vector [S1_hat, S2_hat] or [S1_hat, S2_hat, S3_hat].
    initial_guess_map : ndarray, optional
        Initial guess stress map [H, W, 3].
    n_jobs : int, optional
        Number of parallel jobs. -1 uses all available cores (default: -1).

    Returns
    -------
    stress_map : ndarray
        Array of shape [H, W, 3] containing [sigma_xx, sigma_yy, sigma_xy] in Pa.
    """
    from joblib import Parallel, delayed

    H, W, _, _ = image_stack.shape
    stress_map = np.zeros((H, W, 3), dtype=np.float32)

    # Create list of all pixel coordinates
    pixel_coords = [(y, x) for y in range(H) for x in range(W)]

    # Create arguments for each pixel
    if initial_guess_map is not None:
        pixel_args = [
            (y, x, image_stack, wavelengths, C_values, nu, L, S_i_hat, initial_guess_map[y, x])
            for y, x in pixel_coords
        ]
    else:
        pixel_args = [(y, x, image_stack, wavelengths, C_values, nu, L, S_i_hat) for y, x in pixel_coords]

    # Process pixels in parallel
    results = Parallel(n_jobs=n_jobs)(
        delayed(_process_pixel)(args) for args in tqdm(pixel_args, desc="Processing pixels")
    )

    # Fill in the stress map
    for y, x, stress_tensor in results:
        stress_map[y, x, :] = stress_tensor

    return stress_map