0
# Registration
1

2
Complete image registration framework for aligning images through optimization of similarity metrics. SimpleITK provides a comprehensive registration system that can handle rigid, affine, and deformable transformations with various similarity metrics and optimization methods.
3

4
## Capabilities
5

6
### ImageRegistrationMethod
7

8
The main registration framework that coordinates all registration components.
9

10
```python { .api }
11
class ImageRegistrationMethod:
12
    def __init__(self):
13
        """Initialize image registration method"""
14
    
15
    def SetFixedImage(self, image: Image):
16
        """Set the fixed (reference) image"""
17
    
18
    def GetFixedImage(self) -> Image:
19
        """Get the fixed image"""
20
    
21
    def SetMovingImage(self, image: Image):
22
        """Set the moving (image to be registered) image"""
23
    
24
    def GetMovingImage(self) -> Image:
25
        """Get the moving image"""
26
    
27
    def SetInitialTransform(self, transform: Transform):
28
        """
29
        Set initial transformation.
30
        
31
        Args:
32
            transform: Initial transform (will be optimized)
33
        """
34
    
35
    def GetInitialTransform(self) -> Transform:
36
        """Get initial transform"""
37
    
38
    def SetMovingInitialTransform(self, transform: Transform):
39
        """Set initial transform for moving image"""
40
    
41
    def GetMovingInitialTransform(self) -> Transform:
42
        """Get moving initial transform"""
43
    
44
    def SetFixedInitialTransform(self, transform: Transform):
45
        """Set initial transform for fixed image"""
46
    
47
    def GetFixedInitialTransform(self) -> Transform:
48
        """Get fixed initial transform"""
49
    
50
    def Execute(self) -> Transform:
51
        """
52
        Execute registration and return optimized transform.
53
        
54
        Returns:
55
            Optimized transformation
56
        """
57
    
58
    # Metric configuration
59
    def SetMetricAsMeanSquares(self):
60
        """Use mean squares similarity metric"""
61
    
62
    def SetMetricAsCorrelation(self):
63
        """Use normalized correlation metric"""
64
    
65
    def SetMetricAsANTSNeighborhoodCorrelation(self, radius: int = 5):
66
        """
67
        Use ANTS neighborhood correlation metric.
68
        
69
        Args:
70
            radius: Neighborhood radius
71
        """
72
    
73
    def SetMetricAsJointHistogramMutualInformation(self, numberOfHistogramBins: int = 20,
74
                                                  varianceForJointPDFSmoothing: float = 1.5):
75
        """
76
        Use joint histogram mutual information metric.
77
        
78
        Args:
79
            numberOfHistogramBins: Histogram bins
80
            varianceForJointPDFSmoothing: Smoothing variance
81
        """
82
    
83
    def SetMetricAsMattesMutualInformation(self, numberOfHistogramBins: int = 50):
84
        """
85
        Use Mattes mutual information metric.
86
        
87
        Args:
88
            numberOfHistogramBins: Histogram bins for estimation
89
        """
90
    
91
    def SetMetricFixedMask(self, mask: Image):
92
        """Set mask for fixed image metric computation"""
93
    
94
    def SetMetricMovingMask(self, mask: Image):
95
        """Set mask for moving image metric computation"""
96
    
97
    def SetMetricSamplingStrategy(self, strategy: str):
98
        """
99
        Set sampling strategy for metric computation.
100
        
101
        Args:
102
            strategy: 'NONE', 'REGULAR', or 'RANDOM'
103
        """
104
    
105
    def SetMetricSamplingPercentage(self, percentage: float):
106
        """
107
        Set percentage of pixels to sample for metric computation.
108
        
109
        Args:
110
            percentage: Sampling percentage (0.0 to 1.0)
111
        """
112
    
113
    # Optimizer configuration
114
    def SetOptimizerAsRegularStepGradientDescent(self, learningRate: float = 1.0,
115
                                                minStep: float = 0.001,
116
                                                numberOfIterations: int = 100,
117
                                                gradientMagnitudeTolerance: float = 1e-6):
118
        """
119
        Use regular step gradient descent optimizer.
120
        
121
        Args:
122
            learningRate: Initial learning rate
123
            minStep: Minimum step size for convergence
124
            numberOfIterations: Maximum iterations
125
            gradientMagnitudeTolerance: Gradient tolerance for convergence
126
        """
127
    
128
    def SetOptimizerAsGradientDescent(self, learningRate: float = 1.0,
129
                                     numberOfIterations: int = 100,
130
                                     convergenceMinimumValue: float = 1e-6,
131
                                     convergenceWindowSize: int = 10):
132
        """Use gradient descent optimizer"""
133
    
134
    def SetOptimizerAsGradientDescentLineSearch(self, learningRate: float = 1.0,
135
                                               numberOfIterations: int = 100,
136
                                               convergenceMinimumValue: float = 1e-6,
137
                                               convergenceWindowSize: int = 10):
138
        """Use gradient descent with line search optimizer"""
139
    
140
    def SetOptimizerAsConjugateGradientLineSearch(self, learningRate: float = 1.0,
141
                                                 numberOfIterations: int = 100,
142
                                                 convergenceMinimumValue: float = 1e-6,
143
                                                 convergenceWindowSize: int = 10):
144
        """Use conjugate gradient line search optimizer"""
145
    
146
    def SetOptimizerAsLBFGSB(self, gradientConvergenceTolerance: float = 1e-5,
147
                            numberOfIterations: int = 500,
148
                            maximumNumberOfCorrections: int = 5,
149
                            maximumNumberOfFunctionEvaluations: int = 2000,
150
                            costFunctionConvergenceFactor: float = 1e7):
151
        """Use L-BFGS-B optimizer for bounded optimization"""
152
    
153
    def SetOptimizerAsOnePlusOneEvolutionary(self, numberOfIterations: int = 100,
154
                                           epsilon: float = 1.5e-4,
155
                                           initialRadius: float = 1.01,
156
                                           growthFactor: float = -1.0,
157
                                           shrinkFactor: float = -1.0):
158
        """Use (1+1) evolutionary strategy optimizer"""
159
    
160
    def SetOptimizerAsPowell(self, maximumNumberOfIterations: int = 100,
161
                           maximumNumberOfLineIterations: int = 100,
162
                           stepLength: float = 1.0, stepTolerance: float = 1e-6,
163
                           valueTolerance: float = 1e-6):
164
        """Use Powell optimizer (derivative-free)"""
165
    
166
    def SetOptimizerScales(self, scales: list):
167
        """
168
        Set parameter scaling for optimizer.
169
        
170
        Args:
171
            scales: Scaling factors for each parameter
172
        """
173
    
174
    def SetOptimizerWeights(self, weights: list):
175
        """Set parameter weights for optimizer"""
176
    
177
    # Interpolator configuration
178
    def SetInterpolator(self, interpolator: int):
179
        """
180
        Set interpolation method.
181
        
182
        Args:
183
            interpolator: Interpolation type (e.g., sitk.sitkLinear, sitk.sitkBSpline)
184
        """
185
    
186
    def GetInterpolator(self) -> int:
187
        """Get interpolation method"""
188
    
189
    # Multi-resolution configuration
190
    def SetShrinkFactorsPerLevel(self, factors: list):
191
        """
192
        Set shrink factors for multi-resolution registration.
193
        
194
        Args:
195
            factors: List of shrink factors per resolution level
196
        """
197
    
198
    def GetShrinkFactorsPerLevel(self) -> list:
199
        """Get shrink factors per level"""
200
    
201
    def SetSmoothingSigmasPerLevel(self, sigmas: list):
202
        """
203
        Set smoothing sigmas for multi-resolution registration.
204
        
205
        Args:
206
            sigmas: List of smoothing sigmas per level
207
        """
208
    
209
    def GetSmoothingSigmasPerLevel(self) -> list:
210
        """Get smoothing sigmas per level"""
211
    
212
    def SetSmoothingSigmasAreSpecifiedInPhysicalUnits(self, inPhysicalUnits: bool):
213
        """Specify if smoothing sigmas are in physical units"""
214
    
215
    def GetSmoothingSigmasAreSpecifiedInPhysicalUnits(self) -> bool:
216
        """Check if smoothing sigmas are in physical units"""
217
    
218
    # Results and monitoring
219
    def GetOptimizerIteration(self) -> int:
220
        """Get current optimizer iteration"""
221
    
222
    def GetOptimizerPosition(self) -> tuple:
223
        """Get current parameter values"""
224
    
225
    def GetOptimizerLearningRate(self) -> float:
226
        """Get current learning rate"""
227
    
228
    def GetOptimizerConvergenceValue(self) -> float:
229
        """Get optimizer convergence value"""
230
    
231
    def GetMetricValue(self) -> float:
232
        """Get current metric value"""
233
    
234
    def GetOptimizerStopConditionDescription(self) -> str:
235
        """Get description of why optimizer stopped"""
236
    
237
    # Callbacks for monitoring
238
    def AddCommand(self, event: int, command: object):
239
        """
240
        Add command callback for registration events.
241
        
242
        Args:
243
            event: Event type (e.g., sitk.sitkIterationEvent)
244
            command: Callback function or command object
245
        """
246
    
247
    def RemoveAllCommands(self):
248
        """Remove all command callbacks"""
249
```
250

