Function iterDeconv

src/api/c/deconvolution.cpp:156–196

template<typename InputType, typename RealType = float>
af_array iterDeconv(const af_array in, const af_array ker, const uint iters,
                    const float rfactor, const af_iterative_deconv_algo algo) {
    using T    = RealType;
    using CT   = typename std::conditional<std::is_same<T, double>::value,
                                         cdouble, cfloat>::type;
    auto input = castArray<T>(in);
    auto psf   = castArray<T>(ker);
    const dim4& idims = input.dims();
    const dim4& fdims = psf.dims();
    dim_t nElems      = 1;

    dim4 inUPad, psfUPad, inLPad, psfLPad, odims(1);

    auto index = calcPadInfo(inLPad, psfLPad, inUPad, psfUPad, odims, nElems,
                             idims, fdims);
    auto paddedIn =
        padArrayBorders<T>(input, inLPad, inUPad, AF_PAD_CLAMP_TO_EDGE);
    auto paddedPsf = padArrayBorders<T>(psf, psfLPad, psfUPad, AF_PAD_ZERO);

    const std::array<int, 4> shiftDims = {-int(fdims[0] / 2),
                                          -int(fdims[1] / 2), 0, 0};
    auto shiftedPsf                    = shift(paddedPsf, shiftDims.data());

    auto P  = fft_r2c<CT, T>(shiftedPsf, BASE_DIM);
    auto Pc = conj(P);

    Array<T> currentEstimate = paddedIn;
    const double normFactor  = 1 / static_cast<double>(nElems);

    switch (algo) {
        case AF_ITERATIVE_DECONV_RICHARDSONLUCY:
            richardsonLucy(currentEstimate, paddedIn, P, Pc, iters, normFactor,
                           odims);
            break;
        case AF_ITERATIVE_DECONV_LANDWEBER:
        default:
            landweber(currentEstimate, paddedIn, P, Pc, iters, rfactor,
                      normFactor, odims);
    }
    return getHandle(createSubArray<T>(currentEstimate, index));
}

af_err af_iterative_deconv(af_array* out, const af_array in, const af_array ker,
                           const unsigned iterations, const float relax_factor,