StructureTensor.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <cstdint>
 
#include <opencv2/core.hpp>
 
#include "vc/core/math/Tensor3D.hpp"
#include "vc/core/types/Volume.hpp"
 
namespace volcart
{
using EigenValue = double;
 
using EigenVector = cv::Vec3d;
 
using EigenPairs = std::array<std::pair<EigenValue, EigenVector>, 3>;
 
using StructureTensor = cv::Matx33d;
 
static const auto ZERO_STRUCTURE_TENSOR =
    StructureTensor(0, 0, 0, 0, 0, 0, 0, 0, 0);
 
StructureTensor ComputeVoxelStructureTensor(
    const Volume::Pointer& volume,
    int vx,
    int vy,
    int vz,
    int radius = 1,
    int kernelSize = 3);
 
StructureTensor ComputeVoxelStructureTensor(
    const Volume::Pointer& volume,
    const cv::Vec3i& index,
    int radius = 1,
    int kernelSize = 3);
 
StructureTensor ComputeSubvoxelStructureTensor(
    const Volume::Pointer& volume,
    double vx,
    double vy,
    double vz,
    int radius = 1,
    int kernelSize = 3);
 
StructureTensor ComputeSubvoxelStructureTensor(
    const Volume::Pointer& volume,
    const cv::Vec3d& index,
    int radius = 1,
    int kernelSize = 3);
 
EigenPairs ComputeVoxelEigenPairs(
    const Volume::Pointer& volume,
    int x,
    int y,
    int z,
    int radius = 1,
    int kernelSize = 3);
 
EigenPairs ComputeVoxelEigenPairs(
    const Volume::Pointer& volume,
    const cv::Vec3i& index,
    int radius = 1,
    int kernelSize = 3);
 
EigenPairs ComputeSubvoxelEigenPairs(
    const Volume::Pointer& volume,
    double x,
    double y,
    double z,
    int radius = 1,
    int kernelSize = 3);
 
EigenPairs ComputeSubvoxelEigenPairs(
    const Volume::Pointer& volume,
    const cv::Vec3d& index,
    int radius = 1,
    int kernelSize = 3);
 
template <typename DType>
Tensor3D<DType> ComputeVoxelNeighbors(
    const Volume::Pointer& volume,
    const cv::Vec3i& center,
    std::int32_t rx,
    std::int32_t ry,
    std::int32_t rz)
{
    Tensor3D<DType> v(2 * rx + 1, 2 * ry + 1, 2 * rz + 1);
    for (std::int32_t k = center(2) - rz, c = 0; k <= center(2) + rz;
         ++k, ++c) {
        for (std::int32_t j = center(1) - ry, b = 0; j <= center(1) + ry;
             ++j, ++b) {
            for (std::int32_t i = center(0) - rx, a = 0; i <= center(0) + rx;
                 ++i, ++a) {
                v(a, b, c) = DType(volume->intensityAt(i, j, k));
            }
        }
    }
 
    return v;
}
 
template <typename DType>
Tensor3D<DType> ComputeSubvoxelNeighbors(
    const Volume::Pointer& volume,
    const cv::Vec3d& center,
    int rx,
    int ry,
    int rz,
    const cv::Vec3d& xvec = {1, 0, 0},
    const cv::Vec3d& yvec = {0, 1, 0},
    const cv::Vec3d& zvec = {0, 0, 1})
{
    Tensor3D<DType> v(2 * rx + 1, 2 * ry + 1, 2 * rz + 1);
    for (int c = 0; c < 2 * rz + 1; ++c) {
        for (int b = 0; b < 2 * ry + 1; ++b) {
            for (int a = 0; a < 2 * rx + 1; ++a) {
                auto xOffset = -rx + a;
                auto yOffset = -ry + b;
                auto zOffset = -rz + c;
                auto p = center + (xvec * xOffset) + (yvec * yOffset) +
                         (zvec * zOffset);
                v(a, b, c) = DType(volume->interpolateAt(p));
            }
        }
    }
 
    return v;
}
}  // namespace volcart