251
### Transform Initialization
252

253
Utilities for initializing transforms based on image geometry and landmarks.
254

255
```python { .api }
256
def CenteredTransformInitializer(fixedImage: Image, movingImage: Image, 
257
                               transform: Transform, operationMode: str = "GEOMETRY") -> Transform:
258
    """
259
    Initialize transform based on image centers.
260
    
261
    Args:
262
        fixedImage: Fixed image
263
        movingImage: Moving image
264
        transform: Transform to initialize
265
        operationMode: 'GEOMETRY' or 'MOMENTS'
266
        
267
    Returns:
268
        Initialized transform
269
    """
270

271
def CenteredVersorTransformInitializer(fixedImage: Image, movingImage: Image,
272
                                     transform: VersorRigid3DTransform) -> Transform:
273
    """Initialize versor transform using image centers"""
274

275
def LandmarkBasedTransformInitializer(fixedLandmarks: list, movingLandmarks: list,
276
                                    transform: Transform) -> Transform:
277
    """
278
    Initialize transform using corresponding landmark points.
279
    
280
    Args:
281
        fixedLandmarks: List of points in fixed image
282
        movingLandmarks: List of corresponding points in moving image
283
        transform: Transform to initialize
284
        
285
    Returns:
286
        Initialized transform
287
    """
288

289
class BSplineTransformInitializer:
290
    def __init__(self):
291
        """Initialize B-spline transform initializer"""
292
    
293
    def SetTransformDomainMeshSize(self, meshSize: list):
294
        """Set mesh size for B-spline grid"""
295
    
296
    def SetImage(self, image: Image):
297
        """Set reference image for transform domain"""
298
    
299
    def Execute(self) -> BSplineTransform:
300
        """Create initialized B-spline transform"""
301
```
302

