diff --git a/python/demo/medical_imaging/data/download_cbis_ddsm.py b/python/demo/medical_imaging/data/download_cbis_ddsm.py
new file mode 100644
index 0000000000..bb4f108df2
--- /dev/null
+++ b/python/demo/medical_imaging/data/download_cbis_ddsm.py
@@ -0,0 +1,112 @@
+import os
+import numpy as np
+import zipfile
+from PIL import Image
+from pathlib import Path
+import argparse
+import shutil
+
+def download_dataset(output_dir="data/cbis_ddsm"):
+    output_path = Path(output_dir)
+    output_path.mkdir(parents=True, exist_ok=True)
+    
+    os.system(f"kaggle datasets download -d awsaf49/cbis-ddsm-breast-cancer-image-dataset -p {output_dir}")
+    
+    zip_file = output_path / "cbis-ddsm-breast-cancer-image-dataset.zip"
+    if zip_file.exists():
+        print("Extracting dataset...")
+        with zipfile.ZipFile(zip_file, 'r') as zip_ref:
+            zip_ref.extractall(output_path)
+        zip_file.unlink()
+        return True
+    else:
+        print("Error: Dataset download failed.")
+        return False
+
+def prepare_sample_images(dataset_dir="data/cbis_ddsm", sample_dir="data/samples", n_samples=5):
+    sample_path = Path(sample_dir)
+    sample_path.mkdir(parents=True, exist_ok=True)
+    
+    dataset_path = Path(dataset_dir)
+    image_files = list(dataset_path.rglob("*.png"))[:n_samples]
+    
+    if not image_files:
+        print("No images found. Please download the dataset first.")
+        return
+    
+    print(f"Copying {len(image_files)} sample images...")
+    
+    for i, img_file in enumerate(image_files):
+        img = Image.open(img_file).convert('L')
+        max_size = 1024
+        if max(img.size) > max_size:
+            img.thumbnail((max_size, max_size), Image.Resampling.LANCZOS)
+        output_file = sample_path / f"sample_{i:02d}.png"
+        img.save(output_file)
+        print(f"  Saved: {output_file}")
+    
+    print(f"Sample images saved to {sample_dir}/")
+
+def create_synthetic_test_image(output_file="data/synthetic_test.png", size=(512, 512)):
+    # Create clean image with geometric shapes
+    img = np.zeros(size, dtype=np.float64)
+    
+    y, x = np.ogrid[:size[0], :size[1]]
+    
+    # Circle 1 (top-left) - bright white
+    cx1, cy1, r1 = size[1]//4, size[0]//4, 60
+    mask1 = (x - cx1)**2 + (y - cy1)**2 <= r1**2
+    img[mask1] = 1.0
+    
+    # Circle 2 (bottom-right) - slightly dimmer
+    cx2, cy2, r2 = 3*size[1]//4, 3*size[0]//4, 80
+    mask2 = (x - cx2)**2 + (y - cy2)**2 <= r2**2
+    img[mask2] = 0.9
+    
+    # Rectangle (center)
+    img[size[0]//3:2*size[0]//3, size[1]//2-40:size[1]//2+40] = 0.8
+    
+    # Add light Gaussian noise
+    np.random.seed(42)  # Reproducible
+    noise = np.random.normal(0, 0.03, size)
+    noisy_img = img + noise
+    
+    # Clip to [0, 1]
+    noisy_img = np.clip(noisy_img, 0, 1)
+    
+    # Convert to 8-bit with high quality
+    img_uint8 = (noisy_img * 255).astype(np.uint8)
+    
+    # Save with maximum quality
+    Image.fromarray(img_uint8).save(output_file, quality=100, optimize=False)
+    
+    print(f"Improved synthetic test image saved to {output_file}")
+    print(f"    Image size: {size}")
+    print(f"    Normalized range: [0, 1]")
+    print(f"    Noise level: LIGHT (std=0.03)")
+    print(f"    Features: 2 circles + 1 rectangle")
+    print(f"    Quality: HIGH (no compression)")
+    
+    return output_file
+
+if __name__ == "__main__":   
+    parser = argparse.ArgumentParser(description="Download and prepare CBIS-DDSM dataset")
+    parser.add_argument("--download", action="store_true", help="Download full dataset")
+    parser.add_argument("--samples", action="store_true", help="Create sample subset")
+    parser.add_argument("--synthetic", action="store_true", help="Create synthetic test image")
+    parser.add_argument("--all", action="store_true", help="Do all of the above")
+    
+    args = parser.parse_args()
+    
+    if args.all or args.synthetic:
+        create_synthetic_test_image()
+    
+    if args.all or args.download:
+        success = download_dataset()
+        if success and (args.all or args.samples):
+            prepare_sample_images()
+    elif args.samples:
+        prepare_sample_images()
+    
+    if not any([args.download, args.samples, args.synthetic, args.all]):
+        parser.print_help()
\ No newline at end of file
diff --git a/python/demo/medical_imaging/data/sample_mammogram.jpg b/python/demo/medical_imaging/data/sample_mammogram.jpg
new file mode 100644
index 0000000000..b2fa341e97
Binary files /dev/null and b/python/demo/medical_imaging/data/sample_mammogram.jpg differ
diff --git a/python/demo/medical_imaging/demo/demo_anisotropic_diffusion_medical.py b/python/demo/medical_imaging/demo/demo_anisotropic_diffusion_medical.py
new file mode 100644
index 0000000000..ba7fb0a68e
--- /dev/null
+++ b/python/demo/medical_imaging/demo/demo_anisotropic_diffusion_medical.py
@@ -0,0 +1,205 @@
+import numpy as np
+from mpi4py import MPI
+from dolfinx import mesh, fem, io
+from dolfinx.fem.petsc import assemble_matrix, assemble_vector, LinearProblem
+import ufl
+from ufl import grad, div, dx, dot, sqrt, inner
+from petsc4py import PETSc
+import sys
+import os
+from basix.ufl import element
+
+sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..'))
+from utils.image_processing import (
+    load_image, save_image, create_image_mesh, 
+    image_to_function, function_to_image,
+    compute_psnr, compute_ssim, edge_preservation_index
+)
+from data.download_cbis_ddsm import create_synthetic_test_image
+
+class AnisotropicDiffusion:
+    def __init__(self, image_shape, kappa=30.0, dt=0.1, comm=MPI.COMM_WORLD):
+        self.image_shape = image_shape
+        self.kappa = kappa
+        self.dt = dt
+        self.comm = comm
+
+        self.domain = create_image_mesh(image_shape, comm)
+        
+        P1 = element("Lagrange", self.domain.basix_cell(), 1)
+        self.V = fem.functionspace(self.domain, P1)
+        
+        self.u = ufl.TrialFunction(self.V)
+        self.v = ufl.TestFunction(self.V)
+        
+        self.u_n = fem.Function(self.V)
+        
+    def diffusion_coefficient(self, grad_u):
+        # Compute gradient magnitude with small epsilon to avoid division by zero
+        grad_magnitude = sqrt(dot(grad_u, grad_u) + 1e-10)
+        
+        # Perona-Malik diffusion coefficient
+        c = 1.0 / (1.0 + (grad_magnitude / self.kappa)**2)
+        
+        return c
+    
+    def setup_variational_problem(self):
+        c = self.diffusion_coefficient(grad(self.u_n))
+        
+        F = (self.u - self.u_n) * self.v * dx + \
+            self.dt * c * dot(grad(self.u), grad(self.v)) * dx
+        
+        a = ufl.lhs(F)
+        L = ufl.rhs(F)
+        
+        return a, L
+    
+    def solve_timestep(self, a, L):
+        A = assemble_matrix(fem.form(a))
+        A.assemble()
+        b = assemble_vector(fem.form(L))
+        b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
+        
+        solver = PETSc.KSP().create(self.comm)
+        solver.setOperators(A)
+        solver.setType(PETSc.KSP.Type.CG)
+        
+        pc = solver.getPC()
+        pc.setType(PETSc.PC.Type.JACOBI)
+        
+        u_new = fem.Function(self.V)
+        solver.solve(b, u_new.x.petsc_vec)
+        u_new.x.scatter_forward()
+        
+        return u_new
+    
+    def denoise(self, image, n_iterations=20, verbose=True):
+        self.u_n = image_to_function(image, self.V)
+        
+        a, L = self.setup_variational_problem()
+        
+        for iteration in range(n_iterations):
+            if verbose and self.comm.rank == 0:
+                print(f"Iteration {iteration + 1}/{n_iterations}")
+            
+            u_new = self.solve_timestep(a, L)
+            
+            self.u_n.x.array[:] = u_new.x.array[:]
+            
+            a, L = self.setup_variational_problem()
+        
+        denoised = function_to_image(self.u_n, self.image_shape)
+        
+        return denoised
+
+def demo_synthetic_image():
+    print(f"Demo: Anisotropic Diffusion on Synthetic Image\n")
+    
+    # Load or create synthetic test image
+    test_image_path = "data/synthetic_test.png"
+    
+    if not os.path.exists(test_image_path):
+        print("Creating synthetic test image...")
+        os.makedirs("data", exist_ok=True)
+        create_synthetic_test_image(test_image_path)
+    
+    # Load image
+    print(f"Loading image: {test_image_path}")
+    image = load_image(test_image_path, normalize=True)
+    print(f"Image shape: {image.shape}")
+    
+    solver = AnisotropicDiffusion(
+        image_shape=image.shape,
+        kappa=20.0,
+        dt=0.05
+    )
+    
+    print("\nApplying anisotropic diffusion...")
+    denoised = solver.denoise(image, n_iterations=40, verbose=True)
+    
+    # Compute quality metrics
+    print(f"Quality Metrics:\n")
+    psnr = compute_psnr(image, denoised)
+    ssim = compute_ssim(image, denoised)
+    epi = edge_preservation_index(image, denoised)
+    
+    print(f"PSNR: {psnr:.2f} dB")
+    print(f"SSIM: {ssim:.4f}")
+    print(f"Edge Preservation Index: {epi:.4f}")
+    
+    # Save results
+    os.makedirs("output", exist_ok=True)
+    save_image(denoised, "output/synthetic_denoised.png")
+    print(f"\nDenoised image saved to: output/synthetic_denoised.png")
+    
+    return image, denoised
+
+
+def demo_medical_image(image_path):
+    print(f"Demo: Anisotropic Diffusion on Medical Image\n")
+    
+    # Load image
+    print(f"Loading image: {image_path}")
+    
+    target_size = (512, 512)
+    image = load_image(image_path, normalize=True, target_size=target_size)
+    print(f"Image shape: {image.shape}")
+    
+    print("\nInitializing anisotropic diffusion solver...")
+    solver = AnisotropicDiffusion(
+        image_shape=image.shape,
+        kappa=20.0,
+        dt=0.05
+    )
+    
+    print("\nApplying anisotropic diffusion...")
+    denoised = solver.denoise(image, n_iterations=30, verbose=True)
+    
+    print(f"Quality Metrics:\n")
+    psnr = compute_psnr(image, denoised)
+    ssim = compute_ssim(image, denoised)
+    epi = edge_preservation_index(image, denoised)
+    
+    print(f"PSNR: {psnr:.2f} dB")
+    print(f"SSIM: {ssim:.4f}")
+    print(f"Edge Preservation Index: {epi:.4f}")
+    
+    # Save results
+    os.makedirs("output", exist_ok=True)
+    output_path = "output/medical_denoised.png"
+    save_image(denoised, output_path)
+    print(f"\nDenoised image saved to: {output_path}")
+    
+    return image, denoised
+
+
+if __name__ == "__main__":
+    import argparse
+    
+    parser = argparse.ArgumentParser(
+        description="Anisotropic diffusion for medical image denoising"
+    )
+    parser.add_argument(
+        "--image",
+        type=str,
+        help="Path to input image (if not provided, uses synthetic test image)"
+    )
+    parser.add_argument(
+        "--kappa",
+        type=float,
+        default=30.0,
+        help="Gradient threshold (default: 30.0)"
+    )
+    parser.add_argument(
+        "--iterations",
+        type=int,
+        default=20,
+        help="Number of iterations (default: 20)"
+    )
+    
+    args = parser.parse_args()
+    
+    if args.image:
+        demo_medical_image(args.image)
+    else:
+        demo_synthetic_image()
\ No newline at end of file
diff --git a/python/demo/medical_imaging/requirements.txt b/python/demo/medical_imaging/requirements.txt
new file mode 100644
index 0000000000..677b01f645
--- /dev/null
+++ b/python/demo/medical_imaging/requirements.txt
@@ -0,0 +1,24 @@
+# Core Dependencies 
+# (Based on Last Stable Release: Dolfinx v0.9.0)
+fenics-dolfinx==0.9.0
+numpy>=1.24.0
+matplotlib>=3.7.0
+scipy>=1.10.0
+
+# Image Processing
+pillow>=10.0.0
+scikit-image>=0.21.0
+
+# Visualization
+pyvista>=0.43.0
+
+# Data handling
+pandas>=2.0.0
+kaggle>=1.6.0 # Unique to the CBIS DDSM Dataset
+
+# Testing
+pytest>=7.4.0
+pytest-mpi>=0.6
+
+# For DICOM
+pydicom>=2.4.0
\ No newline at end of file
diff --git a/python/demo/medical_imaging/utils/__init__.py b/python/demo/medical_imaging/utils/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/python/demo/medical_imaging/utils/image_processing.py b/python/demo/medical_imaging/utils/image_processing.py
new file mode 100644
index 0000000000..8bb6245d0d
--- /dev/null
+++ b/python/demo/medical_imaging/utils/image_processing.py
@@ -0,0 +1,141 @@
+import numpy as np
+from PIL import Image
+from dolfinx import mesh, fem
+from mpi4py import MPI
+import ufl
+from scipy import ndimage
+from skimage.metrics import structural_similarity
+from scipy.interpolate import griddata
+from basix.ufl import element
+
+
+def load_image(filepath, normalize=True, target_size=None):
+    img = Image.open(filepath).convert('L')  # Convert to grayscale
+    
+    if target_size is not None:
+        img = img.resize((target_size[1], target_size[0]), Image.Resampling.LANCZOS)
+    
+    image = np.array(img, dtype=np.float64)
+    
+    if normalize:
+        image = image / 255.0
+    
+    return image
+
+
+def save_image(image, filepath, denormalize=True):
+    if denormalize:
+        image = np.clip(image * 255.0, 0, 255)
+    
+    img = Image.fromarray(image.astype(np.uint8))
+    img.save(filepath)
+
+
+def create_image_mesh(image_shape, comm=MPI.COMM_WORLD):
+    height, width = image_shape
+    
+    # Create unit square mesh and scale to image dimensions
+    domain = mesh.create_rectangle(
+        comm,
+        [[0.0, 0.0], [float(width), float(height)]],
+        [width-1, height-1],
+        cell_type=mesh.CellType.triangle
+    )
+    
+    return domain
+
+
+def image_to_function(image, V):
+    u = fem.Function(V)
+    
+    # Get mesh coordinates
+    height, width = image.shape
+    
+    # Interpolate image values onto function space
+    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]
+            
+            # Convert to pixel coordinates (with bounds checking)
+            px = int(np.clip(np.round(x_coord), 0, width - 1))
+            py = int(np.clip(np.round(height - 1 - y_coord), 0, height - 1))
+            
+            values[i] = image[py, px]
+        
+        return values
+    
+    u.interpolate(image_values)
+    
+    return u
+
+
+def function_to_image(u, image_shape): 
+    height, width = image_shape
+    
+    # Get mesh coordinates and function values
+    mesh = u.function_space.mesh
+    coords = mesh.geometry.x[:, :2]  # Get x, y coordinates (drop z)
+    values = u.x.array
+    
+    x_pixels = np.linspace(0, width - 1, width)
+    y_pixels = np.linspace(0, height - 1, height)
+    X, Y = np.meshgrid(x_pixels, y_pixels)
+    
+    Y_flipped = height - 1 - Y
+    
+    # Interpolate function values onto regular pixel grid
+    # Use 'linear' interpolation for smooth results
+    image = griddata(
+        coords,
+        values,
+        (X, Y_flipped),
+        method='linear',
+        fill_value=0.0
+    )
+    
+    image = np.nan_to_num(image, nan=0.0)
+    
+    return image
+
+
+def compute_psnr(original, denoised, max_value=1.0):
+    mse = np.mean((original - denoised) ** 2)
+    if mse == 0:
+        return float('inf')
+    
+    psnr = 20 * np.log10(max_value / np.sqrt(mse))
+    return psnr
+
+
+def compute_ssim(original, denoised):   
+    return structural_similarity(original, denoised, data_range=denoised.max() - denoised.min())
+
+
+def edge_preservation_index(original, denoised):
+    sobel_orig = ndimage.sobel(original)
+    sobel_denoised = ndimage.sobel(denoised)
+    
+    # Normalize
+    sobel_orig = sobel_orig / (np.max(sobel_orig) + 1e-10)
+    sobel_denoised = sobel_denoised / (np.max(sobel_denoised) + 1e-10)
+    
+    # Compute correlation
+    numerator = np.sum(sobel_orig * sobel_denoised)
+    denominator = np.sqrt(np.sum(sobel_orig**2) * np.sum(sobel_denoised**2))
+    
+    epi = numerator / (denominator + 1e-10)
+    
+    return epi
+
+
+def add_gaussian_noise(image, var=0.01, seed=None):
+    if seed is not None:
+        np.random.seed(seed)
+    
+    noise = np.random.normal(0, np.sqrt(var), image.shape)
+    noisy = image + noise
+    
+    return np.clip(noisy, 0, 1)
\ No newline at end of file