MCXAdapterReflectance does not use anisotropy for simulation

[RecurvedBow]

Describe the bug Despite the anisotropy parameter defined in the MCXAdapterReflectance.forward_model(...) method, it is unused. Instead, the parameter _assumed_anisotropy is used for the simulations (optical_forward_model_mcx_reflectance_adapter.py, line 71). Therefore, setting the anisotropy does not affect the simulation.

Specify a priority (low, medium, high) Medium

Potential solution To solve this, remove the method override MCXAdapterReflectance.forward_model(...) as the base class MCXAdapter contains the same implementation, but with the bug already fixed.

To Reproduce Run this script:

from simpa import Tags
import simpa as sp
import numpy as np

# FIXME temporary workaround for newest Intel architectures
import os
os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"

path_manager = sp.PathManager()

VOLUME_TRANSDUCER_DIM_IN_MM = 10
VOLUME_PLANAR_DIM_IN_MM = 10
VOLUME_HEIGHT_IN_MM = 12
SPACING = 0.5
RANDOM_SEED = 471
VOLUME_NAME = "MyVolumeName_"+str(RANDOM_SEED)
SAVE_REFLECTANCE = True
SAVE_PHOTON_DIRECTION = False

# If VISUALIZE is set to True, the simulation result will be plotted
VISUALIZE = True

def create_example_tissue(anisotropy: float):
    background_dictionary = sp.Settings()
    background_dictionary[Tags.MOLECULE_COMPOSITION] = sp.TISSUE_LIBRARY.constant(1e-4, 1e-4, 0.9)
    background_dictionary[Tags.STRUCTURE_TYPE] = Tags.BACKGROUND

    tissue_dict = sp.Settings()
    tissue_dict[Tags.BACKGROUND] = background_dictionary

    tissue_dict["tissue"] = sp.define_horizontal_layer_structure_settings(
            molecular_composition=sp.TISSUE_LIBRARY.constant(1e-6, 1e-1, anisotropy),
            z_start_mm=0,
            thickness_mm=VOLUME_HEIGHT_IN_MM)

    return tissue_dict

np.random.seed(RANDOM_SEED)

general_settings = {
    Tags.RANDOM_SEED: RANDOM_SEED,
    Tags.VOLUME_NAME: VOLUME_NAME,
    Tags.SIMULATION_PATH: path_manager.get_hdf5_file_save_path(),
    Tags.SPACING_MM: SPACING,
    Tags.DIM_VOLUME_Z_MM: VOLUME_HEIGHT_IN_MM,
    Tags.DIM_VOLUME_X_MM: VOLUME_TRANSDUCER_DIM_IN_MM,
    Tags.DIM_VOLUME_Y_MM: VOLUME_PLANAR_DIM_IN_MM,
    Tags.WAVELENGTHS: np.arange(300, 800 + 1, 50),
    Tags.DO_FILE_COMPRESSION: True
}

settings = sp.Settings(general_settings)

settings.set_optical_settings({
    Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 1e7,
    Tags.OPTICAL_MODEL_BINARY_PATH: path_manager.get_mcx_binary_path(),
    Tags.COMPUTE_DIFFUSE_REFLECTANCE: SAVE_REFLECTANCE,
    Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT: SAVE_PHOTON_DIRECTION
})

anisotropies = [0.7, 0.9, 1.1]
reflectance_dicts = []
reflectance_pos_dicts = []
for anisotropy in anisotropies:
    settings.set_volume_creation_settings({
        Tags.STRUCTURES: create_example_tissue(anisotropy)
    })

    pipeline = [
        sp.ModelBasedVolumeCreationAdapter(settings),
        sp.MCXAdapterReflectance(settings),
    ]

    device = sp.PencilBeamIlluminationGeometry(device_position_mm=np.asarray([VOLUME_TRANSDUCER_DIM_IN_MM/2,
                                                                              VOLUME_PLANAR_DIM_IN_MM/2, 0]))

    sp.simulate(pipeline, settings, device)

    refl_by_wavelength = sp.load_data_field(path_manager.get_hdf5_file_save_path() + "/" + VOLUME_NAME + ".hdf5",
                                            sp.Tags.DATA_FIELD_DIFFUSE_REFLECTANCE)
    refl_pos_by_wavelength = sp.load_data_field(path_manager.get_hdf5_file_save_path() + "/" + VOLUME_NAME + ".hdf5",
                                                sp.Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS)
    reflectance_dicts.append(refl_by_wavelength)
    reflectance_pos_dicts.append(refl_pos_by_wavelength)

print(len(reflectance_dicts))
print(reflectance_dicts[0]["300"])

