0
# NumPy Integration
1

2
Seamless integration utilities for converting between ITK images and NumPy arrays, enabling interoperability with the Python scientific computing ecosystem. This integration allows efficient data exchange and processing workflows that combine ITK's medical image processing capabilities with NumPy's numerical computing functions.
3

4
## Capabilities
5

6
### Array Conversion Functions
7

8
Functions for converting between ITK images and NumPy arrays with proper handling of metadata and memory layout.
9

10
```python { .api }
11
def array_from_image(image: Image) -> numpy.ndarray:
12
    """
13
    Convert ITK image to NumPy array.
14
    
15
    Parameters:
16
    - image: ITK image to convert
17
    
18
    Returns:
19
    - numpy.ndarray: Array with pixel data, axes ordered as [z, y, x] for 3D
20
    
21
    Notes:
22
    - Creates a copy of the image data
23
    - Array dtype matches the ITK pixel type
24
    - Spatial metadata (spacing, origin, direction) is not preserved
25
    """
26

27
def image_from_array(array: numpy.ndarray, is_vector: bool = False) -> Image:
28
    """
29
    Create ITK image from NumPy array.
30
    
31
    Parameters:
32
    - array: NumPy array containing image data
33
    - is_vector: Whether to interpret the array as vector-valued pixels
34
    
35
    Returns:
36
    - Image: ITK image with default spacing, origin, and direction
37
    
38
    Notes:
39
    - Creates a copy of the array data
40
    - Array axes are interpreted as [z, y, x] for 3D
41
    - For vector images, last axis represents vector components
42
    """
43

44
def array_view_from_image(image: Image) -> numpy.ndarray:
45
    """
46
    Create NumPy array view of ITK image data.
47
    
48
    Parameters:
49
    - image: ITK image to view
50
    
51
    Returns:
52
    - numpy.ndarray: Array view sharing memory with ITK image
53
    
54
    Notes:
55
    - No data copying - modifications affect original image
56
    - More memory efficient than array_from_image
57
    - Image must remain alive while view is used
58
    """
59

60
def image_view_from_array(array: numpy.ndarray, is_vector: bool = False) -> Image:
61
    """
62
    Create ITK image view of NumPy array data.
63
    
64
    Parameters:
65
    - array: NumPy array to view
66
    - is_vector: Whether to interpret as vector-valued pixels
67
    
68
    Returns:
69
    - Image: ITK image sharing memory with NumPy array
70
    
71
    Notes:
72
    - No data copying - modifications affect original array
73
    - Array must remain alive while image is used
74
    - Default spacing, origin, and direction are set
75
    """
76
```
77

78
### Metadata Handling
79

80
Functions for preserving and transferring spatial metadata between ITK images and processing workflows.
81

82
```python { .api }
83
def spacing_from_image(image: Image) -> numpy.ndarray:
84
    """Extract spacing information from ITK image as NumPy array."""
85

86
def origin_from_image(image: Image) -> numpy.ndarray:
87
    """Extract origin information from ITK image as NumPy array."""
88

89
def direction_from_image(image: Image) -> numpy.ndarray:
90
    """Extract direction matrix from ITK image as NumPy array."""
91

92
def set_spacing(image: Image, spacing: numpy.ndarray) -> None:
93
    """Set ITK image spacing from NumPy array."""
94

95
def set_origin(image: Image, origin: numpy.ndarray) -> None:
96
    """Set ITK image origin from NumPy array."""
97

98
def set_direction(image: Image, direction: numpy.ndarray) -> None:
99
    """Set ITK image direction matrix from NumPy array."""
100
```
101

102
### Type Instantiation Helpers
103

104
Utility functions for working with ITK's template system from Python.
105