303
### Resampling
304

305
Apply transforms to images for registration results or coordinate transformation.
306

307
```python { .api }
308
def Resample(image: Image, transform: Transform, interpolator: int = None,
309
            defaultPixelValue: float = 0.0, outputPixelType: int = None,
310
            outputSpacing: tuple = None, outputOrigin: tuple = None,
311
            outputDirection: tuple = None, outputSize: tuple = None,
312
            referenceImage: Image = None) -> Image:
313
    """
314
    Resample image using transform.
315
    
316
    Args:
317
        image: Input image to resample
318
        transform: Transformation to apply
319
        interpolator: Interpolation method
320
        defaultPixelValue: Value for pixels outside image
321
        outputPixelType: Output pixel type
322
        outputSpacing: Output image spacing
323
        outputOrigin: Output image origin
324
        outputDirection: Output image direction
325
        outputSize: Output image size
326
        referenceImage: Use geometry from reference image
327
        
328
    Returns:
329
        Resampled image
330
    """
331

332
class ResampleImageFilter:
333
    def __init__(self):
334
        """Initialize resampling filter"""
335
    
336
    def SetTransform(self, transform: Transform):
337
        """Set transformation"""
338
    
339
    def GetTransform(self) -> Transform:
340
        """Get transformation"""
341
    
342
    def SetInterpolator(self, interpolator: int):
343
        """Set interpolation method"""
344
    
345
    def GetInterpolator(self) -> int:
346
        """Get interpolation method"""
347
    
348
    def SetDefaultPixelValue(self, value: float):
349
        """Set default pixel value for outside regions"""
350
    
351
    def GetDefaultPixelValue(self) -> float:
352
        """Get default pixel value"""
353
    
354
    def SetSize(self, size: tuple):
355
        """Set output image size"""
356
    
357
    def GetSize(self) -> tuple:
358
        """Get output image size"""
359
    
360
    def SetOutputSpacing(self, spacing: tuple):
361
        """Set output image spacing"""
362
    
363
    def GetOutputSpacing(self) -> tuple:
364
        """Get output image spacing"""
365
    
366
    def SetOutputOrigin(self, origin: tuple):
367
        """Set output image origin"""
368
    
369
    def GetOutputOrigin(self) -> tuple:
370
        """Get output image origin"""
371
    
372
    def SetOutputDirection(self, direction: tuple):
373
        """Set output image direction"""
374
    
375
    def GetOutputDirection(self) -> tuple:
376
        """Get output image direction"""
377
    
378
    def SetReferenceImage(self, image: Image):
379
        """Set reference image for output geometry"""
380
    
381
    def Execute(self, image: Image) -> Image:
382
        """Execute resampling"""
383
```
384

