0
# Feature Detection
1
2
Advanced feature detection and analysis techniques for identifying anatomical landmarks, tissue boundaries, and structural patterns in medical and scientific images using SimpleITK.
3
4
## Capabilities
5
6
### Edge and Boundary Detection
7
8
Sophisticated edge detection algorithms for identifying boundaries and transitions in medical images.
9
10
```python { .api }
11
def CannyEdgeDetection(image: Image, lowerThreshold: float = 0.0,
12
upperThreshold: float = 0.0, variance: list = [2.0, 2.0]) -> Image:
13
"""
14
Canny edge detection with hysteresis thresholding
15
16
Args:
17
image: Input grayscale image
18
lowerThreshold: Lower threshold for edge linking
19
upperThreshold: Upper threshold for strong edges
20
variance: Gaussian smoothing variance per dimension
21
22
Returns:
23
Binary edge map with connected edge pixels
24
"""
25
26
def ZeroCrossingBasedEdgeDetection(image: Image, variance: list = [1.0, 1.0],
27
foregroundValue: float = 1.0,
28
backgroundValue: float = 0.0,
29
maximumError: float = 0.1) -> Image:
30
"""
31
Zero-crossing edge detection (Laplacian of Gaussian)
32
33
Args:
34
image: Input grayscale image
35
variance: Gaussian kernel variance per dimension
36
foregroundValue: Value for detected edges
37
backgroundValue: Value for non-edge pixels
38
maximumError: Maximum approximation error
39
40
Returns:
41
Binary edge image from zero-crossings
42
"""
43
44
def LaplacianImageFilter(image: Image, useImageSpacing: bool = True) -> Image:
45
"""
46
Laplacian edge detection (second derivative)
47
48
Args:
49
image: Input image
50
useImageSpacing: Use image spacing in computation
51
52
Returns:
53
Laplacian filtered image highlighting edges
54
"""
55
56
def RecursiveGaussianImageFilter(image: Image, sigma: float = 1.0,
57
direction: int = 0, order: int = 0,
58
normalizeAcrossScale: bool = False) -> Image:
59
"""
60
Recursive Gaussian filtering with derivative computation
61
62
Args:
63
image: Input image
64
sigma: Standard deviation of Gaussian kernel
65
direction: Direction for derivative (0, 1, 2 for x, y, z)
66
order: Derivative order (0=smoothing, 1=first derivative, 2=second derivative)
67
normalizeAcrossScale: Normalize across different scales
68
69
Returns:
70
Gaussian filtered or derivative image
71
"""
72
```
73
74
**Usage Examples:**
75
76
```python
77
import SimpleITK as sitk
78
import numpy as np
79
80
# Load medical image
81
image = sitk.ReadImage('brain_mri.nii', sitk.sitkFloat32)
82
83
# Canny edge detection
84
canny_edges = sitk.CannyEdgeDetection(image,
85
lowerThreshold=10.0,
86
upperThreshold=50.0,
87
variance=[1.5, 1.5, 1.0])
88
89
# Zero-crossing edge detection (LoG)
90
log_edges = sitk.ZeroCrossingBasedEdgeDetection(image,
91
variance=[2.0, 2.0, 1.5],
92
maximumError=0.05)
93
94
# Multi-scale edge detection
95
scales = [1.0, 2.0, 4.0]
96
multiscale_edges = []
97
98
for sigma in scales:
99
# Compute image gradients
100
grad_x = sitk.RecursiveGaussianImageFilter(image, sigma=sigma,
101
direction=0, order=1)
102
grad_y = sitk.RecursiveGaussianImageFilter(image, sigma=sigma,
103
direction=1, order=1)
104
105
# Gradient magnitude
106
grad_mag = sitk.Sqrt(sitk.Add(sitk.Square(grad_x), sitk.Square(grad_y)))
107
multiscale_edges.append(grad_mag)
108
109
# Combine multi-scale edges
110
combined_edges = multiscale_edges[0]
111
for edge_img in multiscale_edges[1:]:
112
combined_edges = sitk.Maximum(combined_edges, edge_img)
113
114
# Non-maximum suppression for thin edges
115
def non_maximum_suppression(gradient_mag, gradient_direction):
116
"""Simple non-maximum suppression"""
117
# Convert to numpy for processing
118
mag_array = sitk.GetArrayFromImage(gradient_mag)
119
dir_array = sitk.GetArrayFromImage(gradient_direction)
120
121
# Simplified NMS implementation
122
suppressed = np.zeros_like(mag_array)
123
124
for i in range(1, mag_array.shape[0] - 1):
125
for j in range(1, mag_array.shape[1] - 1):
126
angle = dir_array[i, j]
127
128
# Determine neighboring pixels based on gradient direction
129
if (0 <= angle < 22.5) or (157.5 <= angle <= 180):
130
# Horizontal edge
131
neighbors = [mag_array[i, j-1], mag_array[i, j+1]]
132
elif 22.5 <= angle < 67.5:
133
# Diagonal edge (/)
134
neighbors = [mag_array[i-1, j-1], mag_array[i+1, j+1]]
135
elif 67.5 <= angle < 112.5:
136
# Vertical edge
137
neighbors = [mag_array[i-1, j], mag_array[i+1, j]]
138
else:
139
# Diagonal edge (\)
140
neighbors = [mag_array[i-1, j+1], mag_array[i+1, j-1]]
141
142
# Keep pixel if it's local maximum
143
if mag_array[i, j] >= max(neighbors):
144
suppressed[i, j] = mag_array[i, j]
145
146
result = sitk.GetImageFromArray(suppressed)
147
result.CopyInformation(gradient_mag)
148
return result
149
150
# Apply non-maximum suppression
151
gradient_direction = sitk.Atan2(grad_y, grad_x) # Compute gradient direction
152
thin_edges = non_maximum_suppression(combined_edges, gradient_direction)
153
154
# Display results
155
sitk.Show(image, "Original")
156
sitk.Show(canny_edges, "Canny Edges")
157
sitk.Show(combined_edges, "Multi-scale Edges")
158
sitk.Show(thin_edges, "Thinned Edges")
159
```
160
161
### Corner and Interest Point Detection
162
163
Detect and characterize corner points and interest points for landmark identification.
164
165
```python { .api }
166
def HarrisCornerDetection(image: Image, sigma: float = 1.0, alpha: float = 0.04,
167
threshold: float = 0.01, suppressionRadius: int = 3) -> Image:
168
"""
169
Harris corner detection algorithm
170
171
Args:
172
image: Input grayscale image
173
sigma: Gaussian smoothing parameter
174
alpha: Harris detector free parameter (typically 0.04-0.06)
175
threshold: Corner response threshold
176
suppressionRadius: Non-maximum suppression radius
177
178
Returns:
179
Corner response image or binary corner map
180
"""
181
182
def ShiTomasiCornerDetection(image: Image, sigma: float = 1.0,
183
threshold: float = 0.01, suppressionRadius: int = 3) -> Image:
184
"""
185
Shi-Tomasi corner detection (modified Harris)
186
187
Args:
188
image: Input grayscale image
189
sigma: Gaussian smoothing parameter
190
threshold: Corner response threshold
191
suppressionRadius: Non-maximum suppression radius
192
193
Returns:
194
Corner strength image
195
"""
196
197
def NobleCornerDetection(image: Image, sigma: float = 1.0,
198
threshold: float = 0.01, suppressionRadius: int = 3) -> Image:
199
"""
200
Noble corner detection algorithm
201
202
Args:
203
image: Input grayscale image
204
sigma: Gaussian smoothing parameter
205
threshold: Corner response threshold
206
suppressionRadius: Non-maximum suppression radius
207
208
Returns:
209
Corner response image
210
"""
211
```
212
213
**Usage Examples:**
214
215
```python
216
import SimpleITK as sitk
217
import numpy as np
218
219
def extract_corner_points(image, method='harris'):
220
"""Extract corner points from image"""
221
222
# Preprocessing
223
if image.GetPixelID() != sitk.sitkFloat32:
224
image = sitk.Cast(image, sitk.sitkFloat32)
225
226
# Smooth image to reduce noise
227
smoothed = sitk.Gaussian(image, sigma=1.0)
228
229
# Compute image gradients
230
grad_x = sitk.Derivative(smoothed, direction=0, order=1)
231
grad_y = sitk.Derivative(smoothed, direction=1, order=1)
232
233
# Compute structure tensor components
234
Ixx = sitk.Gaussian(sitk.Square(grad_x), sigma=1.5)
235
Iyy = sitk.Gaussian(sitk.Square(grad_y), sigma=1.5)
236
Ixy = sitk.Gaussian(sitk.Multiply(grad_x, grad_y), sigma=1.5)
237
238
if method == 'harris':
239
# Harris corner response
240
det = sitk.Subtract(sitk.Multiply(Ixx, Iyy), sitk.Square(Ixy))
241
trace = sitk.Add(Ixx, Iyy)
242
harris_response = sitk.Subtract(det, sitk.Multiply(0.04, sitk.Square(trace)))
243
corner_response = harris_response
244
245
elif method == 'shi_tomasi':
246
# Shi-Tomasi corner response (minimum eigenvalue)
247
# λ_min = 0.5 * (trace - sqrt(trace² - 4*det))
248
det = sitk.Subtract(sitk.Multiply(Ixx, Iyy), sitk.Square(Ixy))
249
trace = sitk.Add(Ixx, Iyy)
250
discriminant = sitk.Subtract(sitk.Square(trace), sitk.Multiply(4.0, det))
251
sqrt_discriminant = sitk.Sqrt(sitk.Maximum(discriminant, 0.0))
252
min_eigenvalue = sitk.Multiply(0.5, sitk.Subtract(trace, sqrt_discriminant))
253
corner_response = min_eigenvalue
254
255
# Threshold corner response
256
threshold = 0.01 * sitk.GetArrayFromImage(corner_response).max()
257
strong_corners = sitk.BinaryThreshold(corner_response, threshold, 1e6, 1, 0)
258
259
return corner_response, strong_corners
260
261
def non_maximum_suppression_corners(corner_response, radius=3):
262
"""Apply non-maximum suppression to corner points"""
263
264
# Convert to numpy
265
response_array = sitk.GetArrayFromImage(corner_response)
266
suppressed = np.zeros_like(response_array)
267
268
# Find local maxima
269
for i in range(radius, response_array.shape[0] - radius):
270
for j in range(radius, response_array.shape[1] - radius):
271
local_region = response_array[i-radius:i+radius+1, j-radius:j+radius+1]
272
if response_array[i, j] == local_region.max() and response_array[i, j] > 0:
273
suppressed[i, j] = response_array[i, j]
274
275
result = sitk.GetImageFromArray(suppressed)
276
result.CopyInformation(corner_response)
277
return result
278
279
# Load image
280
image = sitk.ReadImage('building.png', sitk.sitkFloat32)
281
282
# Extract corners using different methods
283
harris_response, harris_corners = extract_corner_points(image, method='harris')
284
shitomasi_response, shitomasi_corners = extract_corner_points(image, method='shi_tomasi')
285
286
# Apply non-maximum suppression
287
harris_nms = non_maximum_suppression_corners(harris_response, radius=5)
288
shitomasi_nms = non_maximum_suppression_corners(shitomasi_response, radius=5)
289
290
# Extract corner coordinates
291
def get_corner_coordinates(corner_image, threshold=0.01):
292
"""Extract corner point coordinates"""
293
array = sitk.GetArrayFromImage(corner_image)
294
max_val = array.max()
295
296
if max_val > 0:
297
thresh = threshold * max_val
298
y_coords, x_coords = np.where(array > thresh)
299
300
# Convert to physical coordinates
301
corners = []
302
for y, x in zip(y_coords, x_coords):
303
point = corner_image.TransformIndexToPhysicalPoint([int(x), int(y)])
304
corners.append(point)
305
306
return corners
307
308
return []
309
310
harris_points = get_corner_coordinates(harris_nms)
311
shitomasi_points = get_corner_coordinates(shitomasi_nms)
312
313
print(f"Harris corners found: {len(harris_points)}")
314
print(f"Shi-Tomasi corners found: {len(shitomasi_points)}")
315
316
# Display results
317
sitk.Show(image, "Original")
318
sitk.Show(harris_response, "Harris Response")
319
sitk.Show(harris_nms, "Harris NMS")
320
```
321
322
### Blob and Region Detection
323
324
Detect blob-like structures and regions of interest in medical images.
325
326
```python { .api }
327
def LaplacianOfGaussian(image: Image, sigma: float = 1.0,
328
useImageSpacing: bool = True) -> Image:
329
"""
330
Laplacian of Gaussian blob detection
331
332
Args:
333
image: Input grayscale image
334
sigma: Standard deviation of Gaussian kernel
335
useImageSpacing: Use image spacing in computation
336
337
Returns:
338
LoG filtered image (blob response)
339
"""
340
341
def DifferenceOfGaussians(image: Image, sigma1: float = 1.0, sigma2: float = 2.0) -> Image:
342
"""
343
Difference of Gaussians blob detection
344
345
Args:
346
image: Input image
347
sigma1: Standard deviation of first Gaussian
348
sigma2: Standard deviation of second Gaussian
349
350
Returns:
351
DoG filtered image highlighting blobs
352
"""
353
354
def HessianRecursiveGaussian(image: Image, sigma: float = 1.0,
355
normalizeAcrossScale: bool = False) -> Image:
356
"""
357
Hessian-based feature detection
358
359
Args:
360
image: Input image
361
sigma: Gaussian kernel standard deviation
362
normalizeAcrossScale: Normalize response across scales
363
364
Returns:
365
Hessian determinant image (blob/ridge response)
366
"""
367
368
def MultiScaleHessianBasedObjectness(image: Image, sigmaMin: float = 1.0,
369
sigmaMax: float = 8.0, numberOfSigmaSteps: int = 10,
370
alpha: float = 0.5, beta: float = 0.5,
371
gamma: float = 5.0, brightObject: bool = True,
372
objectDimension: int = 1) -> Image:
373
"""
374
Multi-scale Hessian-based objectness measure
375
376
Args:
377
image: Input image
378
sigmaMin: Minimum scale (sigma)
379
sigmaMax: Maximum scale (sigma)
380
numberOfSigmaSteps: Number of scale steps
381
alpha: Plate-like vs blob-like discrimination (typically 0.5)
382
beta: Blob-like vs line-like discrimination (typically 0.5)
383
gamma: Background suppression (typically 5-25)
384
brightObject: Whether to detect bright objects on dark background
385
objectDimension: Expected object dimension (0=points, 1=lines, 2=surfaces)
386
387
Returns:
388
Objectness measure highlighting target structures
389
"""
390
```
391
392
**Usage Examples:**
393
394
```python
395
import SimpleITK as sitk
396
import numpy as np
397
398
def multiscale_blob_detection(image, scales=None):
399
"""Multi-scale blob detection using LoG and DoG"""
400
401
if scales is None:
402
scales = np.logspace(0, 2, 10) # Log-spaced from 1 to 100
403
404
log_responses = []
405
dog_responses = []
406
407
for sigma in scales:
408
# Laplacian of Gaussian
409
log_response = sitk.LaplacianRecursiveGaussian(image, sigma=sigma)
410
# Normalize by scale
411
log_normalized = sitk.Multiply(log_response, sigma * sigma)
412
log_responses.append(log_normalized)
413
414
# Difference of Gaussians
415
sigma2 = sigma * 1.6 # Standard ratio
416
dog_response = sitk.DifferenceOfGaussians(image, sigma1=sigma, sigma2=sigma2)
417
dog_responses.append(dog_response)
418
419
# Find maximum response across scales
420
max_log_response = log_responses[0]
421
max_dog_response = dog_responses[0]
422
423
for log_resp, dog_resp in zip(log_responses[1:], dog_responses[1:]):
424
max_log_response = sitk.Maximum(max_log_response, sitk.Abs(log_resp))
425
max_dog_response = sitk.Maximum(max_dog_response, sitk.Abs(dog_resp))
426
427
return max_log_response, max_dog_response, log_responses, dog_responses
428
429
def vessel_enhancement_filter(image, scales=None):
430
"""Multi-scale vessel enhancement using Hessian eigenvalues"""
431
432
if scales is None:
433
scales = [1.0, 2.0, 4.0, 8.0]
434
435
vessel_responses = []
436
437
for sigma in scales:
438
# Compute Hessian matrix eigenvalues
439
hessian = sitk.HessianRecursiveGaussian(image, sigma=sigma)
440
441
# For 2D: compute vesselness measure
442
# λ1, λ2 are eigenvalues with |λ1| ≤ |λ2|
443
# Vesselness = exp(-Rb²/2β²) * (1 - exp(-S²/2γ²))
444
# where Rb = |λ1|/|λ2|, S = sqrt(λ1² + λ2²)
445
446
# Simplified vesselness computation
447
vessel_response = sitk.HessianRecursiveGaussian(image, sigma=sigma)
448
vessel_response = sitk.Multiply(vessel_response, sigma * sigma) # Scale normalization
449
450
# Keep only negative values (dark vessels on bright background)
451
vessel_response = sitk.Clamp(vessel_response, upperBound=0.0)
452
vessel_response = sitk.Abs(vessel_response)
453
454
vessel_responses.append(vessel_response)
455
456
# Maximum response across scales
457
max_vessel_response = vessel_responses[0]
458
for resp in vessel_responses[1:]:
459
max_vessel_response = sitk.Maximum(max_vessel_response, resp)
460
461
return max_vessel_response
462
463
def detect_and_characterize_blobs(image, min_sigma=1.0, max_sigma=10.0, threshold=0.01):
464
"""Detect and characterize blob-like structures"""
465
466
# Multi-scale LoG detection
467
scales = np.logspace(np.log10(min_sigma), np.log10(max_sigma), 8)
468
responses = []
469
scale_map = sitk.Image(image.GetSize(), sitk.sitkFloat32)
470
scale_map.CopyInformation(image)
471
472
# Compute LoG at multiple scales
473
for i, sigma in enumerate(scales):
474
log_response = sitk.LaplacianRecursiveGaussian(image, sigma=sigma)
475
# Scale-normalized response
476
normalized_response = sitk.Multiply(log_response, sigma * sigma)
477
responses.append((normalized_response, sigma))
478
479
# Find scale-space maxima
480
max_response = sitk.Image(image.GetSize(), sitk.sitkFloat32)
481
max_response.CopyInformation(image)
482
483
response_arrays = [sitk.GetArrayFromImage(resp[0]) for resp in responses]
484
scale_values = [resp[1] for resp in responses]
485
486
max_array = sitk.GetArrayFromImage(max_response)
487
scale_array = sitk.GetArrayFromImage(scale_map)
488
489
# For each pixel, find the scale with maximum response
490
for i in range(len(response_arrays[0])):
491
for j in range(len(response_arrays[0][0]) if len(response_arrays[0].shape) > 1 else [0]):
492
if len(response_arrays[0].shape) > 2:
493
for k in range(len(response_arrays[0][0][0])):
494
pixel_responses = [arr[i, j, k] for arr in response_arrays]
495
max_idx = np.argmax(np.abs(pixel_responses))
496
max_array[i, j, k] = pixel_responses[max_idx]
497
scale_array[i, j, k] = scale_values[max_idx]
498
else:
499
if len(response_arrays[0].shape) > 1:
500
pixel_responses = [arr[i, j] for arr in response_arrays]
501
max_idx = np.argmax(np.abs(pixel_responses))
502
max_array[i, j] = pixel_responses[max_idx]
503
scale_array[i, j] = scale_values[max_idx]
504
505
max_response = sitk.GetImageFromArray(max_array)
506
max_response.CopyInformation(image)
507
508
scale_map = sitk.GetImageFromArray(scale_array)
509
scale_map.CopyInformation(image)
510
511
# Threshold and extract blob locations
512
abs_response = sitk.Abs(max_response)
513
threshold_value = threshold * sitk.GetArrayFromImage(abs_response).max()
514
blob_mask = sitk.BinaryThreshold(abs_response, threshold_value, 1e6, 1, 0)
515
516
return max_response, scale_map, blob_mask
517
518
# Apply blob detection
519
image = sitk.ReadImage('spots.tif', sitk.sitkFloat32)
520
521
# Multi-scale blob detection
522
log_max, dog_max, log_scales, dog_scales = multiscale_blob_detection(image)
523
524
# Vessel enhancement
525
vessels = vessel_enhancement_filter(image)
526
527
# Detailed blob characterization
528
blob_response, blob_scales, blob_mask = detect_and_characterize_blobs(image,
529
min_sigma=2.0,
530
max_sigma=20.0,
531
threshold=0.1)
532
533
# Display results
534
sitk.Show(image, "Original")
535
sitk.Show(log_max, "LoG Blobs")
536
sitk.Show(vessels, "Vessel Enhancement")
537
sitk.Show(blob_mask, "Detected Blobs")
538
sitk.Show(blob_scales, "Blob Scales")
539
```
540
541
### Texture and Pattern Analysis
542
543
Analyze texture patterns and local image structure for tissue characterization.
544
545
```python { .api }
546
def ScalarImageToTextureFeaturesFilter(image: Image, mask: Image = None,
547
insideValue: int = 1, numberOfBinsPerAxis: int = 256,
548
pixelValueMinimum: float = None,
549
pixelValueMaximum: float = None,
550
distanceValueMinimum: int = 0,
551
distanceValueMaximum: int = 1,
552
requestedFeatures: list = None) -> dict:
553
"""
554
Compute texture features from gray-level co-occurrence matrix
555
556
Args:
557
image: Input grayscale image
558
mask: Binary mask defining region of interest
559
insideValue: Mask value representing valid region
560
numberOfBinsPerAxis: Number of intensity bins for GLCM
561
pixelValueMinimum: Minimum pixel value for binning
562
pixelValueMaximum: Maximum pixel value for binning
563
distanceValueMinimum: Minimum distance for co-occurrence
564
distanceValueMaximum: Maximum distance for co-occurrence
565
requestedFeatures: List of features to compute
566
567
Returns:
568
Dictionary of computed texture features
569
"""
570
571
def LocalBinaryPattern(image: Image, radius: int = 1, points: int = 8,
572
method: str = 'uniform') -> Image:
573
"""
574
Local Binary Pattern texture analysis
575
576
Args:
577
image: Input grayscale image
578
radius: Radius of circular neighborhood
579
points: Number of sample points on circle
580
method: LBP variant ('uniform', 'default', 'ror', 'var')
581
582
Returns:
583
LBP pattern image
584
"""
585
586
def GaborFilter(image: Image, frequency: float = 0.1, theta: float = 0.0,
587
sigma_x: float = 2.0, sigma_y: float = 2.0,
588
offset: float = 0.0) -> Image:
589
"""
590
Gabor filter for texture analysis
591
592
Args:
593
image: Input image
594
frequency: Spatial frequency of harmonic function
595
theta: Orientation angle of parallel stripes
596
sigma_x: Standard deviation in x direction
597
sigma_y: Standard deviation in y direction
598
offset: Phase offset
599
600
Returns:
601
Gabor filtered image
602
"""
603
```
604
605
**Usage Examples:**
606
607
```python
608
import SimpleITK as sitk
609
import numpy as np
610
611
def comprehensive_texture_analysis(image, mask=None):
612
"""Comprehensive texture feature extraction"""
613
614
features = {}
615
616
# GLCM-based features
617
if mask is None:
618
mask = sitk.Image(image.GetSize(), sitk.sitkUInt8)
619
mask.CopyInformation(image)
620
mask = sitk.Add(mask, 1) # All ones mask
621
622
# Compute GLCM features
623
glcm_filter = sitk.ScalarImageToTextureFeaturesFilter()
624
glcm_filter.SetMaskImage(mask)
625
glcm_filter.SetNumberOfBinsPerAxis(64)
626
glcm_filter.Execute(image)
627
628
# Extract key GLCM features
629
features['energy'] = glcm_filter.GetEnergy()
630
features['entropy'] = glcm_filter.GetEntropy()
631
features['correlation'] = glcm_filter.GetCorrelation()
632
features['inertia'] = glcm_filter.GetInertia()
633
features['inverse_difference_moment'] = glcm_filter.GetInverseDifferenceMoment()
634
features['haralick_correlation'] = glcm_filter.GetHaralickCorrelation()
635
features['cluster_shade'] = glcm_filter.GetClusterShade()
636
features['cluster_prominence'] = glcm_filter.GetClusterProminence()
637
638
# Multi-scale Gabor filtering
639
gabor_responses = []
640
frequencies = [0.05, 0.1, 0.2, 0.4]
641
orientations = [0, np.pi/4, np.pi/2, 3*np.pi/4]
642
643
for freq in frequencies:
644
for theta in orientations:
645
try:
646
gabor_real = sitk.GaborImageSource(image.GetSize(),
647
frequency=freq,
648
theta=theta,
649
sigma=[2.0, 2.0])
650
gabor_real.CopyInformation(image)
651
652
# Convolve with image
653
gabor_response = sitk.Convolution(image, gabor_real)
654
gabor_responses.append(gabor_response)
655
except:
656
# Fallback: create simple oriented filter
657
pass
658
659
# Local statistics
660
mean_filter = sitk.Mean(image, radius=[5, 5])
661
variance_filter = sitk.Variance(image, radius=[5, 5])
662
663
features['local_mean'] = mean_filter
664
features['local_variance'] = variance_filter
665
features['gabor_responses'] = gabor_responses
666
667
return features
668
669
def Laws_texture_measures(image):
670
"""Laws texture energy measures"""
671
672
# Laws kernels (simplified 1D versions)
673
L5 = [1, 4, 6, 4, 1] # Level
674
E5 = [-1, -2, 0, 2, 1] # Edge
675
S5 = [-1, 0, 2, 0, -1] # Spot
676
W5 = [-1, 2, 0, -2, 1] # Wave
677
R5 = [1, -4, 6, -4, 1] # Ripple
678
679
kernels_1d = {'L': L5, 'E': E5, 'S': S5, 'W': W5, 'R': R5}
680
681
# Create 2D kernels by outer product
682
kernels_2d = {}
683
for name1, k1 in kernels_1d.items():
684
for name2, k2 in kernels_1d.items():
685
kernel_2d = np.outer(k1, k2) / 25.0 # Normalize
686
kernel_name = name1 + name2
687
688
# Convert to SimpleITK image
689
kernel_image = sitk.GetImageFromArray(kernel_2d.astype(np.float32))
690
kernels_2d[kernel_name] = kernel_image
691
692
# Apply kernels and compute energy
693
texture_measures = {}
694
695
for name, kernel in kernels_2d.items():
696
# Convolve image with kernel
697
filtered = sitk.Convolution(image, kernel)
698
699
# Compute local energy (squared response)
700
energy = sitk.Square(filtered)
701
702
# Average energy in local neighborhood
703
energy_avg = sitk.Mean(energy, radius=[7, 7])
704
705
texture_measures[name] = energy_avg
706
707
return texture_measures
708
709
def fractal_dimension_estimation(binary_image, box_sizes=None):
710
"""Estimate fractal dimension using box-counting method"""
711
712
if box_sizes is None:
713
max_size = min(binary_image.GetSize()) // 4
714
box_sizes = [2**i for i in range(1, int(np.log2(max_size)) + 1)]
715
716
counts = []
717
718
for box_size in box_sizes:
719
# Downsample image to box_size grid
720
shrink_factor = [box_size] * binary_image.GetDimension()
721
downsampled = sitk.Shrink(binary_image, shrink_factor)
722
723
# Count non-zero boxes
724
array = sitk.GetArrayFromImage(downsampled)
725
non_zero_count = np.count_nonzero(array)
726
counts.append(non_zero_count)
727
728
# Linear regression on log-log plot
729
log_box_sizes = np.log(box_sizes)
730
log_counts = np.log(counts)
731
732
# Fit line: log(count) = -D * log(box_size) + constant
733
coeffs = np.polyfit(log_box_sizes, log_counts, 1)
734
fractal_dimension = -coeffs[0] # Negative slope is fractal dimension
735
736
return fractal_dimension, box_sizes, counts
737
738
# Apply texture analysis
739
texture_image = sitk.ReadImage('texture_sample.tif', sitk.sitkFloat32)
740
741
# Comprehensive texture features
742
texture_features = comprehensive_texture_analysis(texture_image)
743
744
print("GLCM Features:")
745
for feature_name in ['energy', 'entropy', 'correlation', 'inertia']:
746
if feature_name in texture_features:
747
print(f" {feature_name}: {texture_features[feature_name]:.4f}")
748
749
# Laws texture measures
750
laws_features = laws_texture_measures(texture_image)
751
752
# Fractal analysis (if binary image available)
753
binary_texture = sitk.BinaryThreshold(texture_image, 128, 255, 1, 0)
754
fractal_dim, box_sizes, counts = fractal_dimension_estimation(binary_texture)
755
print(f"Fractal Dimension: {fractal_dim:.3f}")
756
757
# Display selected features
758
sitk.Show(texture_image, "Original Texture")
759
if 'local_variance' in texture_features:
760
sitk.Show(texture_features['local_variance'], "Local Variance")
761
762
# Display Laws features
763
if 'LL' in laws_features:
764
sitk.Show(laws_features['LL'], "Laws LL (Level-Level)")
765
if 'EE' in laws_features:
766
sitk.Show(laws_features['EE'], "Laws EE (Edge-Edge)")
767
```
768
769
### Ridge and Valley Detection
770
771
Specialized detection of ridge-like and valley-like structures in medical images.
772
773
```python { .api }
774
def MultiScaleHessianBasedRidgeFilter(image: Image, sigmaMin: float = 1.0,
775
sigmaMax: float = 8.0, numberOfSigmaSteps: int = 10,
776
alpha: float = 0.5, beta: float = 0.5,
777
gamma: float = 15.0) -> Image:
778
"""
779
Multi-scale Hessian-based ridge detection
780
781
Args:
782
image: Input grayscale image
783
sigmaMin: Minimum scale parameter
784
sigmaMax: Maximum scale parameter
785
numberOfSigmaSteps: Number of scales to test
786
alpha: Plate-like vs line-like discrimination
787
beta: Line-like vs blob-like discrimination
788
gamma: Background suppression parameter
789
790
Returns:
791
Ridge strength image
792
"""
793
794
def SymmetricEigenAnalysis(image: Image, dimension: int = 3,
795
orderEigenValuesBy: str = 'OrderByValue') -> tuple:
796
"""
797
Compute eigenvalues and eigenvectors of Hessian matrix
798
799
Args:
800
image: Input Hessian matrix image (symmetric)
801
dimension: Spatial dimension (2D or 3D)
802
orderEigenValuesBy: Ordering method for eigenvalues
803
804
Returns:
805
Tuple of (eigenvalue images, eigenvector images)
806
"""
807
808
def StructureTensorAnalysis(image: Image, innerScale: float = 1.0,
809
outerScale: float = 3.0) -> dict:
810
"""
811
Structure tensor analysis for ridge/valley detection
812
813
Args:
814
image: Input grayscale image
815
innerScale: Scale for gradient computation
816
outerScale: Scale for tensor averaging
817
818
Returns:
819
Dictionary with eigenvalues and orientation information
820
"""
821
```
822
823
**Usage Examples:**
824
825
```python
826
import SimpleITK as sitk
827
import numpy as np
828
829
def ridge_valley_analysis(image):
830
"""Comprehensive ridge and valley analysis"""
831
832
results = {}
833
834
# Multi-scale Hessian analysis
835
scales = np.logspace(0, 1.5, 8) # Scales from 1 to ~32
836
837
ridge_responses = []
838
valley_responses = []
839
840
for sigma in scales:
841
# Compute Hessian matrix
842
hxx = sitk.RecursiveGaussianImageFilter(image, sigma=sigma,
843
direction=0, order=2)
844
hyy = sitk.RecursiveGaussianImageFilter(image, sigma=sigma,
845
direction=1, order=2)
846
hxy = sitk.RecursiveGaussianImageFilter(image, sigma=sigma,
847
direction=0, order=1)
848
hxy = sitk.RecursiveGaussianImageFilter(hxy, sigma=sigma,
849
direction=1, order=1)
850
851
# For 2D images, eigenvalues of Hessian
852
# λ1, λ2 where λ1 ≤ λ2
853
trace = sitk.Add(hxx, hyy)
854
det = sitk.Subtract(sitk.Multiply(hxx, hyy), sitk.Square(hxy))
855
856
# Eigenvalues: λ = (trace ± sqrt(trace² - 4*det)) / 2
857
discriminant = sitk.Subtract(sitk.Square(trace), sitk.Multiply(4.0, det))
858
sqrt_disc = sitk.Sqrt(sitk.Maximum(discriminant, 0.0))
859
860
lambda1 = sitk.Multiply(0.5, sitk.Subtract(trace, sqrt_disc)) # Smaller eigenvalue
861
lambda2 = sitk.Multiply(0.5, sitk.Add(trace, sqrt_disc)) # Larger eigenvalue
862
863
# Ridge detection: λ1 < 0, λ2 ≈ 0
864
# Valley detection: λ1 ≈ 0, λ2 > 0
865
866
# Ridge strength (bright ridges on dark background)
867
ridge_condition = sitk.And(sitk.Less(lambda1, -0.001),
868
sitk.Greater(sitk.Abs(lambda2), sitk.Abs(lambda1)))
869
ridge_strength = sitk.Multiply(sitk.Cast(ridge_condition, sitk.sitkFloat32),
870
sitk.Abs(lambda1))
871
ridge_strength = sitk.Multiply(ridge_strength, sigma * sigma) # Scale normalization
872
873
# Valley strength (dark valleys on bright background)
874
valley_condition = sitk.And(sitk.Greater(lambda2, 0.001),
875
sitk.Greater(sitk.Abs(lambda2), sitk.Abs(lambda1)))
876
valley_strength = sitk.Multiply(sitk.Cast(valley_condition, sitk.sitkFloat32),
877
lambda2)
878
valley_strength = sitk.Multiply(valley_strength, sigma * sigma)
879
880
ridge_responses.append(ridge_strength)
881
valley_responses.append(valley_strength)
882
883
# Maximum response across scales
884
max_ridge = ridge_responses[0]
885
max_valley = valley_responses[0]
886
887
for ridge_resp, valley_resp in zip(ridge_responses[1:], valley_responses[1:]):
888
max_ridge = sitk.Maximum(max_ridge, ridge_resp)
889
max_valley = sitk.Maximum(max_valley, valley_resp)
890
891
results['ridge_strength'] = max_ridge
892
results['valley_strength'] = max_valley
893
results['ridge_responses'] = ridge_responses
894
results['valley_responses'] = valley_responses
895
896
return results
897
898
def vesselness_measure(image, scales=None, alpha=0.5, beta=0.5, gamma=15.0):
899
"""Frangi vesselness measure for tube-like structure detection"""
900
901
if scales is None:
902
scales = [1.0, 2.0, 4.0, 8.0]
903
904
vesselness_responses = []
905
906
for sigma in scales:
907
# Compute Hessian eigenvalues
908
hxx = sitk.RecursiveGaussianImageFilter(image, sigma=sigma, direction=0, order=2)
909
hyy = sitk.RecursiveGaussianImageFilter(image, sigma=sigma, direction=1, order=2)
910
hxy = sitk.RecursiveGaussianImageFilter(image, sigma=sigma, direction=0, order=1)
911
hxy = sitk.RecursiveGaussianImageFilter(hxy, sigma=sigma, direction=1, order=1)
912
913
# Eigenvalues
914
trace = sitk.Add(hxx, hyy)
915
det = sitk.Subtract(sitk.Multiply(hxx, hyy), sitk.Square(hxy))
916
discriminant = sitk.Subtract(sitk.Square(trace), sitk.Multiply(4.0, det))
917
sqrt_disc = sitk.Sqrt(sitk.Maximum(discriminant, 0.0))
918
919
lambda1 = sitk.Multiply(0.5, sitk.Subtract(trace, sqrt_disc))
920
lambda2 = sitk.Multiply(0.5, sitk.Add(trace, sqrt_disc))
921
922
# Frangi vesselness measures
923
# Rb = |λ1| / |λ2| (ratio measure)
924
# S = sqrt(λ1² + λ2²) (structureness measure)
925
926
abs_lambda1 = sitk.Abs(lambda1)
927
abs_lambda2 = sitk.Abs(lambda2)
928
929
# Avoid division by zero
930
lambda2_safe = sitk.Maximum(abs_lambda2, 1e-10)
931
Rb = sitk.Divide(abs_lambda1, lambda2_safe)
932
933
S = sitk.Sqrt(sitk.Add(sitk.Square(lambda1), sitk.Square(lambda2)))
934
935
# Vesselness response
936
# V = exp(-Rb²/(2α²)) * (1 - exp(-S²/(2γ²))) * condition
937
term1 = sitk.Exp(sitk.Multiply(-0.5 / (alpha * alpha), sitk.Square(Rb)))
938
term2 = sitk.Subtract(1.0, sitk.Exp(sitk.Multiply(-0.5 / (gamma * gamma),
939
sitk.Square(S))))
940
941
# Only consider pixels where λ2 < 0 (bright structures on dark background)
942
condition = sitk.Cast(sitk.Less(lambda2, 0.0), sitk.sitkFloat32)
943
944
vesselness = sitk.Multiply(sitk.Multiply(term1, term2), condition)
945
vesselness = sitk.Multiply(vesselness, sigma * sigma) # Scale normalization
946
947
vesselness_responses.append(vesselness)
948
949
# Maximum response across scales
950
max_vesselness = vesselness_responses[0]
951
for resp in vesselness_responses[1:]:
952
max_vesselness = sitk.Maximum(max_vesselness, resp)
953
954
return max_vesselness, vesselness_responses
955
956
# Apply ridge/valley detection
957
image = sitk.ReadImage('retinal_vessels.png', sitk.sitkFloat32)
958
959
# Ridge and valley analysis
960
ridge_valley_results = ridge_valley_analysis(image)
961
962
# Vesselness analysis
963
vessel_response, vessel_scales = vesselness_measure(image,
964
scales=[0.5, 1.0, 2.0, 4.0],
965
alpha=0.5, beta=0.5, gamma=15.0)
966
967
# Threshold responses to extract features
968
ridge_threshold = 0.1 * sitk.GetArrayFromImage(ridge_valley_results['ridge_strength']).max()
969
valley_threshold = 0.1 * sitk.GetArrayFromImage(ridge_valley_results['valley_strength']).max()
970
vessel_threshold = 0.1 * sitk.GetArrayFromImage(vessel_response).max()
971
972
ridges = sitk.BinaryThreshold(ridge_valley_results['ridge_strength'],
973
ridge_threshold, 1e6, 1, 0)
974
valleys = sitk.BinaryThreshold(ridge_valley_results['valley_strength'],
975
valley_threshold, 1e6, 1, 0)
976
vessels = sitk.BinaryThreshold(vessel_response, vessel_threshold, 1e6, 1, 0)
977
978
# Display results
979
sitk.Show(image, "Original")
980
sitk.Show(ridge_valley_results['ridge_strength'], "Ridge Strength")
981
sitk.Show(ridge_valley_results['valley_strength'], "Valley Strength")
982
sitk.Show(vessel_response, "Vesselness")
983
sitk.Show(ridges, "Detected Ridges")
984
sitk.Show(vessels, "Detected Vessels")
985
```
986
987
## Type Definitions
988
989
```python { .api }
990
# Core Types
991
Image = sitk.Image
992
"""SimpleITK Image object"""
993
994
FeatureImage = Image
995
"""Image containing feature detection results"""
996
997
ResponseImage = Image
998
"""Feature response/strength image"""
999
1000
# Detection Parameter Types
1001
ThresholdValue = float
1002
"""Threshold for feature detection"""
1003
1004
ScaleValue = float
1005
"""Scale parameter (typically sigma for Gaussian operations)"""
1006
1007
ScaleList = list[float]
1008
"""List of scales for multi-scale analysis"""
1009
1010
# Corner Detection Types
1011
CornerResponse = Image
1012
"""Corner strength/response image"""
1013
1014
CornerPoints = list[tuple[float, ...]]
1015
"""List of detected corner coordinates"""
1016
1017
HarrisAlpha = float
1018
"""Harris corner detector free parameter (typically 0.04-0.06)"""
1019
1020
# Edge Detection Types
1021
EdgeMap = Image
1022
"""Binary edge detection result"""
1023
1024
GradientMagnitude = Image
1025
"""Magnitude of image gradient"""
1026
1027
GradientDirection = Image
1028
"""Direction/orientation of image gradient"""
1029
1030
CannyThresholds = tuple[float, float]
1031
"""Lower and upper thresholds for Canny edge detection"""
1032
1033
# Blob Detection Types
1034
BlobResponse = Image
1035
"""Blob detection response image"""
1036
1037
BlobScale = Image
1038
"""Scale at which each blob was detected"""
1039
1040
BlobMask = Image
1041
"""Binary mask of detected blobs"""
1042
1043
LaplacianResponse = Image
1044
"""Laplacian of Gaussian response"""
1045
1046
HessianResponse = Image
1047
"""Hessian-based feature response"""
1048
1049
# Texture Analysis Types
1050
TextureFeatures = dict[str, float]
1051
"""Dictionary of computed texture features"""
1052
1053
GLCMFeatures = dict[str, float]
1054
"""Gray-Level Co-occurrence Matrix features"""
1055
1056
LawsFeatures = dict[str, Image]
1057
"""Laws texture energy measures"""
1058
1059
GaborResponse = Image
1060
"""Gabor filter response"""
1061
1062
# Ridge/Valley Detection Types
1063
RidgeStrength = Image
1064
"""Ridge detection strength image"""
1065
1066
ValleyStrength = Image
1067
"""Valley detection strength image"""
1068
1069
VesselnessResponse = Image
1070
"""Vesselness measure for tube-like structures"""
1071
1072
EigenValue = Image
1073
"""Eigenvalue image from structure analysis"""
1074
1075
EigenVector = Image
1076
"""Eigenvector image from structure analysis"""
1077
1078
# Multi-scale Analysis Types
1079
ScaleSpaceResponse = list[Image]
1080
"""Responses across multiple scales"""
1081
1082
MaximumResponse = Image
1083
"""Maximum response across scale space"""
1084
1085
OptimalScale = Image
1086
"""Scale producing maximum response at each pixel"""
1087
1088
# Advanced Feature Types
1089
StructureTensor = dict[str, Image]
1090
"""Structure tensor components and derived measures"""
1091
1092
FractalDimension = float
1093
"""Estimated fractal dimension value"""
1094
1095
OrientationField = Image
1096
"""Local orientation/direction field"""
1097
1098
CoherenceMap = Image
1099
"""Local coherence/anisotropy measure"""
1100
```