Segmentation Workflow

import SimpleITK as sitk

# Manual threshold
binary = sitk.BinaryThreshold(
    image,
    lowerThreshold=100,
    upperThreshold=255,
    insideValue=1,
    outsideValue=0
)

# Otsu's method (minimizes within-class variance)
otsu = sitk.OtsuThreshold(image)

# Get computed threshold
otsu_filter = sitk.OtsuThresholdImageFilter()
otsu_filter.Execute(image)
threshold_value = otsu_filter.GetThreshold()
print(f"Otsu threshold: {threshold_value}")

# Other automatic methods
li_threshold = sitk.LiThreshold(image)
triangle_threshold = sitk.TriangleThreshold(image)

def segment_with_seeds(image, seed_points, lower, upper):
    """Segment region from seeds within intensity range."""
    
    segmentation = sitk.ConnectedThreshold(
        image,
        seedList=seed_points,
        lower=lower,
        upper=upper,
        replaceValue=1
    )
    
    return segmentation

# Example usage
seeds = [(128, 128, 64), (130, 130, 64)]
tumor = segment_with_seeds(ct_image, seeds, lower=100, upper=200)

def confidence_segment(image, seed_points):
    """Automatic threshold from seed region statistics."""
    
    segmentation = sitk.ConfidenceConnected(
        image,
        seedList=seed_points,
        numberOfIterations=5,
        multiplier=2.5,  # Threshold = mean ± 2.5 * stddev
        initialNeighborhoodRadius=1
    )
    
    return segmentation

def watershed_segmentation(image):
    """Watershed segmentation from gradient."""
    
    # Compute gradient magnitude
    gradient = sitk.GradientMagnitudeRecursiveGaussian(image, sigma=1.0)
    
    # Watershed
    labels = sitk.MorphologicalWatershed(
        gradient,
        level=10,
        markWatershedLine=False,
        fullyConnected=False
    )
    
    return labels

def marker_watershed(image, markers):
    """Watershed with predefined markers."""
    
    # Compute gradient
    gradient = sitk.GradientMagnitude(image)
    
    # Watershed from markers
    labels = sitk.MorphologicalWatershedFromMarkers(
        gradient,
        markers,
        markWatershedLine=False
    )
    
    return labels

def complete_segmentation_pipeline(image_path, seed_points, output_path):
    """Complete segmentation workflow with preprocessing and postprocessing."""
    
    # 1. Read image
    image = sitk.ReadImage(image_path)
    
    # 2. Preprocessing
    # Denoise
    image = sitk.CurvatureAnisotropicDiffusion(
        image,
        timeStep=0.0625,
        conductanceParameter=3.0,
        numberOfIterations=5
    )
    
    # Normalize
    image = sitk.RescaleIntensity(image, outputMinimum=0, outputMaximum=1000)
    
    # 3. Segmentation
    segmentation = sitk.ConfidenceConnected(
        image,
        seedList=seed_points,
        numberOfIterations=5,
        multiplier=2.5
    )
    
    # 4. Postprocessing
    # Remove small objects
    segmentation = sitk.BinaryMorphologicalOpening(
        segmentation,
        kernelRadius=[2, 2, 2],
        kernelType=sitk.sitkBall
    )
    
    # Fill holes
    segmentation = sitk.BinaryMorphologicalClosing(
        segmentation,
        kernelRadius=[3, 3, 3],
        kernelType=sitk.sitkBall
    )
    
    # 5. Label connected components
    labeled = sitk.ConnectedComponent(segmentation)
    labeled = sitk.RelabelComponent(labeled, minimumObjectSize=100)
    
    # 6. Analyze results
    shape_stats = sitk.LabelShapeStatisticsImageFilter()
    shape_stats.Execute(labeled)
    
    print(f"Found {shape_stats.GetNumberOfLabels()} objects")
    for label in shape_stats.GetLabels():
        size = shape_stats.GetNumberOfPixels(label)
        centroid = shape_stats.GetCentroid(label)
        print(f"  Label {label}: {size} pixels, centroid at {centroid}")
    
    # 7. Save
    sitk.WriteImage(labeled, output_path)
    
    return labeled

def superpixel_segmentation(image):
    """SLIC superpixel segmentation."""
    
    # Generate superpixels
    superpixels = sitk.SLIC(
        image,
        superGridSize=(20, 20, 20),
        spatialProximityWeight=10.0,
        maximumNumberOfIterations=5,
        enforceConnectivity=True
    )
    
    # Analyze superpixel statistics
    intensity_stats = sitk.LabelIntensityStatisticsImageFilter()
    intensity_stats.Execute(image, superpixels)
    
    # Merge superpixels based on intensity
    merged_labels = {}
    for label in intensity_stats.GetLabels():
        mean_intensity = intensity_stats.GetMean(label)
        # Group by intensity ranges
        group = int(mean_intensity / 100)
        merged_labels[label] = group
    
    result = sitk.ChangeLabel(superpixels, merged_labels)
    
    return result

def level_set_segmentation(image, initial_segmentation):
    """Geodesic active contour level set segmentation."""
    
    # Compute feature image (edge map)
    gradient = sitk.GradientMagnitudeRecursiveGaussian(image, sigma=1.0)
    feature = sitk.Sigmoid(
        gradient,
        alpha=-1.0,
        beta=50.0,
        outputMinimum=0,
        outputMaximum=1
    )
    
    # Create initial level set (signed distance)
    initial_distance = sitk.SignedMaurerDistanceMap(
        initial_segmentation,
        insideIsPositive=False,
        useImageSpacing=True
    )
    
    # Geodesic active contour
    level_set = sitk.GeodesicActiveContourLevelSetImageFilter()
    level_set.SetPropagationScaling(1.0)
    level_set.SetCurvatureScaling(1.0)
    level_set.SetAdvectionScaling(1.0)
    level_set.SetMaximumRMSError(0.02)
    level_set.SetNumberOfIterations(100)
    
    final_level_set = level_set.Execute(initial_distance, feature)
    
    # Extract zero level set
    segmentation = sitk.BinaryThreshold(
        final_level_set,
        lowerThreshold=-1000,
        upperThreshold=0,
        insideValue=1,
        outsideValue=0
    )
    
    return segmentation

def assess_segmentation(ground_truth, segmentation):
    """Assess segmentation quality against ground truth."""
    
    # Overlap measures
    overlap = sitk.LabelOverlapMeasuresImageFilter()
    overlap.Execute(ground_truth, segmentation)
    
    dice = overlap.GetDiceCoefficient()
    jaccard = overlap.GetJaccardCoefficient()
    
    print(f"Dice coefficient: {dice:.3f}")
    print(f"Jaccard coefficient: {jaccard:.3f}")
    
    # Hausdorff distance
    hausdorff = sitk.HausdorffDistanceImageFilter()
    hausdorff.Execute(ground_truth, segmentation)
    
    print(f"Hausdorff distance: {hausdorff.GetHausdorffDistance():.2f}")
    print(f"Average Hausdorff: {hausdorff.GetAverageHausdorffDistance():.2f}")
    
    return {
        'dice': dice,
        'jaccard': jaccard,
        'hausdorff': hausdorff.GetHausdorffDistance()
    }