385
### Registration Evaluation
386

387
Metrics and utilities for evaluating registration quality.
388

389
```python { .api }
390
def LabelOverlapMeasures(sourceImage: Image, targetImage: Image) -> dict:
391
    """
392
    Compute overlap measures between labeled images.
393
    
394
    Args:
395
        sourceImage: Source labeled image
396
        targetImage: Target labeled image
397
        
398
    Returns:
399
        Dictionary with overlap measures per label
400
    """
401

402
def HausdorffDistance(image1: Image, image2: Image) -> dict:
403
    """
404
    Compute Hausdorff distance between binary images.
405
    
406
    Args:
407
        image1: First binary image
408
        image2: Second binary image
409
        
410
    Returns:
411
        Dictionary with distance measures
412
    """
413

414
def SimilarityIndex(image1: Image, image2: Image) -> float:
415
    """
416
    Compute similarity index (Dice coefficient) between binary images.
417
    
418
    Args:
419
        image1: First binary image  
420
        image2: Second binary image
421
        
422
    Returns:
423
        Similarity index (0 to 1)
424
    """
425

426
class LabelOverlapMeasuresImageFilter:
427
    def __init__(self):
428
        """Initialize label overlap measures filter"""
429
    
430
    def Execute(self, sourceImage: Image, targetImage: Image):
431
        """Compute overlap measures"""
432
    
433
    def GetJaccardCoefficient(self, label: int = 1) -> float:
434
        """Get Jaccard coefficient for label"""
435
    
436
    def GetDiceCoefficient(self, label: int = 1) -> float:
437
        """Get Dice coefficient for label"""
438
    
439
    def GetVolumeSimilarity(self, label: int = 1) -> float:
440
        """Get volume similarity for label"""
441
    
442
    def GetFalseNegativeError(self, label: int = 1) -> float:
443
        """Get false negative error for label"""
444
    
445
    def GetFalsePositiveError(self, label: int = 1) -> float:
446
        """Get false positive error for label"""
447
    
448
    def GetTargetOverlap(self, label: int = 1) -> float:
449
        """Get target overlap for label"""
450
    
451
    def GetLabels(self) -> list:
452
        """Get list of labels processed"""
453
```
454

455
### Advanced Registration Components
456

457
Specialized registration components for complex scenarios.
458

