Skip to content

3D PA simulation question #426

@Cuuuuu123

Description

@Cuuuuu123

Whenever I use the L-system to generate simulated blood vessels and then perform three-dimensional photoacoustic simulation with SIMPA, errors always occur.
Below are some basic codes that I wrote:
import os
import numpy as np
os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"

import tifffile
from scipy.ndimage import zoom

import simpa as sp
from simpa import Tags

TIFF_PATH = "./raw_data/15.tiff" # Path to your 3D TIFF (mask stack)
OUTPUT_DIR = "./output" # Where HDF5 and logs will be stored

input_spacing_mm = 1
target_spacing_mm = 1

DO_ACOUSTIC = True

RANDOM_SEED = 1234

path_manager = sp.PathManager()

def load_tiff_as_segmentation(path, assume_binary=True):
vol_zyx = tifffile.imread(path) # common shape: (Z, Y, X)
if vol_zyx.ndim != 3:
raise ValueError(f"Expected a 3D TIFF, got shape {vol_zyx.shape}")

seg_xyz = np.transpose(vol_zyx, (2, 1, 0))

if assume_binary and np.unique(seg_xyz).size > 2:
    # Heuristic threshold for binary: Otsu is available in skimage, but to reduce deps, use percentile
    thr = np.percentile(seg_xyz, 50)  # median threshold
    seg_xyz = (seg_xyz > thr).astype(np.uint8)
elif assume_binary and np.unique(seg_xyz).size <= 2:
    mn, mx = int(seg_xyz.min()), int(seg_xyz.max())
    if mn == mx:
        raise ValueError("Segmentation is constant; please check your TIFF data.")
    seg_xyz = (seg_xyz > mn).astype(np.uint8)
else:
    seg_xyz = seg_xyz.astype(np.int32)

return seg_xyz

def resample_segmentation_nearest(seg_xyz, in_spacing, out_spacing):
if in_spacing == out_spacing:
return seg_xyz
scale = in_spacing / out_spacing
# zoom takes per-axis zoom. Our seg is (X, Y, Z).
zoom_factors = (scale, scale, scale)
seg_resampled = sp.round_x5_away_from_zero(zoom(seg_xyz, zoom_factors, order=0)).astype(int)
return seg_resampled

def make_segmentation_class_mapping(unique_labels):
mapping = dict()
for label in unique_labels:
if label == 0:
mapping[label] = sp.TISSUE_LIBRARY.heavy_water()
else:
# Treat any non-zero as blood by default
mapping[label] = sp.TISSUE_LIBRARY.blood()
return mapping

def build_device(volume_transducer_dim_mm, volume_planar_dim_mm):
import numpy as _np

pitch_mm = 0.25
min_elements = 4
margin_mm = 0.5

CENTER_FREQ_HZ = 3.96e6
BANDWIDTH_PERCENT = 55

element_width_u_mm = 0.24
element_width_v_mm = 0.5

fov_half = max((volume_transducer_dim_mm / 2.0) - margin_mm, 1.0)
number_detector_elements = int(_np.floor((2 * fov_half) / pitch_mm)) + 1
number_detector_elements = max(number_detector_elements, min_elements)
number_detector_elements = min(number_detector_elements, 128)

device_pos = _np.array([volume_transducer_dim_mm / 2,
                        volume_planar_dim_mm / 2,
                        1.0])


z_extent = min(5.0, max(1.0, volume_transducer_dim_mm / 4.0))
fov = _np.asarray([-fov_half, fov_half, 0.0, 0.0, 0.0, z_extent])

device = sp.PhotoacousticDevice(
    device_position_mm=device_pos,
    field_of_view_extent_mm=fov,
)

detection_geom = sp.PlanarArrayDetectionGeometry(
    device_position_mm=device_pos,
    field_of_view_extent_mm=fov,
    center_frequency_hz=CENTER_FREQ_HZ,
    bandwidth_percent=BANDWIDTH_PERCENT,
)

device.set_detection_geometry(detection_geom)

slit_length = min(20.0, volume_transducer_dim_mm)
device.add_illumination_geometry(sp.SlitIlluminationGeometry(
    slit_vector_mm=[slit_length, 0.0, 0.0]
))

print(f"[DEBUG] Planar array (as Linear) with {number_detector_elements}x1 elements created.")
print(f"[DEBUG]   Acoustic props: F_center={CENTER_FREQ_HZ / 1e6} MHz, BW={BANDWIDTH_PERCENT}%")
print(f"[DEBUG]   Element size: {element_width_u_mm} mm (U) x {element_width_v_mm} mm (V)")

