import numpy as np  # NumPy: core array ops for CT HU manipulation and Poisson noise generation
from scipy.ndimage import gaussian_filter  # SciPy: unused here but available for post-smoothing low-dose CT

def simulate_low_dose(ct_slice, target_mAs_ratio=0.25):  # Simulate reduced tube current; mAs ratio 0.25 ≈ 75% dose reduction
    hu_image = ct_slice * 1.0  # Cast to float64 — HU values are typically int16 from DICOM; avoid integer division
    mu_image = (hu_image + 1000) / 1000.0  # Convert Hounsfield Units to linear attenuation coefficient μ; μ_water = 1.0 at 0 HU
    mu_image = np.clip(mu_image, 0.001, None)  # Floor at 0.001 to prevent log(0) in Poisson — air (-1000 HU) would give μ=0

    photon_count = mu_image * 100000 * target_mAs_ratio  # Expected photon count per detector element; higher μ = more attenuation = fewer photons
    noisy_photons = np.random.poisson(photon_count).astype(np.float32)  # Poisson noise models quantum mottle — the dominant CT noise source, signal-dependent

    noisy_mu = noisy_photons / (100000 * target_mAs_ratio)  # Invert: measured photons → reconstructed μ (noisy version)
    noisy_hu = noisy_mu * 1000 - 1000  # Convert μ back to Hounsfield Units: HU = 1000 × (μ − μ_water) / μ_water

    return noisy_hu  # Returns simulated low-dose CT slice in HU; train your model on this paired with clean reference

from skimage.restoration import denoise_nl_means  # Non-local means denoising — preserves edges while suppressing quantum mottle

def preprocess_ct(slice_hu):  # Inference-time denoising for low-dose CT; runs per 2D slice
    slice_hu = np.clip(slice_hu, -200, 400)  # Soft-tissue window: -200 HU (air/fat) to 400 HU (bone); discard irrelevant range
    denoised = denoise_nl_means(  # NLM works by averaging similar patches across the whole image — excellent for CT noise
        slice_hu / 600.0,  # Normalize to [≈-0.33, ≈0.67] so patch distances are well-scaled for the h parameter
        patch_size=5,  # 5×5 patch — small enough to capture local texture, large enough for reliable statistics
        patch_distance=9,  # Search window radius of 9 pixels — limits computational cost while finding enough similar patches
        h=0.1,  # Filter strength: 0.1 × standard deviation of noise; too low → under-denoised, too high → oversmoothed
        fast_mode=True  # Uses integral images for faster patch distance computation — acceptable quality, ~4× speedup
    ) * 600.0  # Rescale back to Hounsfield Units
    return denoised  # Cleaned slice ready for your segmentation/detection model

import numpy as np  # NumPy: FFT routines for k-space manipulation and random motion line selection

def simulate_motion_artifact(mri_slice, num_corrupt_lines=20):  # Bulk motion during readout: swap k-space lines to mimic patient movement
    kspace = np.fft.fft2(mri_slice)  # 2D FFT: image space → k-space (frequency/phase encoding domain)
    kspace_shifted = np.fft.fftshift(kspace)  # Shift zero-frequency component to center — DC is at [height/2, width/2]

    height = kspace_shifted.shape[0]  # Number of phase-encoding steps (rows in k-space)
    center = height // 2  # Center row index — central k-space lines encode contrast; peripheral encode fine detail
    corrupted = kspace_shifted.copy()  # Work on a copy so each motion realization is independent; keep original for comparison

    for _ in range(num_corrupt_lines):  # Each iteration simulates one "shot" where the patient moved during readout
        line = np.random.randint(center + 20, height - 10)  # Corrupt only outer k-space lines: central lines would obliterate image
        shift = np.random.randint(-5, 5)  # Random displacement in phase-encoding direction — mimics sudden body motion
        source_line = line + shift  # The line that was actually acquired instead of the intended line
        if 0 <= source_line < height:  # Clamp check: shifting off the k-space grid would cause index errors
            corrupted[line, :] = kspace_shifted[source_line, :]  # Replace corrupted line with the mis-acquired line

    result = np.fft.ifft2(np.fft.ifftshift(corrupted))  # Inverse shift then inverse FFT: corrupted k-space → image with motion ghosting
    return np.abs(result)  # Magnitude image: discard phase; clinical MRI typically displays magnitude only

