Line data    Source code

       1             : #include "FiniteDifferences.h"
       2             : #include "Timer.h"
       3             : #include "VolumeDescriptor.h"
       4             : #include "TypeCasts.hpp"
       5             : 
       6             : namespace elsa
       7             : {
       8             :     template <typename data_t>
       9           0 :     FiniteDifferences<data_t>::FiniteDifferences(const DataDescriptor& domainDescriptor,
      10             :                                                  DiffType type)
      11             :         : FiniteDifferences(domainDescriptor,
      12           0 :                             BooleanVector_t::Ones(domainDescriptor.getNumberOfDimensions()), type)
      13             :     {
      14           0 :     }
      15             : 
      16             :     template <typename data_t>
      17           0 :     FiniteDifferences<data_t>::FiniteDifferences(const DataDescriptor& domainDescriptor,
      18             :                                                  const BooleanVector_t& activeDims, DiffType type)
      19             :         : LinearOperator<data_t>(domainDescriptor,
      20             :                                  domainDescriptor), // setting range in body of constructor
      21           0 :           _type{type},
      22           0 :           _activeDims{activeDims},
      23           0 :           _coordDiff{activeDims.size()},
      24           0 :           _coordDelta{activeDims.size()},
      25           0 :           _dimCounter{activeDims.size()}
      26             :     {
      27             :         // build the range descriptor of appropriate size
      28           0 :         IndexVector_t coefficients(domainDescriptor.getNumberOfDimensions() + 1);
      29           0 :         coefficients << domainDescriptor.getNumberOfCoefficientsPerDimension(),
      30           0 :             activeDims.cast<index_t>().sum();
      31             : 
      32           0 :         RealVector_t spacing(domainDescriptor.getNumberOfDimensions() + 1);
      33           0 :         spacing << domainDescriptor.getSpacingPerDimension(), 1;
      34             : 
      35           0 :         this->_rangeDescriptor = std::make_unique<VolumeDescriptor>(coefficients, spacing);
      36             : 
      37           0 :         precomputeHelpers();
      38           0 :     }
      39             : 
      40             :     template <typename data_t>
      41           0 :     void FiniteDifferences<data_t>::precomputeHelpers()
      42             :     {
      43           0 :         IndexVector_t numberOfCoefficients =
      44             :             this->_rangeDescriptor->getNumberOfCoefficientsPerDimension();
      45             : 
      46           0 :         index_t deltaTmp = 1;
      47           0 :         int count = -1;
      48           0 :         for (index_t ic = 0; ic < this->getDomainDescriptor().getNumberOfDimensions(); ++ic) {
      49           0 :             _coordDiff[ic] = numberOfCoefficients.head(ic).prod();
      50             : 
      51           0 :             deltaTmp *= numberOfCoefficients[ic];
      52           0 :             _coordDelta[ic] = deltaTmp;
      53             : 
      54           0 :             if (_activeDims[ic])
      55           0 :                 ++count;
      56           0 :             _dimCounter[ic] = count;
      57             :         }
      58           0 :     }
      59             : 
      60             :     template <typename data_t>
      61           0 :     void FiniteDifferences<data_t>::applyImpl(const DataContainer<data_t>& x,
      62             :                                               DataContainer<data_t>& Ax) const
      63             :     {
      64           0 :         Timer timeguard("FiniteDifferences", "apply");
      65             : 
      66           0 :         Ax = 0;
      67             : 
      68           0 :         switch (_type) {
      69           0 :             case DiffType::FORWARD:
      70           0 :                 applyHelper(x, Ax, DiffType::FORWARD);
      71           0 :                 break;
      72           0 :             case DiffType::BACKWARD:
      73           0 :                 applyHelper(x, Ax, DiffType::BACKWARD);
      74           0 :                 break;
      75           0 :             case DiffType::CENTRAL:
      76           0 :                 applyHelper(x, Ax, DiffType::CENTRAL);
      77           0 :                 break;
      78           0 :             default:
      79           0 :                 throw LogicError("FiniteDifferences::apply: invalid DiffType");
      80             :         }
      81           0 :     }
      82             : 
      83             :     template <typename data_t>
      84           0 :     void FiniteDifferences<data_t>::applyAdjointImpl(const DataContainer<data_t>& y,
      85             :                                                      DataContainer<data_t>& Aty) const
      86             :     {
      87           0 :         Timer timeguard("FiniteDifferences", "applyAdjoint");
      88             : 
      89           0 :         Aty = 0;
      90             : 
      91           0 :         switch (_type) {
      92           0 :             case DiffType::FORWARD:
      93           0 :                 applyAdjointHelper(y, Aty, DiffType::FORWARD);
      94           0 :                 break;
      95           0 :             case DiffType::BACKWARD:
      96           0 :                 applyAdjointHelper(y, Aty, DiffType::BACKWARD);
      97           0 :                 break;
      98           0 :             case DiffType::CENTRAL:
      99           0 :                 applyAdjointHelper(y, Aty, DiffType::CENTRAL);
     100           0 :                 break;
     101           0 :             default:
     102           0 :                 throw LogicError("FiniteDifferences::applyAdjoint: invalid DiffType");
     103             :         }
     104           0 :     }
     105             : 
     106             :     template <typename data_t>
     107             :     template <typename FDtype>
     108           0 :     void FiniteDifferences<data_t>::applyHelper(const DataContainer<data_t>& x,
     109             :                                                 DataContainer<data_t>& Ax, FDtype type) const
     110             :     {
     111           0 :         index_t sizeOfDomain = this->getDomainDescriptor().getNumberOfCoefficients();
     112           0 :         index_t numDim = this->getDomainDescriptor().getNumberOfDimensions();
     113             : 
     114           0 :         IndexVector_t numberOfCoefficients =
     115           0 :             this->getRangeDescriptor().getNumberOfCoefficientsPerDimension();
     116           0 :         IndexVector_t decrementedCoefficients =
     117             :             numberOfCoefficients
     118           0 :             - IndexVector_t::Ones(this->getRangeDescriptor().getNumberOfDimensions());
     119             : 
     120           0 : #pragma omp parallel
     121             :         for (long currDim = 0; currDim < numDim; ++currDim) {
     122             :             if (!_activeDims[currDim])
     123             :                 continue;
     124             : 
     125             :             index_t modulus = numberOfCoefficients.head(currDim + 1).prod();
     126             :             index_t divisor = numberOfCoefficients.head(currDim).prod();
     127             : 
     128             : #pragma omp for nowait
     129             :             for (index_t id = 0; id < sizeOfDomain; ++id) {
     130             :                 index_t icCount = (id % modulus) / divisor; //_domainDescriptor.index(id, ic);
     131             :                 index_t ir = id + _dimCounter[currDim] * _coordDelta[currDim];
     132             : 
     133             :                 // store result depending on mode
     134             :                 switch (type) {
     135             :                     case DiffType::FORWARD:
     136             :                         Ax[ir] = -x[id];
     137             :                         if (icCount < decrementedCoefficients[currDim])
     138             :                             Ax[ir] += x[id + _coordDiff[currDim]];
     139             :                         break;
     140             :                     case DiffType::BACKWARD:
     141             :                         Ax[ir] = x[id];
     142             :                         if (icCount > 0)
     143             :                             Ax[ir] -= x[id - _coordDiff[currDim]];
     144             :                         break;
     145             :                     case DiffType::CENTRAL:
     146             :                         Ax[ir] = static_cast<data_t>(0.0);
     147             :                         if (icCount < decrementedCoefficients[currDim])
     148             :                             Ax[ir] += static_cast<data_t>(0.5) * x[id + _coordDiff[currDim]];
     149             :                         if (icCount > 0)
     150             :                             Ax[ir] -= static_cast<data_t>(0.5) * x[id - _coordDiff[currDim]];
     151             :                         break;
     152             :                 }
     153             :             }
     154             :         }
     155           0 :     }
     156             : 
     157             :     template <typename data_t>
     158             :     template <typename FDtype>
     159           0 :     void FiniteDifferences<data_t>::applyAdjointHelper(const DataContainer<data_t>& y,
     160             :                                                        DataContainer<data_t>& Aty,
     161             :                                                        FDtype type) const
     162             :     {
     163           0 :         index_t sizeOfDomain = this->getDomainDescriptor().getNumberOfCoefficients();
     164           0 :         index_t numDim = this->getDomainDescriptor().getNumberOfDimensions();
     165             : 
     166           0 :         IndexVector_t numberOfCoefficients =
     167           0 :             this->getDomainDescriptor().getNumberOfCoefficientsPerDimension();
     168           0 :         IndexVector_t decrementedCoefficients = numberOfCoefficients - IndexVector_t::Ones(numDim);
     169             : 
     170           0 : #pragma omp parallel
     171             :         for (long currDim = 0; currDim < numDim; ++currDim) {
     172             :             if (!_activeDims[currDim])
     173             :                 continue;
     174             : 
     175             :             index_t modulus = numberOfCoefficients.head(currDim + 1).prod();
     176             :             index_t divisor = numberOfCoefficients.head(currDim).prod();
     177             : 
     178             : #pragma omp for nowait
     179             :             for (index_t id = 0; id < sizeOfDomain; ++id) {
     180             :                 index_t icCount = (id % modulus) / divisor;
     181             :                 index_t ir = id + _dimCounter[currDim] * _coordDelta[currDim];
     182             : 
     183             :                 switch (type) {
     184             :                     case DiffType::FORWARD:
     185             :                         if (icCount > 0)
     186             :                             Aty[id] += y[ir - _coordDiff[currDim]];
     187             : 
     188             :                         Aty[id] -= y[ir];
     189             :                         break;
     190             :                     case DiffType::BACKWARD:
     191             :                         if (icCount < decrementedCoefficients(currDim))
     192             :                             Aty[id] -= y[ir + _coordDiff[currDim]];
     193             : 
     194             :                         Aty[id] += y[ir];
     195             :                         break;
     196             :                     case DiffType::CENTRAL:
     197             :                         if (icCount > 0)
     198             :                             Aty[ir] += static_cast<data_t>(0.5) * y[ir - _coordDiff[currDim]];
     199             : 
     200             :                         if (icCount < decrementedCoefficients(currDim))
     201             :                             Aty[ir] -= static_cast<data_t>(0.5) * y[ir + _coordDiff[currDim]];
     202             :                         break;
     203             :                 }
     204             :             }
     205             :         }
     206           0 :     }
     207             : 
     208             :     template <typename data_t>
     209           0 :     FiniteDifferences<data_t>* FiniteDifferences<data_t>::cloneImpl() const
     210             :     {
     211           0 :         return new FiniteDifferences(this->getDomainDescriptor(), _activeDims, _type);
     212             :     }
     213             : 
     214             :     template <typename data_t>
     215           0 :     bool FiniteDifferences<data_t>::isEqual(const LinearOperator<data_t>& other) const
     216             :     {
     217           0 :         if (!LinearOperator<data_t>::isEqual(other))
     218           0 :             return false;
     219             : 
     220           0 :         auto otherFD = downcast_safe<FiniteDifferences>(&other);
     221           0 :         if (!otherFD)
     222           0 :             return false;
     223             : 
     224           0 :         if (_type != otherFD->_type || _activeDims != otherFD->_activeDims)
     225           0 :             return false;
     226             : 
     227           0 :         return true;
     228             :     }
     229             : 
     230             :     // ------------------------------------------
     231             :     // explicit template instantiation
     232             :     template class FiniteDifferences<float>;
     233             :     template class FiniteDifferences<double>;
     234             :     template class FiniteDifferences<complex<float>>;
     235             :     template class FiniteDifferences<complex<double>>;
     236             : 
     237             : } // namespace elsa