implement embedding

.gitignore

@@ -0,0 +1 @@
+/build/

CMakeLists.txt

@@ -0,0 +1,30 @@
+cmake_minimum_required(VERSION 3.18)
+
+project(blocky C CXX)
+
+option(BLOCKY_STATIC OFF)
+
+set(CMAKE_CXX_STANDARD          20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+set(BLOCKY_C_FLAGS
+    $<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>,$<CXX_COMPILER_ID:GNU>>:
+      -Wall -Werror -pedantic-errors -Wextra -Wconversion -Wsign-conversion>
+    $<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>>:
+      -Wno-overloaded-virtual>
+    $<$<CXX_COMPILER_ID:MSVC>:
+    /W4 /WX>
+)
+set(BLOCKY_CXX_FLAGS
+  $<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>,$<CXX_COMPILER_ID:GNU>>:
+    -Wall -Werror -pedantic-errors -Wextra -Wconversion -Wsign-conversion>
+  $<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>>:
+    -Wno-overloaded-virtual>
+  $<$<CXX_COMPILER_ID:MSVC>:
+    /W4 /WX>
+)
+
+add_subdirectory(thirdparty)
+add_subdirectory(liblocky)
+add_subdirectory(blocky)
+add_subdirectory(playground)

README.md

@@ -0,0 +1,23 @@
+blocky
+====
+
+## CMake command
+
+When you built openh264 on `/home/user/openh264`:
+```
+cmake -DCMAKE_CXX_FLAGS=-isystem\ /home/user/openh264/codec/api\ -L/home/user/openh264 ..
+```
+
+## ffmpeg useful commands
+
+```
+ffmpeg -i in.mp4 -vframes 300 "%d.png"
+```
+
+```
+ffmpeg -r 30 -i "%d.png" -vcodec libx264 -pix_fmt yuv420p out.mp4
+```
+
+```
+ffmpeg -i src.mp4 -t 10s cut.mp4
+```

blocky/CMakeLists.txt

@@ -0,0 +1,19 @@
+add_executable(blocky)
+target_compile_options(blocky PRIVATE ${BLOCKY_C_FLAGS})
+
+target_sources(blocky
+  PUBLIC
+    main.cc
+
+    bytes.hh
+    common.hh
+    video_encoder.hh
+    video_decoder.hh
+)
+target_link_libraries(blocky
+  PUBLIC
+    args
+    liblocky
+    minimp4
+    openh264
+)

blocky/bytes.hh

@@ -0,0 +1,75 @@
+#pragma once
+
+#include <cstdint>
+#include <tuple>
+#include <vector>
+
+#include "common.hh"
+
+
+namespace blky {
+
+std::vector<uint32_t> BytesEncoder(
+    const std::vector<uint8_t>&     bytes,
+    std::tuple<uint32_t, uint32_t>& block_num,
+    uint8_t                         feat_bits,
+    uint32_t                        block_idx,
+    uint64_t                        seed) {
+  const uint32_t block_num_all = std::get<0>(block_num)*std::get<1>(block_num);
+  if (block_num_all == 0) throw std::runtime_error {"block num is zero"};
+
+  blky_encoder_t enc = {};
+  enc.block_num   = block_num_all;
+  enc.feat_bits   = feat_bits;
+  enc.block_index = block_idx;
+  enc.seed        = seed;
+
+  std::vector<uint32_t> ret;
+  for (auto c : bytes) {
+    blky_encoder_feed(&enc, c);
+
+    uint32_t feat;
+    while (blky_encoder_pop(&enc, &feat, false)) {
+      ret.push_back(feat);
+    }
+  }
+
+  uint32_t feat;
+  if (blky_encoder_pop(&enc, &feat, true)) {
+    ret.push_back(feat);
+  }
+  return ret;
+}
+
+std::vector<uint8_t> BytesDecoder(
+    const std::vector<uint32_t>&    features,
+    std::tuple<uint32_t, uint32_t>& block_num,
+    uint8_t                         feat_bits,
+    uint64_t                        seed) {
+  const uint32_t block_num_all = std::get<0>(block_num)*std::get<1>(block_num);
+  if (block_num_all == 0) throw std::runtime_error {"block num is zero"};
+
+  blky_decoder_t de = {};
+  de.block_num   = block_num_all;
+  de.feat_bits   = feat_bits;
+  de.seed        = seed;
+
+  std::vector<uint8_t> ret;
+  for (auto c : features) {
+    if (!blky_decoder_feed(&de, c)) {
+      throw std::runtime_error {"path corruption"};
+    }
+
+    uint8_t feat;
+    while (blky_decoder_pop(&de, &feat, false)) {
+      ret.push_back(feat);
+    }
+  }
+  uint8_t feat;
+  if (blky_decoder_pop(&de, &feat, true)) {
+    ret.push_back(feat);
+  }
+  return ret;
+}
+
+}  // namespace blky

blocky/common.hh

@@ -0,0 +1,59 @@
+#pragma once
+
+extern "C" {
+# include <liblocky.h>
+}
+
+#include <cctype>
+#include <cstdint>
+#include <iostream>
+#include <tuple>
+#include <unordered_map>
+#include <vector>
+
+
+namespace blky {
+
+enum DataFlow {
+  kBytes,
+  kFeatures,
+  kFeatureProbs,
+  kVideo,
+};
+
+static inline const std::string kDataFlowList = "bytes/features/video";
+static inline const std::unordered_map<std::string, DataFlow> kDataFlowMap = {
+  {"bytes",         kBytes},
+  {"features",      kFeatures},
+  {"feature-probs", kFeatureProbs},
+  {"video",         kVideo},
+};
+
+
+template <typename T>
+std::vector<T> ReadAll(auto& ist) noexcept {
+  std::vector<T> ret;
+  for (;;) {
+    T v;
+    ist >> v;
+    if (ist.eof()) return ret;
+    ret.push_back(v);
+  }
+}
+
+static uint8_t ToHex(char c) {
+  if (!std::isxdigit(c)) throw std::runtime_error {"not xdigit"};
+  return static_cast<uint8_t>(std::isalpha(c)? std::tolower(c)-'a'+0xA: c-'0');
+}
+
+}  // namespace blky
+
+
+namespace args {
+
+auto& operator>>(auto& ist, std::tuple<uint32_t, uint32_t>& v) {
+  ist >> std::get<0>(v) >> std::get<1>(v);
+  return ist;
+}
+
+}  // namespace args

blocky/embed.hh

@@ -0,0 +1,107 @@
+#pragma once
+
+#include <cassert>
+#include <cstring>
+#include <optional>
+
+#include "video_decoder.hh"
+#include "video_encoder.hh"
+
+
+void Embed(const std::vector<uint32_t>& features,
+           const std::string& dst_path,
+           const std::string& video_path,
+           const std::tuple<uint32_t, uint32_t>& div,
+           uint32_t                              utime) {
+  assert(features.size());
+
+  const uint32_t div_x = std::get<0>(div);
+  const uint32_t div_y = std::get<1>(div);
+
+  VideoDecoder dec {video_path};
+  std::optional<VideoEncoder> enc;
+
+  std::vector<uint8_t> Y;
+  std::vector<uint8_t> U;
+  std::vector<uint8_t> V;
+
+  uint32_t last_feat = UINT32_MAX;
+  std::vector<uint8_t> feat_pix;
+  for (uint32_t time = 0; dec.Decode();) {
+    const auto& src = dec.frame();
+    if (src.iBufferStatus != 1) continue;
+    ++time;
+
+    const uint32_t w = static_cast<uint32_t>(src.UsrData.sSystemBuffer.iWidth);
+    const uint32_t h = static_cast<uint32_t>(src.UsrData.sSystemBuffer.iHeight);
+
+    const uint32_t stride_y  = static_cast<uint32_t>(src.UsrData.sSystemBuffer.iStride[0]);
+    const uint32_t stride_uv = static_cast<uint32_t>(src.UsrData.sSystemBuffer.iStride[1]);
+
+    const uint8_t* const* srcp = src.pDst;
+
+    // copy buffer to modify
+    Y.resize(w*h);
+    U.resize(w*h/2/2);
+    V.resize(w*h/2/2);
+    for (uint32_t y = 0; y < h; ++y) {
+      std::memcpy(Y.data()+y*w, srcp[0]+stride_y*y, w);
+    }
+    for (uint32_t y = 0; y < h/2; ++y) {
+      std::memcpy(U.data()+y*(w/2), srcp[1]+stride_uv*y, w/2);
+      std::memcpy(V.data()+y*(w/2), srcp[2]+stride_uv*y, w/2);
+    }
+
+    // embed a feature to the buffer
+    const uint32_t feat = features[(time/utime)%features.size()];
+
+    const uint32_t feat_x = feat%div_x;
+    const uint32_t feat_y = feat/div_x;
+
+    const uint32_t feat_size_x = w/div_x;
+    const uint32_t feat_size_y = h/div_y;
+
+    const uint32_t feat_offset_x = feat_x*feat_size_x;
+    const uint32_t feat_offset_y = feat_y*feat_size_y;
+
+    if (feat != last_feat) {
+      feat_pix.resize(feat_size_x*feat_size_y);
+      for (uint32_t y = 0; y < feat_size_y; ++y) {
+        const uint32_t ay = y+feat_offset_y;
+        std::memcpy(feat_pix.data()+(y/2)*(feat_size_x/2), U.data()+(ay/2)*(w/2)+feat_offset_x/2, feat_size_x/2);
+        std::memcpy(feat_pix.data()+(y/2)*(feat_size_x/2)+feat_pix.size()/2, V.data()+(ay/2)*(w/2)+feat_offset_x/2, feat_size_x/2);
+      }
+      last_feat = feat;
+    }
+    for (uint32_t y = 0; y < feat_size_y; ++y) {
+      const uint32_t ay = y+feat_offset_y;
+      std::memcpy(U.data()+(ay/2)*(w/2)+feat_offset_x/2, feat_pix.data()+(y/2)*(feat_size_x/2), feat_size_x/2);
+      std::memcpy(V.data()+(ay/2)*(w/2)+feat_offset_x/2, feat_pix.data()+(y/2)*(feat_size_x/2)+feat_pix.size()/2, feat_size_x/2);
+    }
+
+    // create an encoder if not yet
+    if (!enc) {
+      SEncParamBase param = {};
+      param.iUsageType     = SCREEN_CONTENT_REAL_TIME;
+      param.iPicWidth      = static_cast<int>(w);
+      param.iPicHeight     = static_cast<int>(h);
+      param.fMaxFrameRate  = 30;
+      param.iTargetBitrate = 5000000;
+      enc.emplace(dst_path, param);
+    }
+
+    // encode
+    SSourcePicture dst = {};
+    dst.iColorFormat = videoFormatI420;
+    dst.pData[0]     = Y.data();
+    dst.pData[1]     = U.data();
+    dst.pData[2]     = V.data();
+    dst.iStride[0]   = static_cast<int>(w);
+    dst.iStride[1]   = static_cast<int>(w/2);
+    dst.iStride[2]   = static_cast<int>(w/2);
+    dst.iPicWidth    = static_cast<int>(w);
+    dst.iPicHeight   = static_cast<int>(h);
+    dst.uiTimeStamp  = static_cast<int64_t>(src.uiOutYuvTimeStamp);
+    enc->Encode(dst);
+  }
+}

blocky/features.hh

@@ -0,0 +1,105 @@
+#pragma once
+
+#include <algorithm>
+#include <iostream>
+#include <tuple>
+
+#include "common.hh"
+
+
+namespace blky {
+
+std::vector<uint32_t> PathfindFeatures(
+    const std::vector<double>&      probs,
+    std::tuple<uint32_t, uint32_t>& block_num,
+    uint8_t                         feat_bits,
+    uint64_t                        seed) {
+  const uint32_t num = std::get<0>(block_num) * std::get<1>(block_num);
+  if (num == 0) throw std::runtime_error {"number of blocks is zero"};
+
+  const uint64_t dur = probs.size()/num;
+  if (probs.size() == 0) {
+    throw std::runtime_error {"size of probability matrix is empty"};
+  }
+  if (dur*num != probs.size()) {
+    throw std::runtime_error {"size of probability matrix is mismatch"};
+  }
+
+  blky_pathfinder_t pf = {};
+  pf.block_num = num;
+  pf.feat_bits = feat_bits;
+  pf.seed      = seed;
+
+  blky_pathfinder_init(&pf);
+
+  for (uint64_t t = 0; t < dur; ++t) {
+    blky_pathfinder_feed(&pf, probs.data()+t*num);
+  }
+  assert(pf.step_last);
+
+  uint32_t max_idx  = 0;
+  double   max_prob = -1;
+  for (uint32_t i = 0; i < num; ++i) {
+    if (max_prob < pf.probs[i]) {
+      max_prob = pf.probs[i];
+      max_idx  = i;
+    }
+  }
+
+  std::vector<uint32_t> ret;
+  blky_pathfinder_step_t* step = pf.step_last;
+  uint32_t idx = max_idx;
+  for (;;) {
+    ret.push_back(idx);
+    if (!step) break;
+    idx  = step->indices[idx];
+    step = step->prev;
+  }
+  std::reverse(ret.begin(), ret.end());
+
+  blky_pathfinder_deinit(&pf);
+  return ret;
+}
+
+std::vector<double> GenerateFeatureProbs(
+    const std::vector<uint32_t>&    features,
+    std::tuple<uint32_t, uint32_t>& block_num,
+    double                          false_positive,
+    double                          false_negative,
+    uint64_t                        seed,
+    bool                            normalize) {
+  const uint32_t num = std::get<0>(block_num) * std::get<1>(block_num);
+  if (num == 0) throw std::runtime_error {"number of blocks is zero"};
+
+  std::vector<double> ret;
+  for (auto c : features) {
+    const auto begin = ret.size();
+
+    double sum = 0;
+    for (uint32_t bi = 0; bi < num; ++bi) {
+      seed = blky_numeric_xorshift64(seed);
+      const double prob1 = static_cast<double>(seed%1000)/1000.;
+
+      seed = blky_numeric_xorshift64(seed);
+      const double prob2 = static_cast<double>(seed%1000)/1000.;
+
+      seed  = blky_numeric_xorshift64(seed);
+      double fprob;
+      if ((c == bi && prob1 >= false_negative) || prob2 < false_positive) {
+        fprob = static_cast<double>(seed%200)/1000. + .8;
+      } else {
+        fprob = static_cast<double>(seed%800)/1000.;
+      }
+      ret.push_back(fprob);
+      sum += fprob;
+    }
+    if (normalize) {
+      for (size_t i = begin; i < ret.size(); ++i) {
+        ret[i] /= sum;
+      }
+    }
+  }
+  return ret;
+}
+
+}  // namespace blky

blocky/main.cc

@@ -0,0 +1,307 @@
+#include <cassert>
+#include <exception>
+#include <iomanip>
+#include <iostream>
+#include <optional>
+#include <tuple>
+
+#include "common.hh"
+#include "bytes.hh"
+#include "embed.hh"
+#include "features.hh"
+
+#include <args.hxx>
+
+#include "video_decoder.hh"
+#include "video_encoder.hh"
+
+
+using namespace blky;
+
+
+args::ArgumentParser parser(
+    "liblocky command line tool",
+    "liblocky allow you to embed a bit array into video data secretly");
+
+args::MapFlag<std::string, DataFlow> from {
+    parser, kDataFlowList, "input layer specifier", {"from"}, kDataFlowMap, args::Options::Required};
+
+args::MapFlag<std::string, DataFlow> to {
+    parser, kDataFlowList, "output layer specifier", {"to"}, kDataFlowMap, args::Options::Required};
+
+args::Group src_group {
+  parser, "source specifier", args::Group::Validators::Xor, args::Options::Required};
+args::Flag src_stdin {src_group, "src-stdin", "read from stdin", {"src-stdin", "stdin"}};
+args::Flag src_stdin_hex {src_group, "src-stdin-hex", "read hex text from stdin", {"src-stdin-hex", "stdin-hex"}};
+args::ValueFlag<std::string> src_video {src_group, "path", "video input", {"src-video"}};
+
+args::Group dst_group  {
+  parser, "destination specifier", args::Group::Validators::Xor, args::Options::Required};
+args::Flag dst_stdout {dst_group, "dst-stdout", "write to stdout", {"dst-stdout", "stdout"}};
+args::Flag dst_stdout_hex {dst_group, "dst-stdout-hex", "write to stdout as hex text", {"dst-stdout-hex", "stdout-hex"}};
+args::ValueFlag<std::string> dst_video {dst_group, "path", "video output", {"dst-video"}};
+
+args::Group param_group {
+  parser, "general parameters", args::Group::Validators::DontCare
+};
+args::ValueFlag<std::tuple<uint32_t, uint32_t>> param_block_num {
+  param_group,
+  "int>0,int>0",
+  "number of features",
+  {"feature-num"},
+  {16, 16}
+};
+args::ValueFlag<uint32_t> param_block_first {
+  param_group,
+  "int>=0",
+  "an index of first block where feature will be embedded. used when encoding",
+  {"feature-first-index"},
+  0
+};
+args::ValueFlag<uint8_t> param_feat_bits {
+  param_group,
+  "int>0",
+  "number of bits that can be represented by a single feature",
+  {"feature-bits"},
+  1
+};
+args::ValueFlag<uint8_t> param_seed {
+  param_group,
+  "int>0",
+  "seed number for hopping randomization",
+  {"seed"},
+  123
+};
+args::ValueFlag<std::string> param_video {
+  param_group,
+  "path",
+  "a video file where information is embed",
+  {"path"},
+};
+args::ValueFlag<uint32_t> param_utime {
+  param_group,
+  "int>0",
+  "a duration (milliseconds) of features",
+  {"utime"},
+  10
+};
+
+args::Group probgen_group {
+  parser, "params for feature probability generator", args::Group::Validators::DontCare
+};
+args::ValueFlag<double> probgen_false_positive {
+  probgen_group,
+  "0<=double<=1",
+  "false positive ratio in feature probability generation",
+  {"probgen-false-positive"},
+  0
+};
+args::ValueFlag<double> probgen_false_negative {
+  probgen_group,
+  "0<=double<=1",
+  "false negative ratio in feature probability generation",
+  {"probgen-false-negative"},
+  0
+};
+args::ValueFlag<uint64_t> probgen_seed {
+  probgen_group,
+  "int>0",
+  "random seed",
+  {"probgen-seed"},
+  1
+};
+args::Flag probgen_normalize {
+  probgen_group,
+  "probgen-normalize",
+  "normalize probabilities",
+  {"probgen-normalize"},
+};
+
+
+std::vector<uint8_t>        bytes;
+std::vector<uint32_t>       features;
+std::vector<double>         feature_probs;
+std::optional<VideoDecoder> decoder;
+
+int main(int argc, char** argv)
+try {
+  parser.ParseCLI(argc, argv);
+
+  // read input
+  switch (args::get(from)) {
+  case kBytes:
+    if (src_stdin) {
+      std::string temp;
+      std::cin >> temp;
+      bytes = {temp.begin(), temp.end()};
+    } else if (src_stdin_hex) {
+      for (;;) {
+        char buf[2];
+        std::cin >> buf[0] >> buf[1];
+        if (std::cin.eof()) break;
+        bytes.push_back(static_cast<uint8_t>((ToHex(buf[0]) << 4) | ToHex(buf[1])));
+      }
+    } else {
+      throw std::runtime_error {"invalid source format for bytes"};
+    }
+    break;
+
+  case kFeatures:
+    if (src_stdin) {
+      features = ReadAll<uint32_t>(std::cin);
+    } else {
+      throw std::runtime_error {"invalid source format for features"};
+    }
+    break;
+
+  case kFeatureProbs:
+    if (src_stdin) {
+      feature_probs = ReadAll<double>(std::cin);
+    } else {
+      throw std::runtime_error {"invalid source format for feature probs"};
+    }
+    break;
+
+  case kVideo:
+    if (src_video) {
+      decoder.emplace(args::get(src_video));
+    } else {
+      throw std::runtime_error {"invalid source format for video"};
+    }
+  }
+
+  if (args::get(from) < args::get(to)) {
+    // execute encoding
+    switch (args::get(from)) {
+    case kBytes:
+      if (args::get(to) == kBytes) break;
+      features = BytesEncoder(
+          bytes,
+          args::get(param_block_num),
+          args::get(param_feat_bits),
+          args::get(param_block_first),
+          args::get(param_seed));
+      /* fallthrough */
+
+    case kFeatures:
+      if (args::get(to) == kFeatures) break;
+      if (args::get(to) == kFeatureProbs) {
+        feature_probs = GenerateFeatureProbs(
+            features,
+            args::get(param_block_num),
+            args::get(probgen_false_positive),
+            args::get(probgen_false_negative),
+            args::get(probgen_seed),
+            probgen_normalize);
+        break;
+      }
+      assert(dst_video);
+      Embed(
+          features,
+          args::get(dst_video),
+          args::get(param_video),
+          args::get(param_block_num),
+          args::get(param_utime));
+      /* fallthrough */
+
+    case kVideo:
+      if (args::get(to) == kVideo) break;
+      assert(false);
+
+    case kFeatureProbs:
+      throw std::runtime_error("couldn't start flow from the data");
+    }
+
+  } else if (args::get(from) > args::get(to)) {
+    // execute decoding
+    switch (args::get(from)) {
+    case kVideo:
+      if (args::get(to) == kVideo) break;
+      // TODO extract feature probs // features = XX
+      assert(false);
+
+    case kFeatureProbs:
+      if (args::get(to) == kFeatureProbs) break;
+      features = PathfindFeatures(
+          feature_probs,
+          args::get(param_block_num),
+          args::get(param_feat_bits),
+          args::get(param_seed));
+      /* fallthrough */
+
+    case kFeatures:
+      if (args::get(to) == kFeatures) break;
+      bytes = BytesDecoder(
+          features,
+          args::get(param_block_num),
+          args::get(param_feat_bits),
+          args::get(param_seed));
+      /* fallthrough */
+
+    case kBytes:
+      if (args::get(to) == kBytes) break;
+      assert(false);
+    }
+  }
+
+  // output
+  switch (args::get(to)) {
+  case kBytes:
+    if (dst_stdout) {
+      std::cout << std::string {bytes.begin(), bytes.end()} << std::endl;
+    } else if (dst_stdout_hex) {
+      for (auto c : bytes) {
+        std::cout << std::hex << (int) c;
+      }
+      std::cout << std::endl;
+    } else {
+      throw std::runtime_error {"invalid destination format for bytes"};
+    }
+    break;
+
+  case kFeatures:
+    if (dst_stdout) {
+      for (auto& f : features) std::cout << f << "\n";
+    } else {
+      throw std::runtime_error {"invalid destination format for features"};
+    }
+    break;
+
+  case kFeatureProbs:
+    if (dst_stdout) {
+      const auto size = args::get(param_block_num);
+      const auto cols = std::get<0>(size) * std::get<1>(size);
+      for (size_t i = 0; i < feature_probs.size();) {
+        for (size_t j = 0; i < feature_probs.size() && j < cols; ++i, ++j) {
+          std::cout << feature_probs[i] << " ";
+        }
+        std::cout << "\n";
+      }
+    } else {
+      throw std::runtime_error {"invalid destination format for feature probs"};
+    }
+    break;
+
+  case kVideo:
+    break;
+  }
+  return 0;
+
+} catch (const args::Help&) {
+  std::cout << parser << std::endl;
+  return 0;
+
+} catch (const args::ParseError& e) {
+  std::cerr << e.what() << std::endl;
+  std::cerr << parser   << std::endl;
+  return 1;
+
+} catch (const args::ValidationError& e) {
+  std::cerr << e.what() << std::endl;
+  std::cerr << parser   << std::endl;
+  return 1;
+
+} catch (const std::runtime_error& e) {
+  std::cerr << "runtime error: " << e.what() << std::endl;
+  return 1;
+}

blocky/video_decoder.hh

@@ -0,0 +1,179 @@
+#pragma once
+
+#include <cassert>
+#include <cstdint>
+#include <cstring>
+#include <exception>
+#include <fstream>
+#include <string_view>
+#include <vector>
+
+#include <minimp4.h>
+
+#include <wels/codec_api.h>
+
+
+class VideoDecoder final {
+ public:
+  VideoDecoder() = delete;
+  VideoDecoder(const std::string& path) :
+      file_(path, std::ifstream::binary | std::ifstream::ate),
+      size_(static_cast<size_t>(file_.tellg())) {
+    // init objects
+    SDecodingParam dparam = {};
+    dparam.sVideoProperty.eVideoBsType = VIDEO_BITSTREAM_DEFAULT;
+    dparam.eEcActiveIdc = ERROR_CON_SLICE_COPY;
+
+    WelsCreateDecoder(&decoder_);
+    decoder_->Initialize(&dparam);
+
+    int lv = WELS_LOG_DEBUG;
+    decoder_->SetOption(DECODER_OPTION_TRACE_LEVEL, &lv);
+
+    demuxer_ = {};
+    MP4D_open(&demuxer_, ReadCallback, this, static_cast<int64_t>(size_));
+
+    // find video track
+    track_ = SIZE_MAX;
+    for (size_t i = 0; i < demuxer_.track_count; ++i) {
+      if (demuxer_.track[i].handler_type == MP4D_HANDLER_TYPE_VIDE) {
+        if (track_ != SIZE_MAX) {
+          throw std::runtime_error {"there are many video tracks"};
+        }
+        track_ = i;
+      }
+    }
+    if (track_ == SIZE_MAX) {
+      throw std::runtime_error {"there is no video track"};
+    }
+
+    // setup decoder
+    std::vector<uint8_t> temp_;
+
+    for (size_t i = 0;; ++i) {
+      int size;
+      auto sps = static_cast<const uint8_t*>(MP4D_read_sps(
+              &demuxer_,
+              static_cast<unsigned int>(track_),
+              static_cast<int>(i),
+              &size));
+      if (!sps) break;
+
+      temp_.resize(static_cast<size_t>(4+size));
+      temp_[0] = 0, temp_[1] = 0, temp_[2] = 0, temp_[3] = 1;
+      std::memcpy(&temp_[4], sps, static_cast<size_t>(size));
+
+      if (decoder_->DecodeFrameNoDelay(temp_.data(), static_cast<int>(temp_.size()), yuv_, &frame_)) {
+        throw std::runtime_error {"failed to decode SPS"};
+      }
+    }
+    for (size_t i = 0;; ++i) {
+      int size;
+      auto pps = static_cast<const uint8_t*>(MP4D_read_pps(
+              &demuxer_,
+              static_cast<unsigned int>(track_),
+              static_cast<int>(i),
+              &size));
+      if (!pps) break;
+
+      temp_.resize(static_cast<size_t>(4+size));
+      temp_[0] = 0, temp_[1] = 0, temp_[2] = 0, temp_[3] = 1;
+      std::memcpy(&temp_[4], pps, static_cast<size_t>(size));
+
+      if (decoder_->DecodeFrameNoDelay(temp_.data(), static_cast<int>(temp_.size()), yuv_, &frame_)) {
+        throw std::runtime_error {"failed to decode SPS"};
+      }
+    }
+    temp_.clear();
+  }
+  ~VideoDecoder() noexcept {
+    decoder_->Uninitialize();
+    WelsDestroyDecoder(decoder_);
+
+    MP4D_close(&demuxer_);
+  }
+  VideoDecoder(const VideoDecoder&) = delete;
+  VideoDecoder(VideoDecoder&&) = delete;
+  VideoDecoder& operator=(const VideoDecoder&) = delete;
+  VideoDecoder& operator=(VideoDecoder&&) = delete;
+
+  bool Decode() {
+    if (temp_consumed_ >= temp_.size()) {
+      if (count_ >= demuxer_.track[track_].sample_count) {
+        return false;
+      }
+
+      unsigned size, time, dur;
+      const auto off = MP4D_frame_offset(
+          &demuxer_,
+          static_cast<unsigned int>(track_),
+          static_cast<unsigned int>(count_),
+          &size, &time, &dur);
+
+      assert(size > 0);
+      temp_.resize(size);
+      temp_consumed_ = 0;
+
+      file_.seekg(static_cast<std::streamoff>(off));
+      assert(file_);
+      file_.read((char*) temp_.data(), size);
+      assert(file_);
+
+      Decode();
+      ++count_;
+      return true;
+
+    } else {
+      auto& i = temp_consumed_;
+
+      const uint32_t nal_size = 4 +
+          static_cast<uint32_t>((temp_[i+0] << 24) |
+                                (temp_[i+1] << 16) |
+                                (temp_[i+2] <<  8) |
+                                (temp_[i+3] <<  0));
+      temp_[i  ] = 0;
+      temp_[i+1] = 0;
+      temp_[i+2] = 0;
+      temp_[i+3] = 1;
+      if (decoder_->DecodeFrameNoDelay(temp_.data()+i, static_cast<int>(nal_size), yuv_, &frame_)) {
+        throw std::runtime_error {"failed to decode a frame"};
+      }
+      i += nal_size;
+      return true;
+    }
+  }
+
+  const SBufferInfo& frame() const noexcept { return frame_; }
+  const uint8_t* const* yuv() const noexcept { return yuv_; }
+
+ private:
+  std::ifstream file_;
+  size_t size_;
+
+  ISVCDecoder* decoder_;
+  MP4D_demux_t demuxer_;
+
+  size_t track_;
+
+  uint8_t*    yuv_[3] = {0};
+  SBufferInfo frame_  = {};
+  size_t      count_  = 0;
+
+  size_t temp_consumed_ = 0;
+  std::vector<uint8_t> temp_;
+
+
+  static int ReadCallback(int64_t off, void* buf, size_t size, void* ptr) noexcept {
+    auto self = (VideoDecoder*) ptr;
+
+    auto n = self->size_ - static_cast<size_t>(off) - size;
+    if (size < n) n = size;
+
+    self->file_.seekg(off);
+    assert(self->file_);
+
+    self->file_.read((char*) buf, static_cast<std::streamsize>(n));
+    assert(self->file_);
+    return 0;
+  }
+};

blocky/video_encoder.hh

@@ -0,0 +1,79 @@
+#pragma once
+
+#include <cassert>
+#include <cstdint>
+#include <cstring>
+#include <exception>
+#include <fstream>
+#include <string>
+
+#include <minimp4.h>
+
+#include <wels/codec_api.h>
+
+
+class VideoEncoder final {
+ public:
+  VideoEncoder() = delete;
+  VideoEncoder(const std::string& path, const SEncParamBase& p) :
+      file_(path, std::ofstream::binary), fps_(p.fMaxFrameRate) {
+    if (WelsCreateSVCEncoder(&encoder_)) {
+      throw std::runtime_error {"failed to init openh264 encoder"};
+    }
+
+    muxer_ = MP4E_open(false, false, this, WriteCallback);
+    if (MP4E_STATUS_OK != mp4_h26x_write_init(&writer_, muxer_, p.iPicWidth, p.iPicHeight, false)) {
+      throw std::runtime_error {"failed to init h26x writer"};
+    }
+
+    encoder_->Initialize(&p);
+
+    int lv = WELS_LOG_DEBUG;
+    encoder_->SetOption(ENCODER_OPTION_TRACE_LEVEL, &lv);
+    int fmt = videoFormatI420;
+    encoder_->SetOption(ENCODER_OPTION_DATAFORMAT, &fmt);
+  }
+  ~VideoEncoder() noexcept {
+    encoder_->Uninitialize();
+    WelsDestroySVCEncoder(encoder_);
+
+    MP4E_close(muxer_);
+    mp4_h26x_write_close(&writer_);
+  }
+
+  void Encode(const SSourcePicture& pic) {
+    SFrameBSInfo info;
+    if (cmResultSuccess != encoder_->EncodeFrame(&pic, &info)) {
+      throw std::runtime_error {"failed to encode a frame"};
+    }
+    if (info.eFrameType == videoFrameTypeSkip) return;
+
+    for (size_t i = 0; i < static_cast<size_t>(info.iLayerNum); ++i) {
+      const auto& lay = info.sLayerInfo[i];
+
+      uint8_t* buf = lay.pBsBuf;
+      for (size_t j = 0; j < static_cast<size_t>(lay.iNalCount); ++j) {
+        mp4_h26x_write_nal(&writer_, buf, lay.pNalLengthInByte[j], static_cast<unsigned int>(90000./fps_));
+        buf += lay.pNalLengthInByte[j];
+      }
+    }
+  }
+
+ private:
+  std::ofstream file_;
+
+  ISVCEncoder*      encoder_;
+  MP4E_mux_t*       muxer_;
+  mp4_h26x_writer_t writer_;
+
+  double fps_;
+
+
+  static int WriteCallback(int64_t off, const void* buf, size_t size, void* ptr) noexcept {
+    auto self = (VideoEncoder*) ptr;
+    self->file_.seekp(off);
+    assert(self->file_);
+    self->file_.write(static_cast<const char*>(buf), static_cast<std::streamsize>(size));
+    return !self->file_;
+  }
+};

exp/pathfinder-plot.sh

@@ -0,0 +1,39 @@
+#!/bin/bash
+
+set -e
+
+BLKY_FLAGS="${BLKY_FLAGS}"
+
+
+echo "this will take 4 hours" > /dev/stderr
+
+function count_bits() {
+  err=$((0))
+  for i in $(seq 0 7); do
+    if [[ $(($1 & (1 << $i))) > 0 ]]; then
+      err=$((err + 1))
+    fi
+  done
+  echo $err
+}
+function proc() {
+  err=$((0))
+  for i in $(seq 0 255); do
+    v=$(printf "%02x" ${i} |  \
+      ${BLKY} ${BLKY_FLAGS} --from bytes --stdin-hex --to feature-probs --stdout   \
+           --probgen-false-positive $(printf "0.%02d" ${1})  \
+           --probgen-false-negative $(printf "0.%02d" ${2}) |  \
+      ${BLKY} ${BLKY_FLAGS}  --from feature-probs --stdin --to bytes --stdout-hex)
+      err=$((err + $(count_bits $(($i ^ 16#$v)))))
+  done
+  printf "0.%02d 0.%02d %d\n" ${fn} ${fp} ${err}
+}
+
+for fn in $(seq 0 50); do
+  for fp in $(seq 0 50); do
+    until [ "$( jobs -lr 2>&1 | wc -l)" -lt 4 ]; do
+      sleep 0.5
+    done
+    proc ${fp} ${fn} &
+  done
+done

liblocky/CMakeLists.txt

@@ -0,0 +1,20 @@
+cmake_minimum_required(VERSION 3.18)
+
+add_library(liblocky)
+set_target_properties(liblocky PROPERTIES PREFIX "")
+target_compile_options(liblocky PRIVATE ${BLOCKY_C_FLAGS})
+
+target_link_libraries(liblocky PUBLIC m)
+target_include_directories(liblocky PUBLIC SYSTEM .)
+target_sources(liblocky
+  PUBLIC
+    liblocky.h
+  PRIVATE
+    block.c
+    decoder.c
+    encoder.c
+    extractor.c
+    image.c
+    pathfinder.c
+    sensor.c
+)

liblocky/block.c

@@ -0,0 +1,24 @@
+#include "liblocky.h"
+
+
+double blky_block_estimate(
+    const blky_sensor_t* sensors, uint64_t n,
+    double max_var, double min_avg) {
+  double sum = 0;
+  for (uint64_t i = 0; i < n; ++i) {
+    sum += blky_abs(sensors[i].correl);
+  }
+  const double avg = sum / (double) n;
+
+  double var = 0;
+  for (uint64_t i = 0; i < n; ++i) {
+    const double diff = blky_abs(sensors[i].correl) - avg;
+    var += diff*diff;
+  }
+  var /= (double) n;
+
+  // FIXME: calculate probability
+  if (var > max_var) return 0;
+  if (avg < min_avg) return 0;
+  return 1;
+}

liblocky/decoder.c

@@ -0,0 +1,51 @@
+#include "liblocky.h"
+
+#include <assert.h>
+
+
+bool blky_decoder_feed(blky_decoder_t* de, uint32_t block_index) {
+  assert(de->block_num > 0);
+  assert(de->feat_bits > 0);
+  assert(de->seed      > 0);
+
+  const uint64_t seed = blky_numeric_xorshift64(de->seed);
+
+  if (de->count++) {
+    const uint32_t feat_max = 1 << de->feat_bits;
+    assert(feat_max < de->block_num);
+
+    uint32_t feat = 0;
+    for (; feat < feat_max; ++feat) {
+      if (blky_numeric_hop(de->block_index, feat, seed)%de->block_num == block_index) {
+        break;
+      }
+    }
+    if (feat >= feat_max) return false;
+
+    assert(de->scrap_bits+8 <= 32);
+    de->scrap      |= feat << de->scrap_bits;
+    de->scrap_bits += de->feat_bits;
+  }
+
+  de->seed        = seed;
+  de->block_index = block_index;
+  return true;
+}
+
+bool blky_decoder_pop(blky_decoder_t* de, uint8_t* b, bool force) {
+  if (force) {
+    if (de->scrap_bits > 0) return false;
+  } else {
+    if (de->scrap_bits < 8) return false;
+  }
+
+  *b = (uint8_t) (de->scrap & 0xFF);
+
+  de->scrap >>= 8;
+  if (de->scrap_bits >= 8) {
+    de->scrap_bits -= 8;
+  } else {
+    de->scrap_bits = 0;
+  }
+  return true;
+}

liblocky/encoder.c

@@ -0,0 +1,41 @@
+#include "liblocky.h"
+
+#include <assert.h>
+#include <stdio.h>
+
+
+void blky_encoder_feed(blky_encoder_t* enc, uint8_t data) {
+  assert(enc);
+  assert(enc->block_num  > 0);
+  assert(enc->feat_bits  > 0);
+  assert(enc->feat_bits <= sizeof(enc->scrap)*8);
+  assert(enc->seed       > 0);
+
+  if (enc->count++ == 0) enc->scrap_bits = 1;
+
+  assert(enc->scrap_bits+8 <= 32);
+  enc->scrap      |= data << enc->scrap_bits;
+  enc->scrap_bits += 8;
+}
+
+bool blky_encoder_pop(blky_encoder_t* enc, uint32_t* feat, bool force) {
+  assert(enc);
+  assert(feat);
+
+  if (force) {
+    if (enc->scrap_bits > 0) return false;
+  } else {
+    if (enc->scrap_bits < enc->feat_bits) return false;
+  }
+
+  const uint32_t feat_max = 1 << enc->feat_bits;
+  assert(feat_max < enc->block_num);
+
+  enc->seed = blky_numeric_xorshift64(enc->seed);
+  *feat     = blky_numeric_hop(enc->block_index, enc->scrap%feat_max, enc->seed)%enc->block_num;
+
+  enc->block_index   = *feat;
+  enc->scrap       >>= enc->feat_bits;
+  enc->scrap_bits   -= enc->feat_bits;
+  return true;
+}

liblocky/extractor.c

@@ -0,0 +1,89 @@
+#include "liblocky.h"
+
+#include <assert.h>
+#include <stdlib.h>
+#include <string.h>
+
+
+void blky_extractor_init(blky_extractor_t* ex) {
+  assert(ex->block_num_x > 0);
+  assert(ex->block_num_y > 0);
+  assert(ex->sensor_num_block_x > 0);
+  assert(ex->sensor_num_block_y > 0);
+  assert(ex->samples_per_pix > 0);
+  assert(ex->utime > 0);
+  assert(ex->utime <= BLKY_SENSOR_MAX_DUR);
+  assert(ex->pix_stride > 0);
+  assert(ex->pix_stride > ex->samples_per_pix);
+  assert(ex->correl_max_var >= 0);
+  assert(ex->correl_min_avg >= 0);
+  assert(ex->correl_min_avg <= 1);
+
+  // set immutable params
+  ex->block_num = ex->block_num_x * ex->block_num_y;
+
+  ex->sensor_num_block =
+    ex->sensor_num_block_x * ex->sensor_num_block_y *
+    ex->samples_per_pix;
+
+  ex->sensor_num_whole = ex->sensor_num_block * ex->block_num;
+
+  ex->block_w = 1./ex->block_num_x;
+  ex->block_h = 1./ex->block_num_y;
+
+  ex->sensor_interval_x = ex->block_w / ex->sensor_num_block_x;
+  ex->sensor_interval_y = ex->block_h / ex->sensor_num_block_y;
+
+  // clear states
+  ex->time = 0;
+
+  const uint64_t sensors_bytes = sizeof(ex->sensors[0]) * ex->sensor_num_whole;
+  const uint64_t probs_bytes   = sizeof(ex->probs[0]) * ex->block_num;
+
+  const uint64_t mem_bytes = sensors_bytes + probs_bytes;
+  ex->mem = calloc(mem_bytes, 1);
+  assert(ex->mem);
+
+  ex->sensors = (blky_sensor_t*) ex->mem;
+  ex->probs   = (double*) (ex->mem + sensors_bytes);
+}
+void blky_extractor_deinit(blky_extractor_t* ex) {
+  free(ex->mem);
+}
+
+bool blky_extractor_feed(
+    blky_extractor_t* ex,
+    const uint8_t* img, uint32_t w, uint32_t h, const double verts[8]) {
+  blky_sensor_t* s = ex->sensors;
+
+  for (uint32_t by = 0; by < ex->block_num_y; ++by) {
+    const double byf = by*ex->block_h + ex->sensor_interval_y/2.;
+    for (uint32_t bx = 0; bx < ex->block_num_x; ++bx) {
+      const double bxf = bx*ex->block_w + ex->sensor_interval_x/2.;
+      for (uint32_t sy = 0; sy < ex->sensor_num_block_y; ++sy) {
+        const double syf = byf + ex->sensor_interval_y*sy;
+        for (uint32_t sx = 0; sx < ex->sensor_num_block_x; ++sx) {
+          const double sxf = bxf + ex->sensor_interval_x*sx;
+
+          const uint8_t* base = img +
+            ex->pix_stride*blky_image_offset(w, h, verts, sxf, syf);
+          for (uint8_t off = 0; off < ex->samples_per_pix; ++off) {
+            const float v = base[off]*1.f / UINT8_MAX;
+            blky_sensor_feed(s++, &v, 1);
+          }
+        }
+      }
+    }
+  }
+
+  ++ex->time;
+  if (ex->time%ex->utime > 0) return false;
+
+  for (uint64_t bi = 0; bi < ex->block_num; ++bi) {
+    blky_sensor_t* s = ex->sensors + (bi*ex->sensor_num_block);
+    ex->probs[bi] = blky_block_estimate(
+        s, ex->sensor_num_block, ex->correl_max_var, ex->correl_min_avg);
+    memset(s, 0, ex->sensor_num_block*sizeof(*s));
+  }
+  return true;
+}

liblocky/image.c

@@ -0,0 +1,32 @@
+#include "liblocky.h"
+
+#include <assert.h>
+
+
+void blky_image_convert_to_normalized_coord(
+    const double verts[8], double* x, double* y) {
+  assert(0 <= *x && *x < 1);
+  assert(0 <= *y && *y < 1);
+
+  const double rx = 1-*x, ry = 1-*y;
+
+  const double top_xf = *x*verts[6] + rx*verts[0];
+  const double bot_xf = *x*verts[4] + rx*verts[2];
+
+  const double lef_yf = *y*verts[3] + ry*verts[1];
+  const double rig_yf = *y*verts[5] + ry*verts[7];
+
+  const double ans_x = *y*bot_xf + ry*top_xf;
+  const double ans_y = *x*rig_yf + rx*lef_yf;
+
+  *x = ans_x;
+  *y = ans_y;
+}
+
+uint64_t blky_image_offset(
+    uint32_t w, uint32_t h, const double verts[8], double x, double y) {
+  blky_image_convert_to_normalized_coord(verts, &x, &y);
+  const uint32_t xi = (uint32_t) (x*w);
+  const uint32_t yi = (uint32_t) (y*h);
+  return yi*w + xi;
+}

liblocky/liblocky.h

@@ -0,0 +1,206 @@
+#pragma once
+
+#include <stdbool.h>
+#include <stdint.h>
+
+
+#define blky_clamp(x, a, b) ((x) < (a)? (a): (x) > (b)? (b): (x))
+#define blky_abs(x) ((x) < 0? -(x): (x))
+
+
+/* ---- Sensor ----
+ * calculates correl from samples sequencially */
+
+typedef struct blky_sensor_t {
+  double sum, avg, var;
+  uint64_t time;
+
+  double cov;
+  double correl;
+  double prev_correl;
+
+# define BLKY_SENSOR_MAX_DUR 64
+  float values[BLKY_SENSOR_MAX_DUR];
+} blky_sensor_t;
+
+void
+blky_sensor_feed(
+    blky_sensor_t*, const float* v, uint64_t n);
+
+void
+blky_sensor_drop(
+    blky_sensor_t*, uint64_t until);
+
+
+/* ---- Block ----
+ * calculates probability representing how likely the block is a feature,
+ * by dealing with multiple sensors */
+
+double  /* 0~1 probability */
+blky_block_estimate(
+    const blky_sensor_t* sensors, uint64_t n,
+    double max_var, double min_avg);
+
+
+/* ---- Extractor ----
+ * extracts all features from a sequence of multiple frames */
+
+typedef struct blky_extractor_t {
+  // must be filled before init()
+  uint32_t block_num_x;
+  uint32_t block_num_y;
+  uint32_t sensor_num_block_x;
+  uint32_t sensor_num_block_y;
+  uint32_t samples_per_pix;
+  uint32_t pix_stride;
+  uint32_t utime;
+
+  double correl_max_var;
+  double correl_min_avg;
+
+  // immutable internal params
+  uint32_t block_num;
+  uint32_t sensor_num_block;
+  uint32_t sensor_num_whole;
+
+  double block_w;
+  double block_h;
+  double sensor_interval_x;
+  double sensor_interval_y;
+
+  // mutable internal state
+  uint64_t time;
+
+  // heap objects
+  uint8_t*       mem;
+  blky_sensor_t* sensors;
+  double*        probs;
+} blky_extractor_t;
+
+typedef struct blky_extractor_feat_t {
+  uint32_t block;
+  uint64_t begin;
+  double   prob;
+} blky_extractor_feat_t;
+
+void
+blky_extractor_init(
+    blky_extractor_t* ex);
+
+void
+blky_extractor_deinit(
+    blky_extractor_t* ex);
+
+bool
+blky_extractor_feed(
+    blky_extractor_t* ex,
+    const uint8_t* img, uint32_t w, uint32_t h, const double verts[8]);
+
+
+/* ---- Pathfinder ---- */
+typedef struct blky_pathfinder_step_t {
+  struct blky_pathfinder_step_t* prev;
+  uint32_t indices[1];
+} blky_pathfinder_step_t;
+
+typedef struct blky_pathfinder_t {
+  // must be filled before init()
+  uint32_t block_num;
+  uint32_t feat_bits;
+  uint64_t seed;
+
+  // internal state
+  blky_pathfinder_step_t* step_last;
+
+  uint64_t steps;
+  uint64_t step_bytes;
+
+  double* probs;
+  double* probs_prev;
+} blky_pathfinder_t;
+
+void
+blky_pathfinder_init(
+    blky_pathfinder_t* pf);
+
+void
+blky_pathfinder_deinit(
+    blky_pathfinder_t* pf);
+
+void
+blky_pathfinder_feed(
+    blky_pathfinder_t* pf,
+    const double*      probs);
+
+
+/* ---- Encoder ----
+ * converts byte to feature */
+
+typedef struct blky_encoder_t {
+  uint32_t block_num;
+  uint32_t block_index;
+  uint8_t  feat_bits;
+  uint64_t seed;
+
+  uint64_t count;
+  uint32_t scrap;
+  uint8_t  scrap_bits;
+} blky_encoder_t;
+
+void
+blky_encoder_feed(
+    blky_encoder_t* enc,
+    uint8_t         data);
+
+bool
+blky_encoder_pop(
+    blky_encoder_t* enc,
+    uint32_t*       feat,
+    bool            force);
+
+
+/* ---- Decoder ----
+ * converts block indices to byte */
+
+typedef struct blky_decoder_t {
+  uint32_t block_num;
+  uint8_t  feat_bits;
+  uint64_t seed;
+
+  uint64_t count;
+  uint32_t block_index;
+  uint32_t scrap;
+  uint8_t  scrap_bits;
+} blky_decoder_t;
+
+bool
+blky_decoder_feed(
+    blky_decoder_t* de,
+    uint32_t        block_index);
+
+bool
+blky_decoder_pop(
+    blky_decoder_t* de,
+    uint8_t*        b,
+    bool            force);
+
+
+/* ---- Image utility ---- */
+void blky_image_convert_to_normalized_coord(
+    const double verts[8], double* x, double* y);
+
+uint64_t blky_image_offset(
+    uint32_t w, uint32_t h, const double verts[8], double x, double y);
+
+
+/* ---- numeric utility ---- */
+static inline uint64_t blky_numeric_xorshift64(uint64_t x) {
+  x ^= x << 13;
+  x ^= x >> 7;
+  x ^= x << 17;
+  return x;
+}
+static inline uint32_t blky_numeric_hop(uint32_t prev, uint32_t offset, uint64_t seed) {
+  seed = (seed^blky_numeric_xorshift64(prev+seed)) + (offset << 4);
+  return (uint32_t) ((seed & 0xFFFFFFFF) ^ (seed >> 32));
+}

liblocky/pathfinder.c

@@ -0,0 +1,74 @@
+#include "liblocky.h"
+
+#include <assert.h>
+#include <stdlib.h>
+
+
+void blky_pathfinder_init(blky_pathfinder_t* pf) {
+  assert(pf->block_num > 0);
+  assert(pf->feat_bits > 0);
+  assert(pf->feat_bits < 32);
+  assert(pf->seed      > 0);
+
+  assert(pf->block_num > (uint32_t) (1 << pf->feat_bits));
+
+  pf->step_last  = NULL;
+
+  pf->step_bytes =
+    sizeof(blky_pathfinder_step_t) +
+    sizeof(uint32_t)*(pf->block_num-1);
+
+  pf->probs = calloc(sizeof(*pf->probs), pf->block_num);
+  assert(pf->probs);
+
+  pf->probs_prev = calloc(sizeof(*pf->probs_prev), pf->block_num);
+  assert(pf->probs_prev);
+}
+
+void blky_pathfinder_deinit(blky_pathfinder_t* pf) {
+  free(pf->probs);
+  free(pf->probs_prev);
+
+  blky_pathfinder_step_t* step = pf->step_last;
+  while (step) {
+    blky_pathfinder_step_t* temp = step->prev;
+    free(step);
+    step = temp;
+  }
+}
+
+void blky_pathfinder_feed(blky_pathfinder_t* pf, const double* probs) {
+  double* temp = pf->probs;
+  pf->probs      = pf->probs_prev;
+  pf->probs_prev = temp;
+
+  blky_pathfinder_step_t* step = NULL;
+  if (++pf->steps > 1) {
+    step = calloc(pf->step_bytes, 1);
+    assert(step);
+  }
+
+  const uint32_t feat_max = 1 << pf->feat_bits;
+  assert(feat_max < pf->block_num);
+
+  pf->seed = blky_numeric_xorshift64(pf->seed);
+  for (uint32_t bi = 0; bi < pf->block_num; ++bi) {
+    const double prob = probs[bi];
+    for (uint32_t pbi = 0; pbi < pf->block_num; ++pbi) {
+      for (uint32_t fi = 0; fi < feat_max; ++fi) {
+        if (blky_numeric_hop(pbi, fi, pf->seed)%pf->block_num != bi) continue;
+
+        const double sum  = pf->probs_prev[pbi] + prob;
+        if (pf->probs[bi] < sum) {
+          pf->probs[bi] = sum;
+          if (step) step->indices[bi] = pbi;
+        }
+      }
+    }
+  }
+
+  if (step) {
+    step->prev    = pf->step_last;
+    pf->step_last = step;
+  }
+}

liblocky/sensor.c

@@ -0,0 +1,53 @@
+#include "liblocky.h"
+
+#include <assert.h>
+#include <math.h>
+#include <string.h>
+
+
+void blky_sensor_feed(blky_sensor_t* b, const float* v, uint64_t n) {
+  assert(b->time+n <= BLKY_SENSOR_MAX_DUR);
+  if (n == 0) return;
+
+  const uint64_t t = b->time;
+  memcpy(b->values+t, v, n*sizeof(float));
+  b->time += n;
+
+  for (uint64_t i = t; i < b->time; ++i) {
+    b->sum += b->values[i];
+  }
+  b->avg = b->sum / (double) b->time;
+
+  b->var = 0;
+  for (uint64_t i = 0; i < b->time; ++i) {
+    const double diff = b->values[i] - b->avg;
+    b->var += diff*diff;
+  }
+  b->var /= (double) b->time;
+
+  b->cov = 0;
+  for (uint64_t i = 0; i < b->time; ++i) {
+    const double diff_v = b->values[i] - b->avg;
+    const double diff_t = (double) i - (double) b->time/2.;
+    b->cov += diff_v * diff_t;
+  }
+  b->cov /= (double) b->time;
+
+  const double tf   = (double) b->time;
+  const double tvar =
+    (tf/6*(tf+1)*(2*tf+2)-tf*tf/2-tf*tf*tf/4)/tf;
+
+  b->prev_correl = b->correl;
+
+  if (b->var == 0) {
+    b->correl = 1;
+  } else {
+    b->correl = b->cov / sqrt(b->var*tvar);
+  }
+}
+
+void blky_sensor_drop(blky_sensor_t* b, uint64_t until) {
+  (void) b;
+  (void) until;
+  // TODO
+}

playground/CMakeLists.txt

@@ -0,0 +1,30 @@
+add_executable(playground)
+target_compile_options(playground PRIVATE ${BLOCKY_CXX_FLAGS})
+
+target_sources(playground
+  PRIVATE
+    app.hh
+    block.cc
+    encoder.cc
+    exporter.cc
+    extractor.cc
+    initiator.hh
+    input.cc
+    input.hh
+    input_gradient.cc
+    input_imgseq.cc
+    input_noise.cc
+    main.cc
+    player.cc
+    sensor.cc
+)
+target_link_libraries(playground
+  PRIVATE
+    liblocky
+    glew
+    glfw
+    imgui
+    implot
+    m
+    stb
+)

playground/app.hh

@@ -0,0 +1,61 @@
+#pragma once
+
+#include <functional>
+#include <map>
+#include <memory>
+#include <string>
+#include <utility>
+
+
+namespace pg {
+
+class App {
+ public:
+  App() = default;
+  virtual ~App() = default;
+  App(const App&) = delete;
+  App(App&&) = delete;
+  App& operator=(const App&) = delete;
+  App& operator=(App&&) = delete;
+
+  virtual void Update() noexcept = 0;
+
+ protected:
+  struct TypeInfo final {
+   public:
+    using Factory = std::function<std::unique_ptr<App>()>;
+
+    template <typename T>
+    static TypeInfo Create(const char* name) noexcept {
+      return {name, []() { return std::make_unique<T>(); }};
+    }
+
+    TypeInfo(const char* name, Factory&& f) noexcept :
+        name_(name), factory_(std::move(f)) {
+      registry_[name_] = this;
+    }
+    ~TypeInfo() noexcept {
+      registry_.erase(name_);
+    }
+
+    std::unique_ptr<App> Create() noexcept {
+      return factory_();
+    }
+
+    const char* name() const noexcept { return name_; }
+
+   private:
+    const char* name_;
+
+    Factory factory_;
+  };
+
+  static const std::map<std::string, TypeInfo*>& registry() noexcept {
+    return registry_;
+  }
+
+ private:
+  static inline std::map<std::string, TypeInfo*> registry_;
+};
+
+}  // namespace pg

playground/block.cc

@@ -0,0 +1,183 @@
+extern "C" {
+# include <liblocky.h>
+}
+
+#include <numeric>
+#include <vector>
+
+#include <imgui.h>
+#include <implot.h>
+
+#include "app.hh"
+#include "input.hh"
+
+
+namespace pg {
+
+class Block final : public App {
+ public:
+  static inline TypeInfo kType = TypeInfo::Create<Block>("Block");
+
+  Block() noexcept {
+  }
+
+  void Update() noexcept {
+    const auto id = "Block | "+
+        std::to_string(reinterpret_cast<uintptr_t>(this));
+
+    const auto em = ImGui::GetFontSize();
+    ImGui::SetNextWindowSize({40*em, 24*em}, ImGuiCond_Once);
+    if (ImGui::Begin(id.c_str())) {
+      UpdateParams();
+      ImGui::SameLine();
+      UpdatePlots();
+      if (msg_.size() > 0) {
+        ImGui::TextUnformatted(msg_.c_str());
+      }
+    }
+    ImGui::End();
+  }
+  void UpdateParams() noexcept {
+    const auto em = ImGui::GetFontSize();
+
+    bool mod = false;
+    ImGui::BeginGroup();
+    ImGui::PushItemWidth(6*em);
+    mod |= ImGui::DragInt("input_slot", &src_, 1, 0, 1024);
+    mod |= ImGui::DragInt("time", &time_, 1, 0, 1024);
+    mod |= ImGui::DragInt("dur", &dur_, 1, 1, BLKY_SENSOR_MAX_DUR);
+    mod |= ImGui::DragInt2("blocks", block_, 1, 1, 16);
+    mod |= ImGui::DragInt2("sensors", sensor_);
+    mod |= ImGui::DragFloat("var thresh", &var_thresh_);
+    ImGui::PopItemWidth();
+    ImGui::EndGroup();
+
+    if (mod) Calc();
+  }
+  void UpdatePlots() noexcept {
+    const auto em = ImGui::GetFontSize();
+
+    ImGui::BeginGroup();
+
+    auto avail = ImGui::GetContentRegionAvail();
+    avail.x -= em;
+    avail.y -= em;
+    if (ImPlot::BeginPlot("correls", avail)) {
+      ImPlot::SetupAxisLimits(ImAxis_X1, 0, 1);
+      ImPlot::SetupAxisLimits(ImAxis_Y1, 0, 1);
+      ImPlot::PlotHeatmap(
+          "correls_avg", correl_avg_.data(), block_[1], block_[0],
+          0, 1, "%.2f");
+      ImPlot::PlotHeatmap(
+          "correls_var", correl_var_.data(), block_[1], block_[0],
+          0, var_thresh_, "%.4f");
+      ImPlot::PlotHeatmap(
+          "correls_avg_filtered", correl_avg_filtered_.data(), block_[1], block_[0],
+          0, 1, "%.2f");
+      ImPlot::EndPlot();
+    }
+    ImGui::EndGroup();
+  }
+
+ private:
+  int src_ = 0;
+  int time_ = 0;
+  int dur_ = 30;
+  int block_[2] = {1, 1};
+  int sensor_[2] = {1, 1};
+
+  float var_thresh_ = .1f;
+
+  std::string msg_;
+
+  std::vector<double> correls_    = {0};
+  std::vector<double> correl_avg_ = {0};
+  std::vector<double> correl_var_ = {0};
+  std::vector<double> correl_avg_filtered_ = {0};
+
+  void Calc() noexcept
+  try {
+    auto data = Input::instance().slots(static_cast<size_t>(src_));
+    if (!data) throw "missing input";
+
+    const auto block_x = static_cast<size_t>(block_[0]);
+    const auto block_y = static_cast<size_t>(block_[1]);
+    const auto block_n = block_x*block_y;
+    const auto block_w = 1.f / static_cast<float>(block_x);
+    const auto block_h = 1.f / static_cast<float>(block_y);
+
+    const auto sensor_x = static_cast<size_t>(sensor_[0]);
+    const auto sensor_y = static_cast<size_t>(sensor_[1]);
+    const auto sensor_n = sensor_x*sensor_y;
+    const auto sensor_interval_x = block_w / static_cast<float>(sensor_x);
+    const auto sensor_interval_y = block_h / static_cast<float>(sensor_y);
+
+    correls_.clear();
+    correls_.reserve(block_n*sensor_n*3);
+    for (size_t by = 0; by < block_y; ++by) {
+      for (size_t bx = 0; bx < block_x; ++bx) {
+        for (size_t sy = 0; sy < sensor_y; ++sy) {
+          for (size_t sx = 0; sx < sensor_x; ++sx) {
+            const auto xf =
+                block_w*static_cast<float>(bx) +
+                sensor_interval_x*(static_cast<float>(sx)+.5f);
+            const auto yf =
+                block_h*static_cast<float>(by) +
+                sensor_interval_y*(static_cast<float>(sy)+.5f);
+
+            for (size_t i = 0; i < 3; ++i) {
+              const auto samp = data->FetchSamples(
+                  static_cast<size_t>(time_),
+                  static_cast<size_t>(dur_),
+                  xf, yf, i);
+
+              blky_sensor_t sensor = {};
+              blky_sensor_feed(&sensor, samp.data(), samp.size());
+              correls_.push_back(sensor.correl);
+            }
+          }
+        }
+      }
+    }
+
+    correl_avg_.clear();
+    correl_avg_.reserve(block_n);
+    correl_var_.clear();
+    correl_var_.reserve(block_n);
+    correl_avg_filtered_.clear();
+    correl_avg_filtered_.reserve(block_n);
+
+    const auto correls_per_block = sensor_n*3;
+    for (size_t bi = 0; bi < block_n; ++bi) {
+      const size_t st = bi*correls_per_block;
+      const size_t ed = st + correls_per_block;
+
+      double sum = 0;
+      for (size_t i = st; i < ed; ++i) {
+        sum += std::abs(correls_[i]);
+      }
+
+      const auto avg = sum / static_cast<double>(correls_per_block);
+      correl_avg_.push_back(avg);
+
+      double var = 0;
+      for (size_t i = st; i < ed; ++i) {
+        const auto diff = std::abs(correls_[i]) - avg;
+        var += diff*diff;
+      }
+      var /= static_cast<double>(correls_per_block);
+      correl_var_.push_back(var);
+
+      if (var < var_thresh_) {
+        correl_avg_filtered_.push_back(avg);
+      } else {
+        correl_avg_filtered_.push_back(-1);
+      }
+    }
+
+  } catch (const char* msg) {
+    msg_ = msg;
+  }
+};
+
+}  // namespace pg

playground/encoder.cc

@@ -0,0 +1,211 @@
+#include "app.hh"
+#include "input.hh"
+
+#include <imgui.h>
+#include <imgui_stdlib.h>
+
+
+namespace pg {
+namespace {
+
+class Encoder final : public App, public Input::Data {
+ public:
+  static inline TypeInfo kType = TypeInfo::Create<Encoder>("Encoder");
+
+  Encoder() noexcept : Data("encoder") {
+  }
+
+  void Update() noexcept override {
+    const auto id = std::to_string(index())+" Encoder | "+
+        std::to_string(reinterpret_cast<uintptr_t>(this));
+
+    if (ImGui::Begin(id.c_str())) {
+      ImGui::DragInt("src", &src_, 1, 0, 1024);
+      ImGui::DragInt("w", &w_, 1, 1, 1024);
+      ImGui::DragInt("h", &h_, 1, 1, 1024);
+      ImGui::DragInt2("block num", block_num_, 1, 1, 1024);
+      ImGui::DragInt("utime", &utime_, 1, 5, 64);
+      ImGui::InputTextMultiline("feats", &feats_);
+      if (ImGui::Button("encode")) {
+        Encode();
+      }
+    }
+    ImGui::End();
+  }
+
+  Frame Fetch(size_t n) noexcept override {
+    if (n >= frames_) return {};
+    const auto w = static_cast<size_t>(w_);
+    const auto h = static_cast<size_t>(h_);
+    return { .w = w, .h = h, .rgba = buf_.data()+4*w*h*n };
+  }
+  size_t frames() noexcept override {
+    return frames_;
+  }
+
+ private:
+  int src_ = 0;
+  int w_ = 128, h_ = 128;
+  int block_num_[2] = {2, 2};
+  int utime_ = 10;
+  std::string feats_;
+
+  std::vector<uint8_t> buf_;
+  size_t frames_ = 0;
+
+
+  struct Feat final {
+   public:
+    size_t block;
+    char type;  // C: const, L: linear
+  };
+
+  void Encode() {
+    const auto src = Input::instance().slots(static_cast<size_t>(src_));
+    if (!src) throw "missing input";
+
+    const auto feats = ParseFeats(feats_);
+
+    const auto w     = static_cast<size_t>(w_);
+    const auto h     = static_cast<size_t>(h_);
+    const auto utime = static_cast<size_t>(utime_);
+    buf_.resize(w*h*4*feats.size()*utime);
+    frames_ = utime * feats.size();
+
+    const auto b_num_x = static_cast<size_t>(block_num_[0]);
+    const auto b_num_y = static_cast<size_t>(block_num_[1]);
+
+    const auto dst_bw   = w/b_num_x;
+    const auto dst_bh   = h/b_num_y;
+    std::vector<std::pair<uint8_t, uint8_t>> block_rgba(dst_bw*dst_bh*4);
+
+    const auto dst_bw_f = static_cast<float>(dst_bw);
+    const auto dst_bh_f = static_cast<float>(dst_bh);
+
+    auto dst_ptr = buf_.data();
+    for (size_t t = 0; t < feats.size(); ++t) {
+      const auto& f = feats[t];
+
+      const auto bx = f.block%b_num_x;
+      const auto by = f.block/b_num_x;
+      if (by < b_num_y) {
+        {
+          // get pixels in start of the feature
+          const auto frame = src->Fetch(t*utime);
+          if (!frame.rgba) throw "got an empty frame";
+
+          const auto bw    = frame.w / b_num_x;
+          const auto bh    = frame.h / b_num_y;
+          const auto bw_f  = static_cast<float>(bw);
+          const auto bh_f  = static_cast<float>(bh);
+          const auto off_x = bw*bx;
+          const auto off_y = bh*by;
+
+          auto itr = block_rgba.begin();
+          for (size_t y = 0; y < dst_bh; ++y) {
+            const auto yf = static_cast<float>(y) / dst_bh_f;
+            for (size_t x = 0; x < dst_bw; ++x) {
+              const auto xf = static_cast<float>(x) / dst_bw_f;
+
+              const auto srcx = static_cast<size_t>(bw_f*xf) + off_x;
+              const auto srcy = static_cast<size_t>(bh_f*yf) + off_y;
+
+              auto ptr = frame.rgba + 4*(frame.w*srcy+srcx);
+              for (size_t i = 0; i < 4; ++i) {
+                itr->first  = *ptr;
+                itr->second = *ptr;
+                ++itr, ++ptr;
+              }
+            }
+          }
+        }
+        if (f.type == 'L') {
+          // get pixels in end frame if it's linear
+          const auto frame = src->Fetch((t+1)*utime);
+          if (!frame.rgba) throw "got an empty frame";
+
+          const auto bw    = frame.w / b_num_x;
+          const auto bh    = frame.h / b_num_y;
+          const auto bw_f  = static_cast<float>(bw);
+          const auto bh_f  = static_cast<float>(bh);
+          const auto off_x = bw*bx;
+          const auto off_y = bh*by;
+
+          auto itr = block_rgba.begin();
+          for (size_t y = 0; y < dst_bh; ++y) {
+            const auto yf = static_cast<float>(y) / dst_bh_f;
+            for (size_t x = 0; x < dst_bw; ++x) {
+              const auto xf = static_cast<float>(x) / dst_bw_f;
+
+              const auto srcx = static_cast<size_t>(bw_f*xf) + off_x;
+              const auto srcy = static_cast<size_t>(bh_f*yf) + off_y;
+
+              auto ptr = frame.rgba + 4*(frame.w*srcy+srcx);
+              for (size_t i = 0; i < 4; ++i) {
+                (itr++)->second = *(ptr++);
+              }
+            }
+          }
+        }
+      }
+
+      const auto bx_off = dst_bw*bx;
+      const auto by_off = dst_bh*by;
+
+      for (size_t bt = 0; bt < utime; ++bt) {
+        const auto bt_f = static_cast<float>(bt) / static_cast<float>(utime);
+
+        const auto frame = src->Fetch(t*utime+bt);
+        if (!frame.rgba) throw "got an empty frame";
+        for (size_t y = 0; y < w; ++y) {
+          for (size_t x = 0; x < h; ++x) {
+            const bool is_target =
+                bx_off <= x && x < bx_off+dst_bw && by_off <= y && y < by_off+dst_bh;
+            if (is_target) {
+              const auto* ptr = &block_rgba[4*((x-bx_off)+(y-by_off)*dst_bw)];
+              for (size_t i = 0; i < 4; ++i) {
+                const auto a = static_cast<float>(ptr->second-ptr->first);
+                const auto b = static_cast<float>(ptr->first);
+                *(dst_ptr++) = static_cast<uint8_t>(bt_f*a+b);
+                ++ptr;
+              }
+            } else {
+              const auto xf = static_cast<float>(x) / static_cast<float>(w);
+              const auto yf = static_cast<float>(y) / static_cast<float>(h);
+
+              const auto src_x = static_cast<size_t>(static_cast<float>(frame.w)*xf);
+              const auto src_y = static_cast<size_t>(static_cast<float>(frame.h)*yf);
+
+              auto ptr = frame.rgba + 4*(src_x+src_y*frame.w);
+              for (size_t i = 0; i < 4; ++i) {
+                *(dst_ptr++) = (*ptr++);
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+  static std::vector<Feat> ParseFeats(const std::string& script) {
+    std::vector<Feat> ret;
+
+    const char* begin = &*script.c_str();
+    const char* end   = begin + script.size();
+
+    size_t time = 0;
+    for (auto itr = begin; itr < end; ++itr, ++time) {
+      char* end;
+      const auto idx = std::strtol(&*itr, &end, 0);
+      if (itr == end) throw "invalid format: expected block index";
+      if (*end != 'C' && *end != 'L') {
+        throw "invalid format: unknown type";
+      }
+      itr = end+1;
+      ret.push_back(Feat { .block = static_cast<size_t>(idx), .type = *end, });
+    }
+    return ret;
+  }
+};
+
+}
+}  // namespace pg

playground/exporter.cc

@@ -0,0 +1,59 @@
+#include <imgui.h>
+#include <imgui_stdlib.h>
+
+#define STB_IMAGE_WRITE_IMPLEMENTATION
+#include <stb_image_write.h>
+
+#include "app.hh"
+#include "input.hh"
+
+
+namespace pg {
+namespace {
+
+class Exporter final : public App {
+ public:
+  static inline TypeInfo kType = TypeInfo::Create<Exporter>("Exporter");
+
+  Exporter() noexcept {
+  }
+
+  void Update() noexcept override {
+    const auto id = " Exporter | "+
+        std::to_string(reinterpret_cast<uintptr_t>(this));
+
+    if (ImGui::Begin(id.c_str())) {
+      ImGui::DragInt("src", &src_);
+      ImGui::InputText("dir", &dir_);
+      if (ImGui::Button("export")) {
+        Export();
+      }
+    }
+    ImGui::End();
+  }
+
+ private:
+  int src_ = 0;
+  std::string dir_;
+
+  void Export() noexcept {
+    const auto src =
+        Input::instance().slots(static_cast<size_t>(src_));
+
+    const auto n = src->frames();
+    for (size_t i = 0; i < n; ++i) {
+      const auto f = src->Fetch(i);
+      if (!f.rgba) break;
+
+      const auto path = dir_+"/"+std::to_string(i)+".png";
+      stbi_write_png(
+          path.c_str(),
+          static_cast<int>(f.w), static_cast<int>(f.h), 4,
+          f.rgba, static_cast<int>(f.w*4));
+    }
+  }
+};
+
+}
+}  // namespace pg
+

playground/extractor.cc

@@ -0,0 +1,130 @@
+extern "C" {
+#include <liblocky.h>
+}
+
+#include <sstream>
+#include <string>
+
+#include <imgui.h>
+#include <imgui_stdlib.h>
+#include <implot.h>
+
+#include "app.hh"
+#include "input.hh"
+
+
+namespace pg {
+
+class Extractor final : public App {
+ public:
+  static inline TypeInfo kType = TypeInfo::Create<Extractor>("Extractor");
+
+  Extractor() noexcept {
+  }
+
+  void Update() noexcept override {
+    const auto em = ImGui::GetFontSize();
+    const auto id = " Extractor | "+
+        std::to_string(reinterpret_cast<uintptr_t>(this));
+
+    ImGui::SetNextWindowSize({32*em, 32*em}, ImGuiCond_Once);
+    if (ImGui::Begin(id.c_str())) {
+      ImGui::DragInt("src", &src_, 1, 0, 1024);
+      ImGui::DragInt2("block num", block_num_, 1, 1, 1024);
+      ImGui::DragInt2("sensor num", sensor_num_, 1, 1, 1024);
+      ImGui::DragInt("utime", &utime_, 1, 5, 64);
+      ImGui::DragInt("len", &len_, 1, 1, 1024);
+      ImGui::DragFloat("correl max var", &correl_max_var_, 1e-3f, 0, 1);
+      ImGui::DragFloat("correl min avg", &correl_min_avg_, 1e-3f, 0, 1);
+      ImGui::DragFloat("feature min probability", &feat_min_prob_, 1e-3f, 0, 1);
+      if (ImGui::Button("extract")) {
+        Extract();
+      }
+
+      time_ = std::clamp(time_, 0, len_-1);
+      ImGui::SliderInt("time", &time_, 0, len_-1);
+
+      const auto offset = static_cast<size_t>(block_num_[0]*block_num_[1]*time_);
+      if (probs_.size() > offset) {
+        auto avail = ImGui::GetContentRegionAvail();
+        avail.x -= 1*em;
+        avail.y -= 1*em;
+        if (ImPlot::BeginPlot("probs", avail)) {
+          ImPlot::SetupAxisLimits(ImAxis_X1, 0, 1);
+          ImPlot::SetupAxisLimits(ImAxis_Y1, 0, 1);
+
+          ImPlot::PlotHeatmap(
+              "probs", probs_.data()+offset, block_num_[1], block_num_[0],
+              0, 1, "%.3f");
+          ImPlot::EndPlot();
+        }
+      } else {
+        ImGui::TextUnformatted("no data");
+      }
+    }
+    ImGui::End();
+  }
+
+ private:
+  int src_ = 0;
+  int block_num_[2] = {2, 2};
+  int sensor_num_[2] = {2, 2};
+  int utime_ = 10;
+  int len_ = 1;
+
+  float correl_max_var_ = 1e-3f;
+  float correl_min_avg_ = 0.8f;
+  float feat_min_prob_  = .5f;
+
+  int time_;
+
+  std::vector<double> probs_;
+
+
+  void Extract() noexcept {
+    const auto src = Input::instance().slots(static_cast<size_t>(src_));
+    if (!src) return;
+
+    const uint64_t dur = static_cast<uint64_t>(utime_ * len_);
+
+    blky_extractor_t ex = {};
+
+    ex.block_num_x = static_cast<uint32_t>(block_num_[0]);
+    ex.block_num_y = static_cast<uint32_t>(block_num_[1]);
+    ex.sensor_num_block_x = static_cast<uint32_t>(sensor_num_[0]);
+    ex.sensor_num_block_y = static_cast<uint32_t>(sensor_num_[1]);
+    ex.samples_per_pix = 3;
+    ex.pix_stride = 4;
+    ex.utime = static_cast<uint32_t>(utime_);
+    ex.correl_max_var = static_cast<double>(correl_max_var_);
+    ex.correl_min_avg = static_cast<double>(correl_min_avg_);
+    blky_extractor_init(&ex);
+
+    static const double verts[] = {
+      0, 0,
+      0, 1,
+      1, 1,
+      1, 0,
+    };
+
+    const uint64_t block_num = ex.block_num_x * ex.block_num_y;
+    probs_.clear();
+    probs_.reserve(block_num*static_cast<uint64_t>(len_));
+    for (uint64_t t = 0; t < dur; ++t) {
+      const auto f = src->Fetch(t);
+      if (!f.rgba) break;
+
+      const bool pop = blky_extractor_feed(
+          &ex, f.rgba,
+          static_cast<uint32_t>(f.w),
+          static_cast<uint32_t>(f.h),
+          verts);
+      if (pop) {
+        probs_.insert(probs_.end(), ex.probs, ex.probs+block_num);
+      }
+    }
+    blky_extractor_deinit(&ex);
+  }
+};
+
+}  // namespace pg

playground/initiator.hh

@@ -0,0 +1,39 @@
+#pragma once
+
+#include <memory>
+#include <vector>
+
+#include <imgui.h>
+
+#include "app.hh"
+
+
+namespace pg {
+
+class Initiator final : App {
+ public:
+  Initiator() noexcept {
+  }
+
+  void Update() noexcept override {
+    if (ImGui::BeginMainMenuBar()) {
+      if (ImGui::BeginMenu("Apps")) {
+        for (const auto& type : registry()) {
+          if (ImGui::MenuItem(type.first.c_str())) {
+            apps_.push_back(type.second->Create());
+          }
+        }
+        ImGui::EndMenu();
+      }
+      ImGui::EndMainMenuBar();
+    }
+    for (auto& app : apps_) {
+      app->Update();
+    }
+  }
+
+ private:
+  std::vector<std::unique_ptr<App>> apps_;
+};
+
+}  // namespace pg

playground/input.cc

@@ -0,0 +1,52 @@
+#include "input.hh"
+
+#include <imgui.h>
+#include <imgui_stdlib.h>
+
+
+namespace pg {
+
+Input Input::instance_;
+
+
+size_t Input::Add(Data* data) noexcept {
+  slots_.push_back(data);
+  return slots_.size()-1;
+}
+
+void Input::Update() noexcept {
+  if (ImGui::BeginMainMenuBar()) {
+    if (ImGui::BeginMenu("Input")) {
+      for (size_t i = 0; i < slots_.size(); ++i) {
+        const auto str =
+            std::to_string(i)+". "+slots_[i]->name();
+        ImGui::MenuItem(str.c_str());
+      }
+      ImGui::EndMenu();
+    }
+    ImGui::EndMainMenuBar();
+  }
+}
+
+
+std::vector<float> Input::Data::FetchSamples(
+    size_t st, size_t dur, float xf, float yf, size_t offset) noexcept {
+  std::vector<float> ret(dur);
+  for (size_t i = 0; i < dur; ++i) {
+    const auto frame = Fetch(i+st);
+    if (!frame.rgba) continue;
+
+    const auto fwf = static_cast<float>(frame.w);
+    const auto fhf = static_cast<float>(frame.h);
+    const auto fhi = static_cast<intmax_t>(frame.h);
+
+    const auto x = static_cast<intmax_t>(xf*fwf);
+    const auto y = static_cast<intmax_t>(yf*fhf);
+
+    const auto v = frame.rgba + 4*(y*fhi+x);
+    ret[i] = static_cast<float>(v[offset]) / UINT8_MAX;
+  }
+  return ret;
+}
+
+}  // namespace pg

playground/input.hh

@@ -0,0 +1,71 @@
+#pragma once
+
+#include <string>
+#include <string_view>
+#include <vector>
+
+#include "app.hh"
+
+
+namespace pg {
+
+class Input final : App {
+ public:
+  class Data;
+
+  static Input& instance() noexcept { return instance_; }
+
+  Input() = default;
+  Input(const Input&) = delete;
+  Input(Input&&) = delete;
+  Input& operator=(const Input&) = delete;
+  Input& operator=(Input&&) = delete;
+
+  size_t Add(Data*) noexcept;
+  void Update() noexcept override;
+
+  Data* slots(size_t idx) noexcept {
+    return idx < slots_.size()? slots_[idx]: nullptr;
+  }
+
+ private:
+  static Input instance_;
+
+  std::vector<Data*> slots_;
+};
+
+class Input::Data {
+ public:
+  struct Frame {
+   public:
+    size_t w, h;
+    const uint8_t* rgba;
+  };
+
+  Data() = delete;
+  Data(std::string_view name) noexcept : name_(name) {
+    idx_ = Input::instance().Add(this);
+  }
+  virtual ~Data() = default;
+  Data(const Data&) = delete;
+  Data(Data&&) = delete;
+  Data& operator=(const Data&) = delete;
+  Data& operator=(Data&&) = delete;
+
+  virtual Frame Fetch(size_t) noexcept = 0;
+  virtual size_t frames() noexcept = 0;
+
+  size_t index() const noexcept { return idx_; }
+
+  std::vector<float> FetchSamples(
+      size_t st, size_t dur, float x, float y, size_t offset) noexcept;
+
+  const std::string& name() const noexcept { return name_; }
+
+ private:
+  std::string name_;
+
+  size_t idx_;
+};
+
+}  // namespace pg

playground/input_gradient.cc

@@ -0,0 +1,52 @@
+#include "app.hh"
+#include "input.hh"
+
+#include <array>
+#include <cstring>
+
+#include <imgui.h>
+
+
+namespace pg {
+namespace {
+
+class Gradient final : public App, public Input::Data {
+ public:
+  static inline TypeInfo kType = TypeInfo::Create<Gradient>("Input_Gradient");
+
+  Gradient() noexcept : Data("gradient") {
+  }
+
+  void Update() noexcept override {
+    const auto id = std::to_string(index())+" Input_Gradient | "+
+        std::to_string(reinterpret_cast<uintptr_t>(this));
+
+    if (ImGui::Begin(id.c_str())) {
+      ImGui::DragInt("w", &w_, 1, 1, 1024);
+      ImGui::DragInt("h", &h_, 1, 1, 1024);
+      ImGui::DragInt("dur", &dur_, 1, 1, 1024);
+    }
+    ImGui::End();
+  }
+
+  Frame Fetch(size_t n) noexcept override {
+    const auto w = static_cast<size_t>(w_);
+    const auto h = static_cast<size_t>(h_);
+    buf_.resize(w*h*4);
+
+    const auto f = static_cast<float>(n)/static_cast<float>(dur_);
+    std::memset(buf_.data(), static_cast<uint8_t>(f*UINT8_MAX), buf_.size());
+    return Frame {.w = w, .h = h, .rgba = buf_.data()};
+  }
+  size_t frames() noexcept override {
+    return static_cast<size_t>(dur_);
+  }
+
+ private:
+  int w_ = 100, h_ = 100;
+  int dur_ = 100;
+  std::vector<uint8_t> buf_;
+};
+
+}
+}  // namespace pg

playground/input_imgseq.cc

@@ -0,0 +1,96 @@
+#include <string>
+#include <numeric>
+#include <vector>
+
+#include <imgui.h>
+#include <imgui_stdlib.h>
+
+#include "app.hh"
+#include "input.hh"
+
+#define STB_IMAGE_IMPLEMENTATION
+#include <stb_image.h>
+
+
+namespace pg {
+namespace {
+
+class ImgSeq final : public App, public Input::Data {
+ public:
+  static inline TypeInfo kType = TypeInfo::Create<ImgSeq>("Input_ImgSeq");
+
+  ImgSeq() noexcept : Data("imgseq") {
+  }
+  ~ImgSeq() noexcept {
+    DropCache();
+  }
+  void DropCache() noexcept {
+    for (auto& frame : frames_) {
+      stbi_image_free(const_cast<uint8_t*>(frame.second.rgba));
+    }
+    frames_.clear();
+  }
+
+  void Update() noexcept override {
+    const auto id = std::to_string(index())+" Input_ImgSeq | "+
+        std::to_string(reinterpret_cast<uintptr_t>(this));
+
+    if (ImGui::Begin(id.c_str())) {
+      bool mod = false;
+      mod |= ImGui::InputText("dir", &path_);
+      mod |= ImGui::DragInt("start", &st_, 1, 0, 1024);
+      mod |= ImGui::DragInt("end",   &ed_, 1, 0, 1024);
+      if (mod) DropCache();
+
+      if (msg_.size() > 0) {
+        ImGui::TextUnformatted(msg_.c_str());
+      }
+    }
+    ImGui::End();
+  }
+
+  Frame Fetch(size_t n) noexcept override
+  try {
+    n += static_cast<size_t>(st_);
+    if (n >= static_cast<size_t>(ed_)) {
+      throw "frame number out of range";
+    }
+
+    auto itr = frames_.find(n);
+    if (itr != frames_.end()) return itr->second;
+
+    const auto fname = path_+"/"+std::to_string(n)+".png";
+
+    int w, h, comp;
+    uint8_t* buf = stbi_load(fname.c_str(), &w, &h, &comp, 4);
+    if (buf == nullptr) {
+      throw stbi_failure_reason();
+    }
+
+    const Frame ret = {
+      .w = static_cast<size_t>(w),
+      .h = static_cast<size_t>(h),
+      .rgba = buf,
+    };
+    frames_[n] = ret;
+    return ret;
+  } catch (const char* msg) {
+    msg_ = msg;
+    return {.w = 0, .h = 0, .rgba = nullptr};
+  }
+  size_t frames() noexcept override {
+    if (ed_ <= st_) return 0;
+    return static_cast<size_t>(ed_-st_);
+  }
+
+ private:
+  std::string path_;
+  int st_ = 1, ed_ = 100;
+
+  std::string msg_;
+
+  std::unordered_map<size_t, Frame> frames_;
+};
+
+}
+}  // namespace pg

playground/input_noise.cc

@@ -0,0 +1,103 @@
+#include "app.hh"
+#include "input.hh"
+
+#include <vector>
+
+#include <imgui.h>
+
+
+namespace pg {
+namespace {
+
+class Noise final : public App, public Input::Data {
+ public:
+  static inline TypeInfo kType = TypeInfo::Create<Noise>("Input_Noise");
+
+  static uint32_t xorshift(uint32_t x) noexcept {
+    x ^= x<<13;
+    x ^= x>>17;
+    x ^= x<<5;
+    return x;
+  }
+
+  Noise() noexcept : Data("random noise") {
+  }
+
+  void Update() noexcept override {
+    const auto id = std::to_string(index())+" Input_Noise | "+
+        std::to_string(reinterpret_cast<uintptr_t>(this));
+
+    if (ImGui::Begin(id.c_str())) {
+      ImGui::DragInt("src", &src_, 1, 0, 1024);
+
+      ImGui::DragInt("w", &w_, 1, 1, 1024);
+      ImGui::DragInt("h", &h_, 1, 1, 1024);
+
+      ImGui::DragFloat("level", &level_, 1e-4f, 0, 1);
+      ImGui::DragInt("seed", &seed_);
+    }
+    ImGui::End();
+  }
+
+  Frame Fetch(size_t n) noexcept override {
+    auto src = Input::instance().slots(static_cast<size_t>(src_));
+    if (!src) return {};
+
+    const auto w = static_cast<size_t>(w_);
+    const auto h = static_cast<size_t>(h_);
+
+    buf_.resize(w*h*4);
+    uint8_t* buf = buf_.data();
+
+    auto srcf = src->Fetch(n);
+    for (size_t y = 0; y < h; ++y) {
+      const auto yf = static_cast<float>(y)/static_cast<float>(h);
+      for (size_t x = 0; x < w; ++x) {
+        const auto xf = static_cast<float>(x)/static_cast<float>(w);
+        const auto srcx  = static_cast<float>(srcf.w) * xf;
+        const auto srcy  = static_cast<float>(srcf.h) * yf;
+        const auto srcxi = static_cast<size_t>(srcx);
+        const auto srcyi = static_cast<size_t>(srcy);
+
+        auto v = srcf.rgba + 4*(srcxi+srcyi*srcf.w);
+        for (size_t i = 0; i < 4; ++i) {
+          *(buf++) = TryAttack(n, x, y, i, *(v++));
+        }
+      }
+    }
+    return Frame { .w = w, .h = h, .rgba = buf_.data(), };
+  }
+  size_t frames() noexcept override {
+    auto src = Input::instance().slots(static_cast<size_t>(src_));
+    if(!src) return 0;
+    return src->frames();
+  }
+
+  uint8_t TryAttack(
+      size_t n, size_t x, size_t y, size_t i, uint8_t v) const noexcept {
+    if (i == 3) return v;
+
+    const auto s = static_cast<uint32_t>(seed_);
+    const auto seed = static_cast<uint32_t>(s*s*s*n + s*s*x + s*y + i);
+
+    const auto rand = xorshift(seed);
+    if (rand%100 < 50) {
+      const auto t = static_cast<float>(rand%100)/100.f*level_;
+      const auto a = static_cast<uint8_t>(t*UINT8_MAX);
+      return v+a;
+    } else {
+      return v;
+    }
+  }
+
+ private:
+  int src_ = 0;
+  int w_ = 100, h_ = 100;
+  float level_ = 0.01f;
+  int seed_ = 1234;
+
+  std::vector<uint8_t> buf_;
+};
+
+}
+}  // namespace pg

playground/main.cc

@@ -0,0 +1,102 @@
+#include <iostream>
+
+#include <GL/glew.h>
+
+#include <imgui.h>
+#include <imgui_impl_glfw.h>
+#include <imgui_impl_opengl3.h>
+#include <implot.h>
+
+// To prevent conflicts caused by fucking windows.h, include GLFW last.
+#include <GLFW/glfw3.h>
+
+#include "initiator.hh"
+#include "input.hh"
+
+
+#if defined(_MSC_VER) && (_MSC_VER >= 1900) && !defined(IMGUI_DISABLE_WIN32_FUNCTIONS)
+# pragma comment(lib, "legacy_stdio_definitions")
+#endif
+
+
+int main(int, char**) {
+  // init display
+  glfwSetErrorCallback(
+      [](int, const char* msg) {
+        std::cout << "GLFW error: " << msg << std::endl;
+      });
+  if (!glfwInit()) return 1;
+
+  GLFWwindow* window;
+  const char* glsl_version;
+  glfwWindowHint(GLFW_VISIBLE, GLFW_FALSE);
+# if defined(__APPLE__)
+    glsl_version = "#version 150";
+    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
+    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 2);
+    glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
+    glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
+# else
+    glsl_version = "#version 130";
+    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
+    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
+    glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
+    glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
+# endif
+  window = glfwCreateWindow(1280, 720, "playground", NULL, NULL);
+  if (window == NULL) return 1;
+  glfwMakeContextCurrent(window);
+  glfwSwapInterval(1);
+  if (glewInit() != GLEW_OK) return 1;
+
+  // init ImGUI
+  IMGUI_CHECKVERSION();
+  ImGui::CreateContext();
+  ImPlot::CreateContext();
+
+  auto& io = ImGui::GetIO();
+  io.IniFilename  = nullptr;
+  io.ConfigFlags |= ImGuiConfigFlags_DockingEnable;
+
+  ImGui::StyleColorsDark();
+  ImGui_ImplGlfw_InitForOpenGL(window, true);
+  ImGui_ImplOpenGL3_Init(glsl_version);
+
+  pg::Initiator initiator;
+  glfwShowWindow(window);
+
+  // main loop
+  while (!glfwWindowShouldClose(window)) {
+    // new frame
+    glfwPollEvents();
+    ImGui_ImplOpenGL3_NewFrame();
+    ImGui_ImplGlfw_NewFrame();
+    ImGui::NewFrame();
+
+    // update GUI state
+    initiator.Update();
+    pg::Input::instance().Update();
+
+    // render windows
+    ImGui::Render();
+
+    int w, h;
+    glfwGetFramebufferSize(window, &w, &h);
+    glViewport(0, 0, w, h);
+
+    glClear(GL_COLOR_BUFFER_BIT);
+    ImGui_ImplOpenGL3_RenderDrawData(ImGui::GetDrawData());
+
+    glfwSwapBuffers(window);
+  }
+  // teardown ImGUI
+  ImGui_ImplOpenGL3_Shutdown();
+  ImGui_ImplGlfw_Shutdown();
+  ImPlot::DestroyContext();
+  ImGui::DestroyContext();
+
+  // teardown display
+  glfwDestroyWindow(window);
+  glfwTerminate();
+  return 0;
+}

playground/player.cc

@@ -0,0 +1,161 @@
+#include <algorithm>
+#include <cassert>
+#include <string>
+#include <vector>
+
+#include <GL/glew.h>
+#include <imgui.h>
+
+#include "app.hh"
+#include "input.hh"
+
+
+namespace pg {
+namespace {
+
+class Player final : public App {
+ public:
+  static inline TypeInfo kType = TypeInfo::Create<Player>("Player");
+
+  static constexpr size_t kTexW = 1024;
+  static constexpr size_t kTexH = 1024;
+
+  Player() noexcept {
+    glGenTextures(1, &tex_);
+
+    glBindTexture(GL_TEXTURE_2D, tex_);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
+    glBindTexture(GL_TEXTURE_2D, 0);
+  }
+  ~Player() noexcept {
+    glDeleteTextures(1, &tex_);
+  }
+
+  void Update() noexcept override {
+    const auto em = ImGui::GetFontSize();
+    const auto id = "Player | "+
+        std::to_string(reinterpret_cast<uintptr_t>(this));
+
+    ImGui::SetNextWindowSize({32*em, 16*em}, ImGuiCond_Once);
+    if (ImGui::Begin(id.c_str())) {
+      bool mod = false;
+
+      ImGui::BeginGroup();
+      {
+        ImGui::PushItemWidth(4*em);
+        mod |= ImGui::DragInt("input_slot", &slot_, 1, 0, 1024);
+        mod |= ImGui::DragInt("time", &t_, 1, 0, 1024);
+
+        ImGui::Checkbox("auto increment", &auto_inc_);
+        if (auto_inc_) {
+          ++t_; mod = true;
+        }
+
+        ImGui::PopItemWidth();
+      }
+      ImGui::EndGroup();
+      ImGui::SameLine();
+
+      ImVec2 img_min = {0, 0};
+      ImVec2 img_max = {0, 0};
+      if (mod) UpdateTex();
+      const auto wf  = static_cast<float>(w_);
+      const auto hf  = static_cast<float>(h_);
+
+      ImGui::BeginGroup();
+      {
+        if (w_ && h_) {
+          const auto tex = (void*) (uintptr_t) tex_;
+          const auto z   = wf / ImGui::GetContentRegionAvail().x;
+          ImGui::Image(tex, {wf/z, hf/z}, {0, 0}, uv_);
+          img_min = ImGui::GetItemRectMin();
+          img_max = ImGui::GetItemRectMax();
+        }
+      }
+      ImGui::EndGroup();
+
+      if (msg_.size() > 0) {
+        ImGui::TextUnformatted(msg_.c_str());
+      } else if (w_ && h_) {
+        const auto m  = ImGui::GetMousePos();
+        const auto xf = (m.x-img_min.x)/(img_max.x-img_min.x);
+        const auto yf = (m.y-img_min.y)/(img_max.y-img_min.y);
+        ImGui::Text("x=%f, y=%f, xf=%f, yf=%f", xf*wf, yf*hf, xf, yf);
+      } else {
+        ImGui::TextUnformatted("no image shown");
+      }
+    }
+    ImGui::End();
+  }
+
+ private:
+  std::string msg_;
+
+  int slot_      = 0;
+  int t_         = 0;
+  bool auto_inc_ = false;
+
+  GLsizei w_ = 0, h_ = 0;
+  ImVec2 uv_;
+
+  GLuint tex_;
+  GLsizei texw_ = 0, texh_ = 0;
+
+
+  void UpdateTex() noexcept
+  try {
+    msg_ = "";
+
+    auto data = Input::instance().slots(static_cast<size_t>(slot_));
+    if (data == nullptr) {
+      throw "missing slot";
+    }
+    if (static_cast<size_t>(t_) >= data->frames()) {
+      throw "time out of range";
+    }
+    const auto frame = data->Fetch(static_cast<size_t>(t_));
+    if (!frame.rgba) {
+      throw "got an empty frame";
+    }
+    w_ = static_cast<GLsizei>(frame.w);
+    h_ = static_cast<GLsizei>(frame.h);
+
+    glBindTexture(GL_TEXTURE_2D, tex_);
+
+    const auto ptexw = texw_;
+    const auto ptexh = texh_;
+    texw_ = std::max(texw_, NextPowerOf2(w_));
+    texh_ = std::max(texh_, NextPowerOf2(h_));
+    if (texw_ != ptexw || texh_ != ptexh) {
+      glTexImage2D(
+          GL_TEXTURE_2D, 0, GL_RGBA, texw_, texh_,
+          0, GL_RGBA, GL_UNSIGNED_BYTE, nullptr);
+    }
+
+    glTexSubImage2D(
+        GL_TEXTURE_2D, 0, 0, 0, w_, h_,
+        GL_RGBA, GL_UNSIGNED_BYTE, frame.rgba);
+    uv_ = {
+      static_cast<float>(w_)/static_cast<float>(texw_),
+      static_cast<float>(h_)/static_cast<float>(texh_)};
+
+    glBindTexture(GL_TEXTURE_2D, 0);
+    assert(glGetError() == 0);
+  } catch (const char* msg) {
+    msg_ = msg;
+  }
+
+
+  template <typename I>
+  static I NextPowerOf2(I x) {
+    I y = 1;
+    while (y < x) y *= 2;
+    return y;
+  }
+};
+
+}
+}  // namespace pg

playground/sensor.cc

@@ -0,0 +1,192 @@
+extern "C" {
+# include <liblocky.h>
+}
+
+#include <algorithm>
+#include <cinttypes>
+#include <cmath>
+#include <numeric>
+#include <string>
+#include <vector>
+
+#include <imgui.h>
+#include <implot.h>
+
+#include "app.hh"
+#include "input.hh"
+
+
+namespace pg {
+namespace {
+
+class Sensor final : public App {
+ public:
+  static inline TypeInfo kType = TypeInfo::Create<Sensor>("Sensor");
+
+  Sensor() noexcept {
+  }
+
+  void Update() noexcept override {
+    const auto em = ImGui::GetFontSize();
+    const auto id = "Sensor | "+
+        std::to_string(reinterpret_cast<uintptr_t>(this));
+
+    ImGui::SetNextWindowSize({64*em, 24*em}, ImGuiCond_Once);
+    if (ImGui::Begin(id.c_str())) {
+      bool mod = false;
+      ImGui::BeginGroup();
+      {
+        ImGui::PushItemWidth(6*em);
+        mod |= ImGui::DragInt("input slot", &src_, 1, 0, 1024);
+        mod |= ImGui::DragInt("start", &start_, 1, 0, 1024);
+        mod |= ImGui::DragInt("dur", &dur_, 1, 1, BLKY_SENSOR_MAX_DUR);
+
+        ImGui::Spacing();
+        ImGui::DragFloat2("pos", pos_, .001f);
+        ImGui::DragInt("offset", &offset_, 1, 0, 3);
+        if (ImGui::Button("add")) {
+          mod = true;
+          data_.emplace_back(pos_[0], pos_[1], offset_);
+        }
+        if (ImGui::BeginListBox("##points")) {
+          for (size_t i = 0; i < data_.size(); ++i) {
+            const auto& data = data_[i];
+            const auto  name =
+                std::to_string(i)+". "+
+                std::to_string(data.pos[0])+","+
+                std::to_string(data.pos[1]);
+            ImGui::Selectable(name.c_str());
+            if (ImGui::BeginPopupContextItem()) {
+              if (ImGui::MenuItem("clear all")) {
+                mod = true;
+                data_.clear();
+              }
+              ImGui::EndPopup();
+            }
+            if (ImGui::IsItemHovered()) {
+              ImGui::SetTooltip(
+                  "pos   : %f,%f + %" PRIiMAX "\n"
+                  "avg   : %f\n"
+                  "var   : %f\n"
+                  "correl: %f",
+                  data.pos[0], data.pos[1], data.offset,
+                  data.avg, data.var, data.correl);
+            }
+          }
+          ImGui::EndListBox();
+        }
+
+        double correl_abs_avg = 0;
+        for (auto& data : data_) {
+          correl_abs_avg += fabs(data.correl);
+        }
+        correl_abs_avg /= static_cast<float>(data_.size());
+        ImGui::Text("avg(abs(r)) = %f", correl_abs_avg);
+
+        double correl_abs_var = 0;
+        for (auto& data : data_) {
+          const auto diff = fabs(data.correl)-correl_abs_avg;
+          correl_abs_var += diff*diff;
+        }
+        correl_abs_var /= static_cast<float>(data_.size());
+        ImGui::Text("var(abs(r)) = %f", correl_abs_var);
+
+        ImGui::PopItemWidth();
+        if (mod) {
+          try {
+            Calc();
+            msg_ = "";
+          } catch (const char* msg) {
+            msg_ = msg;
+          }
+        }
+      }
+      ImGui::EndGroup();
+      ImGui::SameLine();
+      ImGui::BeginGroup();
+      auto avail = ImGui::GetContentRegionAvail();
+      avail.x -= em;
+      avail.y -= em;
+
+      const ImVec2 size = {avail.x/2, avail.y};
+
+      if (ImPlot::BeginPlot("input value", size)) {
+        ImPlot::SetupAxis(ImAxis_X1, nullptr, ImPlotAxisFlags_AutoFit);
+        ImPlot::SetupAxisLimits(ImAxis_Y1, -0.1, 1.1);
+        for (size_t i = 0; i < data_.size(); ++i) {
+          const auto& d = data_[i];
+          ImPlot::PlotLine(
+              std::to_string(i).c_str(),
+              d.values.data(), static_cast<int>(d.values.size()), 1, start_);
+        }
+        ImPlot::EndPlot();
+      }
+      ImGui::SameLine();
+      if (ImPlot::BeginPlot("input histogram", size)) {
+        ImPlot::SetupAxisLimits(ImAxis_X1, 0, 1);
+        ImPlot::SetupAxis(ImAxis_Y1, nullptr, ImPlotAxisFlags_AutoFit);
+        for (size_t i = 0; i < data_.size(); ++i) {
+          const auto& d = data_[i];
+          ImPlot::PlotHistogram(
+              std::to_string(i).c_str(),
+              d.values.data(), static_cast<int>(d.values.size()));
+        }
+        ImPlot::EndPlot();
+      }
+      ImGui::EndGroup();
+      ImGui::TextUnformatted(msg_.c_str());
+    }
+    ImGui::End();
+  }
+
+ private:
+  int src_ = 0, start_ = 0, dur_ = 30;
+  float pos_[2];
+  int offset_;
+
+  struct Data final {
+   public:
+    float pos[2];
+    size_t offset;
+    std::vector<float> values;
+
+    double avg;
+    double var;
+    double cov;
+    double correl;
+
+    Data(float x, float y, int off) noexcept :
+        pos{x, y}, offset(static_cast<size_t>(off)) {
+    }
+  };
+  std::vector<Data> data_;
+  std::string msg_;
+
+  void Calc();
+};
+
+void Sensor::Calc() {
+  auto in = Input::instance().slots(static_cast<size_t>(src_));
+  if (!in) throw "missing slot";
+
+  if (dur_ == 0) throw "invalid time range";
+
+  for (auto& data : data_) {
+    const auto xf = std::clamp(data.pos[0], 0.f, 1.f);
+    const auto yf = std::clamp(data.pos[1], 0.f, 1.f);
+
+    const auto dur = static_cast<size_t>(dur_);
+    const auto st  = static_cast<size_t>(start_);
+    data.values = in->FetchSamples(st, dur, xf, yf, data.offset);
+
+    blky_sensor_t sensor = {};
+    blky_sensor_feed(&sensor, data.values.data(), data.values.size());
+    data.avg    = sensor.avg;
+    data.var    = sensor.var;
+    data.cov    = sensor.cov;
+    data.correl = sensor.correl;
+  }
+}
+
+}
+}  // namespace pg

thirdparty/CMakeLists.txt

@@ -0,0 +1,154 @@
+include(FetchContent)
+
+
+# ---- args ----
+# repository: https://github.com/Taywee/args
+# license   : MIT
+
+FetchContent_Declare(
+  args
+  URL "https://github.com/Taywee/args/archive/refs/tags/6.3.0.zip"
+)
+
+set(ARGS_BUILD_EXAMPLE   OFF)
+set(ARGS_BUILD_UNITTESTS OFF)
+FetchContent_MakeAvailable(args)
+
+
+# ---- GLEW ----
+# repository: https://github.com/Perlmint/glew-cmake
+# license   :  Modified BSD License, the Mesa 3-D License (MIT) and the Khronos License (MIT).
+
+FetchContent_Declare(
+  glew
+  URL "https://github.com/Perlmint/glew-cmake/archive/refs/tags/glew-cmake-2.2.0.zip"
+)
+FetchContent_MakeAvailable(glew)
+
+if (BLOCKY_STATIC)
+  add_library(glew ALIAS libglew_static)
+else()
+  add_library(glew ALIAS libglew_shared)
+endif()
+
+
+# ---- GLFW ----
+# repository: https://github.com/glfw/glfw
+# license   : zlib
+
+FetchContent_Declare(
+  glfw
+  URL "https://github.com/glfw/glfw/archive/refs/tags/3.3.4.zip"
+)
+
+set(GLFW_BUILD_EXAMPLES OFF CACHE BOOL "" FORCE)
+set(GLFW_BUILD_TESTS    OFF CACHE BOOL "" FORCE)
+set(GLFW_BUILD_DOCS     OFF CACHE BOOL "" FORCE)
+set(GLFW_INSTALL        OFF CACHE BOOL "" FORCE)
+
+FetchContent_MakeAvailable(glfw)
+
+
+# ---- ImGUI (docking branch) ----
+# repository: https://github.com/ocornut/imgui/
+# license   : MIT
+
+FetchContent_Declare(
+  imgui
+  URL "https://github.com/ocornut/imgui/archive/9b0c26b0b2adae3ccf66dc9552fae4945d735a0c.zip"
+)
+FetchContent_Populate(imgui)
+
+add_library(imgui)
+target_sources(imgui
+  PRIVATE
+    "${imgui_SOURCE_DIR}/imgui.cpp"
+    "${imgui_SOURCE_DIR}/imgui_demo.cpp"
+    "${imgui_SOURCE_DIR}/imgui_draw.cpp"
+    "${imgui_SOURCE_DIR}/imgui_internal.h"
+    "${imgui_SOURCE_DIR}/imgui_tables.cpp"
+    "${imgui_SOURCE_DIR}/imgui_widgets.cpp"
+    "${imgui_SOURCE_DIR}/backends/imgui_impl_glfw.cpp"
+    "${imgui_SOURCE_DIR}/backends/imgui_impl_opengl3.cpp"
+    "${imgui_SOURCE_DIR}/misc/cpp/imgui_stdlib.cpp"
+  PUBLIC
+    "${imgui_SOURCE_DIR}/imgui.h"
+    "${imgui_SOURCE_DIR}/imstb_rectpack.h"
+    "${imgui_SOURCE_DIR}/imstb_textedit.h"
+    "${imgui_SOURCE_DIR}/imstb_truetype.h"
+    "${imgui_SOURCE_DIR}/backends/imgui_impl_glfw.h"
+    "${imgui_SOURCE_DIR}/backends/imgui_impl_opengl3.h"
+    "${imgui_SOURCE_DIR}/misc/cpp/imgui_stdlib.h"
+)
+target_include_directories(imgui SYSTEM
+  PUBLIC
+    "${imgui_SOURCE_DIR}"
+    "${imgui_SOURCE_DIR}/backends"
+    "${imgui_SOURCE_DIR}/misc/cpp"
+)
+target_link_libraries(imgui
+  PRIVATE glfw
+)
+
+
+# ---- ImPlot ----
+# repository: https://github.com/epezent/implot
+# license   : MIT
+
+FetchContent_Declare(
+  implot
+  URL "https://github.com/epezent/implot/archive/refs/heads/master.zip"
+)
+FetchContent_Populate(implot)
+
+add_library(implot)
+target_link_libraries(implot
+  PRIVATE
+    imgui
+)
+target_include_directories(implot SYSTEM
+  PUBLIC
+    "${implot_SOURCE_DIR}"
+)
+target_sources(implot
+  PUBLIC
+    "${implot_SOURCE_DIR}/implot.h"
+    "${implot_SOURCE_DIR}/implot_internal.h"
+  PRIVATE
+    "${implot_SOURCE_DIR}/implot.cpp"
+    "${implot_SOURCE_DIR}/implot_items.cpp"
+)
+
+
+# ---- minimp4 ----
+# repository: https://github.com/lieff/minimp4
+# license   : CC0
+
+FetchContent_Declare(
+  minimp4
+  URL "https://github.com/lieff/minimp4/archive/4575afb4f69ace25a1a048e25cc86bf8c8d14f2b.zip"
+)
+FetchContent_Populate(minimp4)
+
+add_library(minimp4)
+target_include_directories(minimp4 PUBLIC SYSTEM ${minimp4_SOURCE_DIR})
+target_sources(minimp4
+  PUBLIC
+    "${minimp4_SOURCE_DIR}/minimp4.h"
+  PRIVATE
+    minimp4.c
+)
+
+
+# ---- stb ----
+FetchContent_Declare(
+  stb
+  URL "https://github.com/nothings/stb/archive/af1a5bc352164740c1cc1354942b1c6b72eacb8a.zip"
+)
+FetchContent_Populate(stb)
+
+add_library(stb INTERFACE)
+target_include_directories(stb SYSTEM
+  INTERFACE
+    "${stb_SOURCE_DIR}"
+)

thirdparty/minimp4.c

@@ -0,0 +1,2 @@
+#define MINIMP4_IMPLEMENTATION
+#include <minimp4.h>