wch
/
nh-rep
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
21 Commits
1 Branch
0 Tags
14 MiB
Branch: main
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'main'
${ noResults }
nh-rep/code/IsoSurfacing/App/console_pytorch/GridComp.h

1160 lines
35 KiB
Raw Permalink Normal View History

init
5 months ago
#pragma once
#include "MyImplicitFunc.h"
#include "ThirdParty/IsoSurfGen/ProgressBar.h"
#include "ThirdParty/termcolor/termcolor.hpp"
#include <unordered_map>
#include <unordered_set>
#include <fstream>
#include <cassert>
#include "ThirdParty/DualContour/octree.h"
int edgeTables[256] =
{
0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
0x190, 0x99, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f, 0x435, 0x53c,
0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
0x3a0, 0x2a9, 0x1a3, 0xaa, 0x7a6, 0x6af, 0x5a5, 0x4ac,
0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
0x460, 0x569, 0x663, 0x76a, 0x66, 0x16f, 0x265, 0x36c,
0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5, 0x2fc,
0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55, 0x15c,
0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc,
0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
0xcc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
0x15c, 0x55, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
0x2fc, 0x3f5, 0xff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663, 0x569, 0x460,
0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0,
0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33, 0x339, 0x230,
0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99, 0x190,
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0};
const int cubestatusmap[256] = {
0, 1, 16, 17, 64, 65, 80, 81, 4, 5, 20, 21, 68, 69, 84, 85, 2, 3,
18, 19, 66, 67, 82, 83, 6, 7, 22, 23, 70, 71, 86, 87, 32, 33, 48,
49, 96, 97, 112, 113, 36, 37, 52, 53, 100, 101, 116, 117, 34, 35,
50, 51, 98, 99, 114, 115, 38, 39, 54, 55, 102, 103, 118, 119, 128,
129, 144, 145, 192, 193, 208, 209, 132, 133, 148, 149, 196, 197,
212, 213, 130, 131, 146, 147, 194, 195, 210, 211, 134, 135, 150,
151, 198, 199, 214, 215, 160, 161, 176, 177, 224, 225, 240, 241,
164, 165, 180, 181, 228, 229, 244, 245, 162, 163, 178, 179, 226,
227, 242, 243, 166, 167, 182, 183, 230, 231, 246, 247, 8, 9, 24,
25, 72, 73, 88, 89, 12, 13, 28, 29, 76, 77, 92, 93, 10, 11, 26,
27, 74, 75, 90, 91, 14, 15, 30, 31, 78, 79, 94, 95, 40, 41, 56,
57, 104, 105, 120, 121, 44, 45, 60, 61, 108, 109, 124, 125, 42,
43, 58, 59, 106, 107, 122, 123, 46, 47, 62, 63, 110, 111, 126,
127, 136, 137, 152, 153, 200, 201, 216, 217, 140, 141, 156, 157,
204, 205, 220, 221, 138, 139, 154, 155, 202, 203, 218, 219, 142,
143, 158, 159, 206, 207, 222, 223, 168, 169, 184, 185, 232, 233,
248, 249, 172, 173, 188, 189, 236, 237, 252, 253, 170, 171, 186,
187, 234, 235, 250, 251, 174, 175, 190, 191, 238, 239, 254, 255};
class IntersectionInfo
{
public:
IntersectionInfo(int _start, int _end)
{
vert[0] = _start;
vert[1] = _end;
}
inline bool operator<(const IntersectionInfo &p) const
{
for (int i = 0; i < 2; i++)
{
if (vert[i] < p.vert[i])
{
return true;
}
else if (vert[i] > p.vert[i])
{
return false;
}
}
return false;
}
inline bool operator==(const IntersectionInfo &p) const
{
for (int i = 0; i < 2; i++)
{
if (vert[i] != p.vert[i])
{
return false;
}
}
return true;
}
inline bool operator!=(const IntersectionInfo &p) const
{
for (int i = 0; i < 2; i++)
{
if (vert[i] != p.vert[i])
{
return true;
}
}
return false;
}
public:
int vert[2];
float shift, nx, ny, nz;
};
//////////////////////////////////////
template <typename Real>
class MyOctreeLayer
{
private:
MyImplicitFunc<Real> *myfunc;
int depth, gridsize, offsets_[8];
Real start_x, end_x, start_y, end_y, start_z, end_z;
Real fCellLengthX, fCellLengthY, fCellLengthZ;
std::vector<int> valid_cell_id;
std::vector<int> valid_cell_types;
std::unordered_map<int, Real> point_value_map;
Real isovalue;
bool verbose;
public:
MyOctreeLayer<Real>(
int _depth,
MyImplicitFunc<Real> *_myfunc,
Real _start_x, Real _end_x,
Real _start_y, Real _end_y,
Real _start_z, Real _end_z,
Real _isovalue = (Real)0,
bool _verbose = true)
: depth(_depth), myfunc(_myfunc),
start_x(_start_x), end_x(_end_x),
start_y(_start_y), end_y(_end_y),
start_z(_start_z), end_z(_end_z), isovalue(_isovalue), verbose(_verbose)
{
assert(depth > 0);
gridsize = pow(2, depth) + 1;
fCellLengthX = (end_x - start_x) / (gridsize - 1);
fCellLengthY = (end_y - start_y) / (gridsize - 1);
fCellLengthZ = (end_z - start_z) / (gridsize - 1);
offsets_[0] = 0;
offsets_[1] = 1;
offsets_[2] = 1 + gridsize;
offsets_[3] = gridsize;
offsets_[4] = gridsize * gridsize;
offsets_[5] = 1 + gridsize * gridsize;
offsets_[6] = 1 + gridsize + gridsize * gridsize;
offsets_[7] = gridsize + gridsize * gridsize;
}
/////////////////////////////////////////////////
const std::vector<int> &get_valid_cells() const { return valid_cell_id; }
/////////////////////////////////////////////////
const std::vector<int> &get_valid_cell_types() const { return valid_cell_types; }
/////////////////////////////////////////////////
const int get_grid_size() const { return gridsize; }
/////////////////////////////////////////////////
const int get_depth() const { return depth; }
/////////////////////////////////////////////////
const int get_point_id(int cell_id, int vertex_id) const
{
unsigned int X(gridsize - 1), Y(gridsize - 1);
unsigned int x = cell_id % X;
cell_id /= X;
unsigned int y = cell_id % Y;
cell_id /= Y;
unsigned int z = cell_id;
unsigned int p = x + y * gridsize + z * gridsize * gridsize;
return p + offsets_[vertex_id];
}
/////////////////////////////////////////////////
const void get_point(int globa_vertex_id, TinyVector<Real, 3> &p)
{
int k = globa_vertex_id / (gridsize * gridsize);
int j = (globa_vertex_id - k * gridsize * gridsize) / gridsize;
int i = globa_vertex_id % gridsize;
p[2] = start_z + k * fCellLengthZ;
p[1] = start_y + j * fCellLengthY;
p[0] = start_x + i * fCellLengthX;
}
/////////////////////////////////////////////////
const std::vector<int> &get_valid_cell_id() const { return valid_cell_id; }
/////////////////////////////////////////////////
const int get_gridsize() const { return gridsize; }
/////////////////////////////////////////////////
const std::unordered_map<int, Real> &get_point_value_map() const { return point_value_map; }
/////////////////////////////////////////////////
void scan_octants(MyOctreeLayer<Real> *parent = 0, bool gather_neighbors = true)
{
valid_cell_id.reserve(10 * gridsize * gridsize);
valid_cell_id.resize(0);
valid_cell_types.reserve(10 * gridsize * gridsize);
valid_cell_types.resize(0);
MyUtil<Real> *util = myfunc->util;
std::vector<Real> &values = myfunc->util->values;
std::vector<int> &point_id_total_cache = myfunc->util->point_id_total_cache;
std::vector<int> &point_id_cache = myfunc->util->point_id_cache;
values.resize(0);
point_value_map.clear();
point_id_total_cache.resize(0);
point_id_cache.resize(0);
////////////////////
std::cout << termcolor::yellow << "---------"
<< "process " << termcolor::blue << termcolor::on_white << depth << termcolor::reset << termcolor::yellow << "-th octree layer!" << std::endl;
std::cout << termcolor::cyan;
print_progress_bar(0.0f, verbose);
if (parent == 0)
{
compute_grid_point_values(gridsize, values);
std::vector<int> total_cell_types((gridsize - 1) * (gridsize - 1) * (gridsize - 1));
#pragma omp parallel for
for (int c = 0; c < (gridsize - 1) * (gridsize - 1) * (gridsize - 1); c++)
{
unsigned char cubetype(0);
for (int s = 0; s < 8; ++s)
if (values[get_point_id(c, s)] > isovalue)
cubetype |= (1 << s);
total_cell_types[c] = cubetype;
}
for (int n = 0; n < (gridsize - 1) * (gridsize - 1) * (gridsize - 1); n++)
{
int k = n / ((gridsize - 1) * (gridsize - 1));
int j = (n - k * (gridsize - 1) * (gridsize - 1)) / (gridsize - 1);
int i = n % (gridsize - 1);
int cubetype = total_cell_types[n];
bool find = true;
if ((cubetype == 0 || cubetype == 255))
{
find = false;
if (gather_neighbors)
{
for (int kshift = -1; kshift <= 1; kshift++)
{
int kk = k + kshift;
if (kk < 0 || kk >= gridsize - 1)
continue;
int cid_k = kk * (gridsize - 1) * (gridsize - 1);
for (int jshift = -1; jshift <= 1; jshift++)
{
int jj = j + jshift;
if (jj < 0 || jj >= gridsize - 1)
continue;
int cid_j = cid_k + jj * (gridsize - 1);
for (int ishift = -1; ishift <= 1; ishift++)
{
int ii = i + ishift;
if (ii < 0 || ii >= gridsize - 1)
continue;
int c_id = cid_j + ii;
if (total_cell_types[c_id] != cubetype)
{
find = true;
break;
}
}
if (find)
break;
}
if (find)
break;
}
}
}
if (find)
{
valid_cell_types.push_back(cubetype);
valid_cell_id.push_back(n);
}
}
for (size_t i = 0; i < values.size(); i++)
point_value_map[i] = values[i];
}
else
{
assert(parent->get_depth() == depth - 1);
for (size_t c = 0; c < parent->get_valid_cell_id().size(); c++)
{
int n = parent->get_valid_cell_id()[c];
int k = n / ((parent->get_gridsize() - 1) * (parent->get_gridsize() - 1));
int j = (n - k * (parent->get_gridsize() - 1) * (parent->get_gridsize() - 1)) / (parent->get_gridsize() - 1);
int i = n % (parent->get_gridsize() - 1);
k *= 2, j *= 2, i *= 2;
int mid[8];
mid[0] = k * (gridsize - 1) * (gridsize - 1) + j * (gridsize - 1) + i;
mid[1] = k * (gridsize - 1) * (gridsize - 1) + j * (gridsize - 1) + i + 1;
mid[2] = k * (gridsize - 1) * (gridsize - 1) + (j + 1) * (gridsize - 1) + i;
mid[3] = k * (gridsize - 1) * (gridsize - 1) + (j + 1) * (gridsize - 1) + i + 1;
mid[4] = (k + 1) * (gridsize - 1) * (gridsize - 1) + j * (gridsize - 1) + i;
mid[5] = (k + 1) * (gridsize - 1) * (gridsize - 1) + j * (gridsize - 1) + i + 1;
mid[6] = (k + 1) * (gridsize - 1) * (gridsize - 1) + (j + 1) * (gridsize - 1) + i;
mid[7] = (k + 1) * (gridsize - 1) * (gridsize - 1) + (j + 1) * (gridsize - 1) + i + 1;
for (int s = 0; s < 8; s++)
{
for (int v = 0; v < 8; v++)
point_value_map[get_point_id(mid[s], v)] = 1;
}
}
typename std::unordered_map<int, Real>::iterator iter = point_value_map.begin();
point_id_total_cache.resize(point_value_map.size());
for (int i = 0; iter != point_value_map.end(); iter++, i++)
point_id_total_cache[i] = iter->first;
at::Tensor m;
if (util->dim_type == 3)
m = torch::zeros({1, std::min((int)point_id_total_cache.size(), util->buffer), 3});
else
m = torch::zeros({std::min((int)point_id_total_cache.size(), util->buffer), 3});
auto mpstorage = static_cast<float *>(m.storage().data());
at::Tensor output;
int counter = 0;
for (size_t pi = 0; pi < point_id_total_cache.size(); pi++)
{
int n = point_id_total_cache[pi];
int k = n / (gridsize * gridsize);
int j = (n - k * gridsize * gridsize) / gridsize;
int i = n % gridsize;
mpstorage[counter + 2] = start_z + k * fCellLengthZ;
mpstorage[counter + 1] = start_y + j * fCellLengthY;
mpstorage[counter] = start_x + i * fCellLengthX;
point_id_cache.push_back(n);
counter += 3;
if (counter == 3 * util->buffer || pi + 1 == point_id_total_cache.size())
{
std::vector<torch::jit::IValue> inputs;
if (util->use_GPU)
{
inputs.push_back(m.to(at::kCUDA));
}
else
inputs.push_back(m);
output = myfunc->net->forward(inputs).toTensor().cpu();
auto outputstorage = static_cast<float *>(output.storage().data());
for (size_t s = 0; s < point_id_cache.size(); s++)
{
point_value_map[point_id_cache[s]] = outputstorage[s];
}
print_progress_bar((pi + 1.0f) / point_id_total_cache.size(), verbose);
if (pi + 1 == point_id_total_cache.size())
printf("\n");
point_id_cache.resize(0);
counter = 0;
}
}
valid_cell_types.reserve(8 * parent->get_valid_cell_id().size());
valid_cell_types.resize(0);
// create its own validcells
std::vector<int> total_cell_types(8 * parent->get_valid_cell_id().size());
std::vector<int> total_cells(8 * parent->get_valid_cell_id().size());
#pragma omp parallel for
for (int c = 0; c < (int)parent->get_valid_cell_id().size(); c++)
{
int n = parent->get_valid_cell_id()[c];
int k = n / ((parent->get_gridsize() - 1) * (parent->get_gridsize() - 1));
int j = (n - k * (parent->get_gridsize() - 1) * (parent->get_gridsize() - 1)) / (parent->get_gridsize() - 1);
int i = n % (parent->get_gridsize() - 1);
k *= 2, j *= 2, i *= 2;
int *mid = &total_cells[8 * c];
mid[0] = k * (gridsize - 1) * (gridsize - 1) + j * (gridsize - 1) + i;
mid[1] = k * (gridsize - 1) * (gridsize - 1) + j * (gridsize - 1) + i + 1;
mid[2] = k * (gridsize - 1) * (gridsize - 1) + (j + 1) * (gridsize - 1) + i;
mid[3] = k * (gridsize - 1) * (gridsize - 1) + (j + 1) * (gridsize - 1) + i + 1;
mid[4] = (k + 1) * (gridsize - 1) * (gridsize - 1) + j * (gridsize - 1) + i;
mid[5] = (k + 1) * (gridsize - 1) * (gridsize - 1) + j * (gridsize - 1) + i + 1;
mid[6] = (k + 1) * (gridsize - 1) * (gridsize - 1) + (j + 1) * (gridsize - 1) + i;
mid[7] = (k + 1) * (gridsize - 1) * (gridsize - 1) + (j + 1) * (gridsize - 1) + i + 1;
for (int ci = 0; ci < 8; ci++)
{
unsigned char cubetype(0);
for (int s = 0; s < 8; ++s)
if (point_value_map.find(get_point_id(mid[ci], s))->second > isovalue)
cubetype |= (1 << s);
total_cell_types[8 * c + ci] = cubetype;
}
}
std::unordered_map<int, int> cell_map;
for (size_t c = 0; c < total_cells.size(); c++)
{
cell_map[total_cells[c]] = c;
}
for (size_t c = 0; c < total_cells.size(); c++)
{
int n = total_cells[c];
int k = n / ((gridsize - 1) * (gridsize - 1));
int j = (n - k * (gridsize - 1) * (gridsize - 1)) / (gridsize - 1);
int i = n % (gridsize - 1);
int cubetype = total_cell_types[c];
bool find = true;
if ((cubetype == 0 || cubetype == 255))
{
find = false;
if (gather_neighbors)
{
for (int kshift = -1; kshift <= 1; kshift++)
{
int kk = k + kshift;
if (kk < 0 || kk >= gridsize - 1)
continue;
int cid_k = kk * (gridsize - 1) * (gridsize - 1);
for (int jshift = -1; jshift <= 1; jshift++)
{
int jj = j + jshift;
if (jj < 0 || jj >= gridsize - 1)
continue;
int cid_j = cid_k + jj * (gridsize - 1);
for (int ishift = -1; ishift <= 1; ishift++)
{
int ii = i + ishift;
if (ii < 0 || ii >= gridsize - 1)
continue;
int c_id = cid_j + ii;
typename std::unordered_map<int, int>::iterator miter = cell_map.find(c_id);
if (miter != cell_map.end())
{
if (total_cell_types[miter->second] != cubetype)
{
find = true;
break;
}
}
}
if (find)
break;
}
if (find)
break;
}
}
}
if (find)
{
valid_cell_types.push_back(cubetype);
valid_cell_id.push_back(n);
}
}
}
// std::cout << termcolor::reset << "valid_cells: " << valid_cell_id.size() << std::endl;
}
/////////////////////////////////////////////////
void compute_grid_point_values(int specified_size, std::vector<Real> &values)
{
MyUtil<Real> *util = myfunc->util;
values.resize(specified_size * specified_size * specified_size);
Real _fCellLengthX = (end_x - start_x) / (specified_size - 1);
Real _fCellLengthY = (end_y - start_y) / (specified_size - 1);
Real _fCellLengthZ = (end_z - start_z) / (specified_size - 1);
at::Tensor m;
if (util->dim_type == 3)
m = torch::zeros({1, std::min((int)values.size(), util->buffer), 3});
else // USE_2_DIM_NETWORK
m = torch::zeros({std::min((int)values.size(), util->buffer), 3});
at::Tensor output;
auto mpstorage = static_cast<float *>(m.storage().data());
int counter = 0;
int start_n = 0;
for (int n = 0; n < specified_size * specified_size * specified_size; n++)
{
int k = n / (specified_size * specified_size);
int j = (n - k * specified_size * specified_size) / specified_size;
int i = n % specified_size;
mpstorage[3 * counter + 2] = start_z + k * _fCellLengthZ;
mpstorage[3 * counter + 1] = start_y + j * _fCellLengthY;
mpstorage[3 * counter] = start_x + i * _fCellLengthX;
counter++;
if (counter == util->buffer || n + 1 == specified_size * specified_size * specified_size)
{
std::vector<torch::jit::IValue> inputs;
if (util->use_GPU)
inputs.push_back(m.to(at::kCUDA));
else
inputs.push_back(m);
output = myfunc->net->forward(inputs).toTensor().cpu();
auto outputstorage = static_cast<float *>(output.storage().data());
#pragma omp parallel for
for (int l = start_n; l <= n; l++)
values[l] = outputstorage[l - start_n];
start_n = n + 1;
counter = 0;
print_progress_bar((n + 1.0f) / (specified_size * specified_size * specified_size), verbose);
if (n + 1 == specified_size * specified_size * specified_size)
printf("\n");
}
}
}
/////////////////////////////////////////////////
void save_cubes(const char filename[]) // *.obj format
{
std::ofstream mfile(filename);
int cubfacets[][4] = {{3, 2, 1, 0}, {6, 5, 1, 2}, {4, 0, 1, 5}, {4, 7, 3, 0}, {4, 5, 6, 7}, {2, 3, 7, 6}};
int counter = 1;
for (size_t c = 0; c < valid_cell_id.size(); c++)
{
int n = valid_cell_id[c];
int k = n / ((gridsize - 1) * (gridsize - 1));
int j = (n - k * (gridsize - 1) * (gridsize - 1)) / (gridsize - 1);
int i = n % (gridsize - 1);
Real x = start_x + i * fCellLengthX;
Real y = start_y + j * fCellLengthY;
Real z = start_z + k * fCellLengthZ;
mfile << "v " << x << ' ' << y << ' ' << z << std::endl;
mfile << "v " << x + fCellLengthX << ' ' << y << ' ' << z << std::endl;
mfile << "v " << x + fCellLengthX << ' ' << y + fCellLengthY << ' ' << z << std::endl;
mfile << "v " << x << ' ' << y + fCellLengthY << ' ' << z << std::endl;
mfile << "v " << x << ' ' << y << ' ' << z + fCellLengthZ << std::endl;
mfile << "v " << x + fCellLengthX << ' ' << y << ' ' << z + fCellLengthZ << std::endl;
mfile << "v " << x + fCellLengthX << ' ' << y + fCellLengthY << ' ' << z + fCellLengthZ << std::endl;
mfile << "v " << x << ' ' << y + fCellLengthY << ' ' << z + fCellLengthZ << std::endl;
for (int i = 0; i < 6; i++)
mfile << "f "
<< counter + cubfacets[i][0] << ' '
<< counter + cubfacets[i][1] << ' '
<< counter + cubfacets[i][2] << ' '
<< counter + cubfacets[i][3] << std::endl;
counter += 8;
}
mfile.close();
}
/////////////////////////////////////////////////
void compute_intersection(std::set<IntersectionInfo> &edgeintersections)
{
edgeintersections.clear();
////////////////////////////////////////////////////////////////////
std::vector<int> cache_cube(valid_cell_id.size() * 2);
for (size_t i = 0; i < valid_cell_id.size(); i++)
{
cache_cube[2 * i] = valid_cell_id[i];
cache_cube[2 * i + 1] = valid_cell_types[i];
}
MyUtil<Real> *util = myfunc->util;
/////////////////////////////////////////////////////////////////////
int buffer = myfunc->util->buffer;
std::vector<TinyVector<Real, 3>> edge_start, edge_end, intersection, normals;
std::vector<Real> start_fvalue, end_fvalue;
edge_start.reserve(buffer + 12);
edge_end.reserve(buffer + 12);
intersection.reserve(buffer + 12);
normals.reserve(buffer + 12);
start_fvalue.reserve(buffer + 12);
end_fvalue.reserve(buffer + 12);
print_progress_bar(0.0f, verbose);
int counter = 0;
TinyVector<Real, 3> ipoint[8];
int label[12] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048};
int edgepair[12][2] = {{0, 1}, {1, 2}, {3, 2}, {0, 3}, {4, 5}, {5, 6}, {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
int corner[8];
Real values[8];
std::set<std::pair<int, int>> edgeset;
std::vector<int> edgestore;
edgestore.reserve(2 * (buffer + 12));
Real unit_length = (end_z - start_z) / (gridsize - 1);
for (int s = 0; s < cache_cube.size(); s += 2)
{
int c = cache_cube[s];
for (int i = 0; i < 8; ++i)
{
corner[i] = get_point_id(c, i);
values[i] = point_value_map.find(corner[i])->second;
}
unsigned char cubetype = cache_cube[s + 1];
for (int i = 0; i < 8; i++)
{
get_point(corner[i], ipoint[i]);
}
for (int i = 0; i < 12; i++)
{
if (edgeTables[cubetype] & label[i])
{
std::pair<int, int> mpair(corner[edgepair[i][0]], corner[edgepair[i][1]]);
if (edgeset.find(mpair) != edgeset.end())
{
;
}
else
{
edgeset.insert(mpair);
edgestore.push_back(corner[edgepair[i][0]]);
edgestore.push_back(corner[edgepair[i][1]]);
edge_start.push_back(ipoint[edgepair[i][0]]);
edge_end.push_back(ipoint[edgepair[i][1]]);
start_fvalue.push_back(values[edgepair[i][0]]);
end_fvalue.push_back(values[edgepair[i][1]]);
counter++;
}
}
}
if (counter != 0 && (counter >= buffer || s + 2 == cache_cube.size()))
{
edge_start.resize(counter);
edge_end.resize(counter);
start_fvalue.resize(counter);
end_fvalue.resize(counter);
intersection.resize(counter);
normals.resize(counter);
myfunc->directed_distance(edge_start, edge_end, start_fvalue, end_fvalue, intersection, normals, isovalue);
for (size_t h = 0; h < intersection.size(); h++)
{
IntersectionInfo info(edgestore[2 * h], edgestore[2 * h + 1]);
info.nx = (float)normals[h][0], info.ny = (float)normals[h][1], info.nz = (float)normals[h][2];
info.shift = (float)((intersection[h] - edge_start[h]).Length() / unit_length);
edgeintersections.insert(info);
}
counter = 0;
edgestore.resize(0);
edgeset.clear();
edge_start.resize(0);
edge_end.resize(0);
start_fvalue.resize(0);
end_fvalue.resize(0);
intersection.resize(0);
normals.resize(0);
print_progress_bar((s + 2.0f) / cache_cube.size(), verbose);
if (s + 2 == cache_cube.size())
printf("\n");
}
}
}
/////////////////////////////////////////////////
void save_dcf(const char filename[], bool binary) // modified dcf format
{
if (valid_cell_id.empty())
{
std::cerr << "scan_octants() should be called first!" << std::endl;
return;
}
//--- compute intersection---
std::set<IntersectionInfo> edgeintersections;
compute_intersection(edgeintersections);
//---------------------------
class Node
{
public:
Node(int _id, int _depth, bool _is_leaf = false, int _cell_type = 0)
: id_(_id), depth_(_depth), is_leaf_(_is_leaf)
{
parent = 0;
cell_type_ = cubestatusmap[_cell_type];
for (int i = 0; i < 8; i++)
child[i] = 0;
for (int i = 0; i < 12; i++)
interinfo[i] = 0;
}
~Node()
{
for (int i = 0; i < 12; i++)
if (interinfo[i])
delete interinfo[i];
}
void travel_print(std::ofstream &out, bool binary)
{
bool is_empty = true;
for (int i = 0; i < 8; i++)
{
if (child[i])
{
is_empty = false;
break;
}
}
int type = is_leaf_ ? 2 : (is_empty ? 1 : 0);
if (binary)
out.write((char *)&type, sizeof(int));
else
{
for (int k = 0; k < depth_; k++)
out << "*";
out << "type: " << type << std::endl;
}
if (type == 0)
{
for (int i = 0; i < 8; i++)
{
if (child[i])
child[i]->travel_print(out, binary);
else
{
if (binary)
{
int ctype = 1;
out.write((char *)&ctype, sizeof(int));
}
else
{
out << "type: 1" << std::endl;
}
}
}
}
else if (type == 1)
{
// short sg = 1; // it is not a true sign value, faked for DCF.
// out.write((char*)&sg, sizeof(short));
}
else
{
if (binary)
out.write((char *)&cell_type_, sizeof(int)); // different with orginal dcf format
else
{
for (int k = 0; k < depth_; k++)
out << "*";
out << "c" << cell_type_ << std::endl;
}
const int zero = 0, one = 1;
for (int l = 0; l < 12; l++)
{
if (binary)
{
if (interinfo[l])
{
out.write((char *)&one, sizeof(int));
if (fabs((*interinfo[l])[0] - 0.524244) < 0.00001)
{
int k = 0;
k = 1;
}
float h = (*interinfo[l])[0], x = (*interinfo[l])[1], y = (*interinfo[l])[2], z = (*interinfo[l])[3];
// out.write((char*)&(*interinfo[l]), sizeof(float)*4);
out.write((char *)&h, sizeof(float));
out.write((char *)&x, sizeof(float));
out.write((char *)&y, sizeof(float));
out.write((char *)&z, sizeof(float));
}
else
out.write((char *)&zero, sizeof(int));
}
else
{
for (int k = 0; k < depth_; k++)
out << "*";
out << "i";
if (interinfo[l])
{
out << 1 << ' ' << *interinfo[l] << std::endl;
}
else
out << 0 << std::endl;
}
}
}
}
public:
int id_, depth_, cell_type_;
bool is_leaf_;
Node *parent;
Node *child[8];
TinyVector<float, 4> *interinfo[12];
};
//--------------------
std::vector<std::unordered_map<int, Node *>> layernodes(depth + 1);
// intialize
const int vertmap2dc[8] = {0, 4, 3, 7, 1, 5, 2, 6};
const int edgevmap[12][2] =
{{0, 4}, {1, 5}, {2, 6}, {3, 7}, {0, 2}, {1, 3}, {4, 6}, {5, 7}, {0, 1}, {2, 3}, {4, 5}, {6, 7}};
int corner[8];
for (size_t c = 0; c < valid_cell_id.size(); c++)
{
int n = valid_cell_id[c];
Node *node = new Node(n, depth, true, valid_cell_types[c]);
layernodes[depth][n] = node;
// todo add 12 edges here
for (int i = 0; i < 8; ++i)
{
corner[i] = get_point_id(n, i);
}
for (int i = 0; i < 12; i++)
{
auto iter = edgeintersections.find(
IntersectionInfo(corner[vertmap2dc[edgevmap[i][0]]], corner[vertmap2dc[edgevmap[i][1]]]));
if (iter != edgeintersections.end())
{
node->interinfo[i] = new TinyVector<float, 4>(iter->shift, iter->nx, iter->ny, iter->nz);
}
}
}
int curdepth = depth;
int res = gridsize - 1;
while (curdepth > 0)
{
int half_res = res / 2;
auto iter = layernodes[curdepth].begin();
for (; iter != layernodes[curdepth].end(); iter++)
{
int n = iter->first;
int k = n / (res * res);
int j = (n - k * res * res) / res;
int i = n % res;
int inck = k % 2, incj = j % 2, inci = i % 2;
k = k / 2, j = j / 2, i = i / 2;
int up_n = k * half_res * half_res + j * half_res + i;
Node *node = 0;
auto siter = layernodes[curdepth - 1].find(up_n);
if (siter == layernodes[curdepth - 1].end())
{
node = new Node(up_n, curdepth - 1);
layernodes[curdepth - 1][up_n] = node;
}
else
{
node = siter->second;
}
assert(node->child[inci * 4 + incj * 2 + inck] == 0);
node->child[inci * 4 + incj * 2 + inck] = iter->second;
}
curdepth--;
res /= 2;
}
///--------------------
std::ofstream dcffile;
if (binary)
dcffile.open(filename, std::ofstream::binary);
else
dcffile.open(filename);
const char version[10] = "multisign";
int dim[3];
dim[0] = dim[1] = dim[2] = gridsize - 1;
if (binary)
{
dcffile.write(version, sizeof(char) * 10);
dcffile.write((char *)dim, sizeof(int) * 3);
}
else
{
dcffile << version << std::endl;
dcffile << dim[0] << ' ' << dim[1] << ' ' << dim[2] << std::endl;
}
layernodes[0][0]->travel_print(dcffile, binary);
dcffile.close();
//--------------------
for (int i = 0; i < depth + 1; i++)
{
for (auto iter = layernodes[i].begin(); iter != layernodes[i].end(); iter++)
delete iter->second;
}
}
/////////////////////////////////////////////////
void dcf_meshing(const char plyname[])
{
if (valid_cell_id.empty())
{
std::cerr << "scan_octants() should be called first!" << std::endl;
return;
}
//--- compute intersection---
std::set<IntersectionInfo> edgeintersections;
compute_intersection(edgeintersections);
//---------------------------
class Node
{
public:
Node(int _id, int _depth, bool _is_leaf = false, int _cell_type = 0)
: id_(_id), depth_(_depth), is_leaf_(_is_leaf)
{
parent = 0;
cell_type_ = cubestatusmap[_cell_type];
for (int i = 0; i < 8; i++)
child[i] = 0;
for (int i = 0; i < 12; i++)
interinfo[i] = 0;
}
~Node()
{
for (int i = 0; i < 12; i++)
if (interinfo[i])
delete interinfo[i];
}
OctreeNode *travel(int *st, int len, int ht)
{
OctreeNode *rvalue = 0;
bool is_empty = true;
for (int i = 0; i < 8; i++)
{
if (child[i])
{
is_empty = false;
break;
}
}
int type = is_leaf_ ? 2 : (is_empty ? 1 : 0);
if (type == 0)
{
rvalue = new InternalNode();
int nlen = len / 2;
int nst[3];
for (int i = 0; i < 8; i++)
{
if (child[i])
{
nst[0] = st[0] + vertMap[i][0] * nlen;
nst[1] = st[1] + vertMap[i][1] * nlen;
nst[2] = st[2] + vertMap[i][2] * nlen;
((InternalNode *)rvalue)->child[i] = child[i]->travel(nst, nlen, ht - 1);
}
}
}
else if (type == 1)
{
return 0;
}
else if (type == 2)
{
float inters[12][3], norms[12][3];
int numinters = 0;
for (int l = 0; l < 12; l++)
{
if (interinfo[l])
{
int dir = l / 4;
int base = edgevmap[l][0];
inters[numinters][0] = st[0] + vertMap[base][0] * len;
inters[numinters][1] = st[1] + vertMap[base][1] * len;
inters[numinters][2] = st[2] + vertMap[base][2] * len;
inters[numinters][dir] += (*interinfo[l])[0];
norms[numinters][0] = (*interinfo[l])[1];
norms[numinters][1] = (*interinfo[l])[2];
norms[numinters][2] = (*interinfo[l])[3];
numinters++;
}
}
if (numinters > 0)
rvalue = new LeafNode(ht, (unsigned char)cell_type_, st, len, numinters, inters, norms);
}
return rvalue;
}
public:
int id_, depth_, cell_type_;
bool is_leaf_;
Node *parent;
Node *child[8];
TinyVector<float, 4> *interinfo[12];
};
//--------------------
std::vector<std::unordered_map<int, Node *>> layernodes(depth + 1);
// intialize
const int vertmap2dc[8] = {0, 4, 3, 7, 1, 5, 2, 6};
const int edgevmap[12][2] =
{{0, 4}, {1, 5}, {2, 6}, {3, 7}, {0, 2}, {1, 3}, {4, 6}, {5, 7}, {0, 1}, {2, 3}, {4, 5}, {6, 7}};
int corner[8];
for (size_t c = 0; c < valid_cell_id.size(); c++)
{
int n = valid_cell_id[c];
Node *node = new Node(n, depth, true, valid_cell_types[c]);
layernodes[depth][n] = node;
// todo add 12 edges here
for (int i = 0; i < 8; ++i)
{
corner[i] = get_point_id(n, i);
}
for (int i = 0; i < 12; i++)
{
auto iter = edgeintersections.find(
IntersectionInfo(corner[vertmap2dc[edgevmap[i][0]]], corner[vertmap2dc[edgevmap[i][1]]]));
if (iter != edgeintersections.end())
{
node->interinfo[i] = new TinyVector<float, 4>(iter->shift, iter->nx, iter->ny, iter->nz);
}
}
}
int curdepth = depth;
int res = gridsize - 1;
while (curdepth > 0)
{
int half_res = res / 2;
auto iter = layernodes[curdepth].begin();
for (; iter != layernodes[curdepth].end(); iter++)
{
int n = iter->first;
int k = n / (res * res);
int j = (n - k * res * res) / res;
int i = n % res;
int inck = k % 2, incj = j % 2, inci = i % 2;
k = k / 2, j = j / 2, i = i / 2;
int up_n = k * half_res * half_res + j * half_res + i;
Node *node = 0;
auto siter = layernodes[curdepth - 1].find(up_n);
if (siter == layernodes[curdepth - 1].end())
{
node = new Node(up_n, curdepth - 1);
layernodes[curdepth - 1][up_n] = node;
}
else
{
node = siter->second;
}
assert(node->child[inci * 4 + incj * 2 + inck] == 0);
node->child[inci * 4 + incj * 2 + inck] = iter->second;
}
curdepth--;
res /= 2;
}
///--------------------
Octree myOctree;
myOctree.dimen = gridsize - 1;
myOctree.hasQEF = 0;
myOctree.maxDepth = depth;
myOctree.scale = std::max(std::max(end_x - start_x, end_y - start_y), end_z - start_z) / pow(2, depth);
myOctree.shift[0] = end_x;
myOctree.shift[1] = end_y;
myOctree.shift[2] = end_z;
int st[3] = {0, 0, 0};
int len = myOctree.dimen;
int ht = myOctree.maxDepth;
myOctree.root = layernodes[0][0]->travel(st, len, ht);
// myOctree.genContour(plyname);
myOctree.genContourNoInter2(plyname);
//--------------------
Node *root = layernodes[0].begin()->second;
delete root;
// for (int i = 0; i < depth + 1; i++)
//{
// for (auto iter = layernodes[i].begin(); iter != layernodes[i].end(); iter++)
// delete iter->second;
// }
}
/////////////////////////////////////////////////
void save_vtk(const char filename[])
{
std::ofstream mfile(filename);
mfile << "# vtk DataFile Version 3.0\nImplicit Function\nASCII\nDATASET UNSTRUCTURED_GRID\n";
mfile << "POINTS " << point_value_map.size() << " float\n";
std::map<int, int> idmap;
TinyVector<Real, 3> pos;
int counter = 0;
for (auto iter = point_value_map.begin(); iter != point_value_map.end(); iter++)
{
int id = iter->first;
get_point(iter->first, pos);
mfile << pos << std::endl;
idmap[id] = counter++;
}
mfile << "CELLS " << valid_cell_id.size() << ' ' << valid_cell_id.size() * 9 << std::endl;
for (size_t c = 0; c < valid_cell_id.size(); c++)
{
mfile << 8;
int n = valid_cell_id[c];
int k = n / ((gridsize - 1) * (gridsize - 1));
int j = (n - k * (gridsize - 1) * (gridsize - 1)) / (gridsize - 1);
int i = n % (gridsize - 1);
n = k * gridsize * gridsize + j * gridsize + i;
mfile << ' ' << idmap[n];
mfile << ' ' << idmap[n + offsets_[1]];
mfile << ' ' << idmap[n + offsets_[2]];
mfile << ' ' << idmap[n + offsets_[3]];
mfile << ' ' << idmap[n + offsets_[4]];
mfile << ' ' << idmap[n + offsets_[5]];
mfile << ' ' << idmap[n + offsets_[6]];
mfile << ' ' << idmap[n + offsets_[7]] << std::endl;
}
mfile << "CELL_TYPES " << valid_cell_id.size() << std::endl;
for (size_t c = 0; c < valid_cell_id.size(); c++)
mfile << 12 << std::endl;
mfile << "POINT_DATA " << point_value_map.size() << std::endl;
mfile << "SCALARS funcvalue float 1\nLOOKUP_TABLE default\n";
for (auto iter = point_value_map.begin(); iter != point_value_map.end(); iter++)
{
mfile << iter->second << std::endl;
}
mfile.close();
}
};
//////////////////////////////////////