return device

def run_from_tiff():
os.makedirs(OUTPUT_DIR, exist_ok=True)

# 1) Load TIFF -> (X,Y,Z) labels
seg_xyz = load_tiff_as_segmentation(TIFF_PATH, assume_binary=True)

# 2) Optional resampling to target spacing
seg_xyz = resample_segmentation_nearest(seg_xyz, input_spacing_mm, target_spacing_mm)

# 3) SIMPA settings following official examples
np.random.seed(RANDOM_SEED)

# Dimensions in mm are derived from voxel count * spacing (X,Y,Z)
dim_x_mm = seg_xyz.shape[0] * target_spacing_mm
dim_y_mm = seg_xyz.shape[1] * target_spacing_mm
dim_z_mm = seg_xyz.shape[2] * target_spacing_mm

# General/global settings
general_settings = {
    Tags.RANDOM_SEED: RANDOM_SEED,
    Tags.VOLUME_NAME: "TIFFSegmentation_" + str(RANDOM_SEED),
    Tags.SIMULATION_PATH: path_manager.get_hdf5_file_save_path() if hasattr(path_manager, "get_hdf5_file_save_path") else OUTPUT_DIR,
    Tags.SPACING_MM: target_spacing_mm,
    Tags.DIM_VOLUME_X_MM: dim_x_mm,
    Tags.DIM_VOLUME_Y_MM: dim_y_mm,
    Tags.DIM_VOLUME_Z_MM: dim_z_mm,
    Tags.GPU: True,
    Tags.WAVELENGTHS: [800],
    Tags.DO_FILE_COMPRESSION: True,
    # You can enable IPASC export if desired
    Tags.DO_IPASC_EXPORT: True,
}
settings = sp.Settings(general_settings)

# Volume creation settings via segmentation-based adapter
settings.set_volume_creation_settings({
    Tags.INPUT_SEGMENTATION_VOLUME: seg_xyz,
    Tags.SEGMENTATION_CLASS_MAPPING: make_segmentation_class_mapping(np.unique(seg_xyz)),
    # You can also set: Tags.VOLUME_CREATOR: Tags.VOLUME_CREATOR_SEGMENTATION_BASED,
})

# Optical settings (MCX)
settings.set_optical_settings({
    Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 1e7,
    Tags.OPTICAL_MODEL_BINARY_PATH: path_manager.get_mcx_binary_path(),
    # Example illumination (MSOT Acuity Echo style)
    Tags.ILLUMINATION_TYPE: Tags.ILLUMINATION_TYPE_MSOT_ACUITY_ECHO,
    Tags.LASER_PULSE_ENERGY_IN_MILLIJOULE: 50,
    Tags.MCX_ASSUMED_ANISOTROPY: 0.9,
    # Optional: print GPU info
    # Tags.ADDITIONAL_FLAGS: ['--printgpu']
})

