1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
|
/*
//
// Copyright 2012, 2013 SRI International
//
// This file is part of the Computational Morphometry Toolkit.
//
// http://www.nitrc.org/projects/cmtk/
//
// The Computational Morphometry Toolkit is free software: you can
// redistribute it and/or modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation, either version 3 of
// the License, or (at your option) any later version.
//
// The Computational Morphometry Toolkit is distributed in the hope that it
// will be useful, but WITHOUT ANY WARRANTY; without even the implied
// warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with the Computational Morphometry Toolkit. If not, see
// <http://www.gnu.org/licenses/>.
//
// $Revision: 5436 $
//
// $LastChangedDate: 2018-12-10 19:01:20 -0800 (Mon, 10 Dec 2018) $
//
// $LastChangedBy: torstenrohlfing $
//
*/
#ifndef __cmtkSphereDetectionNormalizedBipolarMatchedFilterFFT_h_included_
#define __cmtkSphereDetectionNormalizedBipolarMatchedFilterFFT_h_included_
#include <cmtkconfig.h>
#include <Base/cmtkUniformVolume.h>
#include <System/cmtkFFTW.h>
namespace
cmtk
{
/** Detect spheres in an image using FFT-based matched bipolar filter.
* An object of this class can be used to detect spheres of different sizes with only two FFTs applied to the test image and its square.
* Each detected sphere size does require a different filter kernel and thus three repeated FFTs of the kernel, its mask, and its square.
*
* The filter kernel is bipolar, i.e., +1 inside the sphere and -1 outside the sphere, each within a user-provided margin inside and outside the
* sphere surface. This makes the filter robust to intensity differences across the images.
*
* Because the filter values are either +1 or -1 or 0, the squared filter is identical to the filter mask. This simplifies the computation and
* saves FT of the squared filter.
*
*\see D. Padfield, "Masked object registration in the Fourier domain," IEEE Transactions on Image Processing, vol. 21, no. 5, pp. 2706-2718, 2012.
* http://dx.doi.org/10.1109/TIP.2011.2181402
* http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6111478&isnumber=4358840
*
*\note This class requires CMTK to be configured with FFTW3 support ("CMTK_USE_FFTW" CMake option).
*\todo The current implementation does not take advantage of the real-valued image and filter data, which could be used to reduce the storage
* requirement of the FT data (and probably the computational cost of the transform) by almost 50%. On the other hand, capitalizing on these
* savings would either require out-of-place, rather than in-place, transforms, or substantially complicate memory layout of the input data.
*
*\todo Currently, the FFT calls within the class do not take advantage of the redundancy in the real-to-complex FT applied to the
* real-valued image and mask. Properly considering this could potentially save about 1/2 of memory and CPU time.
*/
class SphereDetectionNormalizedBipolarMatchedFilterFFT
{
public:
/// Constructor: initialize FFTW plans and compute image FT.
SphereDetectionNormalizedBipolarMatchedFilterFFT( const UniformVolume& image );
/// Destructor: destroy FFTW plans and de-allocate transformed image memories.
virtual ~SphereDetectionNormalizedBipolarMatchedFilterFFT();
/// Get image filtered with spherical matched filter kernel.
cmtk::TypedArray::SmartPtr GetFilteredImageData( const Types::Coordinate sphereRadius /*!< Radius of detected spheres in world coordinate units (e.g., mm) */,
const int marginWidth = 1 /*!< Half width of the filter margin in pixels: positive filter coefficients in a band of this width inside radius, negative coeffiecients outside radius.*/ );
private:
/// Image number of pixels.
size_t m_NumberOfPixels;
/// Image dimensions.
DataGrid::IndexType m_ImageDims;
/// Image pixel size.
UniformVolume::SpaceVectorType m_PixelSize;
/// Store previous sphere radius to avoid unnecessary recomputation.
Types::Coordinate m_SphereRadius;
/// Store previous filter margin to avoid unnecessary recomputation.
int m_MarginWidth;
/// Store computed filter response.
TypedArray::SmartPtr m_FilterResponse;
/// The Fourier-transformed image.
fftw_complex* m_ImageFT;
/// The Fourier-transformed squared image.
fftw_complex* m_ImageSquareFT;
/// The Fourier-transformed matched filter.
fftw_complex* m_FilterFT;
/// The Fourier-transformed matched filter mask.
fftw_complex* m_FilterMaskFT;
/// Copy of the Fourier-transformed matched filter mask for computing a separate product.
fftw_complex* m_FilterMaskFT2;
/// The filter FFT plan.
fftw_plan m_PlanFilter;
/// The filter mask FFT plan.
fftw_plan m_PlanFilterMask;
/// The backward (filtered data to space domain) FFT plan.
fftw_plan m_PlanBackward;
/// The backward FFT plan for the mask.
fftw_plan m_PlanBackwardMask;
/// The backward FFT plan for the copy of the multiplied mask.
fftw_plan m_PlanBackwardMask2;
/// Sum of filter elements.
Types::DataItem m_SumFilter;
/// Sum of filter mask elements (number of non-zero elements).
Types::DataItem m_SumFilterMask;
/** Make the filter kernel.
*/
void MakeFilter( const Types::Coordinate sphereRadius, const int marginWidth );
};
} // namespace cmtk
#endif // #ifndef __cmtkSphereDetectionNormalizedBipolarMatchedFilterFFT_h_included_
|
