Radon2D.h

Go to the documentation of this file.

 
#pragma once
 
#include "Common.h"
 
namespace reg23 {
 
at::Tensor Radon2D_CPU(const at::Tensor &image, const at::Tensor &imageSpacing, const at::Tensor &phiValues,
                       const at::Tensor &rValues, int64_t samplesPerLine);
 
at::Tensor DRadon2DDR_CPU(const at::Tensor &image, const at::Tensor &imageSpacing, const at::Tensor &phiValues,
                          const at::Tensor &rValues, int64_t samplesPerLine);
 
__host__ at::Tensor Radon2D_CUDA(const at::Tensor &image, const at::Tensor &imageSpacing, const at::Tensor &phiValues,
                                 const at::Tensor &rValues, int64_t samplesPerLine);
 
__host__ at::Tensor Radon2D_CUDA_V2(const at::Tensor &image, const at::Tensor &imageSpacing,
                                    const at::Tensor &phiValues, const at::Tensor &rValues, int64_t samplesPerLine);
 
__host__ at::Tensor DRadon2DDR_CUDA(const at::Tensor &image, const at::Tensor &imageSpacing,
                                    const at::Tensor &phiValues, const at::Tensor &rValues, int64_t samplesPerLine);
 
template <typename texture_t> struct Radon2D {
 
    static_assert(texture_t::DIMENSIONALITY == 2);
 
    using IntType = typename texture_t::IntType;
    using FloatType = typename texture_t::FloatType;
    using SizeType = typename texture_t::SizeType;
    using VectorType = typename texture_t::VectorType;
    using AddressModeType = typename texture_t::AddressModeType;
 
    struct CommonData {
        texture_t inputTexture{};
        Linear<VectorType> mappingIndexToOffset{};
        FloatType scaleFactor{};
        at::Tensor flatOutput{};
    };
 
    __host__ static CommonData Common(const at::Tensor &image, const at::Tensor &imageSpacing,
                                      const at::Tensor &phiValues, const at::Tensor &rValues, int64_t samplesPerLine,
                                      at::DeviceType device) {
        // image should be a 2D tensor of floats on the chosen device
        TORCH_CHECK(image.sizes().size() == 2);
        TORCH_CHECK(image.dtype() == at::kFloat);
        TORCH_INTERNAL_ASSERT(image.device().type() == device);
        // imageSpacing should be a 1D tensor of 2 doubles
        TORCH_CHECK(imageSpacing.sizes() == at::IntArrayRef{2});
        TORCH_CHECK(imageSpacing.dtype() == at::kDouble);
        // phiValues and rValues should have the same sizes, should contain floats, and be on the chosen device
        TORCH_CHECK(phiValues.sizes() == rValues.sizes());
        TORCH_CHECK(phiValues.dtype() == at::kFloat);
        TORCH_CHECK(rValues.dtype() == at::kFloat);
        TORCH_INTERNAL_ASSERT(phiValues.device().type() == device);
        TORCH_INTERNAL_ASSERT(rValues.device().type() == device);
 
        CommonData ret{};
        ret.inputTexture = texture_t::FromTensor(image, VectorType::FromTensor(imageSpacing));
        const FloatType lineLength = sqrt(
            ret.inputTexture.SizeWorld().template Apply<FloatType>(&Square<FloatType>).Sum());
        ret.mappingIndexToOffset = GetMappingIToOffset(lineLength, samplesPerLine);
        ret.scaleFactor = lineLength / static_cast<float>(samplesPerLine);
        ret.flatOutput = torch::zeros(at::IntArrayRef({phiValues.numel()}), image.contiguous().options());
        return ret;
    }
 
    __host__ __device__ [[nodiscard]] static Linear<VectorType> GetMappingIToOffset(FloatType lineLength,
        int64_t samplesPerLine) {
        return {VectorType::Full(-.5f * lineLength),
                VectorType::Full(lineLength / static_cast<FloatType>(samplesPerLine - 1))};
    }
 
    __host__ __device__ [[nodiscard]] static Linear<VectorType>
    GetMappingIndexToTexCoord(const texture_t &textureIn, FloatType phi, FloatType r,
                              const Linear<VectorType> &mappingIToOffset) {
        const FloatType s = sin(phi);
        const FloatType c = cos(phi);
        const Linear<VectorType> mappingOffsetToWorld{{r * c, r * s}, {-s, c}};
        return textureIn.MappingWorldToTexCoord()(mappingOffsetToWorld(mappingIToOffset));
    }
 
    __host__ __device__ [[nodiscard]] static VectorType GetDTexCoordDR(const texture_t &textureIn, FloatType phi,
                                                                       FloatType r) {
        const FloatType sign = static_cast<FloatType>(r > 0.) - static_cast<FloatType>(r < 0.);
        const FloatType s = sin(phi);
        const FloatType c = cos(phi);
        return textureIn.MappingWorldToTexCoord().gradient * sign * VectorType{c, s};
    }
 
    __host__ __device__ [[nodiscard]] static float IntegrateLooped(const texture_t &texture,
                                                                   const Linear<VectorType> &mappingIndexToTexCoord,
                                                                   int64_t samplesPerLine) {
        float ret = 0.f;
        for (int64_t i = 0; i < samplesPerLine; ++i) {
            const FloatType iF = static_cast<FloatType>(i);
            ret += texture.Sample(mappingIndexToTexCoord(VectorType::Full(iF)));
        }
        return ret;
    }
 
    __host__ __device__ [[nodiscard]] static float
    DIntegrateLoopedDMappingParameter(const texture_t &texture, const Linear<VectorType> &mappingIndexToTexCoord,
                                      const VectorType &dTexCoordDMappingParameter, int64_t samplesPerLine) {
        float ret = 0.f;
        for (int64_t i = 0; i < samplesPerLine; ++i) {
            const FloatType iF = static_cast<FloatType>(i);
            const VectorType texCoord = mappingIndexToTexCoord(VectorType::Full(iF));
            ret += VecDot(texture.DSampleDTexCoord(texCoord), dTexCoordDMappingParameter);
        }
        return ret;
    }
};
 
} // namespace reg23