Line data    Source code

       1             : /**
       2             :  * @file test_Huber.cpp
       3             :  *
       4             :  * @brief Tests for the Huber class
       5             :  *
       6             :  * @author Matthias Wieczorek - initial code
       7             :  * @author David Frank - rewrite
       8             :  * @author Tobias Lasser - modernization
       9             :  */
      10             : 
      11             : #include <doctest/doctest.h>
      12             : 
      13             : #include "Huber.h"
      14             : #include "LinearResidual.h"
      15             : #include "Identity.h"
      16             : #include "VolumeDescriptor.h"
      17             : #include "TypeCasts.hpp"
      18             : #include "testHelpers.h"
      19             : 
      20             : using namespace elsa;
      21             : using namespace doctest;
      22             : 
      23             : TYPE_TO_STRING(complex<float>);
      24             : TYPE_TO_STRING(complex<double>);
      25             : 
      26             : TEST_SUITE_BEGIN("functionals");
      27             : 
      28             : TEST_CASE_TEMPLATE("Huber: Testing without residual only data", TestType, float, double)
      29          10 : {
      30          10 :     using Vector = Eigen::Matrix<TestType, Eigen::Dynamic, 1>;
      31             : 
      32          10 :     GIVEN("just data (no residual)")
      33          10 :     {
      34          10 :         IndexVector_t numCoeff(3);
      35          10 :         numCoeff << 7, 16, 29;
      36          10 :         VolumeDescriptor dd(numCoeff);
      37             : 
      38          10 :         real_t delta = 10;
      39             : 
      40          10 :         WHEN("instantiating")
      41          10 :         {
      42          10 :             Huber<TestType> func(dd, delta);
      43             : 
      44          10 :             THEN("the functional is as expected")
      45          10 :             {
      46           2 :                 REQUIRE_EQ(func.getDomainDescriptor(), dd);
      47             : 
      48           2 :                 auto* linRes = downcast_safe<LinearResidual<TestType>>(&func.getResidual());
      49           2 :                 REQUIRE_UNARY(linRes);
      50           2 :                 REQUIRE_UNARY_FALSE(linRes->hasDataVector());
      51           2 :                 REQUIRE_UNARY_FALSE(linRes->hasOperator());
      52           2 :             }
      53             : 
      54          10 :             THEN("a clone behaves as expected")
      55          10 :             {
      56           2 :                 auto huberClone = func.clone();
      57             : 
      58           2 :                 REQUIRE_NE(huberClone.get(), &func);
      59           2 :                 REQUIRE_EQ(*huberClone, func);
      60           2 :             }
      61             : 
      62          10 :             Vector dataVec(dd.getNumberOfCoefficients());
      63          10 :             dataVec.setRandom();
      64             : 
      65             :             // fix the first entries to be bigger/smaller than delta
      66          10 :             dataVec[0] = delta + 1;
      67          10 :             dataVec[1] = delta + 2;
      68          10 :             dataVec[2] = delta + 3;
      69          10 :             dataVec[3] = delta - 1;
      70          10 :             dataVec[4] = delta - 2;
      71          10 :             dataVec[5] = delta - 3;
      72             : 
      73          10 :             DataContainer<TestType> x(dd, dataVec);
      74             : 
      75             :             // compute the "true" values
      76          10 :             TestType trueValue = 0;
      77          10 :             Vector trueGrad(dd.getNumberOfCoefficients());
      78       32490 :             for (index_t i = 0; i < dataVec.size(); ++i) {
      79       32480 :                 TestType value = dataVec[i];
      80       32480 :                 if (std::abs(value) <= delta) {
      81       32450 :                     trueValue += 0.5f * value * value;
      82       32450 :                     trueGrad[i] = value;
      83       32450 :                 } else {
      84          30 :                     trueValue += delta * (std::abs(value) - 0.5f * delta);
      85          30 :                     trueGrad[i] = (value > 0) ? delta : -delta;
      86          30 :                 }
      87       32480 :             }
      88             : 
      89          10 :             THEN("the evaluate works as expected")
      90          10 :             {
      91           2 :                 REQUIRE_UNARY(checkApproxEq(func.evaluate(x), trueValue));
      92           2 :             }
      93          10 :             THEN("the gradient works as expected")
      94          10 :             {
      95           2 :                 DataContainer<TestType> dcTrueGrad(dd, trueGrad);
      96           2 :                 REQUIRE_UNARY(checkApproxEq(func.getGradient(x), dcTrueGrad));
      97           2 :             }
      98          10 :             THEN("the Hessian works as expected")
      99          10 :             {
     100           2 :                 auto hessian = func.getHessian(x);
     101           2 :                 auto hx = hessian.apply(x);
     102        6498 :                 for (index_t i = 0; i < hx.getSize(); ++i)
     103           2 :                     REQUIRE_UNARY(
     104           2 :                         checkApproxEq(hx[i], ((std::abs(dataVec[i]) <= delta) ? x[i] : 0)));
     105           2 :             }
     106          10 :         }
     107          10 :     }
     108          10 : }
     109             : 
     110             : TEST_CASE_TEMPLATE("Huber<TestType>: Testing with residual", TestType, float, double)
     111          10 : {
     112          10 :     using Vector = Eigen::Matrix<TestType, Eigen::Dynamic, 1>;
     113             : 
     114          10 :     GIVEN("a residual with data")
     115          10 :     {
     116             :         // linear residual
     117          10 :         IndexVector_t numCoeff(2);
     118          10 :         numCoeff << 47, 11;
     119          10 :         VolumeDescriptor dd(numCoeff);
     120             : 
     121          10 :         Vector randomData(dd.getNumberOfCoefficients());
     122          10 :         randomData.setRandom();
     123          10 :         DataContainer<TestType> b(dd, randomData);
     124             : 
     125          10 :         Identity<TestType> A(dd);
     126             : 
     127          10 :         LinearResidual<TestType> linRes(A, b);
     128             : 
     129          10 :         real_t delta = 20;
     130             : 
     131          10 :         WHEN("instantiating")
     132          10 :         {
     133          10 :             Huber<TestType> func(linRes, delta);
     134             : 
     135          10 :             THEN("the functional is as expected")
     136          10 :             {
     137           2 :                 REQUIRE_EQ(func.getDomainDescriptor(), dd);
     138             : 
     139           2 :                 auto* lRes = downcast_safe<LinearResidual<TestType>>(&func.getResidual());
     140           2 :                 REQUIRE_UNARY(lRes);
     141           2 :                 REQUIRE_EQ(*lRes, linRes);
     142           2 :             }
     143             : 
     144          10 :             THEN("a clone behaves as expected")
     145          10 :             {
     146           2 :                 auto huberClone = func.clone();
     147             : 
     148           2 :                 REQUIRE_NE(huberClone.get(), &func);
     149           2 :                 REQUIRE_EQ(*huberClone, func);
     150           2 :             }
     151             : 
     152          10 :             Vector dataVec(dd.getNumberOfCoefficients());
     153          10 :             dataVec.setRandom();
     154          10 :             DataContainer<TestType> x(dd, dataVec);
     155             : 
     156             :             // compute the "true" values
     157          10 :             TestType trueValue = 0;
     158          10 :             Vector trueGrad(dd.getNumberOfCoefficients());
     159        5180 :             for (index_t i = 0; i < dataVec.size(); ++i) {
     160        5170 :                 TestType value = dataVec[i] - randomData[i];
     161        5170 :                 if (std::abs(value) <= delta) {
     162        5170 :                     trueValue += 0.5f * value * value;
     163        5170 :                     trueGrad[i] = value;
     164        5170 :                 } else {
     165           0 :                     trueValue += delta * (std::abs(value) - 0.5f * delta);
     166           0 :                     trueGrad[i] = (value > 0) ? delta : -delta;
     167           0 :                 }
     168        5170 :             }
     169             : 
     170          10 :             THEN("the evaluate works as expected")
     171          10 :             {
     172           2 :                 REQUIRE_UNARY(checkApproxEq(func.evaluate(x), trueValue));
     173           2 :             }
     174             : 
     175          10 :             THEN("the gradient works as expected")
     176          10 :             {
     177           2 :                 DataContainer<TestType> dcTrueGrad(dd, trueGrad);
     178           2 :                 REQUIRE_UNARY(isApprox(func.getGradient(x), dcTrueGrad));
     179           2 :             }
     180          10 :             THEN("the Hessian works as expected")
     181          10 :             {
     182             : 
     183           2 :                 auto hessian = func.getHessian(x);
     184           2 :                 auto hx = hessian.apply(x);
     185        1036 :                 for (index_t i = 0; i < hx.getSize(); ++i)
     186           2 :                     REQUIRE(hx[i] == ((std::abs(dataVec[i] - randomData[i]) <= delta) ? x[i] : 0));
     187           2 :             }
     188          10 :         }
     189          10 :     }
     190          10 : }
     191             : 
     192             : TEST_SUITE_END();