106
```python { .api }
107
def template(template_class):
108
    """
109
    Helper for template class instantiation.
110
    
111
    Parameters:
112
    - template_class: ITK template class
113
    
114
    Returns:
115
    - Template instantiation helper
116
    
117
    Notes:
118
    - Simplifies working with ITK's C++ template system
119
    - Enables dynamic template parameter selection
120
    """
121

122
def auto_pipeline(*args, **kwargs):
123
    """
124
    Automatically construct processing pipelines.
125
    
126
    Parameters:
127
    - args: Pipeline components and parameters
128
    - kwargs: Named parameters for pipeline configuration
129
    
130
    Returns:
131
    - Configured processing pipeline
132
    
133
    Notes:
134
    - Infers appropriate filter types and connections
135
    - Reduces boilerplate code for common operations
136
    """
137

138
def size(sequence) -> Size:
139
    """
140
    Create ITK Size object from Python sequence.
141
    
142
    Parameters:
143
    - sequence: Python list or tuple with size values
144
    
145
    Returns:
146
    - Size: ITK Size object
147
    """
148

149
def index(sequence) -> Index:
150
    """
151
    Create ITK Index object from Python sequence.
152
    
153
    Parameters:
154
    - sequence: Python list or tuple with index values
155
    
156
    Returns:
157
    - Index: ITK Index object
158
    """
159

160
def point(sequence) -> Point:
161
    """
162
    Create ITK Point object from Python sequence.
163
    
164
    Parameters:
165
    - sequence: Python list or tuple with coordinate values
166
    
167
    Returns:
168
    - Point: ITK Point object
169
    """
170

171
def vector(sequence) -> Vector:
172
    """
173
    Create ITK Vector object from Python sequence.
174
    
175
    Parameters:
176
    - sequence: Python list or tuple with vector components
177
    
178
    Returns:
179
    - Vector: ITK Vector object
180
    """
181
```
182

183
### I/O Convenience Functions
184

185
High-level functions for reading and writing images with automatic format detection.
186

187
```python { .api }
188
def imread(filename: str, pixel_type=None) -> Image:
189
    """
190
    Read image from file with automatic format detection.
191
    
192
    Parameters:
193
    - filename: Path to image file
194
    - pixel_type: Optional pixel type specification
195
    
196
    Returns:
197
    - Image: Loaded ITK image
198
    
199
    Notes:
200
    - Automatically detects file format and pixel type
201
    - Supports most medical imaging formats (DICOM, NIfTI, etc.)
202
    - Returns appropriate image type based on file contents
203
    """
204

205
def imwrite(image: Image, filename: str, compression: bool = False) -> None:
206
    """
207
    Write image to file with automatic format detection.
208
    
209
    Parameters:
210
    - image: ITK image to write
211
    - filename: Output file path
212
    - compression: Whether to use compression if supported
213
    
214
    Notes:
215
    - Format determined by file extension
216
    - Automatically handles pixel type conversions
217
    - Preserves spatial metadata when format supports it
218
    """
219

220
def transformread(filename: str) -> Transform:
221
    """
222
    Read transform from file.
223
    
224
    Parameters:
225
    - filename: Path to transform file
226
    
227
    Returns:
228
    - Transform: Loaded ITK transform
229
    """
230

231
def transformwrite(transform: Transform, filename: str) -> None:
232
    """
233
    Write transform to file.
234
    
235
    Parameters:
236
    - transform: ITK transform to write
237
    - filename: Output file path
238
    """
239

240
def meshread(filename: str) -> Mesh:
241
    """
242
    Read mesh from file.
243
    
244
    Parameters:
245
    - filename: Path to mesh file
246
    
247
    Returns:
248
    - Mesh: Loaded ITK mesh
249
    """
250

251
def meshwrite(mesh: Mesh, filename: str) -> None:
252
    """
253
    Write mesh to file.
254
    
255
    Parameters:
256
    - mesh: ITK mesh to write
257
    - filename: Output file path
258
    """
259
```
260

261
### Type System Integration
262

263
Functions for working with ITK's extensive type system from Python.
264

