diff --git a/cmake/fetch-fftw.cmake b/cmake/fetch-fftw.cmake
index 226be675..5148f75c 100644
--- a/cmake/fetch-fftw.cmake
+++ b/cmake/fetch-fftw.cmake
@@ -42,7 +42,7 @@ function(fetch_fftw)
       "-DENABLE_THREADS:BOOL=ON" # Use pthread
       "-DWITH_COMBINED_THREADS:BOOL=ON" # Don't need to link in fftw3f_threads
   
-      "-DENABLE_FLOAT:BOOL=ON" # <float>
+      # "-DENABLE_FLOAT:BOOL=ON" # <float>
   
       # Use SSE, but not AVX.
       "-DENABLE_SSE:BOOL=ON" 
@@ -102,7 +102,7 @@ function(fetch_fftw)
     )
 
     unset(generator)
-    
+
     if(result)
       file(READ ${fftw_BINARY_DIR}/build_output.log build_log)
       message(FATAL_ERROR "Result = ${result}\nFailed FFTW-${args_VERSION} build, see build log:\n"
@@ -116,9 +116,9 @@ function(fetch_fftw)
   # Confirm that we can find FFTW.
 
   # Ugly work-around for CMake errors in CI builds.
-  set(_cmake_import_check_xcframework_for_FFTW3::fftw3f "")
+  set(_cmake_import_check_xcframework_for_FFTW3::fftw3 "")
 
-  find_package(FFTW3f
+  find_package(FFTW3
     QUIET
     REQUIRED 
     CONFIG
@@ -126,9 +126,9 @@ function(fetch_fftw)
     NO_DEFAULT_PATH 
   )
 
-  unset(_cmake_import_check_xcframework_for_FFTW3::fftw3f)
+  unset(_cmake_import_check_xcframework_for_FFTW3::fftw3)
 
-  if(NOT FFTW3f_FOUND) 
+  if(NOT FFTW3_FOUND) 
     message(FATAL_ERROR "FFTW-${args_VERSION} not found")
   endif()
 endfunction()
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index c615e48e..d413cb9e 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -85,7 +85,7 @@ if (LINUX)
   add_example(
     NAME periodogram
     SRC periodogram.c 
-    DEPS FFTW3::fftw3f nothings::stb)
+    DEPS FFTW3::fftw3 nothings::stb)
 endif()
 
 # ---- End-of-file commands ----
diff --git a/examples/custom_alloc.c b/examples/custom_alloc.c
index 9948dc82..f2e22a08 100644
--- a/examples/custom_alloc.c
+++ b/examples/custom_alloc.c
@@ -8,6 +8,8 @@
 #define TPH_POISSON_IMPLEMENTATION
 #include "thinks/tph_poisson.h"
 
+static_assert(sizeof(tph_poisson_real) == 4, "");
+
 typedef struct my_alloc_ctx_
 {
   ptrdiff_t total_malloc;
diff --git a/examples/periodogram.c b/examples/periodogram.c
index 6ca067ae..0f0b2739 100644
--- a/examples/periodogram.c
+++ b/examples/periodogram.c
@@ -1,328 +1,189 @@
+#include <math.h> /* sqrt, ceil, floor, round */
+#include <string.h> /* memset, memcpy */
 
 #include <fftw3.h>
 
 #define STB_IMAGE_WRITE_IMPLEMENTATION
 #include <stb_image_write.h>
 
-int main(int argc, char *argv[])
-{
-  (void)argc;
-  (void)argv;
-
-  // const int kImageCount = 100;
-  // const int kPixelSize = 2048;
-
-
-  return 0;
-}
+#define TPH_POISSON_REAL_TYPE double
+#define TPH_POISSON_SQRT sqrt
+#define TPH_POISSON_CEIL ceil
+#define TPH_POISSON_FLOOR floor
+#define TPH_POISSON_IMPLEMENTATION
+#include "thinks/tph_poisson.h"
 
-#if 0
-// Copyright(C) Tommy Hinks <tommy.hinks@gmail.com>
-// This file is subject to the license terms in the LICENSE file
-// found in the top-level directory of this distribution.
-
-// clang-format off
-
-#include <algorithm>
-#include <cmath>
-#include <complex>
-#include <cstdlib>
-#include <iostream>
-
-#include "hedley.h"
-#include "thinks/poisson_disk_sampling/poisson_disk_sampling.h"
-
-HEDLEY_DIAGNOSTIC_PUSH
-#if defined(HEDLEY_MSVC_VERSION)
-#define STBI_MSC_SECURE_CRT
-#elif defined(HEDLEY_GCC_VERSION)
-_Pragma("GCC diagnostic ignored \"-Wold-style-cast\"")
-_Pragma("GCC diagnostic ignored \"-Wsign-conversion\"")
-_Pragma("GCC diagnostic ignored \"-Wconversion\"")
-_Pragma("GCC diagnostic ignored \"-Wcast-qual\"")
-_Pragma("GCC diagnostic ignored \"-Wmissing-declarations\"")
-#elif defined(__clang__)
-_Pragma("clang diagnostic ignored \"-Wcast-qual\"")
-_Pragma("clang diagnostic ignored \"-Wmissing-prototypes\"")
-_Pragma("clang diagnostic ignored \"-Wimplicit-fallthrough\"")
-#endif
+static_assert(sizeof(tph_poisson_real) == sizeof(double), "");
 
-#define STB_IMAGE_WRITE_IMPLEMENTATION
-#include "stb_image_write.h"
-HEDLEY_DIAGNOSTIC_POP
-
-HEDLEY_DIAGNOSTIC_PUSH
-#if defined (HEDLEY_GCC_VERSION)
-_Pragma("GCC diagnostic ignored \"-Wconversion\"")
-_Pragma("GCC diagnostic ignored \"-Wsign-conversion\"")
-#endif
-#include "simple_fft/fft.h"
-HEDLEY_DIAGNOSTIC_POP
-
-  HEDLEY_WARN_UNUSED_RESULT static HEDLEY_CONSTEXPR auto Reinterval(const double in_val,
-    const double old_min,
-    const double old_max,
-    const double new_min,
-    const double new_max) noexcept -> double
+static void fft_shift(const int n0, const int n1, double *inout)
 {
-  return (old_max - old_min) == 0.0
-           ? new_max
-           : (((in_val - old_min) * (new_max - new_min)) / (old_max - old_min)) + new_min;
+  /* Simple and slow implementation, not optimized. */
+
+  const int s0 = n0 / 2;
+  const int s1 = n1 / 2;
+
+  /* Shift rows. */
+  double x = 0.0; /* swap */
+  int jj = 0;
+  for (int j = 0; j < n1; ++j) {
+    jj = n0 * j;
+    for (int i = 0; i < s0; ++i) {
+      x = inout[i + jj];
+      inout[i + jj] = inout[i + s0 + jj];
+      inout[i + s0 + jj] = x;
+    }
+  }
+
+  /* Shift columns. */
+  for (int i = 0; i < n0; ++i) {
+    for (int j = 0; j < s1; ++j) {
+      jj = n0 * j;
+      x = inout[i + jj];
+      inout[i + jj] = inout[i + n0 * (s1 + j)];
+      inout[i + n0 * (s1 + j)] = x;
+    }
+  }
 }
 
-template<typename PixelT> class Image
+static void sampling_image(const tph_poisson_args *args,
+  const tph_poisson_sampling *s,
+  const int n0,
+  const int n1,
+  fftw_complex *img)
 {
-public:
-  Image(const std::size_t width, const std::size_t height)
-    : _width(width), _height(height), _pixels(_width * _height, PixelT{ 0 })
-  {}
-
-  // NOLINTNEXTLINE
-  HEDLEY_WARN_UNUSED_RESULT constexpr auto width() const -> std::size_t { return _width; }
-
-  // NOLINTNEXTLINE
-  HEDLEY_WARN_UNUSED_RESULT constexpr auto height() const -> std::size_t { return _height; }
-
-  // NOLINTNEXTLINE
-  HEDLEY_WARN_UNUSED_RESULT auto data() const -> const PixelT * { return _pixels.data(); }
-
-  // NOLINTNEXTLINE
-  HEDLEY_WARN_UNUSED_RESULT auto data() -> PixelT * { return _pixels.data(); }
-
-  HEDLEY_WARN_UNUSED_RESULT
-  auto operator()(const std::size_t i, const std::size_t j) -> PixelT &
-  {
-    return _pixels[_linearIndex(i, j)];
+  /* Reset, start from zero image. Imaginary part will remain zero. */
+  const int sz = n0 * n1;
+  for (int i = 0; i < sz; ++i) {
+    img[i][0] = 0.0;
+    img[i][1] = 0.0;
   }
 
-  HEDLEY_WARN_UNUSED_RESULT
-  auto operator()(const std::size_t i, const std::size_t j) const -> const PixelT &
-  {
-    return _pixels[_linearIndex(i, j)];
-  }
+  const double *p = tph_poisson_get_samples(s);
+  if (p == NULL) { abort(); }
 
-private:
-  HEDLEY_WARN_UNUSED_RESULT
-  // NOLINTNEXTLINE
-  constexpr std::size_t _linearIndex(const std::size_t i, const std::size_t j) const
-  {
-    return i + _width * j;
-  }
+  const ptrdiff_t nsamples = s->nsamples;
+  const ptrdiff_t ndims = s->ndims;
+  const double x_min = args->bounds_min[0];
+  const double y_min = args->bounds_min[1];
+  const double x_max = args->bounds_max[0];
+  const double y_max = args->bounds_max[1];
+  for (ptrdiff_t i = 0; i < nsamples; ++i) {
+    int ix = (int)floor(((p[i * ndims] - x_min) / (x_max - x_min)) * (double)n0);
+    ix = ix < 0 ? 0 : ((n0 - 1) < ix ? (n0 - 1) : ix);
 
-  std::size_t _width;
-  std::size_t _height;
-  std::vector<PixelT> _pixels;
-};
+    int iy = (int)floor(((p[i * ndims + 1] - y_min) / (y_max - y_min)) * (double)n1);
+    iy = iy < 0 ? 0 : ((n1 - 1) < iy ? (n1 - 1) : iy);
 
-static auto AddEq(const Image<double> &rhs, Image<double> *lhs) noexcept -> Image<double> &
-{
-  if (!(rhs.width() == lhs->width() && rhs.height() == lhs->height())) {
-    std::cerr << "AddEq width/height mismatch";
-    std::abort();
+    img[ix + n0 * iy][0] = 1.0;
   }
 
-  const auto w = rhs.width();
-  const auto h = rhs.height();
-  for (auto i = 0U; i < w; ++i) {
-    for (auto j = 0U; j < h; ++j) { (*lhs)(i, j) += rhs(i, j); }
-  }
-  return *lhs;
+  /* Subtract average. */
+  double avg = 0.0;
+  for (int i = 0; i < sz; ++i) { avg += img[i][0]; }
+  avg /= (double)sz;
+  for (int i = 0; i < sz; ++i) { img[i][0] -= avg; }
 }
 
-HEDLEY_WARN_UNUSED_RESULT static auto Scaled(const Image<double> &img,
-  const double scalar) noexcept -> Image<double>
+static void
+  accum_periodogram(const double scale, const int n0, const int n1, fftw_complex *in, double *out)
 {
-  Image<double> out = img;
-  const auto w = out.width();
-  const auto h = out.height();
-  for (auto i = 0U; i < w; ++i) {
-    for (auto j = 0U; j < h; ++j) { out(i, j) *= scalar; }
-  }
-  return out;
+  const int sz = n0 * n1;
+  for (int i = 0; i < sz; ++i) { out[i] += scale * (in[i][0] * in[i][0] + in[i][1] * in[i][1]); }
 }
 
-HEDLEY_WARN_UNUSED_RESULT static auto SubtractAverage(
-  const Image<double> &img) noexcept -> Image<double>
+static bool write_png(const char *filename, const int n0, const int n1, const double *data)
 {
-  auto out = img;
-  const auto w = out.width();
-  const auto h = out.height();
-  double avg = 0.0;
-  for (auto i = 0U; i < w; ++i) {
-    for (auto j = 0U; j < h; ++j) { avg += out(i, j); }
-  }
-  avg /= static_cast<double>(w * h);
+  static const int comp = 1; /* Greyscale. */
 
-  for (auto i = 0U; i < w; ++i) {
-    for (auto j = 0U; j < h; ++j) { out(i, j) -= avg; }
+  const ptrdiff_t sz = (ptrdiff_t)n0 * (ptrdiff_t)n1;
+  double min_val = data[0];
+  double max_val = data[0];
+  for (ptrdiff_t i = 1; i < sz; ++i) {
+    if (data[i] < min_val) { min_val = data[i]; }
+    if (data[i] > max_val) { max_val = data[i]; }
   }
 
-  return out;
-}
+  const size_t buf_size = (size_t)sz * sizeof(uint8_t);
+  uint8_t *buf = (uint8_t *)malloc(buf_size);
+  memset(buf, 0, buf_size);
 
-HEDLEY_WARN_UNUSED_RESULT static auto Fft2d(
-  const Image<double> &img_in) noexcept -> Image<std::complex<double>>
-{
-  // Create a complex image with imaginary components set to zero.
-  auto c_img = Image<std::complex<double>>(img_in.width(), img_in.height());
-  const auto w = c_img.width();
-  const auto h = c_img.height();
-  for (auto i = 0U; i < w; ++i) {
-    for (auto j = 0U; j < h; ++j) { c_img(i, j) = { img_in(i, j), 0.0 }; }
+  for (ptrdiff_t i = 0; i < sz; ++i) {
+    const int iv = (int)round(((data[i] - min_val) / (max_val - min_val)) * 255.0);
+    buf[i] = (uint8_t)(iv < 0 ? 0 : (255 < iv ? 255 : iv));
   }
 
-  const char *error_description = nullptr;
-  if (!simple_fft::FFT(c_img, c_img.width(), c_img.height(), error_description)) {
-    std::cerr << "FFT error: " << error_description;
-    std::abort();
-  }
+  const int ret = stbi_write_png(filename, n0, n1, comp, buf, n0);
+  free(buf);
 
-  return c_img;
+  return ret != 0;
 }
 
-HEDLEY_WARN_UNUSED_RESULT static auto FftShift2d(
-  const Image<std::complex<double>> &fft_img) noexcept -> Image<std::complex<double>>
+int main(int argc, char *argv[])
 {
-  if (!(fft_img.width() % 2 == 0 && fft_img.height() % 2 == 0)) {
-    std::cerr << "FFT shift error: odd dimension";
-    std::abort();
-  }
+  (void)argc;
+  (void)argv;
 
-  auto shift_img = fft_img;
-  const auto w = fft_img.width();
-  const auto h = fft_img.height();
-  const auto half_w = w / 2;
-  const auto half_h = h / 2;
-
-  // Simple and slow implementation, not optimized.
-  //
-  // Shift rows.
-  for (auto j = 0U; j < h; ++j) {
-    for (auto i = 0U; i < half_w; ++i) { std::swap(shift_img(i, j), shift_img(half_w + i, j)); }
-  }
+  /* Periodogram settings. */
+  const int image_count = 100;
+  const int n0 = 2048;
+  const int n1 = 2048;
 
-  // Shift cols.
-  for (auto i = 0U; i < w; ++i) {
-    for (auto j = 0U; j < half_h; ++j) { std::swap(shift_img(i, j), shift_img(i, half_h + j)); }
-  }
+  /* clang-format off */
+  const tph_poisson_real bounds_min[2] = { 
+    (tph_poisson_real)0, (tph_poisson_real)0 };
+  const tph_poisson_real bounds_max[2] = { 
+    (tph_poisson_real)128, (tph_poisson_real)128 };
+  /* clang-format on */
 
-  return shift_img;
-}
+  /* Configure tph_poisson arguments. Set varying seed later. */
+  tph_poisson_args args = { .bounds_min = bounds_min,
+    .bounds_max = bounds_max,
+    .radius = (tph_poisson_real)1,
+    .ndims = INT32_C(2),
+    .max_sample_attempts = UINT32_C(30),
+    .seed = UINT64_C(1981) };
 
-// The periodogram is the squared magnitude of the FFT bins.
-HEDLEY_WARN_UNUSED_RESULT static auto Periodogram(
-  const Image<std::complex<double>> &fft_img) noexcept -> Image<double>
-{
-  auto img = Image<double>(fft_img.width(), fft_img.height());
-  const auto w = img.width();
-  const auto h = img.height();
-  for (auto i = 0U; i < w; ++i) {
-    for (auto j = 0U; j < h; ++j) {
-      const auto p = std::abs(fft_img(i, j));
-      img(i, j) = p * p;
-    }
-  }
-  return img;
-}
+  /* Initialize empty tph_poisson sampling. */
+  tph_poisson_sampling sampling;
+  memset(&sampling, 0, sizeof(tph_poisson_sampling));
 
-HEDLEY_WARN_UNUSED_RESULT static auto CreateSampleImage(const std::size_t width,
-  const std::size_t height,
-  const std::array<double, 2> &sample_min,
-  const std::array<double, 2> &sample_max,
-  const std::vector<std::array<double, 2>> &samples) noexcept -> Image<double>
-{
-  Image<double> img(width, height);
-  for (const auto &sample : samples) {
-    const auto i = static_cast<std::size_t>(
-      Reinterval(sample[0], sample_min[0], sample_max[0], 0.0, static_cast<double>(width - 1)));
-    const auto j = static_cast<std::size_t>(
-      Reinterval(sample[1], sample_min[1], sample_max[1], 0.0, static_cast<double>(height - 1)));
-    img(i, j) = 1.0;
-  }
-  return img;
-}
+  /* Initlialize buffers used to accumulate the average periodogram. */
+  double *periodogram = (double *)malloc((size_t)n0 * (size_t)n1 * sizeof(double));
+  for (int i = 0; i < (n0 * n1); ++i) { periodogram[i] = 0.0; }
 
-HEDLEY_WARN_UNUSED_RESULT static auto To8bits(
-  const Image<double> &img) noexcept -> Image<std::uint8_t>
-{
-  double min_pixel = std::numeric_limits<double>::max();
-  double max_pixel = std::numeric_limits<double>::lowest();
-  const auto w = img.width();
-  const auto h = img.height();
-  for (auto i = 0U; i < w; ++i) {
-    for (auto j = 0U; j < h; ++j) {
-      min_pixel = std::min(min_pixel, img(i, j));
-      max_pixel = std::max(max_pixel, img(i, j));
-    }
-  }
+  /* Initialize FFT buffers and plan. */
+  fftw_complex *in = fftw_alloc_complex((size_t)n0 * (size_t)n1);
+  fftw_complex *out = fftw_alloc_complex((size_t)n0 * (size_t)n1);
+  fftw_plan plan = fftw_plan_dft_2d(n0, n1, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
 
-  constexpr auto kMinMappedPixel = 0.0;
-  constexpr auto kMaxMappedPixel = 255.999;
-  Image<std::uint8_t> out(w, h);
-  for (auto i = 0U; i < w; ++i) {
-    for (auto j = 0U; j < h; ++j) {
-      out(i, j) = static_cast<std::uint8_t>(
-        Reinterval(img(i, j), min_pixel, max_pixel, kMinMappedPixel, kMaxMappedPixel));
-    }
-  }
-  return out;
-}
+  const double scale = 1. / (double)image_count;
+  for (int i = 0; i < image_count; ++i) {
+    /* Vary the seed for each image. */
+    args.seed = (uint64_t)i;
 
-static void WriteImage(const std::string &filename, const Image<double> &img)
-{
-  constexpr auto kComp = 1;// Greyscale.
-
-  const auto w = static_cast<int>(img.width());
-  const auto h = static_cast<int>(img.height());
-
-  // HACK - clang-tidy will analyze this call and find errors in stb.
-  //        If we ignore this line it doesn't seem to analyze the stb code,
-  //        even though it is inlined in this file.
-  //
-  // NOLINTNEXTLINE
-  if (stbi_write_png(filename.c_str(), w, h, kComp, To8bits(img).data(), w) == 0) {
-    std::cerr << "failed writing image";
-    std::abort();
-  }
-}
+    /* Populate sampling with points. Using default allocator (libc malloc). */
+    if (tph_poisson_create(&args, /*alloc=*/NULL, &sampling) != TPH_POISSON_SUCCESS) { abort(); }
 
-int main(int /*argc*/, char * /*argv*/[])
-{// NOLINT
-  try {
-    constexpr auto kImageCount = 100U;
-    constexpr auto kPixelSize = 2048U;
-
-    constexpr auto kRadius = 1.0;
-    constexpr std::array<double, 2> kXMin = { 0.0, 0.0 };
-    constexpr std::array<double, 2> kXMax = { 128.0, 128.0 };
-    constexpr auto kMaxSampleAttempts = std::uint32_t{ 30 };
-
-    Image<double> avg_periodogram_img(kPixelSize, kPixelSize);
-    for (auto i = 0U; i < kImageCount; ++i) {
-      AddEq(Scaled(Periodogram(
-                     FftShift2d(Fft2d(SubtractAverage(CreateSampleImage(avg_periodogram_img.width(),
-                       avg_periodogram_img.height(),
-                       kXMin,
-                       kXMax,
-                       thinks::PoissonDiskSampling(kRadius,
-                         kXMin,
-                         kXMax,
-                         kMaxSampleAttempts,
-                         /* seed */ i)))))),
-              1.0 / kImageCount),
-        &avg_periodogram_img);
-    }
-    WriteImage("avg_periodogram.png", avg_periodogram_img);
-  } catch (std::exception &e) {
-    std::cerr << "std::exception: " << e.what();
-    std::abort();
-  } catch (...) {
-    std::cerr << "unknown exception";
-    std::abort();
+    /* Construct FFT input from sampling. */
+    sampling_image(&args, &sampling, n0, n1, in);
+
+    /* Perform FFT. */
+    fftw_execute(plan);
+
+    /* Accumulate (scaled) results. */
+    accum_periodogram(scale, n0, n1, out, periodogram);
   }
 
+  /* Shift DC bin to the center of the image and write png file. */
+  fft_shift(n0, n1, periodogram);
+  if (!write_png("./tph_poisson_periodogram.png", n0, n1, periodogram)) { abort(); }
+
+  /* Free resources. */
+  fftw_destroy_plan(plan);
+  free(in);
+  free(out);
+  free(periodogram);
+  tph_poisson_destroy(&sampling);
+
   return EXIT_SUCCESS;
 }
-
-// clang-format on
-#endif
diff --git a/include/thinks/tph_poisson.h b/include/thinks/tph_poisson.h
index 8ca50dd2..132c343b 100644
--- a/include/thinks/tph_poisson.h
+++ b/include/thinks/tph_poisson.h
@@ -135,7 +135,7 @@ extern void tph_poisson_destroy(tph_poisson_sampling *sampling);
  * call to the tph_poisson_create function (without being destroy by the tph_poisson_destroy
  * function in between).
  */
-extern const tph_poisson_real *tph_poisson_get_samples(tph_poisson_sampling *sampling);
+extern const tph_poisson_real *tph_poisson_get_samples(const tph_poisson_sampling *sampling);
 
 /* END PUBLIC API ----------------------------------------------------------- */
 
@@ -1151,7 +1151,7 @@ void tph_poisson_destroy(tph_poisson_sampling *sampling)
   }
 }
 
-const tph_poisson_real *tph_poisson_get_samples(tph_poisson_sampling *sampling)
+const tph_poisson_real *tph_poisson_get_samples(const tph_poisson_sampling *sampling)
 {
   /* Make sure that a 'destroyed' sampling does not return any samples. */
   if (sampling != NULL && sampling->internal != NULL) {