459
```python { .api }
460
# Demons registration (deformable)
461
def DemonsRegistration(fixedImage: Image, movingImage: Image,
462
                      numberOfIterations: int = 10, standardDeviations: float = 1.0) -> Image:
463
    """
464
    Demons deformable registration.
465
    
466
    Args:
467
        fixedImage: Fixed image
468
        movingImage: Moving image  
469
        numberOfIterations: Number of iterations
470
        standardDeviations: Smoothing standard deviation
471
        
472
    Returns:
473
        Displacement field image
474
    """
475

476
def DiffeomorphicDemonsRegistration(fixedImage: Image, movingImage: Image,
477
                                   numberOfIterations: int = 10, 
478
                                   standardDeviations: float = 1.0) -> Image:
479
    """Diffeomorphic demons registration for topology preservation"""
480

481
def FastSymmetricForcesDemonsRegistration(fixedImage: Image, movingImage: Image,
482
                                        numberOfIterations: int = 10,
483
                                        standardDeviations: float = 1.0) -> Image:
484
    """Fast symmetric forces demons registration"""
485

486
def SymmetricForcesDemonsRegistration(fixedImage: Image, movingImage: Image,
487
                                    numberOfIterations: int = 10,
488
                                    standardDeviations: float = 1.0) -> Image:
489
    """Symmetric forces demons registration"""
490

491
# Level set motion registration  
492
def LevelSetMotionRegistration(fixedImage: Image, movingImage: Image,
493
                              numberOfIterations: int = 10, 
494
                              standardDeviations: float = 1.0) -> Image:
495
    """Level set motion registration"""
496

497
# Displacement field utilities
498
def InverseDisplacementField(displacementField: Image, 
499
                           maximumNumberOfIterations: int = 10,
500
                           maxErrorToleranceThreshold: float = 0.1,
501
                           meanErrorToleranceThreshold: float = 0.001,
502
                           enforceBoundaryCondition: bool = True) -> Image:
503
    """
504
    Compute inverse of displacement field.
505
    
506
    Args:
507
        displacementField: Forward displacement field
508
        maximumNumberOfIterations: Maximum iterations for inversion
509
        maxErrorToleranceThreshold: Maximum error tolerance
510
        meanErrorToleranceThreshold: Mean error tolerance
511
        enforceBoundaryCondition: Enforce boundary conditions
512
        
513
    Returns:
514
        Inverse displacement field
515
    """
516

517
def IterativeInverseDisplacementField(displacementField: Image,
518
                                    numberOfIterations: int = 5,
519
                                    stopValue: float = 0.0) -> Image:
520
    """Iterative computation of inverse displacement field"""
521

522
def TransformToDisplacementField(transform: Transform, outputSize: tuple,
523
                               outputSpacing: tuple, outputOrigin: tuple,
524
                               outputDirection: tuple = None) -> Image:
525
    """
526
    Convert transform to displacement field.
527
    
528
    Args:
529
        transform: Input transformation
530
        outputSize: Size of output displacement field
531
        outputSpacing: Spacing of output field
532
        outputOrigin: Origin of output field
533
        outputDirection: Direction of output field
534
        
535
    Returns:
536
        Displacement field image
537
    """
538

539
def DisplacementFieldJacobianDeterminant(displacementField: Image,
540
                                        useImageSpacing: bool = True,
541
                                        derivativeWeights: tuple = None) -> Image:
542
    """
543
    Compute Jacobian determinant of displacement field.
544
    
545
    Args:
546
        displacementField: Input displacement field
547
        useImageSpacing: Account for image spacing in derivatives
548
        derivativeWeights: Weights for derivative computation
549
        
550
    Returns:
551
        Jacobian determinant image
552
    """
553
```
554

555
### Usage Examples
556

557
#### Basic Rigid Registration
558

559
```python
560
import SimpleITK as sitk
561

562
# Load images
563
fixed_image = sitk.ReadImage('fixed.nii')
564
moving_image = sitk.ReadImage('moving.nii')
565

566
# Create registration method
567
registration = sitk.ImageRegistrationMethod()
568

569
# Set up similarity metric
570
registration.SetMetricAsMeanSquares()
571

572
# Set up optimizer
573
registration.SetOptimizerAsRegularStepGradientDescent(
574
    learningRate=1.0,
575
    minStep=0.001, 
576
    numberOfIterations=200
577
)
578

579
# Set interpolator
580
registration.SetInterpolator(sitk.sitkLinear)
581

582
# Create and initialize transform
583
initial_transform = sitk.CenteredTransformInitializer(
584
    fixed_image, moving_image, 
585
    sitk.Euler3DTransform(), 
586
    sitk.CenteredTransformInitializerFilter.GEOMETRY
587
)
588
registration.SetInitialTransform(initial_transform)
589

590
# Execute registration
591
final_transform = registration.Execute()
592

593
print(f"Final metric value: {registration.GetMetricValue()}")
594
print(f"Optimizer stop condition: {registration.GetOptimizerStopConditionDescription()}")
595

596
# Apply transform to moving image
597
resampled = sitk.Resample(moving_image, final_transform, sitk.sitkLinear, 0.0,
598
                         moving_image.GetPixelID())
599

600
# Save result
601
sitk.WriteImage(resampled, 'registered_result.nii')
602
```
603