265
```python { .api }
266
# Common type abbreviations for convenience
267
F = float    # 32-bit float
268
D = float    # 64-bit double (Python float)
269
UC = int     # unsigned char
270
US = int     # unsigned short
271
UI = int     # unsigned int
272
UL = int     # unsigned long
273
SC = int     # signed char
274
SS = int     # signed short
275
SI = int     # signed int
276
SL = int     # signed long
277
B = bool     # boolean
278

279
def ctype(typename: str):
280
    """
281
    Get ITK C type from string name.
282
    
283
    Parameters:
284
    - typename: String name of C type ('float', 'double', etc.)
285
    
286
    Returns:
287
    - ITK type object
288
    """
289

290
def output(filter_instance):
291
    """
292
    Get output from ITK filter instance.
293
    
294
    Parameters:
295
    - filter_instance: ITK filter object
296
    
297
    Returns:
298
    - Filter output
299
    
300
    Notes:
301
    - Handles both single and multiple outputs
302
    - Automatically calls Update() if needed
303
    """
304

305
def input(filter_instance, input_data):
306
    """
307
    Set input for ITK filter instance.
308
    
309
    Parameters:
310
    - filter_instance: ITK filter object
311
    - input_data: Input data (image, mesh, etc.)
312
    
313
    Notes:
314
    - Automatically handles type compatibility
315
    - Supports chaining filters
316
    """
317
```
318

319
### Scientific Computing Integration
320

321
Functions for integration with other scientific Python libraries.
322

323
```python { .api }
324
def xarray_from_image(image: Image) -> xarray.DataArray:
325
    """
326
    Convert ITK image to xarray DataArray with coordinate information.
327
    
328
    Parameters:
329
    - image: ITK image to convert
330
    
331
    Returns:
332
    - xarray.DataArray: Array with spatial coordinates and metadata
333
    
334
    Notes:
335
    - Preserves spatial metadata as coordinates
336
    - Enables integration with xarray ecosystem
337
    - Useful for analysis and visualization workflows
338
    """
339

340
def image_from_xarray(data_array: xarray.DataArray) -> Image:
341
    """
342
    Create ITK image from xarray DataArray.
343
    
344
    Parameters:
345
    - data_array: xarray DataArray with spatial coordinates
346
    
347
    Returns:
348
    - Image: ITK image with spatial metadata from coordinates
349
    """
350

351
def scipy_ndimage_compatible(image: Image) -> tuple:
352
    """
353
    Convert ITK image to format compatible with scipy.ndimage.
354
    
355
    Parameters:
356
    - image: ITK image to convert
357
    
358
    Returns:
359
    - tuple: (array, spacing) compatible with scipy.ndimage functions
360
    """
361

362
def pandas_from_pointset(pointset: PointSet) -> pandas.DataFrame:
363
    """
364
    Convert ITK PointSet to pandas DataFrame.
365
    
366
    Parameters:
367
    - pointset: ITK PointSet to convert
368
    
369
    Returns:
370
    - pandas.DataFrame: DataFrame with point coordinates and data
371
    """
372

373
def pointset_from_pandas(dataframe: pandas.DataFrame) -> PointSet:
374
    """
375
    Create ITK PointSet from pandas DataFrame.
376
    
377
    Parameters:
378
    - dataframe: DataFrame with coordinate columns
379
    
380
    Returns:
381
    - PointSet: ITK PointSet object
382
    """
383
```
384

385
## Usage Examples
386

387
### Basic Array Conversion
388

389
```python
390
import itk
391
import numpy as np
392

393
# Create ITK image
394
ImageType = itk.Image[itk.F, 3]
395
image = ImageType.New()
396

397
size = itk.Size[3]([50, 60, 40])
398
region = itk.ImageRegion[3]()
399
region.SetSize(size)
400
image.SetRegions(region)
401
image.SetSpacing([1.5, 1.0, 2.0])
402
image.SetOrigin([10.0, 20.0, 30.0])
403
image.Allocate()
404

405
# Fill with test data
406
for x in range(50):
407
    for y in range(60):
408
        for z in range(40):
409
            idx = itk.Index[3]([x, y, z])
410
            value = float(x + y + z)
411
            image.SetPixel(idx, value)
412

413
print("Created ITK image:")
414
print(f"Size: {image.GetLargestPossibleRegion().GetSize()}")
415
print(f"Spacing: {image.GetSpacing()}")
416
print(f"Origin: {image.GetOrigin()}")
417

418
# Convert to NumPy array
419
array = itk.array_from_image(image)
420

421
print(f"\nNumPy array shape: {array.shape}")
422
print(f"Array dtype: {array.dtype}")
423
print(f"Array min/max: {array.min():.1f}/{array.max():.1f}")
424

425
# Note: ITK uses [x,y,z] indexing, NumPy array is [z,y,x]
426
print(f"ITK pixel at (10,20,15): {image.GetPixel(itk.Index[3]([10, 20, 15]))}")
427
print(f"NumPy value at [15,20,10]: {array[15, 20, 10]}")
428

429
# Convert back to ITK image
430
new_image = itk.image_from_array(array)
431

432
print(f"\nNew ITK image size: {new_image.GetLargestPossibleRegion().GetSize()}")
433
print(f"New image spacing: {new_image.GetSpacing()}")  # Default spacing
434
print(f"New image origin: {new_image.GetOrigin()}")    # Default origin
435

436
# Copy spatial metadata
437
new_image.SetSpacing(image.GetSpacing())
438
new_image.SetOrigin(image.GetOrigin())
439
new_image.SetDirection(image.GetDirection())
440

441
print(f"After copying metadata - spacing: {new_image.GetSpacing()}")
442
```
443

