1 /*=========================================================================
2  *
3  *  Copyright Insight Software Consortium
4  *
5  *  Licensed under the Apache License, Version 2.0 (the "License");
6  *  you may not use this file except in compliance with the License.
7  *  You may obtain a copy of the License at
8  *
9  *         http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  *  Unless required by applicable law or agreed to in writing, software
12  *  distributed under the License is distributed on an "AS IS" BASIS,
13  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  *  See the License for the specific language governing permissions and
15  *  limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef itkBSplineKernelFunction_h
19 #define itkBSplineKernelFunction_h
20 
21 #include "itkKernelFunctionBase.h"
22 #include "itkMath.h"
23 
24 namespace itk
25 {
26 /** \class BSplineKernelFunction
27  * \brief BSpline kernel used for density estimation and nonparameteric
28  *  regression.
29  *
30  * This class enscapsulates BSpline kernel for
31  * density estimation or nonparameteric regression.
32  * See documentation for KernelFunctionBase for more details.
33  *
34  * This class is templated over the spline order.
35  * \warning Evaluate is only implemented for spline order 0 to 3
36  *
37  * \sa KernelFunctionBase
38  *
39  * \ingroup Functions
40  * \ingroup ITKCommon
41  */
42 template< unsigned int VSplineOrder = 3, typename TRealValueType = double >
43 class ITK_TEMPLATE_EXPORT BSplineKernelFunction:public KernelFunctionBase<TRealValueType>
44 {
45 public:
46   ITK_DISALLOW_COPY_AND_ASSIGN(BSplineKernelFunction);
47 
48   /** Standard class type aliases. */
49   using Self = BSplineKernelFunction;
50   using Superclass = KernelFunctionBase<TRealValueType>;
51   using Pointer = SmartPointer< Self >;
52 
53   using RealType = typename Superclass::RealType;
54   /** Method for creation through the object factory. */
55   itkNewMacro(Self);
56 
57   /** Run-time type information (and related methods). */
58   itkTypeMacro(BSplineKernelFunction, KernelFunctionBase);
59 
60   /** Enum of for spline order. */
61   static constexpr unsigned int SplineOrder = VSplineOrder;
62 
63   /** Evaluate the function. */
Evaluate(const TRealValueType & u)64   TRealValueType Evaluate(const TRealValueType & u) const override
65   {
66     return this->Evaluate(Dispatch< VSplineOrder >(), u);
67   }
68 
69 protected:
70   BSplineKernelFunction()= default;
71   ~BSplineKernelFunction() override = default;
PrintSelf(std::ostream & os,Indent indent)72   void PrintSelf(std::ostream & os, Indent indent) const override
73   {
74     Superclass::PrintSelf(os, indent);
75     os << indent  << "Spline Order: " << SplineOrder << std::endl;
76   }
77 
78 private:
79   /** Structures to control overloaded versions of Evaluate */
80   struct DispatchBase {};
81   template< unsigned int >
82   struct Dispatch: public DispatchBase {};
83 
84   /** Zeroth order spline. */
Evaluate(const Dispatch<0> &,const TRealValueType & u)85   inline TRealValueType Evaluate(const Dispatch< 0 > &, const TRealValueType & u) const
86   {
87     const TRealValueType absValue = itk::Math::abs(u);
88     if ( absValue  < static_cast< TRealValueType >(0.5) )
89       {
90       return NumericTraits< TRealValueType >::OneValue();
91       }
92     else if ( Math::ExactlyEquals(absValue, static_cast< TRealValueType >(0.5)) )
93       {
94       return static_cast< TRealValueType >(0.5);
95       }
96     else
97       {
98       return NumericTraits< TRealValueType >::ZeroValue();
99       }
100   }
101 
102   /** First order spline */
Evaluate(const Dispatch<1> &,const TRealValueType & u)103   inline TRealValueType Evaluate(const Dispatch< 1 > &, const TRealValueType & u) const
104   {
105     const TRealValueType absValue = itk::Math::abs(u);
106     if ( absValue  < NumericTraits< TRealValueType >::OneValue() )
107       {
108       return NumericTraits< TRealValueType >::OneValue() - absValue;
109       }
110     else
111       {
112       return NumericTraits< TRealValueType >::ZeroValue();
113       }
114   }
115 
116   /** Second order spline. */
Evaluate(const Dispatch<2> &,const TRealValueType & u)117   inline TRealValueType Evaluate(const Dispatch< 2 > &, const TRealValueType & u) const
118   {
119     const TRealValueType absValue = itk::Math::abs(u);
120     if ( absValue  < static_cast< TRealValueType >(0.5) )
121       {
122       const TRealValueType sqrValue = itk::Math::sqr(absValue);
123       return static_cast< TRealValueType >(0.75) - sqrValue;
124       }
125     else if ( absValue < static_cast< TRealValueType >(1.5) )
126       {
127       const TRealValueType sqrValue = itk::Math::sqr(absValue);
128       // NOTE: 1.0/8.0 == static_cast< TRealValueType >( 0.125 )
129       return ( static_cast< TRealValueType >(9.0) - static_cast< TRealValueType >(12.0) * absValue
130         + static_cast< TRealValueType >(4.0) * sqrValue ) * static_cast< TRealValueType >(0.125);
131       }
132     else
133       {
134       return NumericTraits< TRealValueType >::ZeroValue();
135       }
136   }
137 
138   /**  Third order spline. */
Evaluate(const Dispatch<3> &,const TRealValueType & u)139   inline TRealValueType Evaluate(const Dispatch< 3 > &, const TRealValueType & u) const
140   {
141     const TRealValueType absValue = itk::Math::abs(u);
142     if ( absValue  < NumericTraits< TRealValueType >::OneValue() )
143       {
144       const TRealValueType sqrValue = itk::Math::sqr(absValue);
145       return ( static_cast< TRealValueType >(4.0) - static_cast< TRealValueType >(6.0) * sqrValue
146         + static_cast< TRealValueType >(3.0) * sqrValue * absValue ) / static_cast< TRealValueType >(6.0);
147       }
148     else if ( absValue < static_cast< TRealValueType >(2.0) )
149       {
150       const TRealValueType sqrValue = itk::Math::sqr(absValue);
151       return ( static_cast< TRealValueType >(8.0) - static_cast< TRealValueType >(12.0) * absValue + static_cast< TRealValueType >(6.0) * sqrValue
152                - sqrValue * absValue ) / static_cast< TRealValueType >(6.0);
153       }
154     else
155       {
156       return NumericTraits< TRealValueType >::ZeroValue();
157       }
158   }
159 
160   /** Unimplemented spline order */
Evaluate(const DispatchBase &,const TRealValueType &)161   inline TRealValueType Evaluate(const DispatchBase &, const TRealValueType &) const
162   {
163     itkExceptionMacro( "Evaluate not implemented for spline order "
164       << SplineOrder);
165   }
166 };
167 } // end namespace itk
168 
169 #endif
170