New rendering (VoronoiSplat) -- incomplete.

* Added voronoi-like splatting to points: the same technique from
  Messias et al., (2014)
* It is now possible to change the projection technique during
  runtime (possible, but still requires some work)

1 files changed, 199 insertions, 0 deletions

diff --git a/skelft_core.cpp b/skelft_core.cpp
new file mode 100644
index 0000000..e4cec20
--- /dev/null
+++ b/skelft_core.cpp
@@ -0,0 +1,199 @@
+#include <cstdlib>
+#include <cmath>
+#include <cstring>
+#include <stack>
+#include <vector>
+#include <iostream>
+
+const float MYINFINITY = 1.0e7f;
+
+inline float SQR(float x) { return x * x; }
+inline int SQR(int x) { return x * x; }
+inline float INTERP(float a, float b, float t) { return a * (1 - t) + b * t; }
+
+struct Coord {
+    int i, j;
+    Coord(int i_, int j_)
+        : i(i_)
+        , j(j_){};
+    Coord() {}
+    int operator==(const Coord &c) const { return i == c.i && j == c.j; }
+};
+
+using namespace std;
+
+const float ALIVE = -10, NARROW_BAND = ALIVE + 1,
+            FAR_AWAY = ALIVE + 2; // Classification values for pixels in an
+                                  // image. The values are arbitrary.
+
+int skelft2DSize(int nx, int ny) {
+    int dx = (int) floor(
+        pow(2.0f, int(log(float(nx)) / log(2.0f) +
+                      1))); // Find minimal pow of 2 which fits the input image
+    int dy = (int) floor(pow(2.0f, int(log(float(ny)) / log(2.0f) + 1)));
+    int fboSize = max(dx, dy); // DT/FT/skeleton image size (pow of 2, should be
+                               // able to contain the input image)
+
+    return fboSize;
+}
+
+// Classify 'f' into inside, outside, and boundary pixels according to isovalue
+// 'iso'
+// Inside pixels are those with value > iso (if thr_upper==true), else those
+// with value < iso.
+template <class T>
+void classify(vector<Coord> &s, const T *f, int xm, int ym, int xM, int yM,
+              int size, T iso, bool thr_upper, float *siteParam,
+              bool scan_bot_to_top) {
+    for (int i = 0; i < size * size; ++i)
+        siteParam[i] = ALIVE; // First: all points are ALIVE, i.e., sites
+
+    for (int j = ym; j < yM; ++j)
+        for (int i = xm; i < xM; ++i) {
+            T fptr = f[i + j * size]; // Mark subset of points which are
+                                      // FAR_AWAY, i.e. NOT sites
+            if (((thr_upper && fptr >= iso) || (!thr_upper && fptr <= iso)))
+                siteParam[i + size * j] = FAR_AWAY;
+        }
+
+    int js, je, jd;
+    if (scan_bot_to_top) {
+        js = ym;
+        je = yM;
+        jd = 1;
+    } else {
+        js = yM - 1;
+        je = ym - 1;
+        jd = -1;
+    }
+
+    for (int j = js; j != je; j += jd)
+        for (int i = xm; i < xM; ++i) {
+            float &cv = siteParam[i + size * j];
+            if (cv == FAR_AWAY)
+                if (((i - 1 + size * j) > 0 &&
+                     siteParam[i - 1 + size * j] == ALIVE) ||
+                    siteParam[i + 1 + size * j] == ALIVE ||
+                    ((i + size * (j - 1)) > 0 &&
+                     siteParam[i + size * (j - 1)] == ALIVE) ||
+                    siteParam[i + size * (j + 1)] == ALIVE) {
+                    cv = NARROW_BAND;
+                    s.push_back(Coord(i, j));
+                }
+        }
+}
+
+// Encode input image (classified as inside/outside/boundary into a
+// boundary-length parameterization image, needed for skeleton computations.
+//
+float encodeBoundary(vector<Coord> &s, int xm, int ym, int xM, int yM,
+                     float *siteParam, int size) {
+    const float LARGE_BOUNDARY_VAL =
+        1000; // Value that the boundary-param will jump with between disjoint,
+              // self-connected, boundary pixel chains
+
+    float length = -LARGE_BOUNDARY_VAL + 1; // Compute the boundary length
+
+    while (s.size()) // Process all narrowband points, add their
+                     // arc-length-parameterization in siteParam
+    {
+        Coord pt = s[s.size() - 1];
+        s.pop_back(); // Pick narrowband-point, test if already numbered
+        if (siteParam[pt.i + size * pt.j] > 0)
+            continue;
+
+        int ci, cj;
+        float cc = 0; // Point pt not numbered
+        for (int dj = -1; dj < 2; ++dj) // Find min-numbered-neighbour of pt
+            for (int di = -1; di < 2; ++di) {
+                int ii = pt.i + di, jj = pt.j + dj;
+                if (ii < xm || ii >= xM || jj < ym || jj >= yM)
+                    continue;
+                float val = siteParam[ii + size * jj];
+                if (val < 0)
+                    continue;
+                if (cc && val > cc)
+                    continue;
+                ci = ii;
+                cj = jj;
+                cc = val;
+            }
+
+        float c = (!cc) ? length + LARGE_BOUNDARY_VAL
+                        : cc + sqrt(float((pt.i - ci) * (pt.i - ci) +
+                                          (pt.j - cj) * (pt.j - cj)));
+
+        siteParam[pt.i + size * pt.j] =
+            c; // Also store the sites' parameterization in siteParam[] for CUDA
+
+        if (c > length)
+            length = c; // Determine length of boundary, used for
+                        // wraparound-distance-computations
+
+        for (int dj = -1; dj < 2; ++dj)
+            for (int di = -1; di < 2; ++di) {
+                if (di == 0 && dj == 0)
+                    continue;
+                int ii = pt.i + di, jj = pt.j + dj;
+                if (ii >= xm && ii < xM && jj >= ym && jj < yM &&
+                    siteParam[ii + size * jj] == NARROW_BAND)
+                    s.push_back(Coord(ii, jj));
+            }
+    }
+
+    for (int i = 0, iend = size * size; i < iend;
+         ++i) // All ALIVE points are marked now as non-sites:
+        if (siteParam[i] < 0)
+            siteParam[i] = 0;
+
+    return length;
+}
+
+float skelft2DMakeBoundary(const float *f, int xm, int ym, int xM, int yM,
+                           float *siteParam, int size, float iso,
+                           bool thr_upper) {
+    bool scan_bot_to_top = false;
+
+    vector<Coord> s;
+    s.reserve((xM - xm) * (yM - ym));
+
+    classify(s, f, xm, ym, xM, yM, size, iso, thr_upper, siteParam,
+             scan_bot_to_top);
+    return encodeBoundary(s, xm, ym, xM, yM, siteParam, size);
+}
+
+float skelft2DMakeBoundary(const unsigned char *f, int xm, int ym, int xM,
+                           int yM, float *siteParam, int size, short iso,
+                           bool thr_upper) {
+    bool scan_bot_to_top = false;
+
+    vector<Coord> s;
+    s.reserve((xM - xm) * (yM - ym));
+
+    classify(s, f, xm, ym, xM, yM, size, (unsigned char) iso, thr_upper,
+             siteParam, scan_bot_to_top);
+    return encodeBoundary(s, xm, ym, xM, yM, siteParam, size);
+}
+
+void skelft2DSave(short *outputFT, int dx, int dy, const char *f) {
+    FILE *fp = fopen(f, "wb");
+    fprintf(fp, "P5 %d %d 255\n", dx, dy);
+
+    const int SIZE = 3000;
+    unsigned char buf[SIZE];
+
+    int size = dx * dy;
+    int bb = 0;
+    float range = max(dx, dy) / 255.0f;
+    for (short *v = outputFT, *vend = outputFT + size; v < vend; ++v) {
+        short val = *v;
+        buf[bb++] = (unsigned char) (val / range);
+        if (bb == SIZE) {
+            fwrite(buf, sizeof(unsigned char), SIZE, fp);
+            bb = 0;
+        }
+    }
+    if (bb)
+        fwrite(buf, sizeof(unsigned char), bb, fp);
+    fclose(fp);
+}