604
#### Multi-modal Registration with Mutual Information
605

606
```python
607
import SimpleITK as sitk
608

609
# Load multi-modal images (e.g., T1 and T2 MRI)
610
fixed_image = sitk.ReadImage('t1_image.nii')
611
moving_image = sitk.ReadImage('t2_image.nii')
612

613
# Registration setup
614
registration = sitk.ImageRegistrationMethod()
615

616
# Use mutual information for multi-modal registration
617
registration.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
618
registration.SetMetricSamplingStrategy(registration.RANDOM)
619
registration.SetMetricSamplingPercentage(0.01)  # Sample 1% of pixels
620

621
# Use more robust optimizer
622
registration.SetOptimizerAsGradientDescentLineSearch(
623
    learningRate=1.0,
624
    numberOfIterations=300,
625
    convergenceMinimumValue=1e-6,
626
    convergenceWindowSize=10
627
)
628

629
# Multi-resolution registration
630
registration.SetShrinkFactorsPerLevel([4, 2, 1])
631
registration.SetSmoothingSigmasPerLevel([2, 1, 0])
632
registration.SetSmoothingSigmasAreSpecifiedInPhysicalUnits(True)
633

634
# Initialize with similarity transform
635
initial_transform = sitk.CenteredTransformInitializer(
636
    fixed_image, moving_image,
637
    sitk.Similarity3DTransform(),
638
    sitk.CenteredTransformInitializerFilter.GEOMETRY
639
)
640
registration.SetInitialTransform(initial_transform, inPlace=False)
641

642
# Execute registration
643
final_transform = registration.Execute()
644

645
print(f"Final metric value: {registration.GetMetricValue()}")
646
print(f"Final parameters: {final_transform.GetParameters()}")
647
```
648

649
#### Registration with Monitoring
650

651
```python
652
import SimpleITK as sitk
653

654
def command_iteration(method):
655
    """Callback function to monitor registration progress"""
656
    print(f"Iteration {method.GetOptimizerIteration():3d} = "
657
          f"{method.GetMetricValue():10.5f} : "
658
          f"{method.GetOptimizerPosition()}")
659

660
# Set up registration
661
registration = sitk.ImageRegistrationMethod()
662
registration.SetMetricAsMeanSquares()
663
registration.SetOptimizerAsRegularStepGradientDescent(1.0, 0.001, 200)
664
registration.SetInterpolator(sitk.sitkLinear)
665

666
# Add monitoring callback
667
registration.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration(registration))
668

669
# Set initial transform
670
initial_transform = sitk.Euler3DTransform()
671
registration.SetInitialTransform(initial_transform)
672

673
# Execute with monitoring
674
final_transform = registration.Execute(fixed_image, moving_image)
675
```
676

677
#### Deformable Registration with B-splines
678

