Support multiple polynomial packing orders.

include/lsst/geom/polynomials.h

@@ -90,8 +90,8 @@ namespace lsst { namespace geom {
  *   - BinomialMatrix provides an efficient and stable way to get
  *     binomial coefficients.
  *
- *   - PackedIndexIterator and PackedIndexRange define a packing order that
- *     allows a conceptual 2-d triangular matrix to be implemented on top of
+ *   - PackedIndexIterator and PackedIndexRange implement PackingOrders that
+ *     allow a conceptual 2-d triangular matrix to be implemented on top of
  *     a 1-d array.
  */
 namespace polynomials {

include/lsst/geom/polynomials/Chebyshev1Basis2d.h

@@ -29,11 +29,25 @@
 
 namespace lsst { namespace geom { namespace polynomials {
 
-/// A Basis2d for Chebyshev polynomials of the first kind.
-using Chebyshev1Basis2d = PackedBasis2d<Chebyshev1Basis1d>;
+/// A Basis2d for Chebyshev polynomials of the first kind, templated on packing order.
+template <PackingOrder packing>
+using Chebyshev1Basis2d = PackedBasis2d<Chebyshev1Basis1d, packing>;
 
-/// A Basis2d for scaled Chebyshev polynomials of the first kind.
-using ScaledChebyshev1Basis2d = ScaledBasis2d<Chebyshev1Basis2d>;
+/// A Basis2d for scaled Chebyshev polynomials of the first kind, templated on packing order.
+template <PackingOrder packing>
+using ScaledChebyshev1Basis2d = ScaledBasis2d<Chebyshev1Basis2d<packing>>;
+
+/// A Basis2d for Chebyshev polynomials of the first kind, ordered via PackingOrder::XY.
+using Chebyshev1Basis2dXY = Chebyshev1Basis2d<PackingOrder::XY>;
+
+/// A Basis2d for Chebyshev polynomials of the first kind, ordered via PackingOrder::YX.
+using Chebyshev1Basis2dYX = Chebyshev1Basis2d<PackingOrder::YX>;
+
+/// A Basis2d for scaled Chebyshev polynomials of the first kind, ordered via PackingOrder::XY.
+using ScaledChebyshev1Basis2dXY = ScaledChebyshev1Basis2d<PackingOrder::XY>;
+
+/// A Basis2d for scaled Chebyshev polynomials of the first kind, ordered via PackingOrder::YX.
+using ScaledChebyshev1Basis2dYX = ScaledChebyshev1Basis2d<PackingOrder::YX>;
 
 }}} // namespace lsst::geom::polynomials
 

include/lsst/geom/polynomials/Chebyshev1Function2d.h

@@ -27,11 +27,25 @@
 
 namespace lsst { namespace geom { namespace polynomials {
 
-/// A Function2d for Chebyshev polynomials of the first kind.
-using Chebyshev1Function2d = Function2d<Chebyshev1Basis2d>;
+/// A Function2d for Chebyshev polynomials of the first kind, templated on packing order.
+template <PackingOrder packing>
+using Chebyshev1Function2d = Function2d<Chebyshev1Basis2d<packing>>;
 
-/// A Function2d for scaled Chebyshev polynomials of the first kind.
-using ScaledChebyshev1Function2d = Function2d<ScaledChebyshev1Basis2d>;
+/// A Function2d for scaled Chebyshev polynomials of the first kind, templated on packing order.
+template <PackingOrder packing>
+using ScaledChebyshev1Function2d = Function2d<ScaledChebyshev1Basis2d<packing>>;
+
+/// A Function2d for Chebyshev polynomials of the first kind, ordered via PackingOrder::XY.
+using Chebyshev1Function2dXY = Chebyshev1Function2d<PackingOrder::XY>;
+
+/// A Function2d for Chebyshev polynomials of the first kind, ordered via PackingOrder::YX.
+using Chebyshev1Function2dYX = Chebyshev1Function2d<PackingOrder::YX>;
+
+/// A Function2d for scaled Chebyshev polynomials of the first kind, ordered via PackingOrder::XY.
+using ScaledChebyshev1Function2dXY = ScaledChebyshev1Function2d<PackingOrder::XY>;
+
+/// A Function2d for scaled Chebyshev polynomials of the first kind, ordered via PackingOrder::YX.
+using ScaledChebyshev1Function2dYX = ScaledChebyshev1Function2d<PackingOrder::YX>;
 
 }}} // namespace lsst::geom::polynomials
 

include/lsst/geom/polynomials/PackedBasis2d.h

@@ -28,7 +28,7 @@
 
 namespace lsst { namespace geom { namespace polynomials {
 
-template <typename Basis1d>
+template <typename Basis1d, PackingOrder packing>
 class PackedBasis2d;
 
 
@@ -47,7 +47,7 @@ class PackedBasisWorkspace2d {
 
 private:
 
-    template <typename Recurrence>
+    template <typename Recurrence, PackingOrder packing>
     friend class PackedBasis2d;
 
     Eigen::VectorXd _x;
@@ -70,7 +70,7 @@ class Function2d;
  *  PackedBasis2d::IndexRange. Note that while all PackedBasis2d subclasses
  *  use this ordering, Basis2d objects in general are not required to.
  */
-template <typename Basis1d>
+template <typename Basis1d, PackingOrder packing>
 class PackedBasis2d {
 public:
 
@@ -84,7 +84,7 @@ class PackedBasis2d {
     using Workspace = PackedBasisWorkspace2d;
 
     /// The type returned by getIndices().
-    using IndexRange = PackedIndexRange;
+    using IndexRange = PackedIndexRange<packing>;
 
     /// Return the size of a PackedBasis with the given order.
     static constexpr std::size_t computeSize(std::size_t order) { return IndexRange::computeSize(order); }
@@ -152,7 +152,7 @@ class PackedBasis2d {
      *  See PackedIndexIterator for documentation of the actual order.
      */
     IndexRange getIndices() const noexcept {
-        return IndexRange(IndexRange::iterator(), IndexRange::iterator::makeEnd(getOrder()));
+        return IndexRange(typename IndexRange::iterator(), IndexRange::iterator::makeEnd(getOrder()));
     }
 
     /// Allocate a workspace that can be passed to sumWith() and fill() to avoid repeated memory allocations.

include/lsst/geom/polynomials/PackedIndex.h

@@ -24,6 +24,39 @@
 
 namespace lsst { namespace geom { namespace polynomials {
 
+/// Enum defining the packing orders used to order 2-d polynomial coefficients.
+enum class PackingOrder {
+
+    /**
+     *  A pair of indices @f$(nx, ny)@f$ is mapped to the flattened position
+     *  @f$i = (nx+ny)(nx+ny+1)/2 + nx@f$, which yields the (nx, ny) ordering
+     *  ```
+     *      (0, 0), (0, 1), (1, 0), (0, 2), (1, 1), (2, 0), ...
+     *  ```
+     *
+     *  In other words, if @f$a_i@f$ are the coefficients of a 2-d polynomial
+     *  expansion with order=2, the full polynomial is:
+     *  @f[
+     *     a_0 + a_1 y + a_2 x + a_3 y^2 + a_4 x y + a_5 x^2
+     *  @f]
+     */
+    YX,
+
+    /**
+     *  A pair of indices @f$(nx, ny)@f$ is mapped to the flattened position
+     *  @f$i = (nx+ny)(nx+ny+1)/2 + ny@f$, which yields the (nx, ny) ordering
+     *  ```
+     *      (0, 0), (1, 0), (0, 1), (2, 0), (1, 1), (0, 2), ...
+     *  ```
+     *
+     *  In other words, if @f$a_i@f$ are the coefficients of a 2-d polynomial
+     *  expansion with order=2, the full polynomial is:
+     *  @f[
+     *     a_0 + a_1 x + a_2 y + a_3 x^2 + a_4 x y + a_5 y^2
+     *  @f]
+     */
+    XY
+};
 
 /**
  *  A custom tuple that relates the indices of two 1-d functions for x and y
@@ -58,6 +91,56 @@ struct Index2d {
     std::size_t ny;    ///< Index into the 1-d functoin for ny.
 };
 
+namespace detail {
+
+// Specialization of all PackedIndexIterator/PackedIndexRange logic
+// that isn't common across PackingOrders.
+template <PackingOrder packing> struct PackingOrderTraits;
+
+template <>
+struct PackingOrderTraits<PackingOrder::YX> {
+
+    static std::size_t computeInnerIndex(std::size_t nx, std::size_t ny) { return nx; }
+
+    static void increment(Index2d & index) {
+        if (index.ny == 0) {
+            index.ny = index.nx + 1;
+            index.nx = 0;
+        } else {
+            --index.ny;
+            ++index.nx;
+        }
+    }
+
+    static std::size_t getEndX(std::size_t order) { return 0; }
+
+    static std::size_t getEndY(std::size_t order) { return order + 1; }
+
+};
+
+template <>
+struct PackingOrderTraits<PackingOrder::XY> {
+
+    static std::size_t computeInnerIndex(std::size_t nx, std::size_t ny) { return ny; }
+
+    static void increment(Index2d & index) {
+        if (index.nx == 0) {
+            index.nx = index.ny + 1;
+            index.ny = 0;
+        } else {
+            --index.nx;
+            ++index.ny;
+        }
+    }
+
+    static std::size_t getEndX(std::size_t order) { return order + 1; }
+
+    static std::size_t getEndY(std::size_t order) { return 0; }
+
+};
+
+} // namespace detail
+
 /**
  *  An iterator for traversing "packed" triangular 2-d series expansions,
  *  in which two 1-d expansions are truncated according to the sum of their
@@ -67,18 +150,17 @@ struct Index2d {
  *  PackedIndexIterator dereferences to Index2d.  Typical usage is via a
  *  PackedIndexRange.
  *
- *  A pair of indices @f$(p, q)@f$ is mapped to the flattened position
- *  @f$i = (nx+ny)(nx+ny+1)/2 + nx@f$, which yields the (nx, ny) ordering
- *  ```
- *      (0, 0), (0, 1), (1, 0), (1, 1), (0, 2), (1, 1), (2, 0), ...
- *  ```
  *  This packing ensures that the coefficients for an nth-order expansion are
  *  a contiguous subset of the coefficients for an (n+1)th-order expansion.
  *
+ *  The packing within each order is set by the `packing` template parameter.
+ *
  *  PackedIndexIterator is an STL input iterator, *not* a forward iterator;
  *  incrementing the iterator invalidates all previously-dereferenced values.
  */
+template <PackingOrder packing>
 class PackedIndexIterator {
+    using Traits = detail::PackingOrderTraits<packing>;
 public:
 
     using difference_type = std::ptrdiff_t;
@@ -99,7 +181,7 @@ class PackedIndexIterator {
 
     /// Return the flattened index for the element with the given x and y orders.
     static constexpr std::size_t computeIndex(std::size_t nx, std::size_t ny) noexcept {
-        return computeOffset(nx + ny) + nx;
+        return computeOffset(nx + ny) + Traits::computeInnerIndex(nx, ny);
     }
 
     /// Construct an iterator one past the end of an expansion with the given order.
@@ -124,13 +206,7 @@ class PackedIndexIterator {
     /// Move to the next element in the packed array and return the iterator.
     PackedIndexIterator & operator++() noexcept {
         ++_index.flat;
-        if (_index.ny == 0) {
-            _index.ny = _index.nx + 1;
-            _index.nx = 0;
-        } else {
-            --_index.ny;
-            ++_index.nx;
-        }
+        Traits::increment(_index);
         return *this;
     }
 
@@ -154,7 +230,7 @@ class PackedIndexIterator {
 private:
 
     constexpr PackedIndexIterator(std::size_t order) noexcept :
-        _index(computeOffset(order + 1), 0, order + 1)
+        _index(computeOffset(order + 1), Traits::getEndX(order), Traits::getEndY(order))
     {}
 
     Index2d _index;
@@ -166,15 +242,16 @@ class PackedIndexIterator {
  *
  *  @see PackedIndexIterator for information on the packing algorithm.
  */
+template <PackingOrder packing>
 class PackedIndexRange {
 public:
 
-    using value_type = PackedIndexIterator::value_type;
-    using reference = PackedIndexIterator::reference;
-    using pointer = PackedIndexIterator::pointer;
-    using iterator = PackedIndexIterator;
-    using const_iterator = PackedIndexIterator;
-    using difference_type = PackedIndexIterator::difference_type;
+    using iterator = PackedIndexIterator<packing>;
+    using const_iterator = iterator;
+    using value_type = typename iterator::value_type;
+    using reference = typename iterator::reference;
+    using pointer = typename iterator::pointer;
+    using difference_type = typename iterator::difference_type;
     using size_type = std::size_t;
 
     /// Return the flattened offset to the start of the given order.
@@ -199,16 +276,16 @@ class PackedIndexRange {
     {}
 
     /// Return an iterator to the start of the range.
-    constexpr PackedIndexIterator begin() const noexcept { return _begin; }
+    constexpr iterator begin() const noexcept { return _begin; }
 
     /// Return an iterator to the start of the range.
-    constexpr PackedIndexIterator cbegin() const noexcept { return _begin; }
+    constexpr iterator cbegin() const noexcept { return _begin; }
 
     /// Return an iterator to one past the end of the range.
-    constexpr PackedIndexIterator end() const noexcept { return _end; }
+    constexpr iterator end() const noexcept { return _end; }
 
     /// Return an iterator to one past the end of the range.
-    constexpr PackedIndexIterator cend() const noexcept { return _end; }
+    constexpr iterator cend() const noexcept { return _end; }
 
     /// Return the number of elements in the flattened expansion.
     constexpr std::size_t size() const noexcept { return _end->flat - _begin->flat; }
@@ -227,8 +304,8 @@ class PackedIndexRange {
     }
 
 private:
-    PackedIndexIterator _begin;
-    PackedIndexIterator _end;
+    iterator _begin;
+    iterator _end;
 };
 
 }}} // namespace lsst::geom::polynomials

include/lsst/geom/polynomials/PolynomialBasis2d.h

@@ -29,11 +29,25 @@
 
 namespace lsst { namespace geom { namespace polynomials {
 
-/// A Basis2d for standard polynomials.
-using PolynomialBasis2d = PackedBasis2d<PolynomialBasis1d>;
+/// A Basis2d for standard polynomials, templated on packing order.
+template <PackingOrder packing>
+using PolynomialBasis2d = PackedBasis2d<PolynomialBasis1d, packing>;
 
-/// A Basis2d for standard polynomials.
-using ScaledPolynomialBasis2d = ScaledBasis2d<PolynomialBasis2d>;
+/// A Basis2d for scaled standard polynomials, templated on packing order.
+template <PackingOrder packing>
+using ScaledPolynomialBasis2d = ScaledBasis2d<PolynomialBasis2d<packing>>;
+
+/// A Basis2d for standard polynomials, ordered via PackingOrder::XY.
+using PolynomialBasis2dXY = PolynomialBasis2d<PackingOrder::XY>;
+
+/// A Basis2d for standard polynomials, ordered via PackingOrder::YX.
+using PolynomialBasis2dYX = PolynomialBasis2d<PackingOrder::YX>;
+
+/// A Basis2d for scaled standard polynomials, ordered via PackingOrder::XY.
+using ScaledPolynomialBasis2dXY = ScaledPolynomialBasis2d<PackingOrder::XY>;
+
+/// A Basis2d for scaled standard polynomials, ordered via PackingOrder::YX.
+using ScaledPolynomialBasis2dYX = ScaledPolynomialBasis2d<PackingOrder::YX>;
 
 }}} // namespace lsst::geom::polynomials
 

include/lsst/geom/polynomials/PolynomialFunction2d.h

@@ -28,10 +28,24 @@
 namespace lsst { namespace geom { namespace polynomials {
 
 /// A Function2d for standard polynomials.
-using PolynomialFunction2d = Function2d<PolynomialBasis2d>;
+template <PackingOrder packing>
+using PolynomialFunction2d = Function2d<PolynomialBasis2d<packing>>;
 
 /// A Function2d for scaled standard polynomials.
-using ScaledPolynomialFunction2d = Function2d<ScaledPolynomialBasis2d>;
+template <PackingOrder packing>
+using ScaledPolynomialFunction2d = Function2d<ScaledPolynomialBasis2d<packing>>;
+
+/// A Function2d for standard polynomials, ordered via PackingOrder::XY.
+using PolynomialFunction2dXY = PolynomialFunction2d<PackingOrder::XY>;
+
+/// A Function2d for standard polynomials, ordered via PackingOrder::YX.
+using PolynomialFunction2dYX = PolynomialFunction2d<PackingOrder::YX>;
+
+/// A Function2d for scaled standard polynomials, ordered via PackingOrder::XY.
+using ScaledPolynomialFunction2dXY = ScaledPolynomialFunction2d<PackingOrder::XY>;
+
+/// A Function2d for scaled standard polynomials, ordered via PackingOrder::YX.
+using ScaledPolynomialFunction2dYX = ScaledPolynomialFunction2d<PackingOrder::YX>;
 
 /**
  *  Calculate the standard polynomial function that is equivalent to a scaled
@@ -42,7 +56,8 @@ using ScaledPolynomialFunction2d = Function2d<ScaledPolynomialBasis2d>;
  *  basis and then (if necessary) use `unscaled` to compute the coefficients
  *  of an equivalent unscaled polynomial.
  */
-PolynomialFunction2d unscaled(ScaledPolynomialFunction2d const & f);
+template <PackingOrder packing>
+PolynomialFunction2d<packing> unscaled(ScaledPolynomialFunction2d<packing> const & f);
 
 }}} // namespace lsst::geom::polynomials