Add medical image denoising demo for anisotropic diffusion

[m-wisell]

Summary

This PR adds medical image denoising using anisotropic diffusion.

Addresses: Issue #3973

Implementation

Mathematical Model:
Perona-Malik anisotropic diffusion for edge-preserving denoising

Numerical Methods:

• P1 Lagrange finite elements for spatial discretization
• Backward Euler time-stepping
• Conjugate Gradient solver with Jacobi preconditioner

Features:

• Image I/O with automatic normalization
• Image-to-mesh and mesh-to-image conversion
• Quality metrics: PSNR, SSIM, Edge Preservation Index
• Auto-generated synthetic test images (no dataset download required)

What Changed (Latest Update)

• Fixed CI failures: removed deprecated basix.ufl imports
• Consolidated to single standalone demo file
• Added test suite (all passing)
• Simplified dependencies and improved documentation

Quality Metrics (Real Mammogram 512x512)

• PSNR: 35.52 dB
• SSIM: 0.8739
• Edge Preservation: 0.8042

Usage

python demo/demo_anisotropic_diffusion_medical.py
python demo/demo_anisotropic_diffusion_medical.py --image path/to/image.jpg
python demo/demo_anisotropic_diffusion_medical.py --kappa 25 --iterations 40

Ready for review. CI check should now pass.

[m-wisell]

Updated to fix CI failures. All tests now passing locally. Ready for review.

[jorgensd]

This function will only work in serial, no?

[jorgensd]

Could you vectorize this? Something along the lines of:

from mpi4py import MPI
import dolfinx
import numpy as np

width = 3
height = 2

mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [[0,0],[width, height]], [int(width), int(height)])
V = dolfinx.fem.functionspace(mesh, ("DG", 0))
u = dolfinx.fem.Function(V)


a = np.linspace(0.5, width-0.5, width, endpoint=True)
b = np.linspace(0.5, height-0.5, height, endpoint=True)
mg = np.meshgrid(a, b)

def f(x, y):
    return x**2 + y**2
z = f(*mg)

    
def image_values(x):
    values = np.zeros(x.shape[1])
    
    for i in range(x.shape[1]):
        x_coord = x[0, i]
        y_coord = x[1, i]
        
        px = int(np.clip(np.round(x_coord), 0, width - 1))
        py = int(np.clip(np.round(height - y_coord), 0, height - 1))
        values[i] = z[py,px]
    
    return values

def images_values_vectorized(x):
    px = np.clip(np.round(x[0]), 0, width - 1).astype(np.int64)
    py = np.clip(np.round(height - x[1]), 0, height-1).astype(np.int64)
    return z[py, px]

u.interpolate(image_values)

u2 = dolfinx.fem.Function(V)
u2.interpolate(images_values_vectorized)
np.testing.assert_allclose(u.x.array, u2.x.array)

[jorgensd]

As your mesh aligns with the image grid, using interpolation to move data from a voxel grid (DG-0) to a P1 Lagrange space is not well defined. I would suggest that you for the first time step use u_n in DG-0, and then after solving, replace this function with a P1 function.

[jorgensd]

This function needs to be fully rewritten.

1. you call fem.form repeatedly, even if the form does not change over time. (with my suggestion above, it would change once, but it wouldn't modify the sparsity pattern of point 2, thus those structures can be created once)
2. The same holds for the matrix and vector. They should only be created once.
3. As you are not using any advanced PETSc features here, why not use dolfinx.fem.petsc.LinearProblem. (even if one changes the form once, by overwriting LinearProblem._a and LinearProblem._b with the form switching from DG-0 to P-1 for u_n)

[jhale]

I'm cautiously positive about the underlying idea for this demo (showing the use of a PDE solve on a problem from outside mechanics), but the execution is totally different to our existing demos:

• README.md
• No inline documentation using jupytext markdown.
• Very heavy dependency set which is going to be hard to maintain.
• We do not accept contributions that don't work in parallel.