DiscretisedDensity.inl

Go to the documentation of this file.

//
//
/*
    Copyright (C) 2000 PARAPET partners
    Copyright (C) 2000- 2009-07-08, Hammersmith Imanet Ltd
    Copyright (C) 2011-07-01 - 2011, Kris Thielemans
    Copyright (C) 2018, 2023, University College London
    Copyright (C) 2018, CSIRO
    This file is part of STIR.
 
    SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license
 
    See STIR/LICENSE.txt for details
*/
#include "stir/round.h"
#include <string>
#include <typeinfo>
#include <stdio.h> // for snprintf
 
START_NAMESPACE_STIR
 
template <int num_dimensions, typename elemT>
DiscretisedDensity<num_dimensions, elemT>::DiscretisedDensity()
{}
 
template <int num_dimensions, typename elemT>
DiscretisedDensity<num_dimensions, elemT>::DiscretisedDensity(const IndexRange<num_dimensions>& range_v,
                                                              const CartesianCoordinate3D<float>& origin_v)
    : base_type(range_v),
      origin(origin_v)
{}
 
template <int num_dimensions, typename elemT>
DiscretisedDensity<num_dimensions, elemT>::DiscretisedDensity(const shared_ptr<const ExamInfo>& exam_info_sptr,
                                                              const IndexRange<num_dimensions>& range_v,
                                                              const CartesianCoordinate3D<float>& origin_v)
    : ExamData(exam_info_sptr),
      base_type(range_v),
      origin(origin_v)
{}
 
template <int num_dimensions, typename elemT>
void
DiscretisedDensity<num_dimensions, elemT>::set_origin(const CartesianCoordinate3D<float>& origin_v)
{
  origin = origin_v;
}
 
template <int num_dimensions, typename elemT>
const CartesianCoordinate3D<float>&
DiscretisedDensity<num_dimensions, elemT>::get_origin() const
{
  return origin;
}
 
template <int num_dimensions, typename elemT>
bool
DiscretisedDensity<num_dimensions, elemT>::actual_has_same_characteristics(DiscretisedDensity<num_dimensions, elemT> const& other,
                                                                           std::string& explanation) const
{
  if (typeid(other) != typeid(*this))
    {
      explanation = "Different type of data";
      return false;
    }
 
  if (norm(other.get_origin() - this->get_origin()) > 1.E-2)
    {
      char tmp[2000];
      snprintf(tmp,
               2000,
               "Not the same origin: (%g,%g,%g) and (%g,%g,%g)",
               other.get_origin()[1],
               other.get_origin()[2],
               other.get_origin()[3],
               this->get_origin()[1],
               this->get_origin()[2],
               this->get_origin()[3]);
      explanation = tmp;
      return false;
    }
  if (other.get_index_range() != this->get_index_range())
    {
      explanation = "Not the same index ranges.";
      return false;
    }
  return true;
}
 
template <int num_dimensions, typename elemT>
bool
DiscretisedDensity<num_dimensions, elemT>::has_same_characteristics(DiscretisedDensity<num_dimensions, elemT> const& other,
                                                                    std::string& explanation) const
{
  return this->actual_has_same_characteristics(other, explanation);
}
 
template <int num_dimensions, typename elemT>
bool
DiscretisedDensity<num_dimensions, elemT>::has_same_characteristics(DiscretisedDensity<num_dimensions, elemT> const& other) const
{
  std::string explanation;
  return this->actual_has_same_characteristics(other, explanation);
}
 
template <int num_dimensions, typename elemT>
bool
DiscretisedDensity<num_dimensions, elemT>::operator==(const self_type& that) const
{
  return this->has_same_characteristics(that) && base_type::operator==(that);
}
 
template <int num_dimensions, typename elemT>
bool
DiscretisedDensity<num_dimensions, elemT>::operator!=(const self_type& that) const
{
  return !((*this) == that);
}
 
template <int num_dimensions, typename elemT>
CartesianCoordinate3D<float>
DiscretisedDensity<num_dimensions, elemT>::get_relative_coordinates_for_indices(
    const BasicCoordinate<num_dimensions, int>& indices) const
{
  return this->actual_get_relative_coordinates_for_indices(BasicCoordinate<num_dimensions, float>(indices));
}
 
template <int num_dimensions, typename elemT>
CartesianCoordinate3D<float>
DiscretisedDensity<num_dimensions, elemT>::get_relative_coordinates_for_indices(
    const BasicCoordinate<num_dimensions, float>& indices) const
{
  return this->actual_get_relative_coordinates_for_indices(indices);
}
 
template <int num_dimensions, typename elemT>
BasicCoordinate<num_dimensions, float>
DiscretisedDensity<num_dimensions, elemT>::get_index_coordinates_for_physical_coordinates(
    const CartesianCoordinate3D<float>& coords) const
{
  return this->actual_get_index_coordinates_for_relative_coordinates(coords - this->get_origin());
}
 
template <int num_dimensions, typename elemT>
BasicCoordinate<num_dimensions, float>
DiscretisedDensity<num_dimensions, elemT>::get_index_coordinates_for_relative_coordinates(
    const CartesianCoordinate3D<float>& coords) const
{
  return this->actual_get_index_coordinates_for_relative_coordinates(coords);
}
 
