SimpleITK  
Python/FFTConvolution.py
1#!/usr/bin/env python
2# =========================================================================
3#
4# Copyright NumFOCUS
5#
6# Licensed under the Apache License, Version 2.0 (the "License");
7# you may not use this file except in compliance with the License.
8# You may obtain a copy of the License at
9#
10# http://www.apache.org/licenses/LICENSE-2.0.txt
11#
12# Unless required by applicable law or agreed to in writing, software
13# distributed under the License is distributed on an "AS IS" BASIS,
14# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15# See the License for the specific language governing permissions and
16# limitations under the License.
17#
18# =========================================================================
19
20""" Example script showing Fast Fourier Transform convolution in SimpleITK """
21
22import sys
23import os
24import math
25import SimpleITK as sitk
26
27if len(sys.argv) < 4:
28 print("Usage: FFTConvolution <input> <kernel> <output>")
29 sys.exit(1)
30
31inputFileName = sys.argv[1]
32kernelFileName = sys.argv[2]
33outputFileName = sys.argv[3]
34
35# -- Input Image --
36# read the input image
37img = sitk.ReadImage(inputFileName)
38
39# save the type of pixel the input is, so that we can cast the result
40# back out to the same type
41pixelID = img.GetPixelID()
42
43# pad the image
44img = sitk.MirrorPad(img, [128] * 2, [128] * 2)
45size = img.GetSize()
46
47# perform the FFT
48fftimg = sitk.ForwardFFT(sitk.Cast(img, sitk.sitkFloat32))
49
50# -- Kernel Image --
51# Read the kernel image file
52kernel = sitk.ReadImage(kernelFileName)
53
54# flip kernel about all axis
55kernel = sitk.Flip(kernel, [True] * 2)
56
57# normalize the kernel to sum to ~1
59stats.Execute(kernel)
60kernel = sitk.Cast(kernel / stats.GetSum(), sitk.sitkFloat32)
61
62upadding = [0] * 2
63upadding[0] = int(math.floor((size[0] - kernel.GetSize()[0]) / 2.0))
64upadding[1] = int(math.floor((size[1] - kernel.GetSize()[1]) / 2.0))
65
66lpadding = [0] * 2
67lpadding[0] = int(math.ceil((size[0] - kernel.GetSize()[0]) / 2.0))
68lpadding[1] = int(math.ceil((size[1] - kernel.GetSize()[1]) / 2.0))
69
70# pad the kernel to prevent edge artifacts
71kernel = sitk.ConstantPad(kernel, upadding, lpadding, 0.0)
72
73# perform FFT on kernel
74fftkernel = sitk.ForwardFFT(sitk.FFTShift(kernel))
75
76# meta-data must match for multiplication
77fftkernel.SetSpacing(fftimg.GetSpacing())
78fftkernel.SetOrigin(fftimg.GetOrigin())
79fftkernel.SetDirection(fftimg.GetDirection())
80
81# -- Convolution --
82# Finally perform the convolution in Fourier space by multiplication
83img = sitk.InverseFFT(fftimg * fftkernel)
84
85# remove the padding
86img = sitk.Crop(img, [128] * 2, [128] * 2)
87
88# -- Writing --
89# write the output image the same type as the input
90sitk.WriteImage(sitk.Cast(img, pixelID), outputFileName)
91
92if "SITK_NOSHOW" not in os.environ:
93 sitk.Show(sitk.ReadImage(inputFileName), "original")
94 sitk.Show(sitk.ReadImage(kernelFileName), "kernel")
95 sitk.Show(sitk.Cast(img, pixelID), "FFT_Convolution")
Compute min, max, variance and mean of an Image .
Image Crop(Image &&image1, std::vector< unsigned int > lowerBoundaryCropSize=std::vector< unsigned int >(3, 0), std::vector< unsigned int > upperBoundaryCropSize=std::vector< unsigned int >(3, 0))
Decrease the image size by cropping the image by an itk::Size at both the upper and lower bounds of t...
Image InverseFFT(const Image &image1)
Base class for inverse Fast Fourier Transform .
Image FFTShift(const Image &image1, bool inverse=false)
Shift the zero-frequency components of a Fourier transform to the center of the image.
Image ForwardFFT(const Image &image1)
Base class for forward Fast Fourier Transform .
SITKIO_EXPORT Image ReadImage(const PathType &filename, PixelIDValueEnum outputPixelType=sitkUnknown, const std::string &imageIO="")
ReadImage is a procedural interface to the ImageFileReader class which is convenient for most image r...
Image Flip(const Image &image1, std::vector< bool > flipAxes=std::vector< bool >(3, false), bool flipAboutOrigin=false)
Flips an image across user specified axes.
void SITKIO_EXPORT Show(const Image &image, const std::string &title="", const bool debugOn=ProcessObject::GetGlobalDefaultDebug())
Image MirrorPad(const Image &image1, std::vector< unsigned int > padLowerBound=std::vector< unsigned int >(3, 0), std::vector< unsigned int > padUpperBound=std::vector< unsigned int >(3, 0), double decayBase=1.0)
Increase the image size by padding with replicants of the input image value.
SITKIO_EXPORT void WriteImage(const Image &image, const PathType &fileName, bool useCompression=false, int compressionLevel=-1)
WriteImage is a procedural interface to the ImageFileWriter. class which is convenient for many image...
Image Cast(const Image &image, PixelIDValueEnum pixelID)
Image ConstantPad(const Image &image1, std::vector< unsigned int > padLowerBound=std::vector< unsigned int >(3, 0), std::vector< unsigned int > padUpperBound=std::vector< unsigned int >(3, 0), double constant=0.0)
Increase the image size by padding with a constant value.