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 itkBSplineDerivativeKernelFunction_h
19 #define itkBSplineDerivativeKernelFunction_h
20 
21 #include "itkBSplineKernelFunction.h"
22 #include "itkMath.h"
23 
24 namespace itk
25 {
26 /** \class BSplineDerivativeKernelFunction
27  * \brief Derivative of a BSpline kernel used for density estimation and
28  *  nonparameteric regression.
29  *
30  * This class encapsulates the derivative of a 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 1 to 4
36  *
37  * \sa KernelFunctionBase
38  *
39  * \ingroup Functions
40  * \ingroup ITKCommon
41  */
42 template< unsigned int VSplineOrder = 3, typename TRealValueType = double >
43 class BSplineDerivativeKernelFunction:public KernelFunctionBase<TRealValueType>
44 {
45 public:
46   ITK_DISALLOW_COPY_AND_ASSIGN(BSplineDerivativeKernelFunction);
47 
48   /** Standard class type aliases. */
49   using Self = BSplineDerivativeKernelFunction;
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(BSplineDerivativeKernelFunction, 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   BSplineDerivativeKernelFunction() = default;
71   ~BSplineDerivativeKernelFunction() override = default;
72 
PrintSelf(std::ostream & os,Indent indent)73   void PrintSelf(std::ostream & os, Indent indent) const override
74     {
75     Superclass::PrintSelf(os, indent);
76     os << indent  << "Spline Order: " << SplineOrder << std::endl;
77     }
78 
79 private:
80   /** Structures to control overloaded versions of Evaluate */
81   struct DispatchBase {};
82   template< unsigned int >
83   struct Dispatch: public DispatchBase {};
84 
85   /** Evaluate the function:  zeroth order spline. */
Evaluate(const Dispatch<0> &,const TRealValueType & itkNotUsed (u))86   inline TRealValueType Evaluate( const Dispatch<0>&, const TRealValueType & itkNotUsed( u ) )
87     const
88     {
89     return NumericTraits< TRealValueType >::ZeroValue();
90     }
91 
92   /** Evaluate the function:  first order spline */
Evaluate(const Dispatch<1> &,const TRealValueType & u)93   inline TRealValueType Evaluate( const Dispatch<1>&, const TRealValueType& u ) const
94     {
95     if( Math::ExactlyEquals(u, -NumericTraits< TRealValueType >::OneValue()) )
96       {
97       return static_cast< TRealValueType >(0.5);
98       }
99     else if( ( u > -NumericTraits< TRealValueType >::OneValue() ) && ( u < NumericTraits< TRealValueType >::ZeroValue() ) )
100       {
101       return NumericTraits< TRealValueType >::OneValue();
102       }
103     else if( Math::ExactlyEquals(u, NumericTraits< TRealValueType >::ZeroValue()) )
104       {
105       return NumericTraits< TRealValueType >::ZeroValue();
106       }
107     else if( ( u > NumericTraits< TRealValueType >::ZeroValue() ) && ( u < NumericTraits< TRealValueType >::OneValue() ) )
108       {
109       return -NumericTraits< TRealValueType >::OneValue();
110       }
111     else if( Math::ExactlyEquals(u, NumericTraits< TRealValueType >::OneValue()) )
112       {
113       return static_cast< TRealValueType >(-0.5);
114       }
115     else
116       {
117       return NumericTraits< TRealValueType >::ZeroValue();
118       }
119     }
120 
121   /** Evaluate the function:  second order spline. */
Evaluate(const Dispatch<2> &,const TRealValueType & u)122   inline TRealValueType Evaluate( const Dispatch<2>&, const TRealValueType& u) const
123     {
124     if( ( u > static_cast< TRealValueType >(-0.5) ) && ( u < static_cast< TRealValueType >(0.5) ) )
125       {
126       return ( static_cast< TRealValueType >(-2.0) * u );
127       }
128     else if( ( u >= static_cast< TRealValueType >(0.5) ) && ( u < static_cast< TRealValueType >(1.5) ) )
129       {
130       return ( static_cast< TRealValueType >(-1.5) + u );
131       }
132     else if( ( u > static_cast< TRealValueType >(-1.5) ) && ( u <= static_cast< TRealValueType >(-0.5) ) )
133       {
134       return ( static_cast< TRealValueType >(1.5) + u );
135       }
136     else
137       {
138       return NumericTraits< TRealValueType >::ZeroValue();
139       }
140     }
141 
142   /** Evaluate the function:  third order spline. */
Evaluate(const Dispatch<3> &,const TRealValueType & u)143   inline TRealValueType Evaluate( const Dispatch<3>&, const TRealValueType& u ) const
144     {
145     if( ( u >= NumericTraits< TRealValueType >::ZeroValue() ) && ( u < NumericTraits< TRealValueType >::OneValue() ) )
146       {
147       return ( static_cast< TRealValueType >(-2.0)* u + static_cast< TRealValueType >(1.5) * u * u );
148       }
149     else if( ( u > -NumericTraits< TRealValueType >::OneValue() ) && ( u < NumericTraits< TRealValueType >::ZeroValue() ) )
150       {
151       return ( static_cast< TRealValueType >(-2.0) * u - static_cast< TRealValueType >(1.5) * u * u );
152       }
153     else if( ( u >= NumericTraits< TRealValueType >::OneValue() ) && ( u < static_cast< TRealValueType >(2.0) ) )
154       {
155       return ( static_cast< TRealValueType >(-2.0) + static_cast< TRealValueType >(2.0) * u - static_cast< TRealValueType >(0.5) * u * u );
156       }
157     else if( ( u > static_cast< TRealValueType >(-2.0) ) && ( u <= -NumericTraits< TRealValueType >::OneValue() ) )
158       {
159       return ( static_cast< TRealValueType >(2.0) + static_cast< TRealValueType >(2.0) * u + static_cast< TRealValueType >(0.5) * u * u );
160       }
161     else
162       {
163       return NumericTraits< TRealValueType >::ZeroValue();
164       }
165     }
166 
167   /** Evaluate the function:  unimplemented spline order */
Evaluate(const DispatchBase &,const TRealValueType &)168   inline TRealValueType Evaluate( const DispatchBase&, const TRealValueType& ) const
169     {
170     itkExceptionMacro( "Evaluate not implemented for spline order " << SplineOrder );
171     }
172 };
173 } // end namespace itk
174 
175 #endif
176