affine_map.hpp

Go to the documentation of this file.

/*
 * Copyright © 2025-2026 Akreon LLC
 * SPDX-License-Identifier: Apache-2.0
 */
 
#pragma once
 
#include "affine_traits.hpp"
#include "akreon/linalg/assert.hpp"
#include "akreon/linalg/operator/matrix_inverse.hpp"
 
namespace affine {
 
template <linalg::Field T>
using AffineMap2d = AffineMap<linalg::Matrix<T, 2, 2>>;
 
template <linalg::Field T>
using AffineMap3d = AffineMap<linalg::Matrix<T, 3, 3>>;
 
template <linalg::SquareMatrix Mat>
constexpr AffineMap<Mat> make_map_id() {
  using T = linalg::matrix_scalar_t<Mat>;
  static constexpr std::size_t N = linalg::mat_rows_v<Mat>;
  return AffineMap<Mat>(linalg::make_identity<T, N>(), linalg::Vector<T, N>{});
}
 
template <linalg::MatrixContainer Mat>
  requires(linalg::mat_rows_v<Mat> == 1)
constexpr AffineMap<Mat> make_selector_map(std::size_t idx) {
  AK_ASSERT(idx < linalg::mat_cols_v<Mat> &&
            "make_selector_map index must be within [0, mat_cols)");
 
  using T = linalg::matrix_scalar_t<Mat>;
  Mat A{};
  A(0, idx) = T(1);
  return AffineMap<Mat>(A, linalg::Vector<T, 1>{});
}
 
template <linalg::MatrixContainer Mat>
class AffineMap {
public:
  using T = linalg::matrix_scalar_t<Mat>;
 
  static constexpr std::size_t M = linalg::mat_rows_v<Mat>;
  static constexpr std::size_t N = linalg::mat_cols_v<Mat>;
 
  using DomVec = linalg::Vector<T, N>;
 
  using AmbientVec = linalg::Vector<T, M>;
 
  constexpr explicit AffineMap() {}
 
  constexpr explicit AffineMap(const Mat& U_in) : U(U_in) {}
 
  constexpr explicit AffineMap(const Mat& U_in, const AmbientVec& c_in)
      : U(U_in), c(c_in) {}
 
  constexpr explicit AffineMap(const linalg::Matrix<T, M + 1, N + 1>& H)
      : U(get_block<0, M, 0, N>(H)),
        c(as_vector(get_block<0, M, N, N + 1>(H))) {}
 
  template <linalg::MatrixContainer OtherMat>
    requires(std::same_as<T, linalg::matrix_scalar_t<OtherMat>> &&
             linalg::mat_rows_v<OtherMat> == N)
  constexpr auto compose(const AffineMap<OtherMat>& other) const {
    constexpr std::size_t P = linalg::mat_cols_v<OtherMat>;
    const auto U_comp = U * other.linear_part();
    const auto c_comp = U * other.translation() + c;
    return AffineMap<linalg::Matrix<T, M, P>>(U_comp, c_comp);
  }
 
  constexpr auto inverse() const {
    static_assert(M >= N, "short-fat TBD");
    const auto [U_inv, event] = linalg::MatrixInverse<decltype(U)>()(U);
    static_cast<void>(event);
    return AffineMap(U_inv, T(-1) * U_inv * c);
  }
 
  constexpr auto linear_part() const { return U; }
 
  constexpr auto translation() const { return c; }
 
  constexpr bool operator==(const AffineMap<Mat>& other) const {
    return is_equal(other);
  }
 
  constexpr bool operator!=(const AffineMap<Mat>& other) const {
    return !is_equal(other);
  }
 
  constexpr auto operator()(const DomVec& x) const { return U * x + c; }
 
private:
  Mat U;        //<  valid linear transformation mapping
  AmbientVec c; //< c bias vector that translates the linear part to
                // a valid affine subspace
 
  constexpr bool is_equal(const AffineMap<Mat>& other) const {
    const auto linear_equal = all_equal(other.linear_part(), linear_part());
    const auto translation_equal =
        all_equal(other.translation(), translation());
    return linear_equal && translation_equal;
  }
};
 
} // namespace affine