444
### Memory Views for Efficient Processing
445

446
```python
447
import itk
448
import numpy as np
449

450
# Create large ITK image
451
ImageType = itk.Image[itk.F, 3]
452
large_image = ImageType.New()
453

454
size = itk.Size[3]([200, 200, 100])
455
region = itk.ImageRegion[3]()
456
region.SetSize(size)
457
large_image.SetRegions(region)
458
large_image.Allocate()
459

460
# Fill with gradient pattern
461
for x in range(200):
462
    for y in range(200):
463
        for z in range(100):
464
            idx = itk.Index[3]([x, y, z])
465
            value = float(x * y * z) / 1000.0
466
            large_image.SetPixel(idx, value)
467

468
print(f"Created large image: {size[0]}x{size[1]}x{size[2]}")
469

470
# Create view (no memory copying)
471
array_view = itk.array_view_from_image(large_image)
472
print(f"Array view shape: {array_view.shape}")
473
print(f"Memory usage: {array_view.nbytes / (1024**2):.1f} MB")
474

475
# Modify through view (affects original image)
476
original_value = large_image.GetPixel(itk.Index[3]([50, 50, 25]))
477
print(f"Original value at (50,50,25): {original_value}")
478

479
# Modify through NumPy view (z=25, y=50, x=50)
480
array_view[25, 50, 50] = 9999.0
481

482
modified_value = large_image.GetPixel(itk.Index[3]([50, 50, 25]))
483
print(f"Modified value at (50,50,25): {modified_value}")
484

485
# Apply NumPy operations efficiently
486
print("Applying NumPy operations...")
487

488
# Compute statistics
489
mean_val = np.mean(array_view)
490
std_val = np.std(array_view)
491
print(f"Array statistics: mean={mean_val:.2f}, std={std_val:.2f}")
492

493
# Apply smoothing using NumPy/SciPy
494
from scipy import ndimage
495

496
# Smooth a slice
497
smoothed_slice = ndimage.gaussian_filter(array_view[50, :, :], sigma=2.0)
498
array_view[50, :, :] = smoothed_slice
499

500
print("Applied Gaussian smoothing to slice 50")
501

502
# Verify changes in ITK image
503
smoothed_itk_value = large_image.GetPixel(itk.Index[3]([100, 100, 50]))
504
print(f"Smoothed ITK value at (100,100,50): {smoothed_itk_value:.3f}")
505
```
506

507
### Working with Vector Images
508

