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.
1159 lines
35 KiB
1159 lines
35 KiB
#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();
|
|
}
|
|
};
|
|
//////////////////////////////////////
|
|
|