diff --git a/.ci/docker/lint.Dockerfile b/.ci/docker/lint.Dockerfile
index f5bf56e..dc8c454 100644
--- a/.ci/docker/lint.Dockerfile
+++ b/.ci/docker/lint.Dockerfile
@@ -18,7 +18,9 @@ RUN wget -qO- https://apt.llvm.org/llvm-snapshot.gpg.key | tee /etc/apt/trusted.
 
 FROM base AS deploy
 
-RUN echo deb http://apt.llvm.org/bookworm/ llvm-toolchain-bookworm-19 main \
+ARG clang_format_VERSION 19
+
+RUN echo deb http://apt.llvm.org/bookworm/ llvm-toolchain-bookworm-${clang_format_VERSION} main \
 >/etc/apt/sources.list.d/llvm.list
 
 COPY --from=tools /etc/apt/trusted.gpg.d/apt.llvm.org.asc /etc/apt/trusted.gpg.d/apt.llvm.org.asc
@@ -26,5 +28,5 @@ COPY --from=tools /etc/apt/trusted.gpg.d/apt.llvm.org.asc /etc/apt/trusted.gpg.d
 RUN --mount=type=cache,target=/var/cache/apt \
 apt-get update && DEBIAN_FRONTEND=noninteractive apt-get install -y \
 --no-install-recommends --no-install-suggests \
-clang-format-19 \
+clang-format-${clang_format_VERSION} \
 fd-find
diff --git a/.clang-format b/.clang-format
index 1d250ca..144e202 100644
--- a/.clang-format
+++ b/.clang-format
@@ -95,7 +95,7 @@ SpacesInParentheses: false
 SpacesInSquareBrackets: false
 SpaceAfterTemplateKeyword: false
 FixNamespaceComments: true
-Standard:           c++17
+Standard:           c++20
 TabWidth:        4
 UseTab:          Never
 ...
diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml
index 9721817..04f07d3 100644
--- a/.github/workflows/linux.yml
+++ b/.github/workflows/linux.yml
@@ -5,7 +5,7 @@ on: [push, pull_request]
 jobs:
   lint:
     runs-on: ubuntu-24.04
-    container: sergiud/hogpp:lint
+    container: sergiud/hogpp:lint-18
     defaults:
       run:
         shell: bash -e -o pipefail {0}
@@ -17,12 +17,11 @@ jobs:
 
       - name: Check code style
         run: |
-          fdfind -g '*.[ch]pp' -x clang-format-19 --dry-run --Werror
+          fdfind -g '*.[ch]pp' -x clang-format-18 --dry-run --Werror
 
   build-native:
     name: ${{matrix.os}}-GCC-${{matrix.build_type}}
     runs-on: ${{matrix.os}}
-    needs: lint
     defaults:
       run:
         shell: bash
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 996a7d9..430a33c 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -255,11 +255,11 @@ lint:pypy3.10:
   needs: []
 
 lint:cpp:
-  image: sergiud/hogpp:lint
+  image: sergiud/hogpp:lint-18
   stage: test
   allow_failure: true
   script:
-    - fdfind -g '*.[ch]pp' -x clang-format-19 --dry-run --Werror
+    - fdfind -g '*.[ch]pp' -x clang-format-18 --dry-run --Werror
   needs: []
 
 build:bookworm:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d129e28..4b24a90 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -84,6 +84,7 @@ endif (MSVC AND MSVC_VERSION GREATER_EQUAL 1930)
 check_cxx_symbol_exists (std::execution::par execution HAVE_EXECUTION)
 
 if (WIN32)
+  add_compile_definitions (NOMINMAX)
   add_compile_definitions (WIN32_LEAN_AND_MEAN)
 endif (WIN32)
 
@@ -100,7 +101,7 @@ target_include_directories (hogpp INTERFACE
   $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
 )
 
-target_compile_features (hogpp INTERFACE cxx_std_17)
+target_compile_features (hogpp INTERFACE cxx_std_20)
 target_link_libraries (hogpp INTERFACE
   Eigen3::Eigen
   fmt::fmt
@@ -110,7 +111,6 @@ if (pybind11_FOUND AND Python_Development.Module_FOUND)
   add_library (python_type_caster STATIC
     python/type_caster/bounds.cpp
     python/type_caster/bounds.hpp
-    python/type_caster/cartesianproduct.hpp
     python/type_caster/stride.hpp
     python/type_caster/tensor.hpp
     python/type_caster/typesequence.hpp
@@ -124,7 +124,6 @@ if (pybind11_FOUND AND Python_Development.Module_FOUND)
     target_link_libraries (python_type_caster PUBLIC opencv_core)
   endif (OpenCV_FOUND)
 
-  target_compile_features (python_type_caster PUBLIC cxx_std_17)
   target_link_libraries (python_type_caster PUBLIC
     Eigen3::Eigen
     hogpp::hogpp
@@ -188,6 +187,18 @@ endif (pybind11_FOUND AND Python_Development.Module_FOUND)
 
 if (BUILD_TESTING)
   if (Boost_unit_test_framework_FOUND)
+    add_executable (test_integral_histogram
+      tests/cpp/test_integral_histogram.cpp
+    )
+    target_link_libraries (test_integral_histogram PRIVATE
+      Boost::unit_test_framework
+      Eigen3::Eigen
+      hogpp::hogpp
+    )
+
+    add_test (NAME integral_histogram_full COMMAND test_integral_histogram -t full)
+    add_test (NAME integral_histogram_sub COMMAND test_integral_histogram -t sub)
+
     add_executable (test_binning tests/cpp/test_binning.cpp)
     target_link_libraries (test_binning PRIVATE hogpp::hogpp
       Boost::unit_test_framework)
diff --git a/docs/build.rst b/docs/build.rst
index 03173c0..0afd1d7 100644
--- a/docs/build.rst
+++ b/docs/build.rst
@@ -4,7 +4,7 @@ Compiling from Source
 Requirements
 ------------
 
--  C++17 compiler
+-  C++20 compiler
 -  `CMake <https://gitlab.kitware.com/cmake/cmake>`__ 3.18
 -  `Eigen <https://gitlab.com/libeigen/eigen>`__ 3.3.7
 -  `fmt <https://github.com/fmtlib/fmt>`__ 6.0
diff --git a/python/type_caster/cartesianproduct.hpp b/include/hogpp/cartesianproduct.hpp
similarity index 57%
rename from python/type_caster/cartesianproduct.hpp
rename to include/hogpp/cartesianproduct.hpp
index 3efadeb..58d5210 100644
--- a/python/type_caster/cartesianproduct.hpp
+++ b/include/hogpp/cartesianproduct.hpp
@@ -17,45 +17,49 @@
 // limitations under the License.
 //
 
-#ifndef PYTHON_TYPE_CASTER_CARTESIANPRODUCT_HPP
-#define PYTHON_TYPE_CASTER_CARTESIANPRODUCT_HPP
+#ifndef HOGPP_CARTESIANPRODUCT_HPP
+#define HOGPP_CARTESIANPRODUCT_HPP
 
+#include <concepts>
 #include <tuple>
 #include <type_traits>
 #include <utility>
 
-template<class... Sizes, class Function, class... Loops>
+namespace hogpp {
+
+template<std::integral... Sizes, class Body, std::integral... Loops>
 constexpr void cartesianProduct(
-    const std::tuple<Sizes...>& /*s*/, const std::tuple<Loops...>& l,
-    Function&& func,
-    std::index_sequence<> /*unused*/) noexcept(noexcept(func(l)))
+    [[maybe_unused]] const std::tuple<Sizes...>& s,
+    const std::tuple<Loops...>& l, Body&& body,
+    std::index_sequence<> /*unused*/) noexcept(noexcept(body(l)))
 {
-    func(l);
+    body(l);
 }
 
-template<class... Sizes, class Function, class... Loops, std::size_t Index,
-         std::size_t... Indices>
+template<std::integral... Sizes, class Body, std::integral... Loops,
+         std::size_t Index, std::size_t... Indices>
 constexpr void cartesianProduct(
-    const std::tuple<Sizes...>& s, const std::tuple<Loops...>& l,
-    Function&& func,
+    const std::tuple<Sizes...>& s, const std::tuple<Loops...>& l, Body&& body,
     std::index_sequence<
         Index,
-        Indices...> /*unused*/) noexcept(noexcept(func(std::declval<std::tuple<Sizes...>>())))
+        Indices...> /*unused*/) noexcept(noexcept(body(std::declval<std::tuple<Sizes...>>())))
 {
     for (auto i = 0; i < std::get<Index>(s); ++i) {
         cartesianProduct(s, std::tuple_cat(l, std::make_tuple(i)),
-                         std::forward<Function>(func),
+                         std::forward<Body>(body),
                          std::index_sequence<Indices...>{});
     }
 }
 
-template<class... Sizes, class Function>
-constexpr void
-cartesianProduct(const std::tuple<Sizes...>& s, Function&& func) noexcept(
-    noexcept(func(std::declval<std::tuple<Sizes...>>())))
+template<std::integral... Sizes, class Body>
+constexpr void cartesianProduct(
+    const std::tuple<Sizes...>& s,
+    Body&& body) noexcept(noexcept(body(std::declval<std::tuple<Sizes...>>())))
 {
-    cartesianProduct(s, std::tuple<>{}, std::forward<Function>(func),
+    cartesianProduct(s, std::tuple<>{}, std::forward<Body>(body),
                      std::index_sequence_for<Sizes...>{});
 }
 
-#endif // PYTHON_TYPE_CASTER_CARTESIANPRODUCT_HPP
+} // namespace hogpp
+
+#endif // HOGPP_CARTESIANPRODUCT_HPP
diff --git a/include/hogpp/constants.hpp b/include/hogpp/constants.hpp
index 6aea0ae..036a85d 100644
--- a/include/hogpp/constants.hpp
+++ b/include/hogpp/constants.hpp
@@ -20,18 +20,18 @@
 #ifndef HOGPP_CONSTANTS_HPP
 #define HOGPP_CONSTANTS_HPP
 
+#include <numbers>
+
 namespace hogpp::constants {
 
 template<class T>
-inline constexpr T pi = T(3.1415926535897932384626433832795028841971693993751);
+inline constexpr T pi = std::numbers::pi_v<T>;
 
 template<class T>
-inline constexpr T two_pi =
-    T(6.28318530717958623199592693708837032318115234375);
+inline constexpr T two_pi = 2 * pi<T>;
 
 template<class T>
-inline constexpr T half_pi =
-    T(1.5707963267948965579989817342720925807952880859375);
+inline constexpr T half_pi = pi<T> / 2;
 
 } // namespace hogpp::constants
 
diff --git a/include/hogpp/integralhistogram.hpp b/include/hogpp/integralhistogram.hpp
new file mode 100644
index 0000000..7904036
--- /dev/null
+++ b/include/hogpp/integralhistogram.hpp
@@ -0,0 +1,357 @@
+//
+// HOGpp - Fast histogram of oriented gradients computation using integral
+// histograms
+//
+// Copyright 2021 Sergiu Deitsch <sergiu.deitsch@gmail.com>
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//   http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+//
+
+#ifndef HOGPP_INTEGRALHISTOGRAM_HPP
+#define HOGPP_INTEGRALHISTOGRAM_HPP
+
+#include <bit>
+#include <concepts>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+#include <unsupported/Eigen/CXX11/Tensor>
+
+#include <hogpp/cartesianproduct.hpp>
+
+namespace hogpp {
+
+template<class T>
+concept Arithmetic = std::is_arithmetic_v<T>;
+
+namespace detail {
+
+template<class T>
+struct IsTuple : std::false_type
+{
+};
+
+template<class... T>
+struct IsTuple<std::tuple<T...>> : std::true_type
+{
+};
+
+template<class T>
+inline constexpr bool IsTuple_v = IsTuple<T>::value;
+
+template<class T>
+concept TupleOfIntegrals = IsTuple_v<std::decay_t<T>>;
+
+template<std::size_t K, class... Types, std::size_t... Indices>
+[[nodiscard]] constexpr decltype(auto) neighbor(
+    const std::tuple<Types...>& i,
+    std::index_sequence<Indices...> /*unused*/) noexcept
+{
+    return std::make_tuple(
+        (std::get<Indices>(i) + (((1U << Indices) & K) != 0))...);
+}
+
+template<std::size_t K, class... Types>
+[[nodiscard]] constexpr decltype(auto) neighbor(
+    const std::tuple<Types...>& i) noexcept
+{
+    return neighbor<K>(i, std::index_sequence_for<Types...>{});
+}
+
+template<class Tensor, class... Types>
+[[nodiscard]] constexpr decltype(auto) chip(
+    Tensor&& t, [[maybe_unused]] const std::tuple<Types...>& i,
+    std::integer_sequence<Eigen::DenseIndex> /*unused*/) noexcept
+{
+    return std::forward<Tensor>(t);
+}
+
+template<class Tensor, class... Types, Eigen::DenseIndex Index,
+         Eigen::DenseIndex... Indices>
+[[nodiscard]] constexpr auto chip(
+    Tensor&& t, const std::tuple<Types...>& i,
+    std::integer_sequence<Eigen::DenseIndex, Index,
+                          Indices...> /*unused*/) noexcept
+{
+    return chip(t.template chip<0>(std::get<Index>(i)), i,
+                std::integer_sequence<Eigen::DenseIndex, Indices...>{});
+}
+
+template<class Tensor, class... Types>
+[[nodiscard]] constexpr decltype(auto) chip(
+    Tensor&& t, const std::tuple<Types...>& i) noexcept
+{
+    return chip(std::forward<Tensor>(t), i,
+                std::make_integer_sequence<Eigen::DenseIndex,
+                                           static_cast<Eigen::DenseIndex>(
+                                               sizeof...(Types))>{});
+}
+
+template<Eigen::DenseIndex Offset, class... Types, Eigen::DenseIndex... Indices>
+[[nodiscard]] constexpr decltype(auto) offset(
+    const std::tuple<Types...>& i,
+    std::integer_sequence<Eigen::DenseIndex, Indices...> /*unused*/) noexcept
+{
+    return std::make_tuple((std::get<Indices>(i) + Offset)...);
+}
+
+template<Eigen::DenseIndex Offset, class... Types>
+[[nodiscard]] constexpr decltype(auto) offset(
+    const std::tuple<Types...>& i) noexcept
+{
+    return offset<Offset>(
+        i, std::make_integer_sequence<Eigen::DenseIndex,
+                                      static_cast<Eigen::DenseIndex>(
+                                          sizeof...(Types))>{});
+}
+
+template<bool Sign, class Tensor, std::enable_if_t<Sign>* = nullptr>
+[[nodiscard]] constexpr decltype(auto) negate(Tensor&& t) noexcept(noexcept(-t))
+{
+    return -t;
+}
+
+template<bool Sign, class Tensor, std::enable_if_t<!Sign>* = nullptr>
+[[nodiscard]] constexpr decltype(auto) negate(Tensor&& t) noexcept
+{
+    return std::forward<Tensor>(t);
+}
+
+template<class Tensor, class... Types, std::size_t... Combinations>
+[[nodiscard]] constexpr decltype(auto) propagate(
+    Tensor&& t, const std::tuple<Types...>& i,
+    std::index_sequence<Combinations...> /*unused*/) noexcept
+{
+    return (negate<Combinations == 0>(
+                chip(std::forward<Tensor>(t), neighbor<Combinations>(i))) +
+            ...);
+}
+
+template<class Tensor, class... Types>
+[[nodiscard]] constexpr decltype(auto) propagate(
+    Tensor&& t, const std::tuple<Types...>& i) noexcept
+{
+    constexpr std::size_t N = (1U << sizeof...(Types)) - 1U;
+    return propagate(std::forward<Tensor>(t), i, std::make_index_sequence<N>{});
+}
+
+template<bool Sign, class... Types1, class... Types2,
+         std::enable_if_t<Sign>* = nullptr>
+[[nodiscard]] constexpr decltype(auto) select(
+    const std::tuple<Types1...>& a,
+    [[maybe_unused]] const std::tuple<Types2...>& b)
+{
+    return a;
+}
+
+template<bool Sign, class... Types1, class... Types2,
+         std::enable_if_t<!Sign>* = nullptr>
+[[nodiscard]] constexpr decltype(auto) select(
+    [[maybe_unused]] const std::tuple<Types1...>& a,
+    const std::tuple<Types2...>& b)
+{
+    return b;
+}
+
+template<std::size_t K, std::integral... Types1, std::integral... Types2,
+         std::size_t... Indices>
+[[nodiscard]] constexpr decltype(auto) cross1(
+    const std::tuple<Types1...>& a, const std::tuple<Types2...>& b,
+    std::index_sequence<Indices...> /*unused*/) noexcept
+{
+    return std::make_tuple(
+        std::get<Indices>(select<(((1U << Indices) & K) != 0)>(a, b))...);
+}
+
+template<std::integral... Types1, std::integral... Types2,
+         std::size_t... Combinations>
+[[nodiscard]] constexpr decltype(auto) combinations(
+    const std::tuple<Types1...>& a, const std::tuple<Types2...>& b,
+    std::index_sequence<Combinations...> /*unused*/) noexcept
+{
+    static_assert(sizeof...(Types1) == sizeof...(Types2),
+                  "combinations tuple size does not match");
+
+    return std::make_tuple(
+        cross1<Combinations>(a, b, std::index_sequence_for<Types1...>{})...);
+}
+
+template<std::integral... Types1, std::integral... Types2>
+[[nodiscard]] constexpr decltype(auto) combinations(
+    const std::tuple<Types1...>& a, const std::tuple<Types2...>& b) noexcept
+{
+    static_assert(sizeof...(Types1) == sizeof...(Types2),
+                  "combinations tuple size does not match");
+
+    constexpr std::size_t N = (1U << sizeof...(Types1));
+    return combinations(a, b, std::make_index_sequence<N>{});
+}
+
+template<std::size_t K, class Tensor>
+[[nodiscard]] constexpr decltype(auto) sign(Tensor&& t)
+{
+    // Cannot use std::bitset::count because it's not constexpr.
+    constexpr std::size_t N = std::popcount(K);
+    // Bit set: + otherwise -
+    // If the number of set bits does not match the number of total bits, we
+    // need to negate the term.
+    return negate<(N & 1) != 0>(t);
+}
+
+template<class Tensor, class... Types, std::size_t... Indices>
+[[nodiscard]] constexpr decltype(auto) intersect(
+    Tensor&& t, const std::tuple<Types...>& ab,
+    std::index_sequence<Indices...> /*unused*/)
+{
+    return (
+        sign<Indices>(chip(std::forward<Tensor>(t), std::get<Indices>(ab))) +
+        ...);
+}
+
+template<class Tensor, TupleOfIntegrals... Types>
+[[nodiscard]] constexpr decltype(auto) intersect(Tensor&& t,
+                                                 const std::tuple<Types...>& ab)
+{
+    return intersect(std::forward<Tensor>(t), ab,
+                     std::index_sequence_for<Types...>{});
+}
+
+template<class Tensor, std::integral... Types1, std::integral... Types2>
+[[nodiscard]] constexpr decltype(auto) intersect(Tensor&& t,
+                                                 const std::tuple<Types1...>& a,
+                                                 const std::tuple<Types2...>& b)
+{
+    return intersect(std::forward<Tensor>(t), combinations(a, b));
+}
+
+} // namespace detail
+
+template<Arithmetic T, Eigen::DenseIndex N, Eigen::DenseIndex K = 1,
+         int Options = 0>
+class IntegralHistogram
+{
+public:
+    static_assert(
+        N > 0,
+        "IntegralHistogram number of input space dimensions must be positive");
+    static_assert(K > 0,
+                  "IntegralHistogram number of histogram space dimensions must "
+                  "be positive");
+
+    using Tensor = Eigen::Tensor<T, N + K, Options>;
+
+    template<std::integral... InputDims, std::integral... BinDims>
+    void resize(const std::tuple<InputDims...>& dims,
+                const std::tuple<BinDims...>& bins)
+    {
+        static_assert(sizeof...(InputDims) == N,
+                      "IntegralHistogram resize number of provided input space "
+                      "dimensions does not match the actual number of input "
+                      "space dimensions");
+        static_assert(
+            sizeof...(BinDims) == K,
+            "IntegralHistogram resize number of provided histogram space "
+            "dimensions does not match the actual number of histogram "
+            "space dimensions");
+
+        resize(std::tuple_cat(detail::offset<1>(dims), bins));
+    }
+
+    template<std::integral... InputDims, std::integral Bins,
+             std::enable_if_t<K == 1>* = nullptr>
+    void resize(const std::tuple<InputDims...>& dims, Bins bins)
+    {
+        resize(dims, std::make_tuple(bins));
+    }
+
+    // TODO Tensor of compatible dimension
+    template<class Binning>
+    void scan(Binning binning) // requires std::invocable<Binning,
+    {
+        scan(histogram_.dimensions(),
+             std::make_integer_sequence<Eigen::DenseIndex, N>{},
+             std::move(binning));
+    }
+
+    template<std::integral... Types1, std::integral... Types2>
+    [[nodiscard]] decltype(auto) intersect(const std::tuple<Types1...>& a,
+                                           const std::tuple<Types2...>& b) const
+        noexcept(noexcept(detail::intersect(std::declval<Tensor>(), a, b)))
+    {
+        static_assert(sizeof...(Types1) == sizeof...(Types2),
+                      "IntegralHistogram intersect coordinates dimensionality "
+                      "does not match");
+        return detail::intersect(histogram_, a, b);
+    }
+
+    [[nodiscard]] const Tensor& histogram() const noexcept
+    {
+        return histogram_;
+    }
+
+    template<class Derived>
+    void setHistogram(
+        const Eigen::TensorBase<Derived, Eigen::ReadOnlyAccessors>& value)
+    {
+        histogram_ = value.eval();
+    }
+
+    [[nodiscard]] bool isEmpty() const noexcept
+    {
+        return histogram_.size() == 0;
+    }
+
+private:
+    template<Eigen::DenseIndex... Indices, class Binning>
+    void scan(
+        const std::array<Eigen::DenseIndex, static_cast<std::size_t>(N + K)>&
+            dims,
+        std::integer_sequence<Eigen::DenseIndex, Indices...> /*unused*/,
+        Binning
+            binning) // noexcept(noexcept(binning(std::declval<Eigen::Tensor<T,
+                     // N> >()
+    {
+        histogram_.setZero();
+
+        // Wavefront scan
+        cartesianProduct(std::make_tuple((dims[Indices] - 1)...),
+                         [this, &binning](const auto& i) constexpr {
+                             auto current = detail::offset<1>(i);
+                             auto&& H = detail::chip(histogram_, current);
+                             H = detail::propagate(histogram_, i);
+                             binning(H, i);
+                         });
+    }
+
+    template<std::integral... Types>
+    void resize(const std::tuple<Types...>& dims)
+    {
+        resize(dims, std::make_integer_sequence<Eigen::DenseIndex,
+                                                static_cast<Eigen::DenseIndex>(
+                                                    sizeof...(Types))>{});
+    }
+
+    template<std::integral... Types, Eigen::DenseIndex... Indices>
+    void resize(const std::tuple<Types...>& dims,
+                std::integer_sequence<Eigen::DenseIndex, Indices...> /*unused*/)
+    {
+        histogram_.resize(std::get<Indices>(dims)...);
+    }
+
+    Tensor histogram_;
+};
+
+} // namespace hogpp
+
+#endif // HOGPP_INTEGRALHISTOGRAM_HPP
diff --git a/include/hogpp/integralhogdescriptor.hpp b/include/hogpp/integralhogdescriptor.hpp
index 151cd16..3b65cca 100644
--- a/include/hogpp/integralhogdescriptor.hpp
+++ b/include/hogpp/integralhogdescriptor.hpp
@@ -33,6 +33,7 @@
 #include <hogpp/bounds.hpp>
 #include <hogpp/gradient.hpp>
 #include <hogpp/gradientmagnitude.hpp>
+#include <hogpp/integralhistogram.hpp>
 #include <hogpp/l2hys.hpp>
 #include <hogpp/unsignedgradient.hpp>
 
@@ -264,54 +265,46 @@ class IntegralHOGDescriptor
         }
 
         const Eigen::Array2i dims{mags.dimension(0), mags.dimension(1)};
-        Eigen::Array2i dims1 = dims + 1;
-
-        histogram_.resize(dims1.x(), dims1.y(), bins_);
-        histogram_.setZero();
+        histogram_.resize(std::make_tuple(dims.x(), dims.y()), bins_);
 
         const Eigen::Tensor<Eigen::DenseIndex, 2, Eigen::RowMajor>& k =
             mags.argmax(2);
 
         const auto scale = static_cast<Scalar>(bins_ - 1);
 
-        // Wavefront scan
-        for (Eigen::DenseIndex i = 0; i < dims.x(); ++i) {
-            for (Eigen::DenseIndex j = 0; j < dims.y(); ++j) {
-                // Histogram propagation
-                const auto& a =
-                    histogram_.template chip<0>(i).template chip<0>(j + 1);
-                const auto& b =
-                    histogram_.template chip<0>(i + 1).template chip<0>(j);
-                const auto& c =
-                    histogram_.template chip<0>(i).template chip<0>(j);
-
-                histogram_.template chip<0>(i + 1).template chip<0>(j + 1) =
-                    a + b - c;
-
+        histogram_.scan(
+            [this, &k, &dxs, &dys, &mags, scale, &masked](
+                Eigen::TensorRef<Eigen::Tensor<Scalar, 1, DataLayout>> bins,
+                const auto& ij) {
+                (void)masked; // Avoid error: lambda capture 'masked' is not
+                              // used [-Werror,-Wunused-lambda-capture] on
+                              // AppleClang
                 if constexpr (!std::is_null_pointer_v<Masking>) {
-                    if (masked(i, j)) {
+                    if (std::apply(masked, ij)) {
                         // Skip masked out pixels
-                        continue;
+                        return;
                     }
                 }
 
                 // Select a channel with the maximum magnitude
-                Eigen::DenseIndex kk = k(i, j);
-                Scalar mag = mags(i, j, kk);
+                Eigen::DenseIndex kk = std::apply(k, ij);
+                const auto ijk = std::tuple_cat(ij, std::make_tuple(kk));
+
+                Scalar mag = std::apply(mags, ijk);
 
                 using std::fpclassify;
 
                 if (fpclassify(mag) == FP_ZERO) {
                     // No gradient; take a shortcut
-                    continue;
+                    return;
                 }
 
                 // The gradient magnitude cannot be negative (or zero at this
                 // point)
                 HOGPP_ASSUME(mag > 0);
 
-                Scalar dx = dxs(i, j, kk);
-                Scalar dy = dys(i, j, kk);
+                Scalar dx = std::apply(dxs, ijk);
+                Scalar dy = std::apply(dys, ijk);
 
                 // Gradient binning
                 Scalar weight = this->binning_(dx, dy);
@@ -341,8 +334,6 @@ class IntegralHOGDescriptor
                 HOGPP_ASSUME(bin1 <= bin2);
 
                 // Distribute weighted values to neighboring bins.
-                Eigen::TensorRef<Eigen::Tensor<Scalar, 1, DataLayout>> bins =
-                    histogram_.template chip<0>(i + 1).template chip<0>(j + 1);
                 Scalar& value1 = bins.coeffRef(bin1);
                 Scalar& value2 = bins.coeffRef(bin2);
 
@@ -351,19 +342,18 @@ class IntegralHOGDescriptor
                 // proportionally a higher magnitude.
                 value1 = fma(1 - alpha, mag, value1);
                 value2 = fma(alpha, mag, value2);
-            }
-        }
+            });
     }
 
     [[nodiscard]] Tensor5 features() const
     {
-        if (histogram_.size() == 0) {
+        if (histogram_.isEmpty()) {
             return Tensor5{};
         }
 
         return features(Bounds{0, 0,
-                               static_cast<int>(histogram_.dimension(1) - 1),
-                               static_cast<int>(histogram_.dimension(0) - 1)});
+                               static_cast<int>(histogram().dimension(1) - 1),
+                               static_cast<int>(histogram().dimension(0) - 1)});
     }
 
     [[nodiscard]] Tensor5 features(const Bounds& roi) const