# Acoustic settings (k-Wave) - only used if DO_ACOUSTIC
if DO_ACOUSTIC:
    settings.set_acoustic_settings({
        Tags.ACOUSTIC_SIMULATION_3D: True,  # set True for 3D k-Wave (much slower)
        Tags.ACOUSTIC_MODEL_BINARY_PATH: path_manager.get_matlab_binary_path(),
        Tags.KWAVE_PROPERTY_ALPHA_POWER: 0.00,
        Tags.KWAVE_PROPERTY_SENSOR_RECORD: "p",
        Tags.KWAVE_PROPERTY_PMLInside: False,
        Tags.KWAVE_PROPERTY_PMLSize: [31, 32],
        Tags.KWAVE_PROPERTY_PMLAlpha: 1.5,
        Tags.KWAVE_PROPERTY_PlotPML: False,
        Tags.RECORDMOVIE: False,
        Tags.MOVIENAME: "kwave_log",
        Tags.ACOUSTIC_LOG_SCALE: True,
    })

    # Reconstruction settings (Time Reversal) similar to example
    settings.set_reconstruction_settings({
        Tags.RECONSTRUCTION_PERFORM_BANDPASS_FILTERING: False,
        Tags.ACOUSTIC_MODEL_BINARY_PATH: path_manager.get_matlab_binary_path(),
        Tags.ACOUSTIC_SIMULATION_3D: True,
        Tags.KWAVE_PROPERTY_ALPHA_POWER: 0.00,
        Tags.TUKEY_WINDOW_ALPHA: 0.5,
        Tags.BANDPASS_CUTOFF_LOWPASS_IN_HZ: int(8e6),
        Tags.BANDPASS_CUTOFF_HIGHPASS_IN_HZ: int(0.1e4),
        Tags.RECONSTRUCTION_BMODE_AFTER_RECONSTRUCTION: False,
        Tags.RECONSTRUCTION_BMODE_METHOD: Tags.RECONSTRUCTION_BMODE_METHOD_HILBERT_TRANSFORM,
        Tags.RECONSTRUCTION_APODIZATION_METHOD: Tags.RECONSTRUCTION_APODIZATION_BOX,
        Tags.RECONSTRUCTION_MODE: Tags.RECONSTRUCTION_MODE_PRESSURE,
        Tags.KWAVE_PROPERTY_SENSOR_RECORD: "p",
        Tags.KWAVE_PROPERTY_PMLInside: False,
        Tags.KWAVE_PROPERTY_PMLSize: [31, 32],
        Tags.KWAVE_PROPERTY_PMLAlpha: 1.5,
        Tags.KWAVE_PROPERTY_PlotPML: False,
        Tags.RECORDMOVIE: False,
        Tags.MOVIENAME: "recon_log",
        Tags.ACOUSTIC_LOG_SCALE: True,
        Tags.DATA_FIELD_SPEED_OF_SOUND: 1540,
        Tags.DATA_FIELD_ALPHA_COEFF: 0.01,
        Tags.DATA_FIELD_DENSITY: 1000,
        Tags.SPACING_MM: target_spacing_mm
    })

# Noise settings (optional, like in example)
settings["noise_initial_pressure"] = {
    Tags.NOISE_MEAN: 1,
    Tags.NOISE_STD: 0.01,
    Tags.NOISE_MODE: Tags.NOISE_MODE_MULTIPLICATIVE,
    Tags.DATA_FIELD: Tags.DATA_FIELD_INITIAL_PRESSURE,
    Tags.NOISE_NON_NEGATIVITY_CONSTRAINT: True
}
settings["noise_time_series"] = {
    Tags.NOISE_STD: 1,
    Tags.NOISE_MODE: Tags.NOISE_MODE_ADDITIVE,
    Tags.DATA_FIELD: Tags.DATA_FIELD_TIME_SERIES_DATA
}

# Build a device (geometry + illumination)
device = build_device(volume_transducer_dim_mm=dim_x_mm, volume_planar_dim_mm=dim_y_mm)

if DO_ACOUSTIC:
    pipeline = [
        sp.SegmentationBasedAdapter(settings),
        sp.MCXAdapter(settings),
        sp.GaussianNoise(settings, "noise_initial_pressure"),
        sp.KWaveAdapter(settings),
        sp.GaussianNoise(settings, "noise_time_series"),
        sp.TimeReversalAdapter(settings),
        sp.FieldOfViewCropping(settings)
    ]
else:
    pipeline = [
        sp.SegmentationBasedAdapter(settings),
        sp.MCXAdapter(settings)
    ]

# Run simulation
sp.simulate(pipeline, settings, device)

# Visualisation
# Determine default wavelength to display
if Tags.WAVELENGTH in settings:
    wl = settings[Tags.WAVELENGTH]
else:
    wl = settings[Tags.WAVELENGTHS][0]

# Show results: segmentation map + initial pressure; if acoustic on, also show time series & reconstruction
sp.visualise_data(
    path_to_hdf5_file=settings[Tags.SIMPA_OUTPUT_FILE_PATH],
    wavelength=wl,
    show_segmentation_map=True,
    show_initial_pressure=True,
    show_time_series_data=DO_ACOUSTIC,
    show_reconstructed_data=DO_ACOUSTIC,
    log_scale=False,
    show_xz_only=False
)

if name == "main":
print("[SIMPA] Starting simulation from TIFF segmentation...")
print(f" TIFF_PATH = {TIFF_PATH}")
print(f" input_spacing_mm = {input_spacing_mm} -> target_spacing_mm = {target_spacing_mm}")
print(f" DO_ACOUSTIC = {DO_ACOUSTIC}")
run_from_tiff()
print("[SIMPA] Done.")

Could you give me some example to show how the code build?

Metadata

Metadata

Assignees

No one assigned

    Labels

    Help NeededThe user has a question they need feedback on and help with

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions