0
# Image Functions
1

2
Interpolation, analysis, and processing functions that operate on image data including various interpolation schemes, neighborhood operations, and image analysis functions. These functions provide the fundamental building blocks for image processing pipelines and spatial analysis operations.
3

4
## Capabilities
5

6
### Image Interpolation
7

8
Functions for interpolating pixel values at arbitrary spatial locations within images.
9

10
```python { .api }
11
class InterpolateImageFunction[TInputImage, TCoordRep]:
12
    """
13
    Base class for image interpolation functions.
14
    
15
    Template Parameters:
16
    - TInputImage: Input image type
17
    - TCoordRep: Coordinate representation type (typically float or double)
18
    """
19
    def SetInputImage(self, image: TInputImage) -> None: ...
20
    def GetInputImage(self) -> TInputImage: ...
21
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
22
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
23
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> OutputType: ...
24
    def IsInsideBuffer(self, index: Index[ImageDimension]) -> bool: ...
25
    def IsInsideBuffer(self, point: Point[TCoordRep, ImageDimension]) -> bool: ...
26
    def IsInsideBuffer(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> bool: ...
27

28
class NearestNeighborInterpolateImageFunction[TInputImage, TCoordRep]:
29
    """Nearest neighbor interpolation (no smoothing)."""
30
    def New() -> NearestNeighborInterpolateImageFunction[TInputImage, TCoordRep]: ...
31
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
32
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
33
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> OutputType: ...
34

35
class LinearInterpolateImageFunction[TInputImage, TCoordRep]:
36
    """Linear (bilinear/trilinear) interpolation."""
37
    def New() -> LinearInterpolateImageFunction[TInputImage, TCoordRep]: ...
38
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
39
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
40
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> OutputType: ...
41

42
class BSplineInterpolateImageFunction[TInputImage, TCoordRep, TCoefficientType]:
43
    """B-spline interpolation with configurable order."""
44
    def New() -> BSplineInterpolateImageFunction[TInputImage, TCoordRep, TCoefficientType]: ...
45
    def SetSplineOrder(self, order: int) -> None: ...
46
    def GetSplineOrder(self) -> int: ...
47
    def SetInputImage(self, image: TInputImage) -> None: ...
48
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
49
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
50
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> OutputType: ...
51
    def EvaluateDerivative(self, point: Point[TCoordRep, ImageDimension]) -> CovariantVector[OutputType, ImageDimension]: ...
52
    def EvaluateDerivativeAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> CovariantVector[OutputType, ImageDimension]: ...
53

54
class WindowedSincInterpolateImageFunction[TInputImage, VRadius, TWindowFunction, TBoundaryCondition, TCoordRep]:
55
    """Windowed sinc interpolation for high-quality resampling."""
56
    def New() -> WindowedSincInterpolateImageFunction[TInputImage, VRadius, TWindowFunction, TBoundaryCondition, TCoordRep]: ...
57
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
58
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
59
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> OutputType: ...
60
    def GetRadius(self) -> int: ...
61
```
62

63
### Vector Field Interpolation
64

65
Specialized interpolation functions for vector and tensor fields.
66

67
```python { .api }
68
class VectorInterpolateImageFunction[TInputImage, TCoordRep]:
69
    """Base class for vector field interpolation."""
70
    def SetInputImage(self, image: TInputImage) -> None: ...
71
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
72
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
73
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> OutputType: ...
74

75
class VectorLinearInterpolateImageFunction[TInputImage, TCoordRep]:
76
    """Linear interpolation for vector fields."""
77
    def New() -> VectorLinearInterpolateImageFunction[TInputImage, TCoordRep]: ...
78
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
79
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
80
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> OutputType: ...
81

82
class VectorNearestNeighborInterpolateImageFunction[TInputImage, TCoordRep]:
83
    """Nearest neighbor interpolation for vector fields."""
84
    def New() -> VectorNearestNeighborInterpolateImageFunction[TInputImage, TCoordRep]: ...
85
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
86
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
87
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> OutputType: ...
88
```
89

90
### Image Analysis Functions
91

92
Functions for computing statistical and geometric properties of image regions.
93

94
```python { .api }
95
class ImageFunction[TInputImage, TOutput, TCoordRep]:
96
    """Base class for image analysis functions."""
97
    def SetInputImage(self, image: TInputImage) -> None: ...
98
    def GetInputImage(self) -> TInputImage: ...
99
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> TOutput: ...
100
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> TOutput: ...
101
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> TOutput: ...
102
    def IsInsideBuffer(self, index: Index[ImageDimension]) -> bool: ...
103
    def IsInsideBuffer(self, point: Point[TCoordRep, ImageDimension]) -> bool: ...
104

105
class MeanImageFunction[TInputImage, TCoordRep]:
106
    """Compute local mean values within a neighborhood."""
107
    def New() -> MeanImageFunction[TInputImage, TCoordRep]: ...
108
    def SetNeighborhoodRadius(self, radius: int) -> None: ...
109
    def GetNeighborhoodRadius(self) -> int: ...
110
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> RealType: ...
111
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> RealType: ...
112
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> RealType: ...
113

114
class MedianImageFunction[TInputImage, TCoordRep]:
115
    """Compute local median values within a neighborhood."""
116
    def New() -> MedianImageFunction[TInputImage, TCoordRep]: ...
117
    def SetNeighborhoodRadius(self, radius: int) -> None: ...
118
    def GetNeighborhoodRadius(self) -> int: ...
119
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
120
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
121

122
class VarianceImageFunction[TInputImage, TCoordRep]:
123
    """Compute local variance within a neighborhood."""
124
    def New() -> VarianceImageFunction[TInputImage, TCoordRep]: ...
125
    def SetNeighborhoodRadius(self, radius: int) -> None: ...
126
    def GetNeighborhoodRadius(self) -> int: ...
127
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> RealType: ...
128
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> RealType: ...
129

130
class BinaryThresholdImageFunction[TInputImage, TCoordRep]:
131
    """Test if pixel values fall within a threshold range."""
132
    def New() -> BinaryThresholdImageFunction[TInputImage, TCoordRep]: ...
133
    def SetLower(self, threshold: PixelType) -> None: ...
134
    def GetLower(self) -> PixelType: ...
135
    def SetUpper(self, threshold: PixelType) -> None: ...
136
    def GetUpper(self) -> PixelType: ...
137
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> bool: ...
138
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> bool: ...
139
    def ThresholdAbove(self, threshold: PixelType) -> None: ...
140
    def ThresholdBelow(self, threshold: PixelType) -> None: ...
141
    def ThresholdBetween(self, lower: PixelType, upper: PixelType) -> None: ...
142
```
143

144
### Derivative and Gradient Functions
145

146
Functions for computing spatial derivatives and gradients.
147

148
```python { .api }
149
class CentralDifferenceImageFunction[TInputImage, TCoordRep]:
150
    """Compute image derivatives using central differences."""
151
    def New() -> CentralDifferenceImageFunction[TInputImage, TCoordRep]: ...
152
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
153
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
154
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> OutputType: ...
155
    def SetUseImageSpacing(self, flag: bool) -> None: ...
156
    def GetUseImageSpacing(self) -> bool: ...
157

158
class GradientImageFunction[TInputImage, TCoordRep]:
159
    """Compute image gradients."""
160
    def New() -> GradientImageFunction[TInputImage, TCoordRep]: ...
161
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
162
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
163
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[TCoordRep, ImageDimension]) -> OutputType: ...
164
    def SetUseImageSpacing(self, flag: bool) -> None: ...
165

166
class HessianRecursiveGaussianImageFunction[TInputImage, TCoordRep]:
167
    """Compute Hessian matrix using recursive Gaussian filters."""
168
    def New() -> HessianRecursiveGaussianImageFunction[TInputImage, TCoordRep]: ...
169
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
170
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
171
    def SetSigma(self, sigma: GaussianDerivativeOperator.ScalarRealType) -> None: ...
172
    def GetSigma(self) -> GaussianDerivativeOperator.ScalarRealType: ...
173
    def SetNormalizeAcrossScale(self, flag: bool) -> None: ...
174
```
175

176
### Neighborhood Operations
177

178
Functions that apply operations over image neighborhoods.
179

180
```python { .api }
181
class NeighborhoodOperatorImageFunction[TInputImage, TOutput]:
182
    """Apply neighborhood operators to images."""
183
    def New() -> NeighborhoodOperatorImageFunction[TInputImage, TOutput]: ...
184
    def SetOperator(self, op: NeighborhoodType) -> None: ...
185
    def GetOperator(self) -> NeighborhoodType: ...
186
    def Evaluate(self, point: Point[CoordRepType, ImageDimension]) -> TOutput: ...
187
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> TOutput: ...
188
    def EvaluateAtContinuousIndex(self, index: ContinuousIndex[CoordRepType, ImageDimension]) -> TOutput: ...
189

190
class GaussianBlurImageFunction[TInputImage, TOutput]:
191
    """Apply Gaussian blur at specified locations."""
192
    def New() -> GaussianBlurImageFunction[TInputImage, TOutput]: ...
193
    def SetSigma(self, sigma: ArrayType) -> None: ...
194
    def SetSigma(self, sigma: double) -> None: ...
195
    def GetSigma(self) -> ArrayType: ...
196
    def SetExtent(self, extent: ArrayType) -> None: ...
197
    def GetExtent(self) -> ArrayType: ...
198
    def SetFilterDimensionality(self, dims: int) -> None: ...
199
    def GetFilterDimensionality(self) -> int: ...
200
    def Evaluate(self, point: Point[CoordRepType, ImageDimension]) -> TOutput: ...
201
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> TOutput: ...
202

203
class DiscreteGaussianDerivativeImageFunction[TInputImage, TOutput]:
204
    """Compute Gaussian derivatives at specified locations."""
205
    def New() -> DiscreteGaussianDerivativeImageFunction[TInputImage, TOutput]: ...
206
    def SetVariance(self, variance: VarianceArrayType) -> None: ...
207
    def SetVariance(self, variance: double) -> None: ...
208
    def GetVariance(self) -> VarianceArrayType: ...
209
    def SetOrder(self, orders: OrderArrayType) -> None: ...
210
    def GetOrder(self) -> OrderArrayType: ...
211
    def SetNormalizeAcrossScale(self, flag: bool) -> None: ...
212
    def GetNormalizeAcrossScale(self) -> bool: ...
213
    def SetUseImageSpacing(self, flag: bool) -> None: ...
214
    def GetUseImageSpacing(self) -> bool: ...
215
    def Evaluate(self, point: Point[CoordRepType, ImageDimension]) -> TOutput: ...
216
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> TOutput: ...
217
```
218

219
### Radon Transform Functions
220

221
Functions for computing Radon transforms and related operations.
222

223
```python { .api }
224
class RadonTransformImageFunction[TInputImage, TCoordRep]:
225
    """Compute Radon transform along specified directions."""
226
    def New() -> RadonTransformImageFunction[TInputImage, TCoordRep]: ...
227
    def Evaluate(self, point: Point[TCoordRep, ImageDimension]) -> OutputType: ...
228
    def EvaluateAtIndex(self, index: Index[ImageDimension]) -> OutputType: ...
229
    def SetTransform(self, transform: TransformType) -> None: ...
230
    def GetTransform(self) -> TransformType: ...
231
    def SetInterpolator(self, interpolator: InterpolatorType) -> None: ...
232
    def GetInterpolator(self) -> InterpolatorType: ...
233
```
234

235
## Usage Examples
236

237
### Basic Image Interpolation
238

239
```python
240
import itk
241
import numpy as np
242

243
# Create a test image
244
ImageType = itk.Image[itk.F, 2]
245
image = ImageType.New()
246

247
# Set up image region
248
size = itk.Size[2]([100, 100])
249
region = itk.ImageRegion[2]()
250
region.SetSize(size)
251
image.SetRegions(region)
252
image.SetSpacing([1.0, 1.0])
253
image.SetOrigin([0.0, 0.0])
254
image.Allocate()
255

256
# Fill with a simple pattern
257
for x in range(100):
258
    for y in range(100):
259
        idx = itk.Index[2]([x, y])
260
        value = float(x + y) / 2.0
261
        image.SetPixel(idx, value)
262

263
# Create different interpolators
264
linear_interp = itk.LinearInterpolateImageFunction[ImageType, itk.D].New()
265
linear_interp.SetInputImage(image)
266

267
nn_interp = itk.NearestNeighborInterpolateImageFunction[ImageType, itk.D].New()
268
nn_interp.SetInputImage(image)
269

270
bspline_interp = itk.BSplineInterpolateImageFunction[ImageType, itk.D, itk.D].New()
271
bspline_interp.SetInputImage(image)
272
bspline_interp.SetSplineOrder(3)  # Cubic B-spline
273

274
# Test interpolation at non-integer coordinates
275
test_point = itk.Point[itk.D, 2]([25.5, 30.7])
276

277
linear_value = linear_interp.Evaluate(test_point)
278
nn_value = nn_interp.Evaluate(test_point)
279
bspline_value = bspline_interp.Evaluate(test_point)
280

281
print(f"Point {test_point}:")
282
print(f"Linear interpolation: {linear_value}")
283
print(f"Nearest neighbor: {nn_value}")
284
print(f"B-spline interpolation: {bspline_value}")
285

286
# Test at continuous index
287
cont_index = itk.ContinuousIndex[itk.D, 2]([25.5, 30.7])
288
linear_value_ci = linear_interp.EvaluateAtContinuousIndex(cont_index)
289
print(f"Linear at continuous index: {linear_value_ci}")
290
```
291

292
### Vector Field Interpolation
293

294
```python
295
import itk
296
import numpy as np
297

298
# Create a 2D vector field
299
VectorType = itk.Vector[itk.F, 2]
300
VectorImageType = itk.Image[VectorType, 2]
301
vector_image = VectorImageType.New()
302

303
# Set up region
304
size = itk.Size[2]([50, 50])
305
region = itk.ImageRegion[2]()
306
region.SetSize(size)
307
vector_image.SetRegions(region)
308
vector_image.SetSpacing([1.0, 1.0])
309
vector_image.Allocate()
310

311
# Fill with a rotation field
312
center_x, center_y = 25.0, 25.0
313
for x in range(50):
314
    for y in range(50):
315
        idx = itk.Index[2]([x, y])
316
        
317
        # Create rotation vector field
318
        dx = float(x) - center_x
319
        dy = float(y) - center_y
320
        
321
        vector = VectorType()
322
        vector.SetElement(0, -dy * 0.1)  # x component
323
        vector.SetElement(1, dx * 0.1)   # y component
324
        
325
        vector_image.SetPixel(idx, vector)
326

327
# Create vector interpolator
328
vector_interp = itk.VectorLinearInterpolateImageFunction[VectorImageType, itk.D].New()
329
vector_interp.SetInputImage(vector_image)
330

331
# Test interpolation
332
test_point = itk.Point[itk.D, 2]([25.5, 30.7])
333
interpolated_vector = vector_interp.Evaluate(test_point)
334

335
print(f"Interpolated vector at {test_point}:")
336
print(f"X component: {interpolated_vector.GetElement(0)}")
337
print(f"Y component: {interpolated_vector.GetElement(1)}")
338
```
339

340
### Image Analysis Functions
341

342
```python
343
import itk
344
import numpy as np
345

346
# Create test image with some structure
347
ImageType = itk.Image[itk.F, 2]
348
image = ImageType.New()
349

350
size = itk.Size[2]([100, 100])
351
region = itk.ImageRegion[2]()
352
region.SetSize(size)
353
image.SetRegions(region)
354
image.Allocate()
355

356
# Create a pattern with noise
357
np.random.seed(42)
358
for x in range(100):
359
    for y in range(100):
360
        idx = itk.Index[2]([x, y])
361
        # Circular pattern with noise
362
        dist = np.sqrt((x-50)**2 + (y-50)**2)
363
        base_value = 100.0 * np.exp(-dist/20.0)
364
        noise = np.random.normal(0, 5)
365
        image.SetPixel(idx, base_value + noise)
366

367
# Create analysis functions
368
mean_func = itk.MeanImageFunction[ImageType, itk.D].New()
369
mean_func.SetInputImage(image)
370
mean_func.SetNeighborhoodRadius(3)
371

372
variance_func = itk.VarianceImageFunction[ImageType, itk.D].New()
373
variance_func.SetInputImage(image)
374
variance_func.SetNeighborhoodRadius(3)
375

376
median_func = itk.MedianImageFunction[ImageType, itk.D].New()
377
median_func.SetInputImage(image)
378
median_func.SetNeighborhoodRadius(2)
379

380
# Analyze different regions
381
test_points = [
382
    itk.Point[itk.D, 2]([50.0, 50.0]),  # Center
383
    itk.Point[itk.D, 2]([25.0, 25.0]),  # Off-center
384
    itk.Point[itk.D, 2]([10.0, 10.0])   # Edge
385
]
386

387
for point in test_points:
388
    mean_val = mean_func.Evaluate(point)
389
    var_val = variance_func.Evaluate(point)
390
    median_val = median_func.Evaluate(point)
391
    
392
    print(f"Point {point}:")
393
    print(f"  Mean: {mean_val:.2f}")
394
    print(f"  Variance: {var_val:.2f}")
395
    print(f"  Median: {median_val:.2f}")
396
    print()
397
```
398

399
### Binary Threshold Function
400

401
```python
402
import itk
403

404
# Create test image
405
ImageType = itk.Image[itk.UC, 2]  # Unsigned char image
406
image = ImageType.New()
407

408
size = itk.Size[2]([50, 50])
409
region = itk.ImageRegion[2]()
410
region.SetSize(size)
411
image.SetRegions(region)
412
image.Allocate()
413

414
# Fill with gradient pattern
415
for x in range(50):
416
    for y in range(50):
417
        idx = itk.Index[2]([x, y])
418
        value = int((x + y) * 255 / 98)  # Scale to 0-255
419
        image.SetPixel(idx, value)
420

421
# Create threshold function
422
threshold_func = itk.BinaryThresholdImageFunction[ImageType, itk.D].New()
423
threshold_func.SetInputImage(image)
424
threshold_func.ThresholdBetween(100, 200)  # Threshold between 100 and 200
425

426
# Test threshold at various points
427
test_points = [
428
    itk.Point[itk.D, 2]([10.0, 10.0]),  # Low values
429
    itk.Point[itk.D, 2]([25.0, 25.0]),  # Middle values
430
    itk.Point[itk.D, 2]([40.0, 40.0])   # High values
431
]
432

433
for point in test_points:
434
    # Get actual pixel value
435
    idx = itk.Index[2]([int(point[0]), int(point[1])])
436
    actual_value = image.GetPixel(idx)
437
    
438
    # Test threshold
439
    is_in_range = threshold_func.Evaluate(point)
440
    
441
    print(f"Point {point}: value={actual_value}, in_threshold={is_in_range}")
442
```
443

444
### Gradient and Derivative Functions
445

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

450
# Create test image with interesting structure
451
ImageType = itk.Image[itk.F, 2]
452
image = ImageType.New()
453

454
size = itk.Size[2]([100, 100])
455
region = itk.ImageRegion[2]()
456
region.SetSize(size)
457
image.SetRegions(region)
458
image.SetSpacing([1.0, 1.0])
459
image.Allocate()
460

461
# Create a 2D Gaussian peak
462
center_x, center_y = 50.0, 50.0
463
sigma = 15.0
464
for x in range(100):
465
    for y in range(100):
466
        idx = itk.Index[2]([x, y])
467
        dx = float(x) - center_x
468
        dy = float(y) - center_y
469
        value = 100.0 * np.exp(-(dx*dx + dy*dy) / (2 * sigma*sigma))
470
        image.SetPixel(idx, value)
471

472
# Create gradient function
473
grad_func = itk.CentralDifferenceImageFunction[ImageType, itk.D].New()
474
grad_func.SetInputImage(image)
475
grad_func.SetUseImageSpacing(True)
476

477
# Create Hessian function for second derivatives
478
hessian_func = itk.HessianRecursiveGaussianImageFunction[ImageType, itk.D].New()
479
hessian_func.SetInputImage(image)
480
hessian_func.SetSigma(2.0)
481

482
# Analyze gradients at different points
483
test_points = [
484
    itk.Point[itk.D, 2]([50.0, 50.0]),  # Center (should be near zero gradient)
485
    itk.Point[itk.D, 2]([60.0, 50.0]),  # Right of center
486
    itk.Point[itk.D, 2]([50.0, 60.0])   # Above center
487
]
488

489
for point in test_points:
490
    gradient = grad_func.Evaluate(point)
491
    hessian = hessian_func.Evaluate(point)
492
    
493
    print(f"Point {point}:")
494
    print(f"  Gradient: [{gradient[0]:.3f}, {gradient[1]:.3f}]")
495
    print(f"  Gradient magnitude: {np.sqrt(gradient[0]**2 + gradient[1]**2):.3f}")
496
    print(f"  Hessian determinant: {hessian[0]*hessian[2] - hessian[1]**2:.3f}")
497
    print()
498
```
499

500
### Neighborhood Operator Function
501

502
```python
503
import itk
504

505
# Create test image
506
ImageType = itk.Image[itk.F, 2]
507
image = ImageType.New()
508

509
size = itk.Size[2]([50, 50])
510
region = itk.ImageRegion[2]()
511
region.SetSize(size)
512
image.SetRegions(region)
513
image.Allocate()
514

515
# Fill with step pattern
516
for x in range(50):
517
    for y in range(50):
518
        idx = itk.Index[2]([x, y])
519
        value = 100.0 if x > 25 else 0.0
520
        image.SetPixel(idx, value)
521

522
# Create a simple edge detection operator (Sobel in X direction)
523
OperatorType = itk.Neighborhood[itk.F, 2]
524
operator = OperatorType()
525
operator.SetRadius(1)
526

527
# Sobel X operator: [-1, 0, 1; -2, 0, 2; -1, 0, 1]
528
operator.SetElement(0, -1.0)  # Top-left
529
operator.SetElement(1, -2.0)  # Middle-left  
530
operator.SetElement(2, -1.0)  # Bottom-left
531
operator.SetElement(3, 0.0)   # Top-center
532
operator.SetElement(4, 0.0)   # Center
533
operator.SetElement(5, 0.0)   # Bottom-center
534
operator.SetElement(6, 1.0)   # Top-right
535
operator.SetElement(7, 2.0)   # Middle-right
536
operator.SetElement(8, 1.0)   # Bottom-right
537

538
# Create neighborhood operator function
539
op_func = itk.NeighborhoodOperatorImageFunction[ImageType, itk.F].New()
540
op_func.SetInputImage(image)
541
op_func.SetOperator(operator)
542

543
# Test edge detection at different points
544
test_points = [
545
    itk.Point[itk.D, 2]([20.0, 25.0]),  # Left side of edge
546
    itk.Point[itk.D, 2]([25.0, 25.0]),  # On the edge
547
    itk.Point[itk.D, 2]([30.0, 25.0])   # Right side of edge
548
]
549

550
for point in test_points:
551
    edge_response = op_func.Evaluate(point)
552
    print(f"Point {point}: Edge response = {edge_response:.2f}")
553
```