@@ -388,8 +378,8 @@ class IntegralHOGDescriptor
                             dims.y())};
         }
 
-        const auto dimX = histogram_.dimension(0) - 1;
-        const auto dimY = histogram_.dimension(1) - 1;
+        const auto dimX = histogram().dimension(0) - 1;
+        const auto dimY = histogram().dimension(1) - 1;
         const Eigen::Array2i offset{roi.y, roi.x};
 
         if (offset.x() < 0 || dimX - offset.x() < 0) {
@@ -430,7 +420,7 @@ class IntegralHOGDescriptor
         Tensor5 X;
 
         X.resize(numBlocks.x(), numBlocks.y(), numCells.x(), numCells.y(),
-                 histogram_.dimension(2));
+                 histogram_.histogram().dimension(2));
         X.setZero();
 
         for (int i = 0; i < numBlocks.x(); ++i) {
@@ -445,18 +435,11 @@ class IntegralHOGDescriptor
                             blockOffset + cellOffset * cellSize_;
                         const Eigen::Array2i& offset2 = offset1 + cellSize_;
 
-                        const auto& a = histogram_.template chip<0>(offset2.x())
-                                            .template chip<0>(offset2.y());
-                        const auto& b = histogram_.template chip<0>(offset1.x())
-                                            .template chip<0>(offset2.y());
-                        const auto& c = histogram_.template chip<0>(offset2.x())
-                                            .template chip<0>(offset1.y());
-                        const auto& d = histogram_.template chip<0>(offset1.x())
-                                            .template chip<0>(offset1.y());
-
                         // Extract the histogram from the integral histogram
                         // as an intersection.
-                        const Eigen::Tensor<Scalar, 1> h = a - b - c + d;
+                        const Eigen::Tensor<Scalar, 1> h = histogram_.intersect(
+                            std::make_tuple(offset1.x(), offset1.y()),
+                            std::make_tuple(offset2.x(), offset2.y()));
 
                         X.template chip<0>(i)
                             .template chip<0>(j)
@@ -545,26 +528,26 @@ class IntegralHOGDescriptor
 
     [[nodiscard]] const Eigen::Tensor<Scalar, 3>& histogram() const noexcept
     {
-        return histogram_;
+        return histogram_.histogram();
     }
 
     template<class Derived>
     void setHistogram(
         const Eigen::TensorBase<Derived, Eigen::ReadOnlyAccessors>& value)
     {
-        histogram_ = value.eval();
+        histogram_.setHistogram(value);
     }
 
     [[nodiscard]] bool isEmpty() const noexcept
     {
-        return histogram_.size() == 0;
+        return histogram_.isEmpty();
     }
 
 private:
     using Tensor2 = Eigen::Tensor<Scalar, 2, Eigen::RowMajor>;
     using Tensor3 = Eigen::Tensor<Scalar, 3, Eigen::RowMajor>;
 
-    Eigen::Tensor<Scalar, 3> histogram_;
+    IntegralHistogram<Scalar, 2, 1> histogram_;
     Eigen::DenseIndex bins_ = 9;
     Eigen::Array2i cellSize_{8, 8};
     Eigen::Array2i blockSize_ = cellSize_ * 2;
diff --git a/python/binning.hpp b/python/binning.hpp
index 7a722a4..cc6e4f4 100644
--- a/python/binning.hpp
+++ b/python/binning.hpp
@@ -20,10 +20,12 @@
 #ifndef PYTHON_HOGPP_BINNING_HPP
 #define PYTHON_HOGPP_BINNING_HPP
 
-#include <pybind11/pybind11.h>
-
+// FIXME GCC 14.x workaround for https://github.com/pybind/pybind11/pull/5208
+#include <algorithm>
 #include <variant>
 
+#include <pybind11/pybind11.h>
+
 #include <hogpp/signedgradient.hpp>
 #include <hogpp/unsignedgradient.hpp>
 
diff --git a/python/blocknormalizer.hpp b/python/blocknormalizer.hpp
index 03f0856..84ee118 100644
--- a/python/blocknormalizer.hpp
+++ b/python/blocknormalizer.hpp
@@ -20,12 +20,14 @@
 #ifndef PYTHON_HOGPP_BLOCKNORMALIZER_HPP
 #define PYTHON_HOGPP_BLOCKNORMALIZER_HPP
 
-#include <pybind11/pybind11.h>
-
+// FIXME GCC 14.x workaround for https://github.com/pybind/pybind11/pull/5208
+#include <algorithm>
 #include <optional>
 #include <type_traits>
 #include <variant>
 
+#include <pybind11/pybind11.h>
+
 #include <hogpp/l1hys.hpp>
 #include <hogpp/l1norm.hpp>
 #include <hogpp/l1sqrt.hpp>
diff --git a/python/hogpp.cpp b/python/hogpp.cpp
index de53a65..70f2395 100644
--- a/python/hogpp.cpp
+++ b/python/hogpp.cpp
@@ -17,7 +17,8 @@
 // limitations under the License.
 //
 
-#include <optional>
+// FIXME GCC 14.x workaround for https://github.com/pybind/pybind11/pull/5208
+#include <algorithm>
 #include <variant>
 
 #include <pybind11/eigen.h>
diff --git a/python/magnitude.hpp b/python/magnitude.hpp
index c95cb7a..aabad34 100644
--- a/python/magnitude.hpp
+++ b/python/magnitude.hpp
@@ -20,11 +20,13 @@
 #ifndef PYTHON_HOGPP_MAGNITUDE_HPP
 #define PYTHON_HOGPP_MAGNITUDE_HPP
 
-#include <pybind11/pybind11.h>
-
+// FIXME GCC 14.x workaround for https://github.com/pybind/pybind11/pull/5208
+#include <algorithm>
 #include <type_traits>
 #include <variant>
 
+#include <pybind11/pybind11.h>
+
 #include <hogpp/gradientmagnitude.hpp>
 #include <hogpp/gradientsqrtmagnitude.hpp>
 #include <hogpp/gradientsquaremagnitude.hpp>
diff --git a/python/type_caster/bounds.hpp b/python/type_caster/bounds.hpp
index 88db88a..245704b 100644
--- a/python/type_caster/bounds.hpp
+++ b/python/type_caster/bounds.hpp
@@ -20,6 +20,9 @@
 #ifndef PYTHON_TYPE_CASTER_BOUNDS_HPP
 #define PYTHON_TYPE_CASTER_BOUNDS_HPP
 
+// FIXME GCC 14.x workaround for https://github.com/pybind/pybind11/pull/5208
+#include <algorithm>
+
 #include <pybind11/pybind11.h>
 
 #include <hogpp/bounds.hpp>
diff --git a/python/type_caster/opencv.hpp b/python/type_caster/opencv.hpp
index ec78836..b98a2ff 100644
--- a/python/type_caster/opencv.hpp
+++ b/python/type_caster/opencv.hpp
@@ -34,7 +34,8 @@
 #include <pybind11/cast.h>
 #include <pybind11/numpy.h>
 
-#include "cartesianproduct.hpp"
+#include <hogpp/cartesianproduct.hpp>
+
 #include "stride.hpp"
 #include "typesequence.hpp"
 
@@ -141,7 +142,7 @@ class type_caster<cv::Mat>
                     in.create(static_cast<int>(shape[0]),
                               static_cast<int>(shape[1]), type);
 
-                    cartesianProduct(
+                    hogpp::cartesianProduct(
                         std::make_tuple(shape[0], shape[1]),
                         [&in, p, strides](const auto& i) constexpr {
                             assign<T>(in, p, i, strides);
@@ -177,7 +178,7 @@ class type_caster<cv::Mat>
                         in.create(static_cast<int>(shape[0]),
                                   static_cast<int>(shape[1]), type);
 
-                        cartesianProduct(
+                        hogpp::cartesianProduct(
                             std::make_tuple(shape[0], shape[1], shape[2]),
                             [&in, p, strides](const auto& i) constexpr {
                                 assign<T>(in, p, i, strides);
@@ -197,7 +198,7 @@ class type_caster<cv::Mat>
                                  const std::tuple<Strides...>& strides,
                                  std::index_sequence<Indices...> /*unused*/)
     {
-        in.at<T>(std::get<Indices>(i)...) =
+        in.at<T>(static_cast<int>(std::get<Indices>(i))...) =
             *reinterpret_cast<const T*>(p + reduceIndices(i, strides));
     }
 
@@ -209,8 +210,8 @@ class type_caster<cv::Mat>
         const std::tuple<Strides...>& strides,
         std::index_sequence<Index0, Index1, Index2> /*unused*/)
     {
-        in.ptr<T>(std::get<Index0>(i),
-                  std::get<Index1>(i))[std::get<Index2>(i)] =
+        in.ptr<T>(static_cast<int>(std::get<Index0>(i)),
+                  static_cast<int>(std::get<Index1>(i)))[std::get<Index2>(i)] =
             *reinterpret_cast<const T*>(p + reduceIndices(i, strides));
     }
 
diff --git a/python/type_caster/tensor.hpp b/python/type_caster/tensor.hpp
index 86e3722..4f8078c 100644
--- a/python/type_caster/tensor.hpp
+++ b/python/type_caster/tensor.hpp
@@ -32,7 +32,8 @@
 #include <pybind11/cast.h>
 #include <pybind11/numpy.h>
 
-#include "cartesianproduct.hpp"
+#include <hogpp/cartesianproduct.hpp>
+
 #include "stride.hpp"
 
 template<class T, std::size_t N, std::size_t... Indices>
@@ -398,11 +399,12 @@ class type_caster
     void copy(const std::uint8_t* ptr, const std::tuple<Sizes...>& sizes,
               const std::tuple<Strides...>& strides)
     {
-        cartesianProduct(sizes, [this, ptr, strides](const auto& i) constexpr {
-            (void)this; // Avoid (incorrect) Clang -Wunused-lambda-capture
-                        // warning
-            assign(ptr, i, strides);
-        });
+        hogpp::cartesianProduct(
+            sizes, [this, ptr, strides](const auto& i) constexpr {
+                (void)this; // Avoid (incorrect) Clang -Wunused-lambda-capture
+                            // warning
+                assign(ptr, i, strides);
+            });
     }
 
     template<std::size_t... Indices>
diff --git a/python/type_caster_test.cpp b/python/type_caster_test.cpp
index b327c59..a2c7323 100644
--- a/python/type_caster_test.cpp
+++ b/python/type_caster_test.cpp
@@ -17,6 +17,9 @@
 // limitations under the License.
 //
 
+// FIXME GCC 14.x workaround for https://github.com/pybind/pybind11/pull/5208
+#include <algorithm>
+
 #include <pybind11/cast.h>
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
diff --git a/tests/cpp/test_integral_histogram.cpp b/tests/cpp/test_integral_histogram.cpp
new file mode 100644
index 0000000..45eb3fd
--- /dev/null
+++ b/tests/cpp/test_integral_histogram.cpp
@@ -0,0 +1,111 @@
+//
+// HOGpp - Fast histogram of oriented gradients computation using integral
+// histograms
+//
+// Copyright 2021 Sergiu Deitsch <sergiu.deitsch@gmail.com>
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//   http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+//
+
+#define BOOST_TEST_MODULE hogpp
+
+#include <Eigen/Core>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <random>
+#include <tuple>
+#include <utility>
+
+#include <boost/histogram.hpp>
+#include <boost/test/included/unit_test.hpp>
+
+#include <hogpp/integralhistogram.hpp>
+
+namespace bh = boost::histogram;
+using Axis = bh::axis::integer<int, bh::use_default, bh::axis::option::none_t>;
+
+struct RandomImage
+{
+    static constexpr Eigen::DenseIndex M = 113;
+    static constexpr Eigen::DenseIndex N = 152;
+    static constexpr Eigen::DenseIndex y1 = 45;
+    static constexpr Eigen::DenseIndex x1 = 35;
+    static constexpr Eigen::DenseIndex Rows = 25;
+    static constexpr Eigen::DenseIndex Cols = 30;
+    static constexpr Eigen::DenseIndex y2 = y1 + Rows;
+    static constexpr Eigen::DenseIndex x2 = x1 + Cols;
+    static constexpr int Bins = std::numeric_limits<std::uint8_t>::max() + 1;
+
+    RandomImage()
+    {
+        std::generate_n(image.data(), image.size(), [&x = x, &g = g] {
+            return static_cast<std::uint8_t>(
+                static_cast<int>(std::round(x(g) * 127)) + 128);
+        });
+
+        auto h1 = bh::make_histogram(Axis{0, Bins});
+        std::for_each(image.data(), image.data() + image.size(), std::ref(h1));
+
+        subimage = image.block(y1, x1, Rows, Cols);
+
+        auto h2 = bh::make_histogram(Axis{0, Bins});
+        std::for_each(subimage.data(), subimage.data() + subimage.size(),
+                      std::ref(h2));
+
+        ih.resize(std::make_tuple(image.rows(), image.cols()), Bins);
+
+        ih.scan([&image = image](Eigen::TensorRef<Eigen::Tensor<int, 1>> h,
+                                 const std::tuple<int, int>& i) {
+            std::uint8_t value = std::apply(std::ref(image), i);
+            ++h.coeffRef(value);
+        });
+
+        Eigen::Tensor<int, 1> ih1 =
+            ih.intersect(std::make_tuple(0, 0), std::make_tuple(M, N));
+        Eigen::Tensor<int, 1> ih2 =
+            ih.intersect(std::make_tuple(y1, x1), std::make_tuple(y2, x2));
+
+        std::copy(h1.cbegin(), h1.cend(), ref_hist1.begin());
+        std::copy(h2.cbegin(), h2.cend(), ref_hist2.begin());
+
+        std::copy(ih1.data(), ih1.data() + ih1.size(), computed_hist1.begin());
+        std::copy(ih2.data(), ih2.data() + ih2.size(), computed_hist2.begin());
+    }
+
+    std::normal_distribution<double> x;
+    std::default_random_engine g;
+    Eigen::Matrix<std::uint8_t, M, N> image;
+    Eigen::Matrix<std::uint8_t, Rows, Cols> subimage;
+
+    hogpp::IntegralHistogram<int, 2, 1> ih;
+
+    std::array<int, Bins> ref_hist1;
+    std::array<int, Bins> ref_hist2;
+
+    std::array<int, Bins> computed_hist1;
+    std::array<int, Bins> computed_hist2;
+};
+
+BOOST_FIXTURE_TEST_CASE(full, RandomImage)
+{
+    BOOST_TEST(ref_hist1 == computed_hist1, boost::test_tools::per_element());
+}
+
+BOOST_FIXTURE_TEST_CASE(sub, RandomImage)
+{
+    BOOST_TEST(ref_hist2 == computed_hist2, boost::test_tools::per_element());
+}