Note
Go to the end to download the full example code.
Comparisons of block b closed for dicom and vtu images.
from pathlib import Path
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import numpy as np
import darsia
# ! ---- Constants
cm2m = 1e-2 # conversion cm -> m
h2s = 60**2 # conversion hours -> seconds
ml_per_hour_to_m3_per_s = cm2m**3 / h2s
# ! ---- Model paramaters
porosity_2d = 0.2321
fracture_aperture = 0.1 * cm2m
depth = 1.95 * cm2m
injection_rate = 15 * ml_per_hour_to_m3_per_s
# ! ---- Read DICOM images
def read_dicom_images() -> darsia.Image:
"""Read space-time image from DICOM format for fractip a."""
# Provide folder with dicom images (note: two folders)
root = Path(
"/home/jakub/images/ift/fracflow/01-single-phase-tracer-in-fracture/fractip-b"
)
images = [
root / Path(p) / Path("DICOM/PA1/ST1/SE1")
for p in [
"fractip b closed 1 min rekon start 0 15 frames",
]
]
if False:
# Stack images (need to use such approach due to lack of common reference
dicom_images = []
for img in images:
dicom_image = darsia.imread(img, suffix=".dcm", dim=3, series=True)
dicom_images.append(dicom_image)
uncorrected_image_3d = darsia.stack(dicom_images)
# Use an assistant to tune the rotation correction and subregion selection
test_image = uncorrected_image_3d.time_slice(9)
test_image.show(
surpress_3d=True, mode="matplotlib", side_view="voxel", cmap="turbo"
)
test_image.show(
mode="plotly",
side_view="voxel",
view="voxel",
cmap="turbo",
threshold=0.8,
relative=True,
)
options = {
"threshold": 0.05,
"relative": True,
"verbosity": True,
}
rotation_assistant = darsia.RotationCorrectionAssistant(test_image, **options)
rotation_corrections = rotation_assistant()
rotated_test_image = test_image.copy()
for rotation in rotation_corrections:
rotated_test_image = rotation(rotated_test_image)
subregion_assistant = darsia.SubregionAssistant(rotated_test_image, **options)
coordinates = subregion_assistant()
else:
# Alternative: Fix tailored input parameters for the correction objects.
rotation_corrections = [
darsia.RotationCorrection(
**{
"anchor": np.array([91, 106, 13]),
"rotation_from_isometry": True,
"pts_src": np.array([[91, 106, 13], [97, 106, 86]]),
"pts_dst": np.array([[91, 106, 13], [91, 106, 86]]),
}
),
darsia.RotationCorrection(
**{
"anchor": np.array([50, 96, 117]),
"rotation_from_isometry": True,
"pts_src": np.array([[50, 96, 117], [90, 97, 117]]),
"pts_dst": np.array([[50, 96, 117], [90, 96, 117]]),
}
),
darsia.RotationCorrection(
**{
"anchor": np.array([106, 94, 0]),
"rotation_from_isometry": True,
"pts_src": np.array([[106, 94, 0], [106, 96, 70]]),
"pts_dst": np.array([[106, 94, 0], [106, 95, 69]]),
}
),
]
coordinates = np.array(
[
[42.21969627956989, 42.29322659139785, -195.13738119462366],
[42.30791713978495, 42.36046250537635, -195.1190467860215],
]
)
# Re-read the dicom image.
# Use the uncorrected images further.
image_3d = darsia.imread(
images,
suffix=".dcm",
dim=3,
series=True,
)
image_3d = image_3d.time_slice(9)
for correction in rotation_corrections:
image_3d = correction(image_3d)
image_3d = image_3d.subregion(coordinates)
# Plot side view to check the result
if True:
image_3d.show(
mode="matplotlib",
surpress_3d=True,
side_view="voxel",
threshold=0.05,
relative=True,
)
# ! ---- Precondition
image_3d.img /= np.max(image_3d.img)
return image_3d
def plot_3d_image(image_3d):
# Plot side view to check the result
image_3d.show(
mode="matplotlib",
surpress_3d=True,
side_view="voxel",
threshold=0.05,
relative=True,
)
def plot_sum(image_3d):
if not isinstance(image_3d, list):
image_3d = [image_3d]
max_val = max(
[
np.max(
np.sum(img, axis=0)
if isinstance(img, np.ndarray)
else np.sum(img.img, axis=0)
)
for img in image_3d
]
)
for i, img in enumerate(image_3d):
plt.figure(f"sum {i}")
if isinstance(img, np.ndarray):
plt.imshow(np.sum(img, axis=0), cmap="turbo", vmin=0, vmax=max_val)
else:
plt.imshow(np.sum(img.img, axis=0), cmap="turbo", vmin=0, vmax=max_val)
plt.colorbar()
plt.show()
def plot_slice(image_3d):
if not isinstance(image_3d, list):
image_3d = [image_3d]
max_val = max(
[np.max(img if isinstance(img, np.ndarray) else img.img) for img in image_3d]
)
for i, img in enumerate(image_3d):
plt.figure(f"slice {i}")
if isinstance(img, np.ndarray):
plt.imshow(img[20], cmap="turbo", vmin=0, vmax=max_val)
else:
plt.imshow(img.img[20], cmap="turbo", vmin=0, vmax=max_val)
plt.colorbar()
plt.show()
# Enrich the signal under the assumption of a monotone displacement experiment.
def accumulate_signal(image: darsia.Image) -> darsia.Image:
new_image = image.copy()
for time_id in range(image.time_num):
new_image.img[..., time_id] = np.max(
image.img[..., slice(0, time_id + 1)], axis=-1
)
return new_image
# ! ---- Apply regularization
def apply_tvd(image: darsia.Image) -> darsia.Image:
"""NOTE: Most time-consuming routine."""
pass
def read_regularized_images(image: darsia.Image) -> darsia.Image:
"""Read preprocessed regularized images from numpy format."""
reg_image = image.copy()
for i in range(0, 13):
reg_image.img[..., i] = np.load(
# f"tvd/block-a/heterogeneous_mu_1000_iter/tvd_a_1000_{i}.npy"
# f"tvd/block-a/mu_0.0001_1000_iter/tvd_{i}.npy" #
f"tvd/delta/block-a/mu_0.0001_10000_iter/tvd_{i}.npy" # latest tvd
)
return reg_image
# ! ---- Reduce 3d image to two spatial dimensions
def dimensional_reduction(image_3d: darsia.Image) -> darsia.Image:
"""Flatten 3d image."""
return darsia.reduce_axis(image_3d, axis="z", mode="average")
# ! ---- Rescale images to concentrations
def extract_concentration(
image_2d: darsia.Image, calibration_interval, injection_rate: float
) -> darsia.Image:
# Define default concentration analysis, enabling calibration
class CalibratedConcentrationAnalysis(
darsia.ConcentrationAnalysis, darsia.InjectionRateModelObjectiveMixin
):
pass
model = darsia.ScalingModel(scaling=1)
dicom_concentration_analysis = CalibratedConcentrationAnalysis(
base=None, model=model
)
# Calibrate concentration analysis
shape_meta = image_2d.shape_metadata()
geometry = darsia.ExtrudedGeometry(expansion=depth, **shape_meta)
dicom_concentration_analysis.calibrate_model(
images=image_2d.time_interval(calibration_interval),
options={
"geometry": geometry,
"injection_rate": injection_rate,
"initial_guess": [1.0],
"tol": 1e-5,
"maxiter": 100,
"method": "Nelder-Mead",
"regression_type": "ransac",
},
plot_result=True,
)
# Interpret the calibration results to effectively determine the injection start.
reference_time = dicom_concentration_analysis.model_calibration_postanalysis()
# The results is that the effective injection start occurs at 3.25 minutes. Thus, update
# the reference time, or simply the relative time. Set reference time with respect to
# the first active image. Use seconds.
image_2d.update_reference_time(reference_time)
# Print some of the specs
if True:
print(
f"The identified reference time / injection start is: {reference_time} [s]"
)
print(
f"The dimensions of the space time dicom image are: {image_2d.dimensions}"
)
print(f"Relative times in seconds: {image_2d.time}")
# Only for plotting reasons - same as above but with updated times and with activated plot:
if False:
dicom_concentration_analysis.calibrate_model(
images=image_2d.time_interval(
calibration_interval
), # Only use some of the first images due to cut-off in the geometry
options={
"geometry": geometry,
"injection_rate": injection_rate,
"initial_guess": [1.0],
"tol": 1e-5,
"maxiter": 100,
"method": "Nelder-Mead",
"regression_type": "ransac",
},
plot_result=True,
)
# Produce space time image respresentation of concentrations.
dicom_concentration = dicom_concentration_analysis(image_2d)
# Debug: Plot the concentration profile. Note that values are larger than 1, which is
# due to the calibration based on a non-smooth signal. The expected volume has to be
# distributed among the active pixels. Pixels with strong activity thereby hold
# unphysical amounts of fluid.
if False:
dicom_concentration.show("dicom concentration", 5)
return dicom_concentration
# Determine total concentrations over time
def plot_concentration_evolution(concentration_2d):
shape_meta = concentration_2d.shape_metadata()
geometry = darsia.ExtrudedGeometry(expansion=depth, **shape_meta)
concentration_values = geometry.integrate(concentration_2d)
plt.figure("Experiments - injected volume over time")
plt.plot(concentration_2d.time, concentration_values)
plt.plot(
concentration_2d.time,
[injection_rate * time for time in concentration_2d.time],
color="black",
linestyle=(0, (5, 5)),
)
plt.xlabel("time [s]")
plt.ylabel("volume [m**3]")
plt.show()
# ! ---- VTU images
def read_vtu_images(vtu_ind: int) -> darsia.Image:
"""Read-in all available VTU data."""
# This corresponds approx. to (only temporary - not used further):
vtu_time = 0 # not relevant here.
# Find the corresponding files
vtu_images_2d = Path(f"data/vtu/block-b/data_2_{str(vtu_ind).zfill(6)}.vtu")
vtu_images_1d = Path(f"data/vtu/block-b/data_1_{str(vtu_ind).zfill(6)}.vtu")
# Define vtu images as DarSIA images
vtu_image_2d = darsia.imread(
vtu_images_2d,
time=vtu_time, # relative time in minutes
key="temperature", # key to address concentration data
shape=(400, 400), # dim = 2
vtu_dim=2, # effective dimension of vtu data
)
vtu_image_1d = darsia.imread(
vtu_images_1d,
time=vtu_time, # relative time in minutes
key="temperature", # key to address concentration data
shape=(1001, 51), # dim = 2
vtu_dim=1, # effective dimension of vtu data
width=fracture_aperture, # for embedding in 2d
)
# PorePy/meshio cuts off coordinate values at some point...
# Correct manually - width / 2 - aperature / 2.
vtu_image_1d.origin[0] = (6.98 / 2 - 0.1 / 2) * cm2m
# Equidimensional reconstructions. Superpose 2d and 1d images.
# And directly interpret
porosity_1d = 1.0 - porosity_2d # for volume conservation
vtu_image = darsia.superpose(
[
darsia.weight(vtu_image_2d, porosity_2d),
darsia.weight(vtu_image_1d, porosity_1d),
]
)
# Plot for testing purposes
if False:
vtu_image.show("equi-dimensionsional reconstruction")
# Concentrations - simple since the equi-dimensional vtu images are interpreted as
# volumetric concentration
vtu_concentration = vtu_image.copy()
return vtu_concentration
# ! --- Align DICOM and vtu images
def align_images(dicom_concentration, vtu_concentration):
# Plot two exemplary images to identify suitable src and dst points which will define
# the alignment procedure
if True:
plt.figure("dicom")
plt.imshow(np.sum(dicom_concentration.img, axis=0))
plt.figure("vtu")
plt.imshow(np.sum(vtu_concentration.img, axis=0))
plt.show()
# Pixels of fracture end points
voxels_src = [[0, 85, 0], [0, 85, 86], [17, 85, 0], [17, 85, 86]] # in DICOM image
voxels_dst = [
[0, 237, 86.5],
[0, 148, 86.5],
[17, 237, 86.5],
[17, 148, 86.5],
] # in VTU image
# Define coordinate transform and apply it
coordinatesystem_src = dicom_concentration.coordinatesystem
coordinatesystem_dst = vtu_concentration.coordinatesystem
coordinate_transformation = darsia.CoordinateTransformation(
coordinatesystem_src,
coordinatesystem_dst,
voxels_src,
voxels_dst,
fit_options={
"tol": 1e-5,
"preconditioning": True,
},
isometry=False,
)
transformed_dicom_concentration = coordinate_transformation(dicom_concentration)
# Restrict to intersecting active canvas
intersection = coordinate_transformation.find_intersection()
aligned_dicom_concentration = transformed_dicom_concentration.subregion(
intersection
)
aligned_vtu_concentration = vtu_concentration.subregion(intersection)
return aligned_dicom_concentration, aligned_vtu_concentration
def single_plot(
mode: str,
full_image: darsia.Image,
add_on: str,
image_path: Path,
):
##############################################################################
# Plot the reconstructed vtu data, vtu plus Gaussian noise, and the dicom data.
if mode == "sum":
image = darsia.reduce_axis(full_image, axis="z", mode="average")
else:
image = darsia.reduce_axis(full_image, axis="z", mode="slice", slice_idx=20)
# Define some plotting options (font, x and ylabels)
# matplotlib.rcParams.update({"font.size": 14})
# Define x and y labels in cm (have to convert from m to cm)
shape = image.num_voxels
dimensions = image.dimensions
x_pixel, y_pixel = np.meshgrid(
np.linspace(dimensions[0], 0, shape[0]),
np.linspace(0, dimensions[1], shape[1]),
indexing="ij",
)
vmax = 1.25
vmin = 0
cmap = "turbo"
# Plot the data with contour line.
fig_combination, axs_combination = plt.subplots(nrows=1, ncols=1)
axs_combination.pcolormesh(
y_pixel, x_pixel, image.img, cmap=cmap, vmin=vmin, vmax=vmax
)
axs_combination.text(
0.005, 0.075, "experiment", color="white", alpha=0.5, rotation=0, fontsize=14
)
axs_combination.text(
0.005, 0.070, add_on, color="white", alpha=0.5, rotation=0, fontsize=14
)
axs_combination.set_ylim(top=0.08)
axs_combination.set_aspect("equal")
axs_combination.set_xlabel("x [m]", fontsize=14)
axs_combination.set_ylabel("y [m]", fontsize=14)
fig_combination.colorbar(
cm.ScalarMappable(cmap=cmap),
ax=axs_combination,
label="volumetric concentration",
)
fig_combination.savefig(image_path, dpi=500, transparent=True)
plt.show()
def qualitative_comparison(
mode: str,
full_dicom_image: darsia.Image,
full_vtu_image: darsia.Image,
add_on: str,
image_path: Path,
colorbar: bool = False,
):
##############################################################################
# Plot the reconstructed vtu data, vtu plus Gaussian noise, and the dicom data.
if mode == "sum":
dicom_image = darsia.reduce_axis(full_dicom_image, axis="z", mode="average")
vtu_image = darsia.reduce_axis(full_vtu_image, axis="z", mode="average")
else:
dicom_image = darsia.reduce_axis(
full_dicom_image, axis="z", mode="slice", slice_idx=20
)
vtu_image = darsia.reduce_axis(
full_vtu_image, axis="z", mode="slice", slice_idx=20
)
# Define some plotting options (font, x and ylabels)
# matplotlib.rcParams.update({"font.size": 14})
# Define x and y labels in cm (have to convert from m to cm)
shape = dicom_image.num_voxels
dimensions = dicom_image.dimensions
x_pixel, y_pixel = np.meshgrid(
np.linspace(dimensions[0], 0, shape[0]),
np.linspace(0, dimensions[1], shape[1]),
indexing="ij",
)
vmax = 1.25
vmin = 0
contourlevels = [0.045 * vmax, 0.055 * vmax]
cmap = "turbo"
# Plot the dicom data with contour line.
fig_combination, axs_combination = plt.subplots(nrows=1, ncols=1)
combined_image = dicom_image.img.copy()
mid = 86
combined_image[:, mid:] = vtu_image.img[:, mid:]
axs_combination.pcolormesh(
y_pixel, x_pixel, combined_image, cmap=cmap, vmin=vmin, vmax=vmax
)
axs_combination.contourf(
y_pixel, x_pixel, vtu_image.img, cmap="Reds", levels=contourlevels, alpha=0.5
)
axs_combination.plot(
[3.45 * cm2m, 3.45 * cm2m],
[0.0, 0.08],
color="white",
alpha=0.3,
linestyle="dashed",
)
axs_combination.text(
0.005, 0.075, "experiment", color="white", alpha=0.5, rotation=0, fontsize=14
)
axs_combination.text(
0.005, 0.070, add_on, color="white", alpha=0.5, rotation=0, fontsize=14
)
axs_combination.text(
0.04, 0.075, "simulation", color="white", alpha=0.5, rotation=0, fontsize=14
)
axs_combination.set_ylim(top=0.08)
axs_combination.set_aspect("equal")
axs_combination.set_xlabel("x [m]", fontsize=14)
axs_combination.set_ylabel("y [m]", fontsize=14)
if colorbar:
fig_combination.colorbar(
cm.ScalarMappable(cmap=cmap),
ax=axs_combination,
label="volumetric concentration",
)
fig_combination.savefig(image_path, dpi=500, transparent=True)
plt.show()
Main analysis
# Original PET images
if True:
dicom_image_3d = read_dicom_images()
dicom_image_3d.save("data/npz/block-b/dicom_raw_3d.npz")
else:
dicom_image_3d = darsia.imread("data/npz/block-b/dicom_raw_3d.npz")
# Pick corresponding vtu images.
vtu_ind = 439
vtu_2d_concentration = read_vtu_images(vtu_ind)
# Resize to similar shape as dicom image
dicom_voxel_size = dicom_image_3d.voxel_size
vtu_2d_concentration = darsia.equalize_voxel_size(
vtu_2d_concentration, min(dicom_voxel_size)
)
# Expand vtu image to 3d
dicom_height = dicom_image_3d.dimensions[0]
dicom_shape = dicom_image_3d.img.shape
vtu_concentration_3d = darsia.extrude_along_axis(
vtu_2d_concentration, dicom_height, dicom_shape[0]
)
# Align dicom and vtu
if True:
dicom_concentration = dicom_image_3d.copy()
aligned_dicom_concentration, aligned_vtu_concentration = align_images(
dicom_concentration, vtu_concentration_3d
)
aligned_dicom_concentration.save("data/npz/block-b/aligned_dicom_concentration.npz")
aligned_vtu_concentration.save(
f"data/npz/block-b/aligned_vtu_{vtu_ind}_concentration.npz"
)
else:
aligned_dicom_concentration = darsia.imread(
"data/npz/block-b/aligned_dicom_concentration.npz"
)
aligned_vtu_concentration = darsia.imread(
f"data/npz/block-b/aligned_vtu_{vtu_ind}_concentration.npz"
)
# Define final vtu concentration, and compute its mass (reference)
vtu_concentration_3d = aligned_vtu_concentration.copy()
vtu_3d_shape = vtu_concentration_3d.shape_metadata()
vtu_3d_geometry = darsia.Geometry(**vtu_3d_shape)
vtu_3d_integral = vtu_3d_geometry.integrate(vtu_concentration_3d)
# DICOM without TVD
def rescale_data(image, ref_integral):
shape = image.shape_metadata()
geometry = darsia.Geometry(**shape)
integral = geometry.integrate(image)
image.img *= ref_integral / integral
return image
# Define dicom concentration with same mass
dicom_concentration_3d = aligned_dicom_concentration.copy()
dicom_concentration_3d = rescale_data(dicom_concentration_3d, vtu_3d_integral)
# single_plot("slice", dicom_concentration_3d, "noisy", "plots/block-b/lab_data_slice.png")
# single_plot("sum", dicom_concentration_3d, "noisy", "plots/block-b/lab_data_avg.png")
# Define mask (omega) for trust for regularization
heterogeneous_omega = True
if heterogeneous_omega:
dicom_rescaled = dicom_concentration_3d.copy()
dicom_rescaled.img /= np.max(dicom_rescaled.img)
omega_bound = 0.15
omega = np.minimum(dicom_rescaled.img, omega_bound)
mask_zero = dicom_rescaled.img < 1e-4
omega[mask_zero] = 10
# plot_slice(omega)
# Save plot
reduced_dicom_concentration_2d = darsia.reduce_axis(
dicom_concentration_3d, axis="z", mode="slice", slice_idx=20
)
shape = reduced_dicom_concentration_2d.num_voxels
dimensions = reduced_dicom_concentration_2d.dimensions
x_pixel, y_pixel = np.meshgrid(
np.linspace(dimensions[0], 0, shape[0]),
np.linspace(0, dimensions[1], shape[1]),
indexing="ij",
)
fig_mask, axs_mask = plt.subplots(nrows=1, ncols=1)
axs_mask.pcolormesh(y_pixel, x_pixel, np.log(omega[20]), cmap="turbo")
axs_mask.set_ylim(top=0.08)
axs_mask.set_aspect("equal")
axs_mask.set_xlabel("x [m]") # , fontsize=14)
axs_mask.set_ylabel("y [m]") # , fontsize=14)
fig_mask.colorbar(
cm.ScalarMappable(cmap="turbo"), ax=axs_mask, label="log($\omega_f$)"
)
fig_mask.savefig("plots/block-b/omega.png", dpi=500, transparent=True)
plt.close()
else:
omega = 0.015
# Choose file suffix depeneding on omega
suffix = "heter" if heterogeneous_omega else "hom"
# DICOM concentration with H1 regularization
if True:
h1_reg_dicom_concentration_3d = darsia.H1_regularization(
dicom_concentration_3d,
mu=0.1,
omega=omega,
dim=3,
solver=darsia.CG(maxiter=10000, tol=1e-5),
)
h1_reg_dicom_concentration_3d.save(
f"data/npz/block-b/h1_reg_dicom_concentration_{suffix}.npz"
)
else:
h1_reg_dicom_concentration_3d = darsia.imread(
f"data/npz/block-b/h1_reg_dicom_concentration_{suffix}.npz"
)
h1_reg_dicom_concentration_3d = rescale_data(
h1_reg_dicom_concentration_3d, vtu_3d_integral
)
# DICOM concentration with TVD regularization
isotropic = True
isotropic_suffix = "iso" if isotropic else "aniso"
if True:
tvd_reg_dicom_concentration_3d = darsia.tvd(
dicom_concentration_3d,
method="heterogeneous bregman",
isotropic=isotropic,
weight=0.02,
omega=omega,
dim=3,
max_num_iter=100,
eps=1e-5,
verbose=True,
solver=darsia.Jacobi(maxiter=20),
)
tvd_reg_dicom_concentration_3d.save(
f"data/npz/block-b/tvd_{isotropic_suffix}_reg_dicom_concentration_{suffix}.npz"
)
else:
tvd_reg_dicom_concentration_3d = darsia.imread(
f"data/npz/block-b/tvd_{isotropic_suffix}_reg_dicom_concentration_{suffix}.npz"
)
tvd_reg_dicom_concentration_3d = rescale_data(
tvd_reg_dicom_concentration_3d, vtu_3d_integral
)
# Make qualitative comparisons
if True:
# plot_sum(
# [
# vtu_concentration_3d,
# dicom_concentration_3d,
# h1_reg_dicom_concentration_3d,
# tvd_reg_dicom_concentration_3d,
# ]
# )
# plot_slice(
# [
# vtu_concentration_3d,
# dicom_concentration_3d,
# h1_reg_dicom_concentration_3d,
# tvd_reg_dicom_concentration_3d,
# ]
# )
# qualitative_comparison(
# "sum",
# dicom_concentration_3d,
# vtu_concentration_3d,
# "noisy",
# "plots/block-b/pure_dicom_avg.png",
# )
qualitative_comparison(
"slice",
dicom_concentration_3d,
vtu_concentration_3d,
"noisy",
"plots/block-b/pure_dicom_slice.png",
)
# qualitative_comparison(
# "sum",
# h1_reg_dicom_concentration_3d,
# vtu_concentration_3d,
# "H1",
# f"plots/block-b/h1_reg_dicom_avg_{suffix}.png",
# )
qualitative_comparison(
"slice",
h1_reg_dicom_concentration_3d,
vtu_concentration_3d,
"H1",
f"plots/block-b/h1_reg_dicom_slice_{suffix}.png",
)
qualitative_comparison(
"sum",
tvd_reg_dicom_concentration_3d,
vtu_concentration_3d,
"TVD",
f"plots/block-b/tvd_{isotropic_suffix}_reg_dicom_avg_{suffix}.png",
colorbar=True,
)
qualitative_comparison(
"slice",
tvd_reg_dicom_concentration_3d,
vtu_concentration_3d,
"TVD",
f"plots/block-b/tvd_{isotropic_suffix}_reg_dicom_slice_{suffix}.png",
)
# Quantitative comparison - Wasserstein distance in 2d (only for illustration purposes)
slice_idx = 20
dicom_slice = darsia.reduce_axis(
dicom_concentration_3d, axis="z", mode="slice", slice_idx=slice_idx
)
vtu_slice = darsia.reduce_axis(
vtu_concentration_3d, axis="z", mode="slice", slice_idx=slice_idx
)
h1_reg_dicom_slice = darsia.reduce_axis(
h1_reg_dicom_concentration_3d, axis="z", mode="slice", slice_idx=slice_idx
)
tvd_reg_dicom_slice = darsia.reduce_axis(
tvd_reg_dicom_concentration_3d, axis="z", mode="slice", slice_idx=slice_idx
)
# Rescale (unphysical but necessary for comparison)
def rescale_slice(image: darsia.Image, ref_integral) -> darsia.Image:
shape = image.shape_metadata()
geometry = darsia.Geometry(**shape)
integral = geometry.integrate(image)
image.img *= ref_integral / integral
return image
# Determine reference value
shape = vtu_slice.shape_metadata()
geometry = darsia.Geometry(**shape)
ref_integral = geometry.integrate(vtu_slice)
dicom_slice = rescale_slice(dicom_slice, ref_integral)
h1_reg_dicom_slice = rescale_slice(h1_reg_dicom_slice, ref_integral)
tvd_reg_dicom_slice = rescale_slice(tvd_reg_dicom_slice, ref_integral)
# Wasserstein for 2d images
options = {
"method": "newton",
"L": 1e-2,
"num_iter": 500,
"tol": 1e-11,
"tol_distance": 1e-12,
"regularization": 1e-16,
"depth": 10,
"lumping": True,
"verbose": True,
"scaling": 10,
"plot_frequency": 10,
}
kwargs = {
"method": "newton",
"options": options,
"plot_solution": True,
}
distance_dicom_vtu = darsia.wasserstein_distance(
dicom_slice, vtu_slice, name="noisy vs. simulation", **kwargs
)
distance_h1_dicom_vtu = darsia.wasserstein_distance(
h1_reg_dicom_slice, vtu_slice, name="H1 vs simulation", **kwargs
)
distance_tvd_dicom_vtu = darsia.wasserstein_distance(
tvd_reg_dicom_slice, vtu_slice, name="TVD vs simulation", **kwargs
)
print("The distances:")
print("Pure DICOM: ", distance_dicom_vtu)
print("H1 reg DICOM: ", distance_h1_dicom_vtu)
print("TVD reg DICOM: ", distance_tvd_dicom_vtu)
assert False, "3d Wasserstein distance computations not sufficiently efficient."
# Quantitative comparison - Wasserstein distance in 3d
options = {
"method": "newton",
"L": 1e-2,
"num_iter": 100,
"tol": 1e-10,
"tol_distance": 1e-11,
"regularization": 1e-16,
"depth": 10,
"lumping": True,
"verbose": True,
}
kwargs = {
"method": "newton",
"options": options,
"plot_solution": False,
}
distance_dicom_vtu_3d = darsia.wasserstein_distance_3d(
dicom_concentration_3d, vtu_concentration_3d, name="", **kwargs
)
distance_h1_dicom_vtu_3d = darsia.wasserstein_distance_3d(
h1_reg_dicom_concentration_3d,
vtu_concentration_3d,
name="",
**kwargs,
)
distance_tvd_dicom_vtu_3d = darsia.wasserstein_distance_3d(
tvd_reg_dicom_concentration_3d,
vtu_concentration_3d,
name="",
**kwargs,
)
print("The distances 3d:")
print("Pure DICOM: ", distance_dicom_vtu_3d)
print("H1 reg DICOM: ", distance_h1_dicom_vtu_3d)
print("TVD reg DICOM: ", distance_tvd_dicom_vtu_3d)