for refl_by_wavelength in reflectance_dicts:
    assert len(set(refl_by_wavelength.keys()) & set(reflectance_dicts[0].keys())) == len(refl_by_wavelength.keys()) == len(reflectance_dicts[0].keys())

    for key, value in refl_by_wavelength.items():
        assert np.max(np.abs(value - reflectance_dicts[0][key])) <= 1e-8

for refl_pos_by_wavelength in reflectance_pos_dicts:
    assert len(set(refl_pos_by_wavelength.keys()) & set(reflectance_pos_dicts[0].keys())) == len(refl_pos_by_wavelength.keys()) == len(
        reflectance_pos_dicts[0].keys())

    for key, value in refl_pos_by_wavelength.items():
        assert (value == reflectance_pos_dicts[0][key]).all()

Current Behavior The simulated reflectance values are the same.

Expected behavior The simulated reflectance values should be different.

[leoyala]

@kdreher do you think that this is because MCX does not support the use of different g values per voxel? I just noticed that there are two different anisotropies used in the base class, do you know why? One comes from the component settings and the other from the forward model parameter. What is the difference between _assummed_anisotropy an anisotropy?

https://github.com/IMSY-DKFZ/simpa/blob/22b4c8e223dd4d5ad32ecccc5e583ff59fdae67f/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_adapter.py#L59-L71

[jgroehl]

Exactly! Because MCX used to only support a global anisotropy, but we defined a per-voxel anisotropy in SIMPA, we had also to define a separate anisotropy that MCX would use. The naming convention is also not ideal here... We convert all scattering coefficients during the MCX volume generation step with the two anisotropies:

def volumes_to_mm(**kwargs) -> Tuple:
    """
    transforms volumes into `mm` units

    :param kwargs: dictionary containing at least the keys `scattering_cm, absorption_cm, anisotropy` and
        `assumed_anisotropy`
    :return: `Tuple` of volumes after transformation
    """
    scattering_cm = kwargs.get('scattering_cm')
    absorption_cm = kwargs.get('absorption_cm')
    absorption_mm = absorption_cm / 10
    scattering_mm = scattering_cm / 10

    # FIXME Currently, mcx only accepts a single value for the anisotropy.
    #   In order to use the correct reduced scattering coefficient throughout the simulation,
    #   we adjust the scattering parameter to be more accurate in the diffuse regime.
    #   This will lead to errors, especially in the quasi-ballistic regime.

    given_reduced_scattering = (scattering_mm * (1 - kwargs.get('anisotropy')))

    # If the anisotropy is 1, all scattering is forward scattering which is equal to no scattering at all
    if kwargs.get("assumed_anisotropy") == 1:
        scattering_mm = given_reduced_scattering * 0
    else:
        scattering_mm = given_reduced_scattering / (1 - kwargs.get('assumed_anisotropy'))
    scattering_mm[scattering_mm < 1e-10] = 1e-10
    return absorption_mm, scattering_mm

There have been multiple bugfixes in the past showing that this practice is not adequate. Since MCX now supports mua/mus/g/n volumes, it might be a good idea to get rid of this and swap completely.

[leoyala]

In theory, the diffuse reflectance adapter should also support now the per-voxel anisotropy definitions. So I would say it is better to only accept now that option as default and remove all traces of the old anisotropy implementation. Which should be the default @jgroehl? The anisotropy parsed to the forward_model method?

[jgroehl]

Yeah, I would agree with you there! The "mcx assumed anisotropy" is artificially introduced because of a legacy shortcoming of the downstream tool. As you said, we should eliminate it from the pipeline and only use the anisotropy provided by the preceding step:

Lines 41-45: simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_adapter

def forward_model(self,
                  absorption_cm: np.ndarray,
                  scattering_cm: np.ndarray,
                  anisotropy: np.ndarray,
                  illumination_geometry: IlluminationGeometryBase) -> Dict:

[leoyala]

ok, will work on it.

[leoyala]

Is the feature to support anisotropy per voxel already in the develop branch? From what I can see, the scattering is still being transformed to reduced scattering, and then the volume is defined only based on scattering and absorption.

[lkeegan]

@leoyala (cc @kdreher) there is an initial implementation of this here:

https://github.com/lkeegan/simpa/tree/gn_mcx_based_on_T265 https://github.com/IMSY-DKFZ/simpa/compare/main...lkeegan:simpa:gn_mcx_based_on_T265

However I didn't make a PR yet because these changes break the reflectance, not yet clear if this is a bug in this code or in mcx:

• to enable refraction between voxels in mcx we need to set -b 1
  • but this also results in refraction at the external boundaries
• to avoid refracting at the external boundaries we set -bc aaaaaa000010
  • but then the minimal_optical_simulation with SAVE_REFLECTANCE = True no longer shows anything for the reflectance

[leoyala]

@lkeegan this actually sounds like a bug that we should probably raise to QianQian if it does turn out that the problem comes from mcx. This issue will have to wait until your PR is integrated then.