from scipy.optimize import minimize  # SciPy optimizer: L-BFGS-B for entropy minimization over k-space degrees of freedom
import numpy as np  # NumPy: FFT for transform loops; flatten/reshape for optimizer interface

def correct_motion(kspace_data):  # Self-nav entropy minimization: motion-free images have lower entropy than motion-corrupted ones
    def image_entropy(kspace_flat):  # Objective function: returns negative entropy of the image reconstructed from flat k-space
        kspace = kspace_flat.reshape(kspace_data.shape)  # Reshape flat parameter vector back to 2D k-space matrix
        image = np.abs(np.fft.ifft2(np.fft.ifftshift(kspace)))  # Full reconstruction pipeline: shift → iFFT → magnitude
        image = image / image.sum()  # Normalize to probability distribution (sum = 1) for Shannon entropy calculation
        return -np.sum(image * np.log(image + 1e-10))  # Negative Shannon entropy: minimizer finds lowest-entropy (sharpest) image

    result = minimize(  # Quasi-Newton optimizer: searches k-space values to minimize image entropy
        image_entropy,  # Objective: entropy as proxy for motion artifact severity
        kspace_data.flatten(),  # Initial guess: the acquired (potentially motion-corrupted) k-space
        method="L-BFGS-B",  # Limited-memory BFGS with bounds — handles high-dimensional k-space without full Hessian
        options={"maxiter": 50}  # 50 iterations max — each iteration requires a full FFT; trade-off between quality and runtime
    )
    corrected = result.x.reshape(kspace_data.shape)  # Reshape the optimal k-space back to 2D
    return np.abs(np.fft.ifft2(np.fft.ifftshift(corrected)))  # Reconstruct the motion-corrected image from optimized k-space

def detect_acoustic_shadow(ultrasound_slice, threshold=0.1):  # Identify regions below a rib/bone/gas where signal is absent
    normalized = ultrasound_slice / (ultrasound_slice.max() + 1e-6)  # Scale to [0, 1] — the 1e-6 prevents division by zero for dark frames
    column_intensity = normalized.mean(axis=1)  # Mean intensity across each row: shadows appear as vertically dim columns
    shadow_mask = np.zeros_like(normalized)  # Initialize binary mask — 0 = visible tissue, 1 = acoustic shadow

    for col in range(normalized.shape[1]):  # Iterate over each column (A-line) independently — shadows are column-localized
        if column_intensity[col] < threshold:  # If this column's mean intensity is suspiciously low, it may be shadowed
            shadow_start = np.argmax(normalized[:, col] > threshold * 2)  # Find the first pixel exceeding 2× threshold — marks the shadow boundary
            shadow_mask[shadow_start:, col] = 1.0  # Everything below the shadow boundary in this column is marked as missing signal

    return shadow_mask  # Binary mask: 1 where ultrasound cannot penetrate; feed as extra channel to your CNN

shadow_mask = detect_acoustic_shadow(us_frame)  # Example: generate shadow mask for one B-mode ultrasound frame
model_input = np.stack([us_frame, shadow_mask], axis=-1)  # Stack intensity + mask as 2-channel input — model learns to ignore shadowed regions

def simulate_low_count_pet(pet_slice, count_fraction=0.25):  # Simulate reduced injected dose or shorter acquisition: fewer counts = noisier
    max_val = pet_slice.max()  # Peak tracer uptake in this slice — used to normalize activity to a relative scale
    activity = pet_slice / max_val  # Normalize to [0, 1]: relative tracer concentration (SUV-like proportional)
    low_count = np.random.poisson(activity * 1000 * count_fraction)  # Poisson noise on counts: fewer total counts = higher relative noise (lower SNR)
    low_count = low_count / (1000 * count_fraction) * max_val  # Rescale back to original intensity range — preserves SUV quantification
    return low_count.astype(pet_slice.dtype)  # Return noisy slice in original dtype (typically float32 from DICOM)