template <int num_dimensions, typename elemT>
BasicCoordinate<num_dimensions, int>
DiscretisedDensity<num_dimensions, elemT>::get_indices_closest_to_relative_coordinates(
    const CartesianCoordinate3D<float>& coords) const
{
  return round(this->actual_get_index_coordinates_for_relative_coordinates(coords));
}
 
template <int num_dimensions, typename elemT>
CartesianCoordinate3D<float>
DiscretisedDensity<num_dimensions, elemT>::get_physical_coordinates_for_indices(
    const BasicCoordinate<num_dimensions, int>& indices) const
{
  return this->get_relative_coordinates_for_indices(indices) + this->get_origin();
}
 
template <int num_dimensions, typename elemT>
CartesianCoordinate3D<float>
DiscretisedDensity<num_dimensions, elemT>::get_physical_coordinates_for_indices(
    const BasicCoordinate<num_dimensions, float>& indices) const
{
  return this->get_relative_coordinates_for_indices(indices) + this->get_origin();
}
 
template <int num_dimensions, typename elemT>
BasicCoordinate<num_dimensions, int>
DiscretisedDensity<num_dimensions, elemT>::get_indices_closest_to_physical_coordinates(
    const CartesianCoordinate3D<float>& coords) const
{
  return this->get_indices_closest_to_relative_coordinates(coords - this->get_origin());
}
 
template <int num_dimensions, typename elemT>
CartesianCoordinate3D<float>
DiscretisedDensity<num_dimensions, elemT>::swap_axes_based_on_orientation(const CartesianCoordinate3D<float>& coords,
                                                                          const PatientPosition patient_position)
{
  CartesianCoordinate3D<float> flip_coords = coords;
  // STIR coordinates run:
  // x: left to right as you face the scanner (as seen from the bed side)
  // y: top to bottom of the scanner
  // z: running from deep in the scanner out along the bed (as seen from the bed
  //    side)
  // ITK coordinates are defined w.r.t LPS
  switch (patient_position.get_position())
    {
    case PatientPosition::unknown_position: // If unknown, assume HFS
    case PatientPosition::HFS:              // HFS means currently in patient LPI
      flip_coords.z() *= -1;
      break;
 
    case PatientPosition::HFP: // HFP means currently in patient RAI
      flip_coords.x() *= -1;
      flip_coords.y() *= -1;
      flip_coords.z() *= -1;
      break;
 
    case PatientPosition::FFS: // FFS means currently in patient RPS
      flip_coords.x() *= -1;
      break;
 
    case PatientPosition::FFP: // FFP means currently in patient LAS
      flip_coords.y() *= -1;
      break;
 
    default:
      throw std::runtime_error("Unsupported patient position, can't convert to LPS.");
    }
  return flip_coords;
}
 
template <int num_dimensions, typename elemT>
CartesianCoordinate3D<float>
DiscretisedDensity<num_dimensions, elemT>::get_LPS_coordinates_for_physical_coordinates(
    const CartesianCoordinate3D<float>& coords) const
{
  return swap_axes_based_on_orientation(coords, this->get_exam_info().patient_position);
}
 
template <int num_dimensions, typename elemT>
CartesianCoordinate3D<float>
DiscretisedDensity<num_dimensions, elemT>::get_LPS_coordinates_for_indices(
    const BasicCoordinate<num_dimensions, float>& indices) const
{
  return this->get_LPS_coordinates_for_physical_coordinates(this->get_physical_coordinates_for_indices(indices));
}
 
template <int num_dimensions, typename elemT>
CartesianCoordinate3D<float>
DiscretisedDensity<num_dimensions, elemT>::get_LPS_coordinates_for_indices(
    const BasicCoordinate<num_dimensions, int>& indices) const
{
  return get_LPS_coordinates_for_indices(BasicCoordinate<num_dimensions, float>(indices));
}
 
template <int num_dimensions, typename elemT>
CartesianCoordinate3D<float>
DiscretisedDensity<num_dimensions, elemT>::get_physical_coordinates_for_LPS_coordinates(
    const CartesianCoordinate3D<float>& coords) const
{
  // operation is symmetric
  return this->get_LPS_coordinates_for_physical_coordinates(coords);
}
 
template <int num_dimensions, typename elemT>
BasicCoordinate<num_dimensions, float>
DiscretisedDensity<num_dimensions, elemT>::get_index_coordinates_for_LPS_coordinates(
    const CartesianCoordinate3D<float>& coords) const
{
  return this->get_index_coordinates_for_physical_coordinates(this->get_physical_coordinates_for_LPS_coordinates(coords));
}
 
template <int num_dimensions, typename elemT>
BasicCoordinate<num_dimensions, int>
DiscretisedDensity<num_dimensions, elemT>::get_indices_closest_to_LPS_coordinates(
    const CartesianCoordinate3D<float>& coords) const
{
  return round(this->get_index_coordinates_for_LPS_coordinates(coords));
}
 
END_NAMESPACE_STIR