PETSIRDInfo_helpers.h

Go to the documentation of this file.

/*
     Copyright 2026, University Medical Center Groningen
     Copyright 2025, MGH / HST A. Martinos Center for Biomedical Imaging
 
    SPDX-License-Identifier: Apache-2.0
    See STIR/LICENSE.txt for details.
 
 */

 
#ifndef __stir_IO_PETSIRDInfo_helpers_H__
#define __stir_IO_PETSIRDInfo_helpers_H__
 
#include "stir/CartesianCoordinate3D.h"
#include <array>
#include <set>
#include "stir/format.h"
#include "stir/info.h"
#include <vector>

 
namespace matrix
{
 
using Mat3 = std::array<std::array<float, 3>, 3>;
using Vec3 = std::array<float, 3>;

inline Mat3
transpose(const Mat3& mat)
{
  std::array<std::array<float, 3>, 3> result{};
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      result[j][i] = mat[i][j];
  return result;
}

inline Mat3
subtract(const Mat3& A, const Mat3& B)
{
  std::array<std::array<float, 3>, 3> result{};
  for (size_t i = 0; i < 3; ++i)
    for (size_t j = 0; j < 3; ++j)
      result[i][j] = A[i][j] - B[i][j];
  return result;
}

inline Vec3
getAxisFromSkew(const Mat3& S)
{
  return {
    0.5f * (S[2][1] - S[1][2]), // x
    0.5f * (S[0][2] - S[2][0]), // y
    0.5f * (S[1][0] - S[0][1])  // z
  };
}
 
} // namespace matrix

namespace vector_utils
{
inline bool
get_spacing_uniform(std::vector<float>& spacing, const std::set<float>& unsorted_block_poss, double epsilon = 1e-4)
{
  std::vector<float> sorted_z(unsorted_block_poss.begin(), unsorted_block_poss.end());
  for (size_t i = 1; i < sorted_z.size(); ++i)
    {
      spacing.push_back(std::abs(sorted_z[i] - sorted_z[i - 1]));
    }
 
  return std::all_of(spacing.begin(), spacing.end(), [&](float s) { return std::abs(s - spacing.front()) <= epsilon; });
}

const std::set<float>&
get_largest_vector(const std::set<float>& x, const std::set<float>& y, const std::set<float>& z)
{
  const std::set<float>* largest = &x;
  if (y.size() > largest->size())
    {
      largest = &y;
    }
  if (z.size() > largest->size())
    {
      largest = &z;
    }
  return *largest;
}

void
find_unique_values_1D(std::set<float>& values, const std::vector<float>& input)
{
  for (float val : input)
    {
      // std::cout << val << std::endl;
      values.insert(val);
    }
}

void
find_unique_values_2D(std::set<float>& values, const std::vector<std::vector<float>>& input)
{
  for (size_t row = 0; row < input.size(); ++row)
    for (size_t col = 0; col < input[row].size(); ++col)
      values.insert(input[row][col]);
}
 
} // namespace vector_utils
 
bool
almostEqual(double a, double b, double tol = 1e-6)
{
  return std::fabs(a - b) <= tol;
}
 
// Detect groupSize along dim2 (y) and loops along dim3 (z)
// Returns true on success, false if pattern doesn't match the assumed structure.
bool
inferGroupSizes_dim2_dim3(const std::vector<stir::CartesianCoordinate3D<float>>& pts,
                          std::size_t& groupSize_dim2,
                          std::size_t& groupSize_dim3,
                          float tol = 1e-5f)
{
  const std::size_t n = pts.size();
  groupSize_dim2 = groupSize_dim3 = 0;
 
  if (n == 0)
    return false;
 
  if (n == 1)
    {
      groupSize_dim2 = 1;
      groupSize_dim3 = 1;
      return true;
    }
 
  // STIR CartesianCoordinate3D<T> is (z, y, x)
  const float x0 = pts[0].x();
  const float z0 = pts[0].z();
 
  // 1) Find how many initial points keep x and z the same
  std::size_t runLen = 1;
  while (runLen < n && almostEqual(pts[runLen].x(), x0, tol) && almostEqual(pts[runLen].z(), z0, tol))
    {
      ++runLen;
    }
 
  groupSize_dim2 = runLen;
 
  // Must tile the full array
  if (groupSize_dim2 == 0 || n % groupSize_dim2 != 0)
    return false;
 
  groupSize_dim3 = n / groupSize_dim2;
 
  // 2) Check each block of size groupSize_dim2 has constant x,z
  for (std::size_t b = 0; b < groupSize_dim3; ++b)
    {
      std::size_t start = b * groupSize_dim2;
      float xb = pts[start].x();
      float zb = pts[start].z();
 
      for (std::size_t i = 1; i < groupSize_dim2; ++i)
        {
          const auto& p = pts[start + i];
          if (!almostEqual(p.x(), xb, tol) || !almostEqual(p.z(), zb, tol))
            {
              return false; // pattern breaks inside a block
            }
        }
    }
 
  // 3) Optionally check that z actually changes between blocks
  for (std::size_t b = 1; b < groupSize_dim3; ++b)
    {
      float z_prev = pts[(b - 1) * groupSize_dim2].z();
      float z_curr = pts[b * groupSize_dim2].z();
      if (almostEqual(z_prev, z_curr, tol))
        {
          return false; // outer loop didn't move in z
        }
    }
 
  return true;
}
 
#endif