def gaussian_post_filter(pet_slice, sigma=1.5):  # Optional post-reconstruction smoothing: Gaussian kernel matches clinical PET reconstruction
    from scipy.ndimage import gaussian_filter  # Import inside function to keep it self-contained
    return gaussian_filter(pet_slice, sigma=sigma)  # sigma=1.5 mm: matches typical clinical PET reconstruction filter (e.g., 4-6 mm FWHM)

import numpy as np  # NumPy: random augmentation, meshgrid for bias field, basic statistics
from skimage.exposure import match_histograms  # skimage: match cumulative distribution of intensities across scanner vendors

def normalize_across_scanners(image, target_histogram=None):  # Z-score normalization with optional histogram matching to a reference scanner
    image = np.clip(image, -100, 300)  # Clip to soft-tissue window: -100 HU to 300 HU — ignores air (-1024) and dense bone (1000+)
    mean = image.mean()  # Mean intensity of the clipped region — varies between vendors due to different reconstruction kernels
    std = image.std() + 1e-6  # Standard deviation with epsilon — prevents division by zero for uniform noise-only regions
    normalized = (image - mean) / std  # Z-score normalization: centers intensity distribution at 0, scales to unit variance

    if target_histogram is not None:  # If a reference scanner histogram is provided, adjust this slice's distribution to match
        normalized = match_histograms(normalized, target_histogram)  # Per-pixel CDF matching — aligns dynamic range across Siemens/GE/Philips

    return normalized  # Normalized slice: closer distribution to the reference scanner; improves cross-vendor model generalisation

def random_intensity_augment(image):  # Domain randomization: applies random gamma and bias field to simulate scanner variability
    gamma = np.random.uniform(0.8, 1.2)  # Random gamma correction: 1.0 = linear; < 1 brightens dark tissue, > 1 darkens bright tissue
    image = image ** gamma  # Power-law transform: models nonlinear reconstruction transfer functions across scanner vendors

    x = np.linspace(-1, 1, image.shape[1])  # Horizontal coordinate normalized to [-1, 1] — centre of image is (0, 0)
    y = np.linspace(-1, 1, image.shape[0])  # Vertical coordinate normalized to [-1, 1]
    xx, yy = np.meshgrid(x, y)  # 2D grid of coordinates for per-pixel bias computation
    bias = 1 + 0.1 * np.random.randn() * (xx**2 + yy**2)  # Quadratic bias field: simulates slow-varying intensity inhomogeneity from B1 field non-uniformity
    image = image * bias  # Apply bias field multiplicatively — mimics receive-coil sensitivity profile variation

    return image  # Augmented slice: your model sees different "scanner feels" during training and becomes vendor-agnostic

import SimpleITK as sitk  # SimpleITK: NiftyReg-powered registration framework with medical image I/O

