Array.h

Go to the documentation of this file.

/*
    Copyright (C) 2000 PARAPET partners
    Copyright (C) 2000 - 2011-10-14, Hammersmith Imanet Ltd
    Copyright (C) 2011-07-01 - 2012, Kris Thielemans
    Copyright (C) 2023 - 2025, University College London
    This file is part of STIR.

    SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license

    See STIR/LICENSE.txt for details
*/

#ifndef __stir_Array_H__
#define __stir_Array_H__

#ifndef ARRAY_FULL
#  define ARRAY_FULL
#endif

#include "stir/NumericVectorWithOffset.h"
#include "stir/ByteOrder.h"
#include "stir/IndexRange.h"
#include "stir/deprecated.h"
#include "stir/shared_ptr.h"
// include forward declaration to ensure consistency, as well as use of
// default parameters (if any)
#include "stir/ArrayFwd.h"

START_NAMESPACE_STIR
class NumericType;

#ifdef ARRAY_FULL
#  ifndef ARRAY_FULL2
template <typename topleveliterT, typename restiterT, typename elemT, typename _Ref, typename _Ptr>
class FullArrayIterator;
#  else
template <int num_dimensions, typename elemT, typename _Ref, typename _Ptr>
class FullArrayIterator;
template <int num_dimensions, typename elemT, typename _Ref, typename _Ptr>
class FullArrayConstIterator;
#  endif

#endif

template <int num_dimensions, typename elemT>
class Array : public NumericVectorWithOffset<Array<num_dimensions - 1, elemT>, elemT>
{
#ifdef STIR_COMPILING_SWIG_WRAPPER
  // work-around swig problem. It gets confused when using a private (or protected)
  // typedef in a definition of a public typedef/member
public:
#endif
  typedef Array<num_dimensions, elemT> self;
  typedef NumericVectorWithOffset<Array<num_dimensions - 1, elemT>, elemT> base_type;

public:
  typedef typename base_type::value_type value_type;
  typedef typename base_type::reference reference;
  typedef typename base_type::const_reference const_reference;
  typedef typename base_type::difference_type difference_type;
  typedef typename base_type::size_type size_type;
  typedef typename base_type::iterator iterator;
  typedef typename base_type::const_iterator const_iterator;
#ifdef ARRAY_FULL

  typedef elemT full_value_type;
  typedef full_value_type* full_pointer;
  typedef const full_value_type* const_full_pointer;
  typedef full_value_type& full_reference;
  typedef const full_value_type& const_full_reference;
#  ifndef ARRAY_FULL2
  typedef FullArrayIterator<typename base_type::iterator,
                            typename Array<num_dimensions - 1, elemT>::full_iterator,
                            elemT,
                            full_reference,
                            full_pointer>
      full_iterator;

  typedef FullArrayIterator<typename base_type::const_iterator,
                            typename Array<num_dimensions - 1, elemT>::const_full_iterator,
                            elemT,
                            const_full_reference,
                            const_full_pointer>
      const_full_iterator;
#  else // ARRAY_FULL2
  typedef FullArrayIterator<num_dimensions, elemT, full_reference, full_pointer> full_iterator;

  typedef FullArrayConstIterator<num_dimensions, elemT, const_full_reference, const_full_pointer> const_full_iterator;

#  endif

#endif
public:
  inline Array();

  inline explicit Array(const IndexRange<num_dimensions>&);


  inline Array(const IndexRange<num_dimensions>& range, shared_ptr<elemT[]> data_sptr);

#ifndef SWIG
  // swig 2.0.4 gets confused by base_type (due to numeric template arguments)
  // therefore, we declare this constructor using the "self" type,
  // i.e. it's just a copy-constructor.
  // This is less powerful as in C++, but swig-generated interfaces don't need to know about the base_type anyway
  inline Array(const base_type& t);
#endif

  // implementation needed as the above doesn't disable the auto-generated copy-constructor
  inline Array(const self& t);

  inline ~Array() override;

  // implementation in .h because of templates/friends/whatever, see https://stackoverflow.com/a/61020224
  friend inline void swap(Array& first, Array& second) // nothrow
  {
    using std::swap;
    //  swap the members of two objects
    swap(static_cast<base_type&>(first), static_cast<base_type&>(second));
    swap(first._allocated_full_data_ptr, second._allocated_full_data_ptr);
  }


  Array(Array&& other) noexcept;


  Array& operator=(Array other);

  inline full_iterator begin_all();
  inline const_full_iterator begin_all() const;
  inline const_full_iterator begin_all_const() const;
  inline full_iterator end_all();
  inline const_full_iterator end_all() const;
  inline const_full_iterator end_all_const() const;

