Line data    Source code

       1             : #include "PseudoHuber.h"
       2             : #include "DataContainer.h"
       3             : #include "Scaling.h"
       4             : #include "TypeCasts.hpp"
       5             : 
       6             : #include <cmath>
       7             : #include <stdexcept>
       8             : 
       9             : namespace elsa
      10             : {
      11             :     template <typename data_t>
      12             :     PseudoHuber<data_t>::PseudoHuber(const DataDescriptor& domainDescriptor, real_t delta)
      13             :         : Functional<data_t>(domainDescriptor), delta_{delta}
      14           8 :     {
      15             :         // sanity check delta
      16           8 :         if (delta <= static_cast<real_t>(0.0))
      17           0 :             throw InvalidArgumentError("PseudoHuber: delta has to be positive.");
      18           8 :     }
      19             : 
      20             :     template <typename data_t>
      21             :     bool PseudoHuber<data_t>::isDifferentiable() const
      22           0 :     {
      23           0 :         return true;
      24           0 :     }
      25             : 
      26             :     template <typename data_t>
      27             :     data_t PseudoHuber<data_t>::evaluateImpl(const DataContainer<data_t>& Rx) const
      28           2 :     {
      29             :         // note: this is currently not a reduction in DataContainer, but implemented here "manually"
      30             : 
      31           2 :         auto result = static_cast<data_t>(0.0);
      32             : 
      33         886 :         for (index_t i = 0; i < Rx.getSize(); ++i) {
      34         884 :             data_t temp = Rx[i] / delta_;
      35         884 :             result += delta_ * delta_
      36         884 :                       * (sqrt(static_cast<data_t>(1.0) + temp * temp) - static_cast<data_t>(1.0));
      37         884 :         }
      38             : 
      39           2 :         return result;
      40           2 :     }
      41             : 
      42             :     template <typename data_t>
      43             :     void PseudoHuber<data_t>::getGradientImpl(const DataContainer<data_t>& Rx,
      44             :                                               DataContainer<data_t>& out) const
      45           2 :     {
      46         886 :         for (index_t i = 0; i < Rx.getSize(); ++i) {
      47         884 :             data_t temp = Rx[i] / delta_;
      48         884 :             out[i] = Rx[i] / sqrt(static_cast<data_t>(1.0) + temp * temp);
      49         884 :         }
      50           2 :     }
      51             : 
      52             :     template <typename data_t>
      53             :     LinearOperator<data_t>
      54             :         PseudoHuber<data_t>::getHessianImpl(const DataContainer<data_t>& Rx) const
      55           2 :     {
      56           2 :         DataContainer<data_t> scaleFactors(Rx.getDataDescriptor());
      57         886 :         for (index_t i = 0; i < Rx.getSize(); ++i) {
      58         884 :             data_t temp = Rx[i] / delta_;
      59         884 :             data_t tempSq = temp * temp;
      60         884 :             data_t sqrtOnePTempSq = sqrt(static_cast<data_t>(1.0) + tempSq);
      61         884 :             scaleFactors[i] =
      62         884 :                 (sqrtOnePTempSq - tempSq / sqrtOnePTempSq) / (static_cast<data_t>(1.0) + tempSq);
      63         884 :         }
      64             : 
      65           2 :         return leaf(Scaling<data_t>(scaleFactors));
      66           2 :     }
      67             : 
      68             :     template <typename data_t>
      69             :     PseudoHuber<data_t>* PseudoHuber<data_t>::cloneImpl() const
      70           2 :     {
      71           2 :         return new PseudoHuber(this->getDomainDescriptor(), delta_);
      72           2 :     }
      73             : 
      74             :     template <typename data_t>
      75             :     bool PseudoHuber<data_t>::isEqual(const Functional<data_t>& other) const
      76           2 :     {
      77           2 :         if (!Functional<data_t>::isEqual(other))
      78           0 :             return false;
      79             : 
      80           2 :         auto otherPHuber = downcast_safe<PseudoHuber>(&other);
      81           2 :         if (!otherPHuber)
      82           0 :             return false;
      83             : 
      84           2 :         return delta_ == otherPHuber->delta_;
      85           2 :     }
      86             : 
      87             :     // ------------------------------------------
      88             :     // explicit template instantiation
      89             :     template class PseudoHuber<float>;
      90             :     template class PseudoHuber<double>;
      91             :     // no complex-number instantiations for PseudoHuber! (they would not really be useful)
      92             : 
      93             : } // namespace elsa