def register_pet_to_ct(pet_path, ct_path):  # Register low-resolution PET (functional) to high-resolution CT (anatomical) — standard clinical fusion
    pet = sitk.ReadImage(pet_path)  # Read PET DICOM/NIfTI: typically 128×128×? with large (3-5 mm) voxels, limited FOV
    ct = sitk.ReadImage(ct_path)  # Read CT DICOM/NIfTI: typically 512×512×? with small (0.5-1 mm) voxels, full anatomical coverage

    resampler = sitk.ResampleImageFilter()  # Filter that will warp PET to CT physical space using the final transform
    resampler.SetReferenceImage(ct)  # Set CT as reference: resampled PET will have same origin, spacing, direction, and size as CT
    resampler.SetInterpolator(sitk.sitkLinear)  # Linear interpolation for PET counts — nearest would destroy SUV quantification; cubic is overkill
    resampler.SetDefaultPixelValue(0)  # Voxels outside PET FOV after transform are set to zero (no tracer activity)

    transform = sitk.CenteredTransformInitializer(  # Compute initial rigid alignment from image centroids (center-of-mass)
        ct, pet,  # Fixed (CT) and moving (PET) images for initializer
        sitk.Euler3DTransform(),  # Rigid transform with 3 rotation + 3 translation DOF — patient position difference, not deformation
        sitk.CenteredTransformInitializerFilter.GEOMETRY  # Use geometric centres (centroids of the image grids), not intensity moments
    )

    optimizer = sitk.RegularStepGradientDescentOptimizerv4()  # Gradient descent with adaptive step size — robust to local minima in mutual information
    optimizer.SetLearningRate(1.0)  # Initial step size of 1.0 (in transform parameter space) — aggressive but decreased on plateau
    optimizer.SetMinimumStepLength(0.001)  # Stopping criterion: step size below 0.001 means the optimizer converged

    registration = sitk.ImageRegistrationMethod()  # Main registration object: ties together metric, optimizer, transform, interpolator
    registration.SetOptimizer(optimizer)  # Attach the gradient descent optimizer
    registration.SetInterpolator(sitk.sitkLinear)  # Interpolation for metric evaluation — linear is fast and sufficiently accurate
    registration.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)  # Mattes MI: measures statistical dependence between CT (HU) and PET (counts) — both modalities have no direct intensity relationship
    registration.SetInitialTransform(transform)  # Start from the centroid-initialized rigid transform

    final_transform = registration.Execute(  # Run the registration: iteratively adjusts the Euler transform to maximise MI
        sitk.Cast(ct, sitk.sitkFloat32),  # Cast CT to float32 — SimpleITK expects matching pixel types for metric computation
        sitk.Cast(pet, sitk.sitkFloat32)  # Cast PET to float32
    )

    resampler.SetTransform(final_transform)  # Apply the optimized transform to the resampler
    registered_pet = resampler.Execute(pet)  # Resample PET into CT space: the output PET image matches CT geometry

    return sitk.GetArrayFromImage(registered_pet)  # Return as NumPy array — shape matches CT, voxel values represent registered SUV

class ClinicalPreprocessor:  # End-to-end pipeline that bridges raw DICOM physics to normalized input tensors
    def __init__(self, modality="CT", target_spacing=(1.0, 1.0, 1.0)):  # Target spacing: 1 mm isotropic for most deep learning pipelines
        self.modality = modality  # Modality tag: controls modality-specific clipping and default augmentations
        self.target_spacing = target_spacing  # Resample target: CT slices are often 0.5-3 mm; models expect consistent voxel spacing

    def __call__(self, volume, metadata):  # Callable interface: accepts 3D volume + dict with "spacing" key (e.g., (0.7, 0.7, 1.25))
        volume = self.resample(volume, metadata["spacing"])  # Resample from acquisition spacing → target spacing: preserves anatomical scale for CNNs

        if self.modality == "CT":  # CT-specific: clip to diagnostic soft-tissue window — excludes air and metal/bone extremes
            volume = np.clip(volume, -200, 400)  # [-200, 400] HU: covers fat (-100), water (0), muscle (40), and contrast-enhanced tissue (100-300)
        elif self.modality == "PET":  # PET-specific: clip negative SUV values to zero — negative counts are reconstruction artifacts, not real uptake
            volume = np.clip(volume, 0, volume.max())  # SUV range: typically 0-15; max() preserves individual scan range

        volume = (volume - volume.mean()) / (volume.std() + 1e-6)  # Z-score normalisation across the whole volume: centres to zero mean, unit variance

        if np.random.random() < 0.3:  # 30% chance: apply Gaussian noise to simulate small coil/respiratory artifacts during training
            volume += np.random.normal(0, 0.05, volume.shape)  # Additive Gaussian noise: std 0.05 after z-score normalisation

        if np.random.random() < 0.2:  # 20% chance: apply smooth bias field to simulate MRI inhomogeneity or CT beam hardening
            volume = self.apply_bias_field(volume)  # Low-frequency intensity modulation; helps model generalise across field strengths (1.5T vs 3.0T)

        return volume.astype(np.float32)  # Cast to float32 — standard dtype for PyTorch/TensorFlow medical imaging pipelines