509
```python
510
import itk
511
import numpy as np
512

513
# Create vector image (e.g., RGB or velocity field)
514
VectorType = itk.Vector[itk.F, 3]
515
VectorImageType = itk.Image[VectorType, 2]
516
vector_image = VectorImageType.New()
517

518
size = itk.Size[2]([100, 100])
519
region = itk.ImageRegion[2]()
520
region.SetSize(size)
521
vector_image.SetRegions(region)
522
vector_image.Allocate()
523

524
# Fill with vector field (e.g., gradient field)
525
for x in range(100):
526
    for y in range(100):
527
        idx = itk.Index[2]([x, y])
528
        
529
        # Create gradient vector
530
        vector = VectorType()
531
        vector.SetElement(0, float(x - 50))     # X component
532
        vector.SetElement(1, float(y - 50))     # Y component  
533
        vector.SetElement(2, float(x + y))      # Z component
534
        
535
        vector_image.SetPixel(idx, vector)
536

537
print("Created vector image")
538

539
# Convert to NumPy array (is_vector=True for proper handling)
540
vector_array = itk.array_from_image(vector_image)
541
print(f"Vector array shape: {vector_array.shape}")  # Should be (100, 100, 3)
542

543
# Access individual components
544
x_component = vector_array[:, :, 0]
545
y_component = vector_array[:, :, 1]
546
z_component = vector_array[:, :, 2]
547

548
print(f"Component shapes: {x_component.shape}")
549

550
# Compute vector magnitude
551
magnitude = np.sqrt(x_component**2 + y_component**2 + z_component**2)
552
print(f"Magnitude range: {magnitude.min():.1f} to {magnitude.max():.1f}")
553

554
# Normalize vectors
555
magnitude_3d = magnitude[:, :, np.newaxis]  # Add dimension for broadcasting
556
normalized_vectors = vector_array / (magnitude_3d + 1e-8)  # Avoid division by zero
557

558
# Create new vector image from normalized data
559
normalized_image = itk.image_from_array(normalized_vectors, is_vector=True)
560

561
# Copy spatial metadata
562
normalized_image.SetSpacing(vector_image.GetSpacing())
563
normalized_image.SetOrigin(vector_image.GetOrigin())
564
normalized_image.SetDirection(vector_image.GetDirection())
565

566
print("Created normalized vector image")
567

568
# Verify normalization
569
test_idx = itk.Index[2]([25, 75])
570
original_vector = vector_image.GetPixel(test_idx)
571
normalized_vector = normalized_image.GetPixel(test_idx)
572

573
orig_mag = np.sqrt(sum(original_vector.GetElement(i)**2 for i in range(3)))
574
norm_mag = np.sqrt(sum(normalized_vector.GetElement(i)**2 for i in range(3)))
575

576
print(f"Original vector magnitude: {orig_mag:.3f}")
577
print(f"Normalized vector magnitude: {norm_mag:.3f}")
578
```
579

580
### High-Level I/O Operations
581

582
```python
583
import itk
584
import numpy as np
585
import tempfile
586
import os
587

588
# Create test image
589
ImageType = itk.Image[itk.US, 3]  # 16-bit unsigned
590
image = ImageType.New()
591

592
size = itk.Size[3]([64, 64, 32])
593
region = itk.ImageRegion[3]()
594
region.SetSize(size)
595
image.SetRegions(region)
596
image.SetSpacing([0.5, 0.5, 1.0])
597
image.SetOrigin([-16.0, -16.0, 0.0])
598
image.Allocate()
599

600
# Fill with synthetic data
601
np.random.seed(42)
602
for x in range(64):
603
    for y in range(64):
604
        for z in range(32):
605
            idx = itk.Index[3]([x, y, z])
606
            # Create pattern with noise
607
            base = int(1000 + 500 * np.sin(x/10.0) * np.cos(y/10.0) * np.sin(z/5.0))
608
            noise = int(np.random.normal(0, 50))
609
            value = max(0, min(65535, base + noise))
610
            image.SetPixel(idx, value)
611

612
print("Created test image for I/O")
613

614
# Create temporary directory for test files
615
with tempfile.TemporaryDirectory() as temp_dir:
616
    
617
    # Test high-level image I/O
618
    nifti_path = os.path.join(temp_dir, "test_image.nii.gz")
619
    
620
    # Write using high-level function
621
    itk.imwrite(image, nifti_path, compression=True)
622
    print(f"Wrote image to {nifti_path}")
623
    
624
    # Read using high-level function
625
    loaded_image = itk.imread(nifti_path)
626
    print(f"Loaded image type: {type(loaded_image)}")
627
    print(f"Loaded image size: {loaded_image.GetLargestPossibleRegion().GetSize()}")
628
    print(f"Loaded image spacing: {loaded_image.GetSpacing()}")
629
    print(f"Loaded image origin: {loaded_image.GetOrigin()}")
630
    
631
    # Verify data integrity
632
    test_idx = itk.Index[3]([32, 32, 16])
633
    original_value = image.GetPixel(test_idx)
634
    loaded_value = loaded_image.GetPixel(test_idx)
635
    print(f"Original value: {original_value}, Loaded value: {loaded_value}")
636
    print(f"Values match: {original_value == loaded_value}")
637
    
638
    # Convert to NumPy for analysis
639
    array = itk.array_from_image(loaded_image)
640
    print(f"Array statistics: min={array.min()}, max={array.max()}, mean={array.mean():.1f}")
641
    
642
    # Test transform I/O
643
    transform = itk.AffineTransform[itk.D, 3].New()
644
    transform.SetIdentity()
645
    transform.Scale(2.0)
646
    transform.SetTranslation(itk.Vector[itk.D, 3]([10.0, 20.0, 30.0]))
647
    
648
    transform_path = os.path.join(temp_dir, "test_transform.tfm")
649
    itk.transformwrite(transform, transform_path)
650
    print(f"Wrote transform to {transform_path}")
651
    
652
    loaded_transform = itk.transformread(transform_path)
653
    print(f"Loaded transform type: {type(loaded_transform)}")
654
    
655
    # Test transform equivalence
656
    test_point = itk.Point[itk.D, 3]([1.0, 2.0, 3.0])
657
    original_result = transform.TransformPoint(test_point)
658
    loaded_result = loaded_transform.TransformPoint(test_point)
659
    
660
    print(f"Original transform result: {original_result}")
661
    print(f"Loaded transform result: {loaded_result}")
662
    
663
    # Check if results are approximately equal
664
    diff = sum((original_result[i] - loaded_result[i])**2 for i in range(3))
665
    print(f"Transform difference: {diff:.10f}")
666
```
667

668
### Integration with Scientific Python Libraries
669

670
```python
671
import itk
672
import numpy as np
673
import matplotlib.pyplot as plt
674
from scipy import ndimage
675
import pandas as pd
676

677
# Create test image
678
ImageType = itk.Image[itk.F, 2]
679
image = ImageType.New()
680

681
size = itk.Size[2]([128, 128])
682
region = itk.ImageRegion[2]()
683
region.SetSize(size)
684
image.SetRegions(region)
685
image.SetSpacing([0.5, 0.5])
686
image.Allocate()
687

688
# Create synthetic data with structures
689
center_x, center_y = 64, 64
690
for x in range(128):
691
    for y in range(128):
692
        idx = itk.Index[2]([x, y])
693
        
694
        # Create multiple circular structures
695
        dist1 = np.sqrt((x - 40)**2 + (y - 40)**2)
696
        dist2 = np.sqrt((x - 90)**2 + (y - 90)**2)
697
        dist3 = np.sqrt((x - 40)**2 + (y - 90)**2)
698
        
699
        value = 0.0
700
        if dist1 < 15:
701
            value += 100.0 * np.exp(-dist1/5.0)
702
        if dist2 < 20:
703
            value += 150.0 * np.exp(-dist2/8.0)
704
        if dist3 < 12:
705
            value += 80.0 * np.exp(-dist3/4.0)
706
            
707
        # Add noise
708
        value += np.random.normal(0, 5)
709
        image.SetPixel(idx, max(0.0, value))
710

711
print("Created synthetic image with multiple structures")
712

713
# Convert to NumPy for processing
714
array = itk.array_from_image(image)
715
print(f"Array shape: {array.shape}, dtype: {array.dtype}")
716

717
# Apply scipy.ndimage operations
718
print("Applying scipy.ndimage operations...")
719

720
# Gaussian smoothing
721
smoothed = ndimage.gaussian_filter(array, sigma=2.0)
722

723
# Edge detection using gradient magnitude
724
grad_x = ndimage.sobel(smoothed, axis=1)
725
grad_y = ndimage.sobel(smoothed, axis=0)
726
gradient_magnitude = np.sqrt(grad_x**2 + grad_y**2)
727

728
# Morphological operations
729
from scipy.ndimage import binary_erosion, binary_dilation
730

731
# Create binary mask
732
threshold = np.mean(smoothed) + 2 * np.std(smoothed)
733
binary_mask = smoothed > threshold
734

735
# Apply morphological operations
736
eroded = binary_erosion(binary_mask, iterations=2)
737
dilated = binary_dilation(eroded, iterations=3)
738

739
print("Applied smoothing, edge detection, and morphological operations")
740

741
# Convert results back to ITK images
742
smoothed_image = itk.image_from_array(smoothed)
743
smoothed_image.CopyInformation(image)
744

745
gradient_image = itk.image_from_array(gradient_magnitude)
746
gradient_image.CopyInformation(image)
747

748
binary_image = itk.image_from_array(dilated.astype(np.float32))
749
binary_image.CopyInformation(image)
750

751
# Extract connected components and create point data
752
from scipy.ndimage import label
753

754
labeled, num_features = label(dilated)
755
print(f"Found {num_features} connected components")
756

757
# Create pandas DataFrame with component analysis
758
components_data = []
759
for label_id in range(1, num_features + 1):
760
    component_mask = labeled == label_id
761
    indices = np.where(component_mask)
762
    
763
    # Calculate properties
764
    area = np.sum(component_mask)
765
    centroid_y = np.mean(indices[0])  # NumPy y-axis
766
    centroid_x = np.mean(indices[1])  # NumPy x-axis
767
    
768
    # Convert to physical coordinates using image spacing and origin
769
    spacing = image.GetSpacing()
770
    origin = image.GetOrigin()
771
    
772
    phys_x = centroid_x * spacing[0] + origin[0]
773
    phys_y = centroid_y * spacing[1] + origin[1]
774
    
775
    # Calculate intensity statistics
776
    component_intensities = smoothed[component_mask]
777
    mean_intensity = np.mean(component_intensities)
778
    max_intensity = np.max(component_intensities)
779
    
780
    components_data.append({
781
        'label': label_id,
782
        'area_pixels': area,
783
        'area_physical': area * spacing[0] * spacing[1],
784
        'centroid_x': phys_x,
785
        'centroid_y': phys_y,
786
        'mean_intensity': mean_intensity,
787
        'max_intensity': max_intensity
788
    })
789

790
# Create pandas DataFrame
791
df = pd.DataFrame(components_data)
792
print("\nComponent analysis results:")
793
print(df)
794

795
# Create ITK PointSet from centroids
796
PointSetType = itk.PointSet[itk.F, 2]
797
pointset = PointSetType.New()
798

799
points = pointset.GetPoints()
800
point_data = pointset.GetPointData()
801

802
for i, row in df.iterrows():
803
    point = itk.Point[itk.F, 2]([row['centroid_x'], row['centroid_y']])
804
    points.InsertElement(i, point)
805
    point_data.InsertElement(i, row['mean_intensity'])
806

807
print(f"Created PointSet with {pointset.GetNumberOfPoints()} points")
808

809
# Display summary statistics
810
print(f"\nSummary statistics:")
811
print(f"Total area: {df['area_physical'].sum():.2f} square units")
812
print(f"Average component area: {df['area_physical'].mean():.2f} ± {df['area_physical'].std():.2f}")
813
print(f"Intensity range: {df['mean_intensity'].min():.1f} to {df['mean_intensity'].max():.1f}")
814

815
# Optional: Create matplotlib visualization if available
816
try:
817
    fig, axes = plt.subplots(2, 2, figsize=(12, 12))
818
    
819
    axes[0, 0].imshow(array, cmap='gray')
820
    axes[0, 0].set_title('Original Image')
821
    
822
    axes[0, 1].imshow(smoothed, cmap='gray')
823
    axes[0, 1].set_title('Smoothed Image')
824
    
825
    axes[1, 0].imshow(gradient_magnitude, cmap='hot')
826
    axes[1, 0].set_title('Gradient Magnitude')
827
    
828
    axes[1, 1].imshow(labeled, cmap='tab10')
829
    axes[1, 1].set_title(f'Connected Components ({num_features})')
830
    
831
    # Add centroids to component plot
832
    axes[1, 1].scatter(df['centroid_x']/spacing[0] - origin[0]/spacing[0], 
833
                      df['centroid_y']/spacing[1] - origin[1]/spacing[1], 
834
                      c='white', s=50, marker='x')
835
    
836
    plt.tight_layout()
837
    plt.savefig('itk_numpy_integration_demo.png', dpi=150, bbox_inches='tight')
838
    print("Saved visualization to 'itk_numpy_integration_demo.png'")
839
    
840
except ImportError:
841
    print("Matplotlib not available for visualization")
842

843
print("ITK-NumPy integration example completed successfully")
844
```
845

846
### Working with Template Classes
847

848
```python
849
import itk
850

851
# Demonstrate template class usage with helper functions
852
print("ITK Template System Integration Examples")
853
print("=" * 50)
854

855
# Using template helper for dynamic type selection
856
pixel_types = [itk.UC, itk.US, itk.F, itk.D]
857
dimensions = [2, 3]
858

859
for pixel_type in pixel_types:
860
    for dim in dimensions:
861
        # Use template helper to create image type
862
        ImageType = itk.Image[pixel_type, dim]
863
        
864
        # Create instance
865
        image = ImageType.New()
866
        
867
        # Set basic properties
868
        size_list = [64] * dim
869
        size = itk.size(size_list)  # Helper function
870
        
871
        region = itk.ImageRegion[dim]()
872
        region.SetSize(size)
873
        image.SetRegions(region)
874
        
875
        spacing_list = [1.0] * dim
876
        spacing = itk.vector(spacing_list)  # Helper function
877
        image.SetSpacing(spacing)
878
        
879
        origin_list = [0.0] * dim  
880
        origin = itk.point(origin_list)  # Helper function
881
        image.SetOrigin(origin)
882
        
883
        image.Allocate()
884
        
885
        print(f"Created {dim}D image with pixel type {pixel_type.__name__}")
886
        print(f"  Size: {image.GetLargestPossibleRegion().GetSize()}")
887
        print(f"  Spacing: {image.GetSpacing()}")
888
        print(f"  Memory size: {image.GetLargestPossibleRegion().GetNumberOfPixels() * image.GetLargestPossibleRegion().GetSize().GetSizeDimension()} elements")
889
        print()
890

891
# Demonstrate auto_pipeline functionality
892
print("Auto-pipeline example:")
893

894
# Create source image
895
ImageType = itk.Image[itk.F, 2]
896
source_image = ImageType.New()
897

898
size = itk.size([100, 100])
899
region = itk.ImageRegion[2]()
900
region.SetSize(size)
901
source_image.SetRegions(region)
902
source_image.Allocate()
903

904
# Fill with test pattern
905
array = itk.array_view_from_image(source_image)
906
x, y = np.meshgrid(range(100), range(100), indexing='ij')
907
array[:, :] = np.sin(x/10.0) * np.cos(y/10.0) * 100 + 100
908

909
# Create processing pipeline using template system
910
print("Building processing pipeline...")
911

912
# Gaussian smoothing
913
GaussianFilterType = itk.SmoothingRecursiveGaussianImageFilter[ImageType, ImageType]
914
gaussian_filter = GaussianFilterType.New()
915
gaussian_filter.SetInput(source_image)
916
gaussian_filter.SetSigma(2.0)
917

918
# Gradient magnitude
919
GradientFilterType = itk.GradientMagnitudeImageFilter[ImageType, ImageType]
920
gradient_filter = GradientFilterType.New()
921
gradient_filter.SetInput(gaussian_filter.GetOutput())
922

923
# Binary threshold
924
ThresholdFilterType = itk.BinaryThresholdImageFilter[ImageType, itk.Image[itk.UC, 2]]
925
threshold_filter = ThresholdFilterType.New()
926
threshold_filter.SetInput(gradient_filter.GetOutput())
927
threshold_filter.SetLowerThreshold(50.0)
928
threshold_filter.SetUpperThreshold(1000.0)
929
threshold_filter.SetInsideValue(255)
930
threshold_filter.SetOutsideValue(0)
931

932
# Execute pipeline
933
threshold_filter.Update()
934
result = threshold_filter.GetOutput()
935

936
print(f"Pipeline executed successfully")
937
print(f"Result image size: {result.GetLargestPossibleRegion().GetSize()}")
938
print(f"Result pixel type: {type(result).__name__}")
939

940
# Convert result to NumPy for analysis
941
result_array = itk.array_from_image(result)
942
unique_values = np.unique(result_array)
943
print(f"Unique values in result: {unique_values}")
944
print(f"Foreground pixels: {np.sum(result_array == 255)}")
945
print(f"Background pixels: {np.sum(result_array == 0)}")
946

947
print("\nTemplate system integration examples completed")
948
```