  inline IndexRange<num_dimensions> get_index_range() const;

  inline size_t size_all() const;


  inline virtual void resize(const IndexRange<num_dimensions>& range);

  virtual inline void grow(const IndexRange<num_dimensions>& range);

  inline elemT sum() const;

  inline elemT sum_positive() const;

  inline elemT find_max() const;

  inline elemT find_min() const;


  inline void fill(const elemT& n);
  inline void apply_lower_threshold(const elemT& l);

  inline void apply_upper_threshold(const elemT& u);


  inline bool is_regular() const;


  bool get_regular_range(BasicCoordinate<num_dimensions, int>& min, BasicCoordinate<num_dimensions, int>& max) const;

  inline Array<num_dimensions - 1, elemT>& operator[](int i);

  inline const Array<num_dimensions - 1, elemT>& operator[](int i) const;

  inline elemT& operator[](const BasicCoordinate<num_dimensions, int>& c);

  // TODO alternative return value: elemT
  inline const elemT& operator[](const BasicCoordinate<num_dimensions, int>& c) const;


  inline Array<num_dimensions - 1, elemT>& at(int i);

  inline const Array<num_dimensions - 1, elemT>& at(int i) const;

  inline elemT& at(const BasicCoordinate<num_dimensions, int>& c);

  inline const elemT& at(const BasicCoordinate<num_dimensions, int>& c) const;


  // tedious reimplementation to fix return types. This could be avoided by using boost::operators.
  // However, reimplementing them explicitly helps SWIG.
  inline self& operator+=(const self& x)
  {
    base_type::operator+=(x);
    return *this;
  }
  inline self& operator-=(const self& x)
  {
    base_type::operator-=(x);
    return *this;
  }
  inline self& operator*=(const self& x)
  {
    base_type::operator*=(x);
    return *this;
  }
  inline self& operator/=(const self& x)
  {
    base_type::operator/=(x);
    return *this;
  }
  inline self& operator+=(const elemT x)
  {
    base_type::operator+=(x);
    return *this;
  }
  inline self& operator-=(const elemT x)
  {
    base_type::operator-=(x);
    return *this;
  }
  inline self& operator*=(const elemT x)
  {
    base_type::operator*=(x);
    return *this;
  }
  inline self& operator/=(const elemT x)
  {
    base_type::operator/=(x);
    return *this;
  }
  inline self operator+(const self& x) const
  {
    self c(*this);
    return c += x;
  }
  inline self operator+(const elemT x) const
  {
    self c(*this);
    return c += x;
  }
  inline self operator-(const self& x) const
  {
    self c(*this);
    return c -= x;
  }
  inline self operator-(const elemT x) const
  {
    self c(*this);
    return c -= x;
  }
  inline self operator*(const self& x) const
  {
    self c(*this);
    return c *= x;
  }
  inline self operator*(const elemT x) const
  {
    self c(*this);
    return c *= x;
  }
  inline self operator/(const self& x) const
  {
    self c(*this);
    return c /= x;
  }
  inline self operator/(const elemT x) const
  {
    self c(*this);
    return c /= x;
  }

  template <typename elemT2>
  STIR_DEPRECATED inline void axpby(const elemT2 a, const Array& x, const elemT2 b, const Array& y);

  inline void xapyb(const Array& x, const elemT a, const Array& y, const elemT b);

  inline void xapyb(const Array& x, const Array& a, const Array& y, const Array& b);

  template <class T>
  inline void sapyb(const T& a, const Array& y, const T& b);


  bool is_contiguous() const;

  inline elemT* get_full_data_ptr();

  inline const elemT* get_const_full_data_ptr() const;

  inline void release_full_data_ptr();

  inline void release_const_full_data_ptr() const;

private:
  // This variable is declared mutable such that get_const_full_data_ptr() can change it.
  mutable bool _full_pointer_access;

  shared_ptr<elemT[]> _allocated_full_data_ptr;


  inline void init_with_copy(const IndexRange<num_dimensions>& range, elemT const* const data_ptr);

  inline void init(const IndexRange<num_dimensions>& range, elemT* const data_ptr, bool copy_data);
  // Make sure that we can access init() recursively
  template <int num_dimensions2, class elemT2>
  friend class Array;

  using base_type::grow;
  using base_type::resize;
};

/**************************************************
 (partial) specialisation for 1 dimensional arrays
 **************************************************/