679
```python
680
import SimpleITK as sitk
681

682
# Load images
683
fixed_image = sitk.ReadImage('fixed.nii')
684
moving_image = sitk.ReadImage('moving.nii')
685

686
# First perform rigid registration
687
rigid_registration = sitk.ImageRegistrationMethod()
688
rigid_registration.SetMetricAsMeanSquares()
689
rigid_registration.SetOptimizerAsRegularStepGradientDescent(1.0, 0.001, 100)
690
rigid_registration.SetInterpolator(sitk.sitkLinear)
691

692
rigid_transform = sitk.CenteredTransformInitializer(
693
    fixed_image, moving_image, sitk.Euler3DTransform(),
694
    sitk.CenteredTransformInitializerFilter.GEOMETRY
695
)
696
rigid_registration.SetInitialTransform(rigid_transform)
697
rigid_result = rigid_registration.Execute()
698

699
# Apply rigid transform to moving image
700
rigidly_registered = sitk.Resample(moving_image, rigid_result, sitk.sitkLinear, 0.0)
701

702
# Set up B-spline deformable registration
703
bspline_registration = sitk.ImageRegistrationMethod()
704
bspline_registration.SetMetricAsMeanSquares()
705
bspline_registration.SetOptimizerAsLBFGSB(
706
    gradientConvergenceTolerance=1e-5,
707
    numberOfIterations=100
708
)
709
bspline_registration.SetInterpolator(sitk.sitkLinear)
710

711
# Create B-spline transform
712
mesh_size = [8, 8, 8]  # Control point grid size
713
bspline_transform = sitk.BSplineTransformInitializer(fixed_image, mesh_size)
714
bspline_registration.SetInitialTransform(bspline_transform, inPlace=False)
715

716
# Multi-resolution deformable registration
717
bspline_registration.SetShrinkFactorsPerLevel([4, 2, 1])
718
bspline_registration.SetSmoothingSigmasPerLevel([2, 1, 0])
719

720
# Execute deformable registration
721
final_bspline = bspline_registration.Execute(fixed_image, rigidly_registered)
722

723
# Apply final transform
724
final_result = sitk.Resample(moving_image, 
725
                           sitk.CompositeTransform([rigid_result, final_bspline]),
726
                           sitk.sitkLinear, 0.0)
727

728
sitk.WriteImage(final_result, 'deformable_registration_result.nii')
729
```
730

731
#### Registration Evaluation
732

733
```python
734
import SimpleITK as sitk
735

736
# Load ground truth and registered segmentations
737
ground_truth = sitk.ReadImage('ground_truth_segmentation.nii')
738
registered_seg = sitk.ReadImage('registered_segmentation.nii')
739

740
# Compute overlap measures
741
overlap_filter = sitk.LabelOverlapMeasuresImageFilter()
742
overlap_filter.Execute(ground_truth, registered_seg)
743

744
# Get overlap statistics for each label
745
labels = overlap_filter.GetLabels()
746
print("Label\tDice\tJaccard\tVolume Similarity")
747
for label in labels:
748
    dice = overlap_filter.GetDiceCoefficient(label)
749
    jaccard = overlap_filter.GetJaccardCoefficient(label)
750
    vol_sim = overlap_filter.GetVolumeSimilarity(label)
751
    print(f"{label}\t{dice:.3f}\t{jaccard:.3f}\t{vol_sim:.3f}")
752

753
# Compute Hausdorff distance for binary images
754
if len(labels) == 2:  # Binary case
755
    hausdorff_filter = sitk.HausdorffDistanceImageFilter()
756
    hausdorff_filter.Execute(ground_truth, registered_seg)
757
    
758
    print(f"Hausdorff distance: {hausdorff_filter.GetHausdorffDistance():.2f}")
759
    print(f"Average Hausdorff distance: {hausdorff_filter.GetAverageHausdorffDistance():.2f}")
760
```
761

762
#### Displacement Field Analysis
763

764
```python
765
import SimpleITK as sitk
766

767
# Load displacement field from registration
768
displacement_field = sitk.ReadImage('displacement_field.nii')
769

770
# Compute Jacobian determinant
771
jacobian = sitk.DisplacementFieldJacobianDeterminant(displacement_field)
772

773
# Analyze deformation
774
stats = sitk.StatisticsImageFilter()
775
stats.Execute(jacobian)
776

777
print(f"Jacobian statistics:")
778
print(f"  Mean: {stats.GetMean():.3f}")
779
print(f"  Std: {stats.GetSigma():.3f}")
780
print(f"  Min: {stats.GetMinimum():.3f}")
781
print(f"  Max: {stats.GetMaximum():.3f}")
782

783
# Check for folding (negative Jacobian)
784
negative_jacobian = sitk.BinaryThreshold(jacobian, -1e10, 0, 1, 0)
785
negative_count = sitk.StatisticsImageFilter()
786
negative_count.Execute(negative_jacobian)
787

788
print(f"Pixels with folding: {int(negative_count.GetSum())}")
789

790
# Save Jacobian for visualization
791
sitk.WriteImage(jacobian, 'jacobian_determinant.nii')
792
```