| 1 |
DEF |
#ifndef _itkWindowedSincInterpolateImageFunction_h |
| 2 |
HRD,DEF |
#define _itkWindowedSincInterpolateImageFunction_h |
| 3 |
|
|
| 4 |
|
#include "itkConstNeighborhoodIterator.h" |
| 5 |
|
#include "itkConstantBoundaryCondition.h" |
| 6 |
|
#include "itkInterpolateImageFunction.h" |
| 7 |
|
|
| 8 |
HRD |
namespace itk |
| 9 |
|
{ |
| 10 |
|
|
| 11 |
HRD |
namespace Function { |
| 12 |
|
|
| 13 |
|
/** |
| 14 |
|
* \class CosineWindowFunction |
| 15 |
|
* \brief Window function for sinc interpolation. |
| 16 |
|
* \f[ w(x) = cos ( \frac{\pi x}{2 m} ) \f] |
| 17 |
|
* \sa WindowedSincInterpolateImageFunction |
| 18 |
|
*/ |
| 19 |
|
template< unsigned int VRadius, |
| 20 |
|
class TInput=double, class TOutput=double> |
| 21 |
|
class CosineWindowFunction |
| 22 |
|
{ |
| 23 |
|
public: |
| 24 |
|
inline TOutput operator()( const TInput & A ) const |
| 25 |
|
{ return (TOutput) cos( A * m_Factor ); } |
| 26 |
|
private: |
| 27 |
IND |
**/** Equal to \f$ \frac{\pi}{2 m} \f$ */ |
| 28 |
IND |
**static const double m_Factor; |
| 29 |
IND |
}; |
| 30 |
|
|
| 31 |
IND |
/** |
| 32 |
IND |
** \class HammingWindowFunction |
| 33 |
IND |
** \brief Window function for sinc interpolation. |
| 34 |
IND |
** \f[ w(x) = 0.54 + 0.46 cos ( \frac{\pi x}{m} ) \f] |
| 35 |
IND |
** \sa WindowedSincInterpolateImageFunction |
| 36 |
IND |
**/ |
| 37 |
IND |
template< unsigned int VRadius, |
| 38 |
|
class TInput=double, class TOutput=double> |
| 39 |
IND |
class HammingWindowFunction |
| 40 |
IND |
{ |
| 41 |
IND |
public: |
| 42 |
IND |
**inline TOutput operator()( const TInput & A ) const |
| 43 |
IND |
***{ return (TOutput) 0.54 + 0.46 * cos( A * m_Factor ); } |
| 44 |
IND |
private: |
| 45 |
IND |
**/** Equal to \f$ \frac{\pi}{m} \f$ */ |
| 46 |
IND |
**static const double m_Factor; |
| 47 |
IND |
}; |
| 48 |
|
|
| 49 |
IND |
/** |
| 50 |
IND |
** \class WelchWindowFunction |
| 51 |
IND |
** \brief Window function for sinc interpolation. |
| 52 |
IND |
** \f[ w(x) = 1 - ( \frac{x^2}{m^2} ) \f] |
| 53 |
IND |
** \sa WindowedSincInterpolateImageFunction |
| 54 |
IND |
**/ |
| 55 |
IND |
template< unsigned int VRadius, |
| 56 |
|
class TInput=double, class TOutput=double> |
| 57 |
IND |
class WelchWindowFunction |
| 58 |
IND |
{ |
| 59 |
IND |
public: |
| 60 |
IND |
**inline TOutput operator()( const TInput & A ) const |
| 61 |
IND |
***{ return (TOutput) (1.0 - A * m_Factor * A); } |
| 62 |
IND |
private: |
| 63 |
IND |
**/** Equal to \f$ \frac{1}{m^2} \f$ */ |
| 64 |
IND |
**static const double m_Factor; |
| 65 |
IND |
}; |
| 66 |
|
|
| 67 |
IND |
/** |
| 68 |
IND |
** \class LancosWindowFunction |
| 69 |
IND |
** \brief Window function for sinc interpolation. |
| 70 |
IND |
** \f[ w(x) = \textrm{sinc} ( \frac{x}{m} ) \f] |
| 71 |
IND |
** Note: Paper referenced in WindowedSincInterpolateImageFunction gives |
| 72 |
IND |
** an incorrect definition of this window function. |
| 73 |
IND |
** \sa WindowedSincInterpolateImageFunction |
| 74 |
IND |
**/ |
| 75 |
IND |
template< unsigned int VRadius, |
| 76 |
|
class TInput=double, class TOutput=double> |
| 77 |
IND |
class LancosWindowFunction |
| 78 |
IND |
{ |
| 79 |
IND |
public: |
| 80 |
IND |
**inline TOutput operator()( const TInput & A ) const |
| 81 |
IND |
****{ |
| 82 |
IND |
****if(A == 0.0) return (TOutput) 1.0; |
| 83 |
IND |
****double z = m_Factor * A; |
| 84 |
IND |
****return (TOutput) ( sin(z) / z ); |
| 85 |
IND |
****} |
| 86 |
IND |
private: |
| 87 |
IND |
**/** Equal to \f$ \frac{\pi}{m} \f$ */ |
| 88 |
IND |
**static const double m_Factor; |
| 89 |
IND |
}; |
| 90 |
|
|
| 91 |
IND |
/** |
| 92 |
IND |
** \class BlackmanWindowFunction |
| 93 |
IND |
** \brief Window function for sinc interpolation. |
| 94 |
IND |
** \f[ w(x) = 0.42 + 0.5 cos(\frac{\pi x}{m}) + 0.08 cos(\frac{2 \pi x}{m}) \f] |
| 95 |
IND |
** \sa WindowedSincInterpolateImageFunction |
| 96 |
IND |
**/ |
| 97 |
IND |
template< unsigned int VRadius, |
| 98 |
|
class TInput=double, class TOutput=double> |
| 99 |
IND |
class BlackmanWindowFunction |
| 100 |
IND |
{ |
| 101 |
IND |
public: |
| 102 |
IND |
**inline TOutput operator()( const TInput & A ) const |
| 103 |
IND |
****{ |
| 104 |
IND |
****return (TOutput) |
| 105 |
IND |
******(0.42 + 0.5 * cos(A * m_Factor1) + 0.08 * cos(A * m_Factor2)); |
| 106 |
IND |
****} |
| 107 |
IND |
private: |
| 108 |
IND |
**/** Equal to \f$ \frac{\pi}{m} \f$ */ |
| 109 |
IND |
**static const double m_Factor1; |
| 110 |
|
|
| 111 |
IND |
**/** Equal to \f$ \frac{2 \pi}{m} \f$ */ |
| 112 |
IND |
**static const double m_Factor2; |
| 113 |
IND |
}; |
| 114 |
|
|
| 115 |
IND |
} // namespace Function |
| 116 |
|
|
| 117 |
IND |
/** |
| 118 |
IND |
**\class WindowedSincInterpolateImageFunction |
| 119 |
IND |
**\brief Use the windowed sinc function to interpolate |
| 120 |
IND |
**\author Paul A. Yushkevich |
| 121 |
|
|
| 122 |
IND |
**\par THEORY |
| 123 |
|
|
| 124 |
IND |
**This function is intended to provide an interpolation function that |
| 125 |
IND |
**has minimum aliasing artifacts, in contrast to linear interpolation. |
| 126 |
IND |
**According to sampling theory, the infinite-support sinc filter, |
| 127 |
IND |
**whose Fourier transform is the box filter, is optimal for resampling |
| 128 |
IND |
**a function. In practice, the infinite support sinc filter is |
| 129 |
IND |
**approximated using a limited support 'windowed' sinc filter. |
| 130 |
|
|
| 131 |
IND |
**\par |
| 132 |
IND |
**This function is based on the following publication: |
| 133 |
|
|
| 134 |
IND |
**\par |
| 135 |
IND |
**Erik H. W. Meijering, Wiro J. Niessen, Josien P. W. Pluim, |
| 136 |
IND |
**Max A. Viergever: Quantitative Comparison of Sinc-Approximating |
| 137 |
IND |
**Kernels for Medical Image Interpolation. MICCAI 1999, pp. 210-217 |
| 138 |
|
|
| 139 |
IND |
**\par |
| 140 |
IND |
**In this work, several 'windows' are estimated. In two dimensions, the |
| 141 |
IND |
**interpolation at a position (x,y) is given by the following |
| 142 |
IND |
**expression: |
| 143 |
|
|
| 144 |
IND |
**\par |
| 145 |
IND |
**\f[ |
| 146 |
IND |
****I(x,y) = |
| 147 |
IND |
******\sum_{i = \lfloor x \rfloor + 1 - m}^{\lfloor x \rfloor + m} |
| 148 |
IND |
******\sum_{j = \lfloor y \rfloor + 1 - m}^{\lfloor y \rfloor + m} |
| 149 |
IND |
******I_{i,j} K(x-i) K(y-j), |
| 150 |
IND |
**\f] |
| 151 |
|
|
| 152 |
IND |
**\par |
| 153 |
IND |
**where m is the 'radius' of the window, (3,4 are reasonable numbers), |
| 154 |
IND |
**and K(t) is the kernel function, composed of the sinc function and |
| 155 |
IND |
**one of several possible window functions: |
| 156 |
|
|
| 157 |
IND |
**\par |
| 158 |
IND |
**\f[ |
| 159 |
IND |
****K(t) = w(t) \textrm{sinc}(t) = w(t) \frac{\sin(\pi t)}{\pi t} |
| 160 |
IND |
**\f] |
| 161 |
|
|
| 162 |
IND |
**\par |
| 163 |
IND |
**Several window functions are provided here in the itk::Function |
| 164 |
IND |
**namespace. The conclusions of the referenced paper suggest to use the |
| 165 |
IND |
**Welch, Cosine, Kaiser, and Lancos windows for m = 4,5. These are based |
| 166 |
IND |
**on error in rotating medical images w.r.t. the linear interpolation |
| 167 |
IND |
**method. In some cases the results achieve a 20-fold improvement in |
| 168 |
IND |
**accuracy. |
| 169 |
|
|
| 170 |
IND |
**\par USING THIS FILTER |
| 171 |
|
|
| 172 |
IND |
**Use this filter the way you would use any ImageInterpolationFunction, |
| 173 |
IND |
**so for instance, you can plug it into the ResampleImageFilter class. |
| 174 |
IND |
**In order to initialize the filter you must choose several template |
| 175 |
IND |
**parameters. |
| 176 |
|
|
| 177 |
IND |
**\par |
| 178 |
IND |
**The first (TInputImage) is the image type, that's standard. |
| 179 |
|
|
| 180 |
IND |
**\par |
| 181 |
IND |
**The second (VRadius) is the radius of the kernel, i.e., the |
| 182 |
IND |
**\f$ m \f$ from the formula above. |
| 183 |
|
|
| 184 |
IND |
**\par |
| 185 |
IND |
**The third (TWindowFunction) is the window function object, which you |
| 186 |
IND |
**can choose from about five different functions defined in this |
| 187 |
IND |
**header. The default is the Hamming window, which is commonly used |
| 188 |
IND |
**but not optimal according to the cited paper. |
| 189 |
|
|
| 190 |
IND |
**\par |
| 191 |
IND |
**The fourth (TBoundaryCondition) is the boundary condition class used |
| 192 |
IND |
**to determine the values of pixels that fall off the image boundary. |
| 193 |
IND |
**This class has the same meaning here as in the NeighborhoodItetator |
| 194 |
IND |
**classes. |
| 195 |
|
|
| 196 |
IND |
**\par |
| 197 |
IND |
**The fifth (TCoordRep) is again standard for interpolating functions, |
| 198 |
IND |
**and should be float or double. |
| 199 |
|
|
| 200 |
IND |
**\par CAVEATS |
| 201 |
|
|
| 202 |
IND |
**There are a few improvements that an enthusiasting ITK developer |
| 203 |
IND |
**could make to this filter. One issue is with the way that the kernel |
| 204 |
IND |
**is applied. The computational expense comes from two sources: |
| 205 |
IND |
**computing the kernel weights K(t) and multiplying the pixels in the |
| 206 |
IND |
**window by the kernel weights. The first is done more or less |
| 207 |
IND |
**efficiently in \f$ 2 m d \f$ operations (where d is the |
| 208 |
|
dimensionality of the image). The second can be done |
| 209 |
IND |
**better. Presently, each pixel \f$ I(i,j,k) \f$ is multiplied by the |
| 210 |
IND |
**weights \f$ K(x-i), K(y-j), K(z-k) \f$ and added to the running |
| 211 |
IND |
**total. This results in \f$ d (2m)^d \f$ multiplication |
| 212 |
IND |
**operations. However, by keeping intermediate sums, it would be |
| 213 |
IND |
**possible to do the operation in \f$ O ( (2m)^d ) \f$ |
| 214 |
IND |
**operations. This would require some creative coding. In addition, in |
| 215 |
IND |
**the case when one of the coordinates is integer, the computation |
| 216 |
IND |
**could be reduced by an order of magnitude. |
| 217 |
|
|
| 218 |
IND |
**\sa LinearInterpolateImageFunction ResampleImageFilter |
| 219 |
IND |
**\sa Function::HammingWindowFunction |
| 220 |
IND |
**\sa Function::CosineWindowFunction |
| 221 |
IND |
**\sa Function::WelchWindowFunction |
| 222 |
IND |
**\sa Function::LancosWindowFunction |
| 223 |
IND |
**\sa Function::BlackmanWindowFunction |
| 224 |
IND |
**\ingroup ImageFunctions ImageInterpolators |
| 225 |
IND |
**/ |
| 226 |
IND |
template < |
| 227 |
|
class TInputImage, |
| 228 |
|
unsigned int VRadius, |
| 229 |
|
class TWindowFunction = Function::HammingWindowFunction<VRadius>, |
| 230 |
IND |
**class TBoundaryCondition = ConstantBoundaryCondition<TInputImage>, |
| 231 |
IND |
**class TCoordRep=double > |
| 232 |
IND |
class ITK_EXPORT WindowedSincInterpolateImageFunction : |
| 233 |
IND |
**public InterpolateImageFunction<TInputImage, TCoordRep> |
| 234 |
IND |
{ |
| 235 |
IND |
public: |
| 236 |
IND |
**/** Standard class typedefs. */ |
| 237 |
IND |
**typedef WindowedSincInterpolateImageFunction Self; |
| 238 |
TDA,IND |
**typedef InterpolateImageFunction<TInputImage,TCoordRep> Superclass; |
| 239 |
TDA,IND |
**typedef SmartPointer<Self> Pointer; |
| 240 |
TDA,IND |
**typedef SmartPointer<const Self> ConstPointer; |
| 241 |
|
|
| 242 |
IND |
**/** Run-time type information (and related methods). */ |
| 243 |
IND |
**itkTypeMacro(WindowedSincInterpolateImageFunction, |
| 244 |
|
InterpolateImageFunction); |
| 245 |
|
|
| 246 |
IND |
**/** Method for creation through the object factory. */ |
| 247 |
IND |
**itkNewMacro(Self); |
| 248 |
|
|
| 249 |
IND |
**/** OutputType typedef support. */ |
| 250 |
IND |
**typedef typename Superclass::OutputType OutputType; |
| 251 |
|
|
| 252 |
IND |
**/** InputImageType typedef support. */ |
| 253 |
IND |
**typedef typename Superclass::InputImageType InputImageType; |
| 254 |
|
|
| 255 |
IND |
**/** RealType typedef support. */ |
| 256 |
IND |
**typedef typename Superclass::RealType RealType; |
| 257 |
|
|
| 258 |
IND |
**/** Dimension underlying input image. */ |
| 259 |
IND |
**itkStaticConstMacro(ImageDimension, unsigned int,Superclass::ImageDimension); |
| 260 |
|
|
| 261 |
IND |
**/** Index typedef support. */ |
| 262 |
IND |
**typedef typename Superclass::IndexType IndexType; |
| 263 |
|
|
| 264 |
IND |
**/** Image type definition */ |
| 265 |
IND |
**typedef TInputImage ImageType; |
| 266 |
|
|
| 267 |
IND |
**/** ContinuousIndex typedef support. */ |
| 268 |
IND |
**typedef typename Superclass::ContinuousIndexType ContinuousIndexType; |
| 269 |
|
|
| 270 |
IND |
**virtual void SetInputImage(const ImageType *image); |
| 271 |
|
|
| 272 |
IND |
**/** Evaluate the function at a ContinuousIndex position |
| 273 |
IND |
**** |
| 274 |
IND |
**** Returns the interpolated image intensity at a |
| 275 |
IND |
**** specified point position. Bounds checking is based on the |
| 276 |
IND |
**** type of the TBoundaryCondition specified. |
| 277 |
IND |
****/ |
| 278 |
IND |
**virtual OutputType EvaluateAtContinuousIndex( |
| 279 |
|
const ContinuousIndexType & index ) const; |
| 280 |
|
|
| 281 |
IND |
protected: |
| 282 |
IND |
**WindowedSincInterpolateImageFunction(); |
| 283 |
IND |
**virtual ~WindowedSincInterpolateImageFunction(); |
| 284 |
IND |
**void PrintSelf(std::ostream& os, Indent indent) const; |
| 285 |
|
|
| 286 |
IND |
private: |
| 287 |
LEN,IND |
**WindowedSincInterpolateImageFunction( const Self& ); //purposely not implemented |
| 288 |
IND |
**void operator=( const Self& ); //purposely not implemented |
| 289 |
|
|
| 290 |
IND |
**// Internal typedefs |
| 291 |
IND |
**typedef ConstNeighborhoodIterator< |
| 292 |
|
ImageType, TBoundaryCondition> IteratorType; |
| 293 |
|
|
| 294 |
IND |
**// Constant to store twice the radius |
| 295 |
IND |
**static const unsigned int m_WindowSize; |
| 296 |
|
|
| 297 |
IND |
**/** The function object, used to compute window */ |
| 298 |
IND |
**TWindowFunction m_WindowFunction; |
| 299 |
|
|
| 300 |
IND |
**/** The offset array, used to keep a list of relevant |
| 301 |
IND |
**** offsets in the neihborhoodIterator */ |
| 302 |
IND |
**unsigned int *m_OffsetTable; |
| 303 |
|
|
| 304 |
IND |
**/** Size of the offset table */ |
| 305 |
IND |
**unsigned int m_OffsetTableSize; |
| 306 |
|
|
| 307 |
IND |
**/** Index into the weights array for each offset */ |
| 308 |
IND |
**unsigned int **m_WeightOffsetTable; |
| 309 |
|
|
| 310 |
IND |
**/** The sinc function */ |
| 311 |
IND |
**inline double Sinc(double x) const |
| 312 |
IND |
****{ |
| 313 |
IND |
****double px = vnl_math::pi * x; |
| 314 |
IND |
****return (x == 0.0) ? 1.0 : sin(px) / px; |
| 315 |
IND |
****} |
| 316 |
IND |
}; |
| 317 |
|
|
| 318 |
IND |
} // namespace itk |
| 319 |
|
|
| 320 |
|
#ifndef ITK_MANUAL_INSTANTIATION |
| 321 |
|
#include "itkWindowedSincInterpolateImageFunction.txx" |
| 322 |
|
#endif |
| 323 |
|
|
| 324 |
|
#endif // _itkWindowedSincInterpolateImageFunction_h |
| 325 |
|
|