template <class elemT>
class Array<1, elemT> : public NumericVectorWithOffset<elemT, elemT>
{
#ifdef STIR_COMPILING_SWIG_WRAPPER
  // work-around swig problem. It gets confused when using a private (or protected)
  // typedef in a definition of a public typedef/member
public:
#endif
  typedef NumericVectorWithOffset<elemT, elemT> base_type;
  typedef Array<1, elemT> self;

public:

  typedef typename base_type::value_type value_type;
  typedef typename base_type::reference reference;
  typedef typename base_type::const_reference const_reference;
  typedef typename base_type::difference_type difference_type;
  typedef typename base_type::size_type size_type;
  typedef typename base_type::iterator iterator;
  typedef typename base_type::const_iterator const_iterator;

  typedef iterator full_iterator;

  typedef const_iterator const_full_iterator;

public:
  inline Array();

  inline explicit Array(const IndexRange<1>& range);

  inline Array(const int min_index, const int max_index);


  inline Array(const IndexRange<1>& range, shared_ptr<elemT[]> data_sptr);


  inline Array(const IndexRange<1>& range, const elemT* const data_ptr);

  inline Array(const NumericVectorWithOffset<elemT, elemT>& il);

  // implementation needed as the above doesn't replace the normal copy-constructor
  // and the auto-generated is disabled because of the move constructor
  inline Array(const self& t);


  Array(Array&& other) noexcept;

  inline ~Array() override;

  // implementation in .h because of templates/friends/whatever, see https://stackoverflow.com/a/61020224
  friend inline void swap(Array& first, Array& second) // nothrow
  {
    swap(static_cast<base_type&>(first), static_cast<base_type&>(second));
  }

  inline Array& operator=(const Array& other);

  inline full_iterator begin_all();
  inline const_full_iterator begin_all() const;
  inline const_full_iterator begin_all_const() const;
  inline full_iterator end_all();
  inline const_full_iterator end_all() const;
  inline const_full_iterator end_all_const() const;

  inline IndexRange<1> get_index_range() const;

  inline size_t size_all() const;

  inline virtual void grow(const IndexRange<1>& range);

  // Array::grow initialises new elements to 0
  inline void grow(const int min_index, const int max_index) override;

  inline virtual void resize(const IndexRange<1>& range);

  inline void resize(const int min_index, const int max_index, bool initialise_with_0);
  inline void resize(const int min_index, const int max_index) override;


  bool is_contiguous() const
  {
    return true;
  }
  inline elemT* get_full_data_ptr()
  {
    return this->get_data_ptr();
  }

  inline const elemT* get_const_full_data_ptr() const
  {
    return this->get_const_data_ptr();
  }

  inline void release_full_data_ptr()
  {
    this->release_data_ptr();
  }

  inline void release_const_full_data_ptr() const
  {
    this->release_const_data_ptr();
  }

  inline elemT sum() const;

  inline elemT sum_positive() const;

  inline elemT find_max() const;

  inline elemT find_min() const;

  inline bool is_regular() const;

  bool get_regular_range(BasicCoordinate<1, int>& min, BasicCoordinate<1, int>& max) const;

  /* Add numerical operators with correct return value, as opposed to those from the base class
   */
  inline self operator+(const base_type& iv) const;

  inline self operator-(const base_type& iv) const;

  inline self operator*(const base_type& iv) const;

  inline self operator/(const base_type& iv) const;

  inline self operator+(const elemT a) const;

  inline self operator-(const elemT a) const;

  inline self operator*(const elemT a) const;

  inline self operator/(const elemT a) const;

  inline elemT& operator[](int i);

  inline const elemT& operator[](int i) const;

  inline const elemT& operator[](const BasicCoordinate<1, int>& c) const;

  inline elemT& operator[](const BasicCoordinate<1, int>& c);


  inline elemT& at(int i);

  inline const elemT& at(int i) const;

  inline elemT& at(const BasicCoordinate<1, int>& c);

  inline const elemT& at(const BasicCoordinate<1, int>& c) const;
private:
  // Make sure we can call init() recursively.
  template <int num_dimensions2, class elemT2>
  friend class Array;


  inline void init_with_copy(const IndexRange<1>& range, elemT const* const data_ptr);

  inline void init(const IndexRange<1>& range, elemT* const data_ptr, bool copy_data);
};

END_NAMESPACE_STIR

#ifdef ARRAY_FULL
#  ifndef ARRAY_FULL2
#    include "FullArrayIterator.h"
#  else
#    include "FullArrayIterator2.h"
#    include "FullArrayConstIterator.h"
#  endif
#endif

#include "stir/Array.inl"

#endif // __Array_H__