Skip to content

Instantly share code, notes, and snippets.

@oxyflour

oxyflour/inc-cuda.h

Created March 7, 2021 15:52
Show Gist options
  • Save oxyflour/4e9abd4a19f5d667cc98dbfc1697046f to your computer and use it in GitHub Desktop.
Save oxyflour/4e9abd4a19f5d667cc98dbfc1697046f to your computer and use it in GitHub Desktop.
Download ZIP
Raw
inc-cuda.h
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cuda_runtime.h>
#ifndef CUDA_H
#define CUDA_H
// https://stackoverflow.com/a/27992604
#ifdef __INTELLISENSE__
dim3 blockIdx;
dim3 blockDim;
dim3 threadIdx;
dim3 gridDim;
#define CU_DIM(grid, block)
#define CU_DIM_MEM(grid, block, bytes)
#define CU_DIM_MEM_STREAM(grid, block, bytes, stream)
extern void __syncthreads();
extern int atomicAdd(int* address, int val);
#else
#define CU_DIM(grid, block) <<<grid, block>>>
#define CU_DIM_MEM(grid, block, bytes) <<<grid, block, bytes>>>
#define CU_DIM_MEM_STREAM(grid, block, bytes, stream) <<<grid, block, bytes, stream>>>
#endif
#define CU_ASSERT(err) do { _cuda_assert((err), __FILE__, __LINE__); } while (0);
inline void _cuda_assert(cudaError_t err, const char *file, int line, bool abort=true)
{
if (err != cudaSuccess) {
fprintf(stderr, "CUDA FATAL: %s at %s:%d\n", cudaGetErrorString(err), file, line);
if (abort) {
exit(1);
}
}
}
#define cuIdx(D) (threadIdx.D + blockIdx.D * blockDim.D)
#define cuDim(D) (blockDim.D * gridDim.D)
template <typename T> T *malloc_device(size_t sz) {
T* out = NULL;
if (sz > 0) {
CU_ASSERT(cudaMalloc(&out, sizeof(T) * sz));
}
return out;
}
template <typename T> T *to_device(const T *in, size_t sz, T *out = NULL) {
if (out == NULL) {
CU_ASSERT(cudaMalloc(&out, sizeof(T) * sz));
}
if (sz > 0) {
CU_ASSERT(cudaMemcpy(out, in, sizeof(T) * sz, cudaMemcpyDefault));
}
return out;
}
template <typename T> T *from_device(const T *in, size_t sz, T *out = NULL) {
if (out == NULL) {
CU_ASSERT(cudaHostAlloc(&out, sizeof(T) * sz, 0));
}
if (sz > 0) {
CU_ASSERT(cudaMemcpy(out, in, sizeof(T) * sz, cudaMemcpyDefault));
}
return out;
}
#endif
Raw
inc-stl.h
#include <string>
#include <map>
#include <iomanip>
#include <iostream>
#include <cuda_runtime.h>
#include "deps/stl_reader/stl_reader.h"
#include "deps/cuda_math/cutil_math.h"
#include "deps/clipper/clipper.hpp"
#ifndef STL_H
#define STL_H
using namespace std;
auto round(const float3 &pos, float tol = 1e-6) {
return float3 { floor(pos.x / tol) * tol, floor(pos.y / tol) * tol, floor(pos.z / tol) * tol };
}
bool operator<(const float3 &l, const float3 &r) {
return (l.x < r.x || (l.x == r.x && l.y < r.y) || (l.x == r.x && l.y == r.y && l.z < r.z));
}
struct joint {
int2 e;
float u;
};
struct mesh_t {
vector<float3> vertices;
vector<int3> indices;
vector<float3> normals;
float3 min = { 1e9, 1e9, 1e9 }, max = { -1e9, -1e9, -1e9 };
auto get(joint &item) {
auto a = vertices[item.e.x], b = vertices[item.e.y];
return lerp(a, b, item.u);
}
};
auto load_stl(string file, float tol = 1e-6) {
map<float3, int> mapping;
mesh_t mesh;
auto stl = stl_reader::StlMesh<float, unsigned int>(file);
for (auto i = 0; i < stl.num_solids(); i ++) {
for (auto j = stl.solid_tris_begin(i), n = stl.solid_tris_end(i); j < n; j ++) {
int idx[3];
for (auto k = 0; k < 3; k ++) {
auto ptr = stl.tri_corner_coords(j, k);
auto pos = float3 { ptr[0], ptr[1], ptr[2] };
auto rounded = round(pos, tol);
if (mapping.count(rounded) == 0) {
idx[k] = mesh.vertices.size();
mapping[rounded] = idx[k];
mesh.vertices.push_back(rounded);
} else {
idx[k] = mapping[rounded];
}
mesh.min = fminf(mesh.min, pos);
mesh.max = fmaxf(mesh.max, pos);
}
if (idx[0] != idx[1] && idx[1] != idx[2] && idx[2] != idx[0]) {
auto norm = stl.tri_normal(j);
mesh.normals.push_back(float3 { norm[0], norm[1], norm[2] });
mesh.indices.push_back(int3 { idx[0], idx[1], idx[2] });
}
}
}
return mesh;
}
auto save_stl(string file, vector<float3> &vertices, vector<int3> &faces) {
ofstream out(file);
out << "solid one\n";
for (auto &face : faces) {
auto &a = vertices[face.x], &b = vertices[face.y], &c = vertices[face.z];
out << "facet\n";
out << " outer loop\n";
out << " vertex " << fixed << setprecision(6) << a.x << " " << a.y << " " << a.z << "\n";
out << " vertex " << fixed << setprecision(6) << b.x << " " << b.y << " " << b.z << "\n";
out << " vertex " << fixed << setprecision(6) << c.x << " " << c.y << " " << c.z << "\n";
out << " endloop\n";
out << "endfacet\n";
}
out << "endsolid\n";
out.close();
}
auto save_loop(string file, vector<float3> &vertices, float3 offset) {
vector<float3> pos;
vector<int3> idx;
int max = vertices.size() * 2;
for (auto &vert : vertices) {
auto start = (int) pos.size();
pos.push_back(vert);
pos.push_back(vert + offset);
idx.push_back(int3 { (start + 0) % max, (start + 1) % max, (start + 2) % max });
idx.push_back(int3 { (start + 2) % max, (start + 1) % max, (start + 3) % max });
}
save_stl(file, pos, idx);
}
struct bucket_t {
float val;
int *faces;
int faceNum;
};
struct splited_t {
int f;
joint from, to;
};
auto idx_of(int2 &i, int n) {
return i.x < i.y ? i.x * n + i.y : i.y * n + i.x;
}
const int DIR_X = 0,
DIR_Y = 1,
DIR_Z = 2;
struct ed_t {
int i, f;
};
struct conn_t {
ed_t a, b;
};
auto get_loops(mesh_t &mesh, splited_t *splited, int num) {
auto maxVert = mesh.vertices.size() + 1;
auto conns = map<int, conn_t>();
auto coord = map<int, float3>();
for (int i = 0; i < num; i ++) {
auto &edge = splited[i];
auto f = edge.f;
auto from = idx_of(edge.from.e, maxVert),
to = idx_of(edge.to.e, maxVert);
conns.count(from) ? (conns[from].b = { to, f }) : (conns[from].a = { to, f });
conns.count(to) ? (conns[to].b = { from, f }) : (conns[to].a = { from, f });
coord[from] = mesh.get(edge.from);
coord[to] = mesh.get(edge.to);
}
auto loops = vector<vector<float3>>();
while (conns.size()) {
auto loop = vector<float3>();
auto &begin = *conns.begin();
auto current = begin.first, prev = begin.second.a.i;
while (conns.count(current)) {
loop.push_back(coord[current]);
auto &pair = conns[current];
auto next = pair.a.i == prev ? pair.b.i : pair.a.i;
conns.erase(current);
prev = current;
current = next;
}
cout << "loop " << loops.size() << ": " << loop.size() << " vertices" << endl;
// save_loop(string("lp") + to_string(bucket.val) + "-" + to_string(loops.size()) + ".stl", loop, float3 { 1, 0, 0 });
loops.push_back(loop);
}
return loops;
}
#endif
Raw
inc-svg.h
#include <cmath>
#include <algorithm>
#include <ctime>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <iomanip>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include "deps/clipper/clipper.hpp"
#ifndef SVG_H
#define SVG_H
using namespace std;
using namespace ClipperLib;
class SVGBuilder
{
static string ColorToHtml(unsigned clr)
{
stringstream ss;
ss << '#' << hex << std::setfill('0') << setw(6) << (clr & 0xFFFFFF);
return ss.str();
}
//------------------------------------------------------------------------------
static float GetAlphaAsFrac(unsigned clr)
{
return ((float)(clr >> 24) / 255);
}
//------------------------------------------------------------------------------
class StyleInfo
{
public:
PolyFillType pft;
unsigned brushClr;
unsigned penClr;
double penWidth;
bool showCoords;
StyleInfo()
{
pft = pftNonZero;
brushClr = 0xFFFFFFCC;
penClr = 0xFF000000;
penWidth = 0.8;
showCoords = false;
}
};
class PolyInfo
{
public:
Paths paths;
StyleInfo si;
PolyInfo(Paths paths, StyleInfo style)
{
this->paths = paths;
this->si = style;
}
};
typedef std::vector<PolyInfo> PolyInfoList;
private:
PolyInfoList polyInfos;
static const std::string svg_xml_start[];
static const std::string poly_end[];
public:
StyleInfo style;
void AddPaths(Paths& poly)
{
if (poly.size() == 0) return;
polyInfos.push_back(PolyInfo(poly, style));
}
bool SaveToFile(const string& filename, double scale = 1.0, int margin = 10)
{
//calculate the bounding rect ...
PolyInfoList::size_type i = 0;
Paths::size_type j;
while (i < polyInfos.size())
{
j = 0;
while (j < polyInfos[i].paths.size() &&
polyInfos[i].paths[j].size() == 0) j++;
if (j < polyInfos[i].paths.size()) break;
i++;
}
if (i == polyInfos.size()) return false;
IntRect rec;
rec.left = polyInfos[i].paths[j][0].X;
rec.right = rec.left;
rec.top = polyInfos[i].paths[j][0].Y;
rec.bottom = rec.top;
for ( ; i < polyInfos.size(); ++i)
for (Paths::size_type j = 0; j < polyInfos[i].paths.size(); ++j)
for (Path::size_type k = 0; k < polyInfos[i].paths[j].size(); ++k)
{
IntPoint ip = polyInfos[i].paths[j][k];
if (ip.X < rec.left) rec.left = ip.X;
else if (ip.X > rec.right) rec.right = ip.X;
if (ip.Y < rec.top) rec.top = ip.Y;
else if (ip.Y > rec.bottom) rec.bottom = ip.Y;
}
if (scale == 0) scale = 1.0;
if (margin < 0) margin = 0;
rec.left = (cInt)((double)rec.left * scale);
rec.top = (cInt)((double)rec.top * scale);
rec.right = (cInt)((double)rec.right * scale);
rec.bottom = (cInt)((double)rec.bottom * scale);
cInt offsetX = -rec.left + margin;
cInt offsetY = -rec.top + margin;
ofstream file;
file.open(filename);
if (!file.is_open()) return false;
file.setf(ios::fixed);
file.precision(0);
file << svg_xml_start[0] <<
((rec.right - rec.left) + margin*2) << "px" << svg_xml_start[1] <<
((rec.bottom - rec.top) + margin*2) << "px" << svg_xml_start[2] <<
((rec.right - rec.left) + margin*2) << " " <<
((rec.bottom - rec.top) + margin*2) << svg_xml_start[3];
setlocale(LC_NUMERIC, "C");
file.precision(2);
for (PolyInfoList::size_type i = 0; i < polyInfos.size(); ++i)
{
file << " <path d=\"";
for (Paths::size_type j = 0; j < polyInfos[i].paths.size(); ++j)
{
if (polyInfos[i].paths[j].size() < 3) continue;
file << " M " << ((double)polyInfos[i].paths[j][0].X * scale + offsetX) <<
" " << ((double)polyInfos[i].paths[j][0].Y * scale + offsetY);
for (Path::size_type k = 1; k < polyInfos[i].paths[j].size(); ++k)
{
IntPoint ip = polyInfos[i].paths[j][k];
double x = (double)ip.X * scale;
double y = (double)ip.Y * scale;
file << " L " << (x + offsetX) << " " << (y + offsetY);
}
file << " z";
}
file << poly_end[0] << ColorToHtml(polyInfos[i].si.brushClr) <<
poly_end[1] << GetAlphaAsFrac(polyInfos[i].si.brushClr) <<
poly_end[2] <<
(polyInfos[i].si.pft == pftEvenOdd ? "evenodd" : "nonzero") <<
poly_end[3] << ColorToHtml(polyInfos[i].si.penClr) <<
poly_end[4] << GetAlphaAsFrac(polyInfos[i].si.penClr) <<
poly_end[5] << polyInfos[i].si.penWidth << poly_end[6];
if (polyInfos[i].si.showCoords)
{
file << "<g font-family=\"Verdana\" font-size=\"11\" fill=\"black\">\n\n";
for (Paths::size_type j = 0; j < polyInfos[i].paths.size(); ++j)
{
if (polyInfos[i].paths[j].size() < 3) continue;
for (Path::size_type k = 0; k < polyInfos[i].paths[j].size(); ++k)
{
IntPoint ip = polyInfos[i].paths[j][k];
file << "<text x=\"" << (int)(ip.X * scale + offsetX) <<
"\" y=\"" << (int)(ip.Y * scale + offsetY) << "\">" <<
ip.X << "," << ip.Y << "</text>\n";
file << "\n";
}
}
file << "</g>\n";
}
}
file << "</svg>\n";
file.close();
setlocale(LC_NUMERIC, "");
return true;
}
}; //SVGBuilder
//------------------------------------------------------------------------------
const std::string SVGBuilder::svg_xml_start [] =
{"<?xml version=\"1.0\" standalone=\"no\"?>\n"
"<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.0//EN\"\n"
"\"http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd\">\n\n"
"<svg width=\"",
"\" height=\"",
"\" viewBox=\"0 0 ",
"\" version=\"1.0\" xmlns=\"http://www.w3.org/2000/svg\">\n\n"
};
const std::string SVGBuilder::poly_end [] =
{"\"\n style=\"fill:",
"; fill-opacity:",
"; fill-rule:",
"; stroke:",
"; stroke-opacity:",
"; stroke-width:",
";\"/>\n\n"
};
#endif
Raw
main.cu
#include <map>
#include <vector>
#include <iostream>
#include <string>
#include "deps/clipper/clipper.hpp"
#include "inc/cuda.h"
#include "inc/svg.h"
#include "inc/stl.h"
using namespace std;
struct bound_t {
float3 min, max;
};
__global__ void kernel_floor(float3 *vertices, int vertNum, float3 *out, float tol, int dir, float shift = 0) {
for (int i = cuIdx(x); i < vertNum; i += cuDim(x)) {
if (dir == DIR_X) {
out[i].x = (floor(vertices[i].x / tol) + shift) * tol;
} else if (dir == DIR_Y) {
out[i].y = (floor(vertices[i].y / tol) + shift) * tol;
} else {
out[i].z = (floor(vertices[i].z / tol) + shift) * tol;
}
}
}
__global__ void kernel_prepare(
float3 *vertices, int3 *faces, int faceNum, bound_t *bounds, float3 *norm) {
for (int i = cuIdx(x); i < faceNum; i += cuDim(x)) {
auto &f = faces[i];
auto &a = vertices[f.x], &b = vertices[f.y], &c = vertices[f.z];
bounds[i].min = fminf(fminf(a, b), c);
bounds[i].max = fmaxf(fmaxf(a, b), c);
norm[i] = normalize(cross(a - b, b - c));
}
}
__global__ void kernel_group(
float3 *vertices, int3 *faces, int faceNum, bound_t *bounds, float val, int dir,
int *out, int *idx) {
for (int i = cuIdx(x); i < faceNum; i += cuDim(x)) {
auto &b = bounds[i];
auto min = dir == DIR_X ? b.min.x : dir == DIR_Y ? b.min.y : b.min.z,
max = dir == DIR_X ? b.max.x : dir == DIR_Y ? b.max.y : b.max.z;
if (min < val && val < max) {
out[atomicAdd(idx, 1)] = i;
}
}
}
__global__ void kernel_split(
float3 *vertices, int3 *faces, int* bucket, int binLen, float val, int dir,
splited_t *out, int *idx) {
for (int i = cuIdx(x); i < binLen; i += cuDim(x)) {
auto &f = faces[bucket[i]];
auto &a = vertices[f.x], &b = vertices[f.y], &c = vertices[f.z];
float m, n, l;
if (dir == DIR_X) {
m = a.x; n = b.x; l = c.x;
} else if (dir == DIR_Y) {
m = a.y; n = b.y; l = c.y;
} else if (dir == DIR_Z) {
m = a.z; n = b.z; l = c.z;
}
float u = (m - val) / (m - n),
v = (n - val) / (n - l),
w = (l - val) / (l - m);
out[atomicAdd(idx, 1)] =
u > 0 && u < 1 && v > 0 && v < 1 ? splited_t { bucket[i], f.x, f.y, u, f.y, f.z, v } :
v > 0 && v < 1 && w > 0 && w < 1 ? splited_t { bucket[i], f.y, f.z, v, f.z, f.x, w } :
splited_t { bucket[i], f.z, f.x, w, f.x, f.y, u };
}
}
template <typename T> auto range(T from, T to, T step = 1) {
vector<T> ret;
for (T i = from; i < to; i += step) {
ret.push_back(i);
}
return ret;
}
int main() {
auto mesh = load_stl("model/tube.stl");
auto xs = range(0.0f, 10.0f, 1.0f),
ys = range(0.0f, 10.0f, 1.0f),
zs = range(0.0f, 10.0f, 1.0f);
auto oriVerts = to_device(mesh.vertices.data(), mesh.vertices.size());
auto vertNum = (int) mesh.vertices.size();
auto faces = to_device(mesh.indices.data(), mesh.indices.size());
auto faceNum = (int) mesh.indices.size();
printf("bound: %f, %f, %f / %f, %f, %f\n",
mesh.min.x, mesh.min.y, mesh.min.z, mesh.max.x, mesh.max.y, mesh.max.z);
printf("faces: %d\n", faceNum);
float tol = 1e-3;
for (int i = 0; i < xs.size(); i ++) {
// make sure that no vertice is at x
xs[i] = (floor(xs[i] / tol) + 0.5) * tol;
}
auto vertices = malloc_device<float3>(vertNum);
auto output = map<float, ClipperLib::Paths>();
for (int shift = 0; shift < 2; shift ++) {
kernel_floor CU_DIM(512, 128) (oriVerts, vertNum, vertices, tol, DIR_X, shift * 1.);
//kernel_floor CU_DIM(512, 128) (oriVerts, vertNum, vertices, tol, DIR_Y, shift * 1.);
//kernel_floor CU_DIM(512, 128) (oriVerts, vertNum, vertices, tol, DIR_Z, shift * 1.);
auto bounds = malloc_device<bound_t>(faceNum);
auto norm = malloc_device<float3>(faceNum);
kernel_prepare CU_DIM(512, 128) (vertices, faces, faceNum, bounds, norm);
CU_ASSERT(cudaGetLastError());
//save_stl("dragon.saved.stl", mesh.vertices, mesh.indices);
for (auto x : xs) {
int binLen = 0;
auto len = to_device(&binLen, 1);
auto bin = malloc_device<int>(faceNum);
kernel_group CU_DIM(512, 128) (vertices, faces, faceNum, bounds, x, DIR_X, bin, len);
CU_ASSERT(cudaGetLastError());
from_device(len, 1, &binLen);
printf("x: %f, bin %d\n", x, binLen);
if (!binLen) {
continue;
}
int outLen = 0;
to_device(&outLen, 1, len);
auto splited = malloc_device<splited_t>(binLen);
kernel_split CU_DIM(512, 128) (vertices, faces, bin, binLen, x, DIR_X, splited, len);
CU_ASSERT(cudaGetLastError());
auto ret = from_device(splited, binLen);
auto loops = get_loops(mesh, ret, binLen);
auto clipper = ClipperLib::Clipper();
for (auto loop : loops) {
auto path = ClipperLib::Path();
for (auto pt : loop) {
path.push_back(IntPoint { (long long) (pt.y * 50), (long long) (pt.z * 50) });
}
clipper.AddPath(path, ptSubject, true);
}
if (output.count(x)) {
clipper.AddPaths(output[x], ptSubject, true);
}
auto paths = ClipperLib::Paths();
clipper.Execute(ctUnion, paths);
output[x] = paths;
}
}
for (auto &pair : output) {
auto svg = SVGBuilder();
svg.AddPaths(pair.second);
svg.SaveToFile("svg/poly-x-" + to_string(pair.first) + ".svg");
}
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment