aboutsummaryrefslogtreecommitdiff
path: root/core/sdf.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'core/sdf.cpp')
-rw-r--r--core/sdf.cpp360
1 files changed, 360 insertions, 0 deletions
diff --git a/core/sdf.cpp b/core/sdf.cpp
new file mode 100644
index 0000000..29d6f13
--- /dev/null
+++ b/core/sdf.cpp
@@ -0,0 +1,360 @@
+// This code contains NVIDIA Confidential Information and is disclosed to you
+// under a form of NVIDIA software license agreement provided separately to you.
+//
+// Notice
+// NVIDIA Corporation and its licensors retain all intellectual property and
+// proprietary rights in and to this software and related documentation and
+// any modifications thereto. Any use, reproduction, disclosure, or
+// distribution of this software and related documentation without an express
+// license agreement from NVIDIA Corporation is strictly prohibited.
+//
+// ALL NVIDIA DESIGN SPECIFICATIONS, CODE ARE PROVIDED "AS IS.". NVIDIA MAKES
+// NO WARRANTIES, EXPRESSED, IMPLIED, STATUTORY, OR OTHERWISE WITH RESPECT TO
+// THE MATERIALS, AND EXPRESSLY DISCLAIMS ALL IMPLIED WARRANTIES OF NONINFRINGEMENT,
+// MERCHANTABILITY, AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// Information and code furnished is believed to be accurate and reliable.
+// However, NVIDIA Corporation assumes no responsibility for the consequences of use of such
+// information or for any infringement of patents or other rights of third parties that may
+// result from its use. No license is granted by implication or otherwise under any patent
+// or patent rights of NVIDIA Corporation. Details are subject to change without notice.
+// This code supersedes and replaces all information previously supplied.
+// NVIDIA Corporation products are not authorized for use as critical
+// components in life support devices or systems without express written approval of
+// NVIDIA Corporation.
+//
+// Copyright (c) 2013-2016 NVIDIA Corporation. All rights reserved.
+
+#include "sdf.h"
+
+#include <vector>
+#include <float.h>
+#include <math.h>
+
+using namespace std;
+
+namespace
+{
+ inline float Sqr(float x) { return x*x; }
+ inline int Clamp(int x, int lower, int upper) { return min(max(lower, x), upper); }
+
+ uint32_t Sample(const uint32_t* image, uint32_t w, uint32_t h, int x, int y)
+ {
+ return image[Clamp(y, 0, h-1)*w + Clamp(x, 0, w-1)];
+ }
+
+ uint32_t Sample(const uint32_t* image, uint32_t w, uint32_t h, uint32_t d, int x, int y, int z)
+ {
+ return image[Clamp(z, 0, d-1)*w*h + Clamp(y, 0, h-1)*w + Clamp(x, 0, w-1)];
+ }
+
+ // returns true if point is on the surface
+ bool EdgeDetect(const uint32_t* img, uint32_t w, uint32_t h, int x, int y)
+ {
+ bool center = Sample(img, w, h, x, y) != 0;
+
+ for (int j=y-1; j <= y+1; ++j)
+ {
+ for (int i=x-1; i <= x+1; ++i)
+ {
+ if ((0 != Sample(img, w, h, i, j)) != center)
+ {
+ return true;
+ }
+ }
+ }
+
+ return false;
+ }
+
+
+ // returns true if point is on the surface
+ bool EdgeDetect(const uint32_t* img, uint32_t w, uint32_t h, uint32_t d, int x, int y, int z, float& dist)
+ {
+ bool center = Sample(img, w, h, d, x, y, z) != 0;
+
+ float minDist = FLT_MAX;
+
+ for (int k=z-1; k <= z+1; ++k)
+ {
+ for (int j=y-1; j <= y+1; ++j)
+ {
+ for (int i=x-1; i <= x+1; ++i)
+ {
+ if ((0 != Sample(img, w, h, d, i, j, k)) != center)
+ {
+ int dx = x-i;
+ int dy = y-j;
+ int dz = z-k;
+
+ minDist = min(sqrtf(float(dx*dx + dy*dy + dz*dz))*0.5f, minDist);
+ }
+ }
+ }
+ }
+
+ dist = minDist;
+
+ return minDist != FLT_MAX;
+ }
+}
+
+// 2D fast marching method (FMM J. Sethian. A fast marching level set method for monotonically advancing fronts. Proc. Natl. Acad. Sci., 93:1591�1595, 1996.)
+namespace
+{
+ struct Coord2D
+ {
+ int i, j;
+
+ float d;
+ int si, sj;
+
+ bool operator < (const Coord2D& c) const { return d > c.d; }
+ };
+}
+
+void MakeSDF(const uint32_t* img, uint32_t w, uint32_t h, float* output)
+{
+ const float scale = 1.0f / max(w, h);
+
+ std::vector<Coord2D> queue;
+
+ // find surface points
+ for (uint32_t y=0; y < h; ++y)
+ {
+ for (uint32_t x=0; x < w; ++x)
+ {
+ if (EdgeDetect(img, w, h, x, y))
+ {
+ Coord2D c = {(int)x, (int)y, 0.0f, (int)x, (int)y};
+ queue.push_back(c);
+ }
+
+ output[y*w + x] = FLT_MAX;
+ }
+ }
+
+ std::make_heap(queue.begin(), queue.end());
+
+ while (!queue.empty())
+ {
+ std::pop_heap(queue.begin(), queue.end());
+
+ Coord2D c = queue.back();
+ queue.pop_back();
+
+ // freeze coord if not already frozen
+ if (output[c.j*w + c.i] == FLT_MAX)
+ {
+ output[c.j*w + c.i] = c.d;
+
+ // update neighbours
+ int xmin = max(c.i-1, 0), xmax = min(c.i+1, int(w-1));
+ int ymin = max(c.j-1, 0), ymax = min(c.j+1, int(h-1));
+
+ for (int y=ymin; y <= ymax; ++y)
+ {
+ for (int x=xmin; x <= xmax; ++x)
+ {
+ if (c.i != x || c.j != y)
+ {
+ int dx = x-c.si;
+ int dy = y-c.sj;
+
+ // calculate distance to source coord
+ float d = sqrtf(float(dx*dx + dy*dy));
+
+ Coord2D newc = {x, y, d, c.si, c.sj};
+ queue.push_back(newc);
+ std::push_heap(queue.begin(), queue.end());
+ }
+ }
+ }
+ }
+ }
+
+ for (uint32_t y=0; y < h; ++y)
+ {
+ for (uint32_t x=0; x < w; ++x)
+ {
+ assert(output[y*w + x] < FLT_MAX);
+
+ // flip sign for interior
+ output[y*w + x] *= (img[y*w + x]?-1.0f:1.0f)*scale;
+ }
+ }
+}
+
+// 3D fast marching method (FMM J. Sethian. A fast marching level set method for monotonically advancing fronts. Proc. Natl. Acad. Sci., 93:1591�1595, 1996.)
+namespace
+{
+ struct Coord3D
+ {
+ int i, j, k;
+
+ float d;
+ int si, sj, sk;
+
+ bool operator < (const Coord3D& c) const { return d > c.d; }
+ };
+}
+
+void MakeSDF(const uint32_t* img, uint32_t w, uint32_t h, uint32_t d, float* output)
+{
+ const float scale = 1.0f / max(max(w, h), d);
+
+ std::vector<Coord3D> queue;
+
+ // find surface points
+ for (uint32_t z=0; z < d; ++z)
+ {
+ for (uint32_t y=0; y < h; ++y)
+ {
+ for (uint32_t x=0; x < w; ++x)
+ {
+ float dist;
+ if (EdgeDetect(img, w, h, d, x, y, z, dist))
+ {
+ Coord3D c = {(int)x, (int)y, (int)z, dist, (int)x, (int)y, (int)z};
+ queue.push_back(c);
+ }
+
+ output[z*w*h + y*w + x] = FLT_MAX;
+ }
+ }
+ }
+
+ // no occupied voxels so quit
+ if (queue.empty())
+ return;
+
+ std::make_heap(queue.begin(), queue.end());
+
+ while (!queue.empty())
+ {
+ std::pop_heap(queue.begin(), queue.end());
+
+ Coord3D c = queue.back();
+ queue.pop_back();
+
+ // freeze coord if not already frozen
+ if (output[c.k*w*h + c.j*w + c.i] == FLT_MAX)
+ {
+ output[c.k*w*h + c.j*w + c.i] = c.d;
+
+ // update neighbours
+ int xmin = max(c.i-1, 0), xmax = min(c.i+1, int(w-1));
+ int ymin = max(c.j-1, 0), ymax = min(c.j+1, int(h-1));
+ int zmin = max(c.k-1, 0), zmax = min(c.k+1, int(d-1));
+
+ for (int z=zmin; z <= zmax; ++z)
+ {
+ for (int y=ymin; y <= ymax; ++y)
+ {
+ for (int x=xmin; x <= xmax; ++x)
+ {
+ if ((c.i != x || c.j != y || c.k != z) && output[z*w*h + y*w + x] == FLT_MAX)
+ {
+ int dx = x-c.si;
+ int dy = y-c.sj;
+ int dz = z-c.sk;
+
+ // calculate distance to source coord
+ float d = sqrtf(float(dx*dx + dy*dy + dz*dz)) + output[c.sk*w*h + c.sj*w + c.si];
+
+ assert(d > 0.0f);
+
+ Coord3D newc = {x, y, z, d, c.si, c.sj, c.sk};
+
+ queue.push_back(newc);
+ std::push_heap(queue.begin(), queue.end());
+ }
+ }
+ }
+ }
+ }
+ }
+
+ for (uint32_t z=0; z < d; ++z)
+ {
+ for (uint32_t y=0; y < h; ++y)
+ {
+ for (uint32_t x=0; x < w; ++x)
+ {
+ assert(output[z*w*h + y*w + x] < FLT_MAX);
+
+ // flip sign for interior
+ output[z*w*h + y*w + x] *= (img[z*w*h + y*w + x]?-1.0f:1.0f)*scale;
+ }
+ }
+ }
+}
+
+
+
+
+/*
+// Brute-force 2D SDF generation
+
+void FindNeighbour(const uint32_t* image, uint32_t w, uint32_t h, uint32_t cx, uint32_t cy, uint32_t& i, uint32_t& j, float& d)
+ {
+ float minDistSq=FLT_MAX;
+
+ float fx = float(cx);
+ float fy = float(cy);
+
+ for (uint32_t y=0; y < h; ++y)
+ {
+ for (uint32_t x=0; x < w; ++x)
+ {
+ if ((x != cx || y != cy) && image[y*w + x])
+ {
+ float dSq = Sqr(fx-float(x)) + Sqr(fy-float(y));
+ if (dSq < minDistSq)
+ {
+ minDistSq = dSq;
+ i = x;
+ j = y;
+ }
+ }
+ }
+ }
+
+ d = sqrtf(minDistSq);
+ }
+
+
+// brute force
+void MakeSDF(const uint32_t* img, uint32_t w, uint32_t h, float* output)
+{
+ // find surface points
+ vector<uint32_t> surface(w*h);
+
+ for (uint32_t y=0; y < h; ++y)
+ {
+ for (uint32_t x=0; x < w; ++x)
+ {
+ if (EdgeDetect(img, w, h, x, y))
+ {
+ surface[y*w + x] = 1;
+ }
+ }
+ }
+
+ // brute force search
+ for (uint32_t y=0; y < h; ++y)
+ {
+ for (uint32_t x=0; x < w; ++x)
+ {
+ uint32_t i, j;
+ float d;
+ FindNeighbour(&surface[0], w, h, x, y, i, j, d);
+
+ // flip sign for pixels inside the shape
+ float sign = (img[y*w + x])?-1.0f:1.0f;
+
+ output[y*w + x] = d*sign/w;
+ }
+ }
+}
+*/
generated by cgit v1.2.3 (git 2.25.1) at 2026-03-11 07:22:44 +0000