create new project

liblocky/CMakeLists.txt

@@ -0,0 +1,17 @@
+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
+    extractor.c
+    sensor.c
+    image.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/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,103 @@
+#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]);
+
+
+/* ---- 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);

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
+}