diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp
index 527474f..5ced7e7 100644
--- a/include/aare/ClusterFinder.hpp
+++ b/include/aare/ClusterFinder.hpp
@@ -37,7 +37,7 @@ class ClusterFinder {
     NDArray<PEDESTAL_TYPE, 2> m_threshold;
     NDArray<PEDESTAL_TYPE, 2> m_pd_corrected_frame;
 
-    double all=0;
+    double all = 0;
 
   public:
     /**
diff --git a/include/aare/ClusterFinderMT.hpp b/include/aare/ClusterFinderMT.hpp
index be92e1a..cf9f6a9 100644
--- a/include/aare/ClusterFinderMT.hpp
+++ b/include/aare/ClusterFinderMT.hpp
@@ -328,7 +328,8 @@ class ClusterFinderMT {
         (*m_frame_pools[tid])[ref.slot].copy_from(frame);
         // auto [wall_ns, cpu_ns] = dt.elapsed();
         // std::cerr << "ClusterFinderMT: find_clusters: copied frame "
-        //           << frame_number << " took " << wall_ns/1000.0 << " wall_us and "
+        //           << frame_number << " took " << wall_ns/1000.0 << " wall_us
+        //           and "
         //           << cpu_ns/1000.0 << " cpu_us" << std::endl;
 
         // Cannot fail, the free list is what limits how many frames are in
diff --git a/include/aare/FastPedestal.hpp b/include/aare/FastPedestal.hpp
index a6c200c..32da6cb 100644
--- a/include/aare/FastPedestal.hpp
+++ b/include/aare/FastPedestal.hpp
@@ -1,5 +1,6 @@
 // SPDX-License-Identifier: MPL-2.0
 #pragma once
+#include "aare/File.hpp"
 #include "aare/Frame.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/NDView.hpp"
@@ -52,7 +53,12 @@ template <typename PEDESTAL_TYPE> class FastPedestal {
         : m_rows(rows), m_cols(cols), m_samples(n_samples),
           m_inv_samples(1.0 / n_samples), m_sum({rows, cols}, Entry{0, 0}),
           m_mean({rows, cols}, PEDESTAL_TYPE(0)) {
-        assert(rows > 0 && cols > 0 && n_samples > 0);
+        if (!(rows > 0 && cols > 0 && n_samples > 0)) {
+            throw std::runtime_error(
+                fmt::format("Invalid parameters for FastPedestal: rows={}, "
+                            "cols={}, n_samples={} need to be positive",
+                            rows, cols, n_samples));
+        }
     }
     ~FastPedestal() = default;
 
@@ -79,9 +85,9 @@ template <typename PEDESTAL_TYPE> class FastPedestal {
     }
 
     PEDESTAL_TYPE variance(ssize_t index) const {
-        auto &entry = m_sum[index];
-        auto m2 = entry.sum * m_inv_samples * entry.sum * m_inv_samples;
-        return entry.sum2 * m_inv_samples - m2;
+        const auto &entry = m_sum[index];
+        const auto m = entry.sum * m_inv_samples;
+        return std::fma(-m, m, entry.sum2 * m_inv_samples);
     }
 
     NDArray<PEDESTAL_TYPE, 2> std() {
@@ -110,18 +116,88 @@ template <typename PEDESTAL_TYPE> class FastPedestal {
         m_ready = false;
     }
 
+    /**
+     * @brief Update the pedestal with the values of the frame. The weight of
+     the update depends on the number of samples. Bounds checks are performed on
+     the frame shape.
+     * @param frame The frame to update the pedestal with.
+     */
     template <typename T> void push(NDView<T, 2> frame) {
         if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
             throw std::runtime_error(
                 "Frame shape does not match pedestal shape");
         }
 
+        if (!ready()) {
+            throw std::runtime_error("Pedestal is not ready, cannot push");
+        }
+
+        // TODO! update with push_fast
         for (size_t row = 0; row < m_rows; row++) {
             for (size_t col = 0; col < m_cols; col++) {
                 push<T>(row, col, frame(row, col));
             }
         }
     }
+
+    /**
+     * @brief Overload for Frame. Bounds checks are performed on the frame shape
+     * in the view.
+     * @param frame The Frame to update the pedestal with.
+     */
+    template <typename T> void push(Frame &frame) { push<T>(frame.view<T>()); }
+
+    /**
+     * @brief Update one pixel using its row and column indices. Checks if the
+     * pedestal is ready.
+     * @param row The row index of the pixel to update.
+     * @param col The column index of the pixel to update.
+     * @param val_ The value of the pixel to update the pedestal with.
+     */
+    template <typename T>
+    void push(const uint32_t row, const uint32_t col, const T val) {
+        if (!ready()) {
+            throw std::runtime_error("Pedestal is not ready, cannot push");
+        }
+
+        push_fast(rc_to_index(row, col), val);
+    }
+
+    /**
+     * @brief Update one pixel using its flat index.
+     *
+     * WARNING: This steady-state fast path assumes the pedestal is ready and
+     * the index is valid. Assertions check those preconditions in debug builds.
+     */
+    template <typename T>
+    void push_fast(const std::size_t index, const T value) noexcept {
+        assert(m_ready);
+        assert(index < static_cast<std::size_t>(m_sum.size()));
+
+        const auto val = static_cast<double>(value);
+        auto &entry = m_sum[index];
+        entry.sum += val - entry.sum * m_inv_samples;
+        entry.sum2 += val * val - entry.sum2 * m_inv_samples;
+        m_mean[index] = static_cast<PEDESTAL_TYPE>(entry.sum * m_inv_samples);
+    }
+
+    // TODO! Do we need fast variants of push_no_update?
+    template <typename T>
+    void push_no_update(const uint32_t row, const uint32_t col, const T val_) {
+        if (!ready()) {
+            throw std::runtime_error("Pedestal is not ready, cannot push");
+        }
+        auto val = static_cast<double>(val_);
+        auto &entry = m_sum(row, col);
+        entry.sum += val - entry.sum * m_inv_samples;
+        entry.sum2 += val * val - entry.sum2 * m_inv_samples;
+    }
+
+    /**
+     * @brief Push a frame to the pedestal to initialize it. We need to push
+     * n_samples frames to be ready.
+     * @param frame The frame to update the pedestal with.
+     */
     template <typename T> void push_init(NDView<T, 2> frame) {
         if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
             throw std::runtime_error(
@@ -149,10 +225,44 @@ template <typename PEDESTAL_TYPE> class FastPedestal {
         }
     }
 
-    template <typename T> void push(Frame &frame) {
-        assert(frame.rows() == static_cast<size_t>(m_rows) &&
-               frame.cols() == static_cast<size_t>(m_cols));
-        push<T>(frame.view<T>());
+    /**
+     * @brief Create a FastPedestal initialized from the first n_samples frames
+     * of a file. Image size is taken from the file.
+     */
+    template <typename T>
+    static FastPedestal from_file(const std::filesystem::path &filename,
+                                  uint32_t n_samples = 1000,
+                                  uint32_t skip_first = 0) {
+        File f(filename);
+
+        if ((f.total_frames() - skip_first) < n_samples) {
+            throw std::runtime_error(
+                "File has less frames than the number of samples needed to "
+                "initialize the pedestal");
+        }
+
+        if (skip_first > 0) {
+            f.seek(static_cast<size_t>(skip_first));
+        }
+        const auto rows = static_cast<uint32_t>(f.rows());
+        const auto cols = static_cast<uint32_t>(f.cols());
+        FastPedestal pedestal(rows, cols, n_samples);
+        NDArray<T, 2> frame({rows, cols});
+
+        uint32_t frame_index = skip_first;
+        while (frame_index < skip_first + n_samples) {
+            f.read_into(frame.buffer());
+            pedestal.template push_init<T>(frame.view());
+            frame_index++;
+        }
+
+        // read the rest of the file
+        while (frame_index < f.total_frames()) {
+            f.read_into(frame.buffer());
+            pedestal.template push<T>(frame.view());
+            frame_index++;
+        }
+        return pedestal;
     }
 
     // getter functions
@@ -160,44 +270,6 @@ template <typename PEDESTAL_TYPE> class FastPedestal {
     uint32_t cols() const { return m_cols; }
     uint32_t n_samples() const { return m_samples; }
 
-    /**
-     * @brief Update one pixel using its flat index.
-     *
-     * This steady-state fast path assumes the pedestal is ready and the index
-     * is valid. Assertions check those preconditions in debug builds.
-     */
-    template <typename T>
-    void push_fast(const std::size_t index, const T value) noexcept {
-        assert(m_ready);
-        assert(index < static_cast<std::size_t>(m_sum.size()));
-
-        const auto val = static_cast<double>(value);
-        auto &entry = m_sum[index];
-        entry.sum += val - entry.sum * m_inv_samples;
-        entry.sum2 += val * val - entry.sum2 * m_inv_samples;
-        m_mean[index] = static_cast<PEDESTAL_TYPE>(entry.sum * m_inv_samples);
-    }
-
-    template <typename T>
-    void push(const uint32_t row, const uint32_t col, const T val_) {
-        if (!ready()) {
-            throw std::runtime_error("Pedestal is not ready, cannot push");
-        }
-        const auto index = (static_cast<std::size_t>(row) * m_cols) + col;
-        push_fast(index, val_);
-    }
-
-    template <typename T>
-    void push_no_update(const uint32_t row, const uint32_t col, const T val_) {
-        if (!ready()) {
-            throw std::runtime_error("Pedestal is not ready, cannot push");
-        }
-        auto val = static_cast<double>(val_);
-        auto &entry = m_sum(row, col);
-        entry.sum += val - entry.sum * m_inv_samples;
-        entry.sum2 += val * val - entry.sum2 * m_inv_samples;
-    }
-
     /**
      * @brief Update the mean of the pedestal. This is used after having done
      * push_no_update. It is not necessary to call this function after push.
diff --git a/python/src/bind_FastPedestal.hpp b/python/src/bind_FastPedestal.hpp
new file mode 100644
index 0000000..a701d8d
--- /dev/null
+++ b/python/src/bind_FastPedestal.hpp
@@ -0,0 +1,131 @@
+// SPDX-License-Identifier: MPL-2.0
+
+#include "aare/FastPedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+template <typename SUM_TYPE>
+void define_fast_pedestal_bindings(py::module &m, const std::string &name) {
+
+    py::class_<FastPedestal<SUM_TYPE>>(m, name.c_str(), py::buffer_protocol())
+        .def(py::init<uint32_t, uint32_t, uint32_t>())
+        .def(py::init<uint32_t, uint32_t>())
+
+        .def(
+            "mean",
+            [](FastPedestal<SUM_TYPE> &self) {
+                auto mean = new NDArray<SUM_TYPE, 2>{};
+                *mean = self.mean();
+                return return_image_data(mean);
+            },
+            "Return a copy of the mean of the pedestal as a NumPy array")
+        .def(
+            "var",
+            [](FastPedestal<SUM_TYPE> &self) {
+                auto variance = new NDArray<SUM_TYPE, 2>{};
+                *variance = self.variance();
+                return return_image_data(variance);
+            },
+            "Return a copy of the variance of the pedestal as a NumPy array")
+        .def(
+            "std",
+            [](FastPedestal<SUM_TYPE> &self) {
+                auto standard_deviation = new NDArray<SUM_TYPE, 2>{};
+                *standard_deviation = self.std();
+                return return_image_data(standard_deviation);
+            },
+            "Return a copy of the standard deviation of the pedestal as a "
+            "NumPy array")
+        .def(
+            "view",
+            [](py::object self_py) {
+                return py::module_::import("numpy").attr("asarray")(self_py);
+            },
+            "Return non-owning, non-writable view of the pedestal as a NumPy "
+            "array")
+
+        // We need to buffer protocol to allow for numpy operations using the
+        // pedestal mean
+        .def_buffer([](FastPedestal<SUM_TYPE> &self) {
+            auto mean = self.view();
+            return py::buffer_info(
+                const_cast<SUM_TYPE *>(mean.data()), sizeof(SUM_TYPE),
+                py::format_descriptor<SUM_TYPE>::format(), 2,
+                {static_cast<py::ssize_t>(mean.shape(0)),
+                 static_cast<py::ssize_t>(mean.shape(1))},
+                {static_cast<py::ssize_t>(mean.strides()[0] * sizeof(SUM_TYPE)),
+                 static_cast<py::ssize_t>(mean.strides()[1] *
+                                          sizeof(SUM_TYPE))},
+                true);
+        })
+        // Subtracting a FastPedestal from a NumPy array
+        .def(
+            "__array_ufunc__",
+            [](py::object self, py::object ufunc, const std::string &method,
+               py::args inputs, py::kwargs kwargs) -> py::object {
+                if (method != "__call__" || inputs.size() != 2 ||
+                    inputs[1].ptr() != self.ptr() ||
+                    py::cast<std::string>(ufunc.attr("__name__")) !=
+                        "subtract") {
+                    return py::reinterpret_borrow<py::object>(
+                        Py_NotImplemented);
+                }
+
+                auto mean =
+                    py::module_::import("builtins").attr("memoryview")(self);
+                return ufunc(inputs[0], mean, **kwargs);
+            },
+            "Support subtracting a FastPedestal from a NumPy array.")
+        .def("clear", py::overload_cast<>(&FastPedestal<SUM_TYPE>::clear))
+        .def_property_readonly("rows", &FastPedestal<SUM_TYPE>::rows)
+        .def_property_readonly("cols", &FastPedestal<SUM_TYPE>::cols)
+        .def_property_readonly("cur_samples",
+                               &FastPedestal<SUM_TYPE>::cur_samples,
+                               "Return the number of samples pushed if < "
+                               "n_samples in the pedestal")
+        .def_property_readonly(
+            "ready", &FastPedestal<SUM_TYPE>::ready,
+            "Return true if the pedestal is ready to be used (i.e. we have "
+            "pushed at least n_samples samples)")
+        .def_property_readonly(
+            "n_samples", &FastPedestal<SUM_TYPE>::n_samples,
+            "Return the number of samples to push to the pedestal to be ready")
+        .def("clone",
+             [](FastPedestal<SUM_TYPE> &pedestal) {
+                 return FastPedestal<SUM_TYPE>(pedestal);
+             })
+        .def(
+            "push",
+            [](FastPedestal<SUM_TYPE> &pedestal,
+               py::array_t<uint16_t, py::array::c_style> &frame) {
+                pedestal.push(make_view_2d(frame));
+            },
+            py::arg("frame").noconvert())
+        .def(
+            "push_init",
+            [](FastPedestal<SUM_TYPE> &pedestal,
+               py::array_t<uint16_t, py::array::c_style> &frame) {
+                pedestal.push_init(make_view_2d(frame));
+            },
+            py::arg("frame").noconvert(),
+            "Push a frame to the pedestal to initialize it. Needs to be called "
+            "n_samples times to be ready")
+        .def_static(
+            "from_file",
+            [](const std::filesystem::path &filename, uint32_t n_samples,
+               uint32_t skip_first) {
+                return FastPedestal<SUM_TYPE>::template from_file<uint16_t>(
+                    filename, n_samples, skip_first);
+            },
+            py::arg("filename"), py::arg("n_samples") = 1000,
+            py::arg("skip_first") = 0,
+            "Create a FastPedestal from a file. Uses n_samples frames for "
+            "initialization (after skip_first), then pushes any remaining "
+            "frames in steady state.");
+}
diff --git a/python/src/fast_pedestal.hpp b/python/src/fast_pedestal.hpp
deleted file mode 100644
index 7596455..0000000
--- a/python/src/fast_pedestal.hpp
+++ /dev/null
@@ -1,111 +0,0 @@
-// SPDX-License-Identifier: MPL-2.0
-
-#include "aare/FastPedestal.hpp"
-#include "np_helper.hpp"
-
-#include <cstdint>
-#include <pybind11/numpy.h>
-#include <pybind11/pybind11.h>
-#include <pybind11/stl.h>
-
-namespace py = pybind11;
-
-template <typename SUM_TYPE>
-void define_fast_pedestal_bindings(py::module &m, const std::string &name) {
-
-    py::class_<FastPedestal<SUM_TYPE>>(m, name.c_str(), py::buffer_protocol())
-        .def(py::init<int, int, int>())
-        .def(py::init<int, int>())
-        .def("mean",
-             [](FastPedestal<SUM_TYPE> &self) {
-                 auto mean = new NDArray<SUM_TYPE, 2>{};
-                 *mean = self.mean();
-                 return return_image_data(mean);
-             })
-        .def("view",
-             [](py::object self_py) {
-                 auto &self = self_py.cast<FastPedestal<SUM_TYPE> &>();
-                 auto view = self.view();
-                 std::array<py::ssize_t, 2> shape{
-                     static_cast<py::ssize_t>(view.shape(0)),
-                     static_cast<py::ssize_t>(view.shape(1))};
-                 std::array<py::ssize_t, 2> byte_strides{
-                     static_cast<py::ssize_t>(view.strides()[0]) *
-                         static_cast<py::ssize_t>(sizeof(SUM_TYPE)),
-                     static_cast<py::ssize_t>(view.strides()[1]) *
-                         static_cast<py::ssize_t>(sizeof(SUM_TYPE))};
-                 auto array = py::array_t<SUM_TYPE>(shape, byte_strides,
-                                                    view.data(), self_py);
-                 array.attr("setflags")(py::arg("write") = false);
-                 return array;
-             })
-        .def("variance",
-             [](FastPedestal<SUM_TYPE> &self) {
-                 auto variance = new NDArray<SUM_TYPE, 2>{};
-                 *variance = self.variance();
-                 return return_image_data(variance);
-             })
-        .def("std",
-             [](FastPedestal<SUM_TYPE> &self) {
-                 auto standard_deviation = new NDArray<SUM_TYPE, 2>{};
-                 *standard_deviation = self.std();
-                 return return_image_data(standard_deviation);
-             })
-        .def(
-            "__array_ufunc__",
-            [](py::object self, py::object ufunc, const std::string &method,
-               py::args inputs, py::kwargs kwargs) -> py::object {
-                if (method != "__call__" || inputs.size() != 2 ||
-                    inputs[1].ptr() != self.ptr() ||
-                    py::cast<std::string>(ufunc.attr("__name__")) !=
-                        "subtract") {
-                    return py::reinterpret_borrow<py::object>(
-                        Py_NotImplemented);
-                }
-
-                auto mean =
-                    py::module_::import("builtins").attr("memoryview")(self);
-                return ufunc(inputs[0], mean, **kwargs);
-            },
-            "Support subtracting a FastPedestal from a NumPy array.")
-        .def("clear", py::overload_cast<>(&FastPedestal<SUM_TYPE>::clear))
-        .def_property_readonly("rows", &FastPedestal<SUM_TYPE>::rows)
-        .def_property_readonly("cols", &FastPedestal<SUM_TYPE>::cols)
-        .def_property_readonly("cur_samples",
-                               &FastPedestal<SUM_TYPE>::cur_samples)
-        .def_property_readonly("ready", &FastPedestal<SUM_TYPE>::ready)
-        .def_property_readonly("n_samples", &FastPedestal<SUM_TYPE>::n_samples)
-        // .def_property_readonly("sum", &FastPedestal<SUM_TYPE>::get_sum)
-        .def("clone",
-             [](FastPedestal<SUM_TYPE> &pedestal) {
-                 return FastPedestal<SUM_TYPE>(pedestal);
-             })
-        .def(
-            "push",
-            [](FastPedestal<SUM_TYPE> &pedestal,
-               py::array_t<uint16_t, py::array::c_style> &frame) {
-                pedestal.push(make_view_2d(frame));
-            },
-            py::arg("frame").noconvert())
-        .def(
-            "push_init",
-            [](FastPedestal<SUM_TYPE> &pedestal,
-               py::array_t<uint16_t, py::array::c_style> &frame) {
-                pedestal.push_init(make_view_2d(frame));
-            },
-            py::arg("frame").noconvert())
-
-        .def("update_mean", &FastPedestal<SUM_TYPE>::update_mean)
-        .def_buffer([](FastPedestal<SUM_TYPE> &self) {
-            auto mean = self.view();
-            return py::buffer_info(
-                const_cast<SUM_TYPE *>(mean.data()), sizeof(SUM_TYPE),
-                py::format_descriptor<SUM_TYPE>::format(), 2,
-                {static_cast<py::ssize_t>(mean.shape(0)),
-                 static_cast<py::ssize_t>(mean.shape(1))},
-                {static_cast<py::ssize_t>(mean.strides()[0] * sizeof(SUM_TYPE)),
-                 static_cast<py::ssize_t>(mean.strides()[1] *
-                                          sizeof(SUM_TYPE))},
-                true);
-        });
-}
diff --git a/python/src/module.cpp b/python/src/module.cpp
index b668ba8..bcd5b7b 100644
--- a/python/src/module.cpp
+++ b/python/src/module.cpp
@@ -13,6 +13,7 @@
 #include "bind_ClusterVector.hpp"
 #include "bind_Defs.hpp"
 #include "bind_Eta.hpp"
+#include "bind_FastPedestal.hpp"
 #include "bind_Interpolator.hpp"
 #include "bind_PedestalTrackingPixelHistogram.hpp"
 #include "bind_PixelHistogram.hpp"
@@ -22,7 +23,6 @@
 
 // TODO! migrate the other names
 #include "ctb_raw_file.hpp"
-#include "fast_pedestal.hpp"
 #include "file.hpp"
 #include "fit.hpp"
 #include "jungfrau_data_file.hpp"