ProjectDRR.h

Go to the documentation of this file.

 
#pragma once
 
#include "Common.h"
 
namespace reg23 {
 
at::Tensor ProjectDRR_CPU(const at::Tensor &volume, const at::Tensor &voxelSpacing,
                          const at::Tensor &homographyMatrixInverse, double sourceDistance, int64_t outputWidth,
                          int64_t outputHeight, const at::Tensor &outputOffset, const at::Tensor &detectorSpacing);
 
__host__ at::Tensor ProjectDRR_CUDA(const at::Tensor &volume, const at::Tensor &voxelSpacing,
                                    const at::Tensor &homographyMatrixInverse, double sourceDistance,
                                    int64_t outputWidth, int64_t outputHeight, const at::Tensor &outputOffset,
                                    const at::Tensor &detectorSpacing);
 
at::Tensor ProjectDRR_backward_CPU(const at::Tensor &volume, const at::Tensor &voxelSpacing,
                                   const at::Tensor &homographyMatrixInverse, double sourceDistance,
                                   int64_t outputWidth, int64_t outputHeight, const at::Tensor &outputOffset,
                                   const at::Tensor &detectorSpacing, const at::Tensor &dLossDDRR);
 
__host__ at::Tensor ProjectDRR_backward_CUDA(const at::Tensor &volume, const at::Tensor &voxelSpacing,
                                             const at::Tensor &homographyMatrixInverse, double sourceDistance,
                                             int64_t outputWidth, int64_t outputHeight, const at::Tensor &outputOffset,
                                             const at::Tensor &detectorSpacing, const at::Tensor &dLossDDRR);
 
__host__ at::Tensor ProjectDRRsBatched_CUDA(const at::Tensor &volume, const at::Tensor &voxelSpacing,
                                    const at::Tensor &invHMatrices, double sourceDistance, int64_t outputWidth,
                                    int64_t outputHeight,  const at::Tensor &outputOffset,
                                    const at::Tensor &detectorSpacing);
 
template <typename texture_t> struct ProjectDRR {
 
    static_assert(texture_t::DIMENSIONALITY == 3);
 
    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{};
        Vec<Vec<double, 4>, 4> homographyMatrixInverse{};
        Vec<double, 2> outputOffset{};
        Vec<double, 2> detectorSpacing{};
        double lambdaStart{};
        FloatType stepSize{};
        int64_t samplesPerRay{};
    };
 
    __host__ static CommonData Common(const at::Tensor &volume, const at::Tensor &voxelSpacing,
                                      const at::Tensor &homographyMatrixInverse, double sourceDistance,
                                      const at::Tensor &outputOffset, const at::Tensor &detectorSpacing,
                                      at::DeviceType device, std::optional<int64_t> samplesPerRay = std::nullopt) {
        // volume should be a 3D tensor of floats on the chosen device
        TORCH_CHECK(volume.sizes().size() == 3);
        TORCH_CHECK(volume.dtype() == at::kFloat);
        TORCH_INTERNAL_ASSERT(volume.device().type() == device);
        // voxelSpacing should be a 1D tensor of 3 doubles
        TORCH_CHECK(voxelSpacing.sizes() == at::IntArrayRef{3});
        TORCH_CHECK(voxelSpacing.dtype() == at::kDouble);
        // homographyMatrixInverse should be of size (4, 4), contain doubles and be on the chosen device
        TORCH_CHECK(homographyMatrixInverse.sizes() == at::IntArrayRef({4, 4}));
        TORCH_CHECK(homographyMatrixInverse.dtype() == at::kDouble);
        TORCH_INTERNAL_ASSERT(homographyMatrixInverse.device().type() == device);
        // outputOffset should be a 1D tensor of 2 doubles
        TORCH_CHECK(outputOffset.sizes() == at::IntArrayRef{2});
        TORCH_CHECK(outputOffset.dtype() == at::kDouble);
        // detectorSpacing should be a 1D tensor of 2 doubles
        TORCH_CHECK(detectorSpacing.sizes() == at::IntArrayRef{2});
        TORCH_CHECK(detectorSpacing.dtype() == at::kDouble);
 
        const int64_t samplesPerRayValue = samplesPerRay.value_or(
            Vec<int64_t, 3>::FromIntArrayRef(volume.sizes()).Max());
 
        CommonData ret{};
        ret.inputTexture = texture_t::FromTensor(volume, VectorType::FromTensor(voxelSpacing));
        ret.homographyMatrixInverse = Vec<Vec<double, 4>, 4>::FromTensor2D(homographyMatrixInverse);
 
        const Vec<int64_t, 3> inputSize = Vec<int64_t, 3>::FromIntArrayRef(volume.sizes()).Flipped();
        const VectorType volumeDiagonal = inputSize.StaticCast<FloatType>() * ret.inputTexture.Spacing();
        const FloatType volumeDiagLength = volumeDiagonal.Length();
        const VectorType sourcePosition = {0.0, 0.0, sourceDistance};
        ret.lambdaStart = MatMul(ret.homographyMatrixInverse, VecCat(sourcePosition, 1.0)).XYZ().Length() - 0.5 *
                          volumeDiagLength;
        ret.stepSize = volumeDiagLength / static_cast<FloatType>(samplesPerRayValue);
        ret.samplesPerRay = samplesPerRayValue;
        ret.outputOffset = Vec<double, 2>::FromTensor(outputOffset);
        ret.detectorSpacing = Vec<double, 2>::FromTensor(detectorSpacing);
        return ret;
    }
};
 
} // namespace reg23