Commit e3c727f7 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Restructure OldDrawGrid example

Separating routines and classes in Examples/OldDrawGrid into files for future reuse. Making streamline calculator universal with respect to definition of data on various element types and adding option to visualize data as streamlines.
parent 61de3bb7
......@@ -3,12 +3,33 @@ set(SOURCE main.cpp
rotate.cpp
clipboard.cpp
rotate.h
my_glut.h
inc_glut.h
clipboard.h
color.h
coord.h
coord.cpp
octree.h)
octree.h
streamline.h
streamline.cpp
svg_line.h
svg_line.cpp
face2gl.h
face2gl.cpp
color_bar.h
color_bar.cpp
printtext.cpp
tga.h
tga.cpp
screenshot.h
screenshot.cpp
volumetric.h
volumetric.cpp
input.h
input.cpp
picker.h
picker.cpp
clipper.h
clipper.cpp)
find_package(OpenGL)
find_package(GLUT)
......
#include "clipper.h"
#include "inc_glut.h"
#include "color_bar.h"
#define CLIP_NONE 0
#define CLIP_NODE 1
#define CLIP_FULL 2
#define CLIP_ENDP 3
#define CLIP_FACE_NONE 0
#define CLIP_FACE_INSIDE 1
#define CLIP_FACE_OUTSIDE 2
#define CLIP_FACE_INTERSECT 3
namespace INMOST
{
kdtree::entry::entry(const entry & other)
{
e = other.e;
xyz[0] = other.xyz[0];
xyz[1] = other.xyz[1];
xyz[2] = other.xyz[2];
}
kdtree::entry & kdtree::entry::operator =(const entry & other)
{
e = other.e;
xyz[0] = other.xyz[0];
xyz[1] = other.xyz[1];
xyz[2] = other.xyz[2];
return *this;
}
static int cmpElements0(const void * a, const void * b)
{
const kdtree::entry * ea = ((const kdtree::entry *)a);
const kdtree::entry * eb = ((const kdtree::entry *)b);
float ad = ea->xyz[0];
float bd = eb->xyz[0];
return (ad > bd) - (ad < bd);
}
static int cmpElements1(const void * a, const void * b)
{
const kdtree::entry * ea = ((const kdtree::entry *)a);
const kdtree::entry * eb = ((const kdtree::entry *)b);
float ad = ea->xyz[1];
float bd = eb->xyz[1];
return (ad > bd) - (ad < bd);
}
static int cmpElements2(const void * a, const void * b)
{
const kdtree::entry * ea = ((const kdtree::entry *)a);
const kdtree::entry * eb = ((const kdtree::entry *)b);
float ad = ea->xyz[2];
float bd = eb->xyz[2];
return (ad > bd) - (ad < bd);
}
inline static unsigned int flip(const unsigned int * fp) { unsigned int mask = -((int)(*fp >> 31)) | 0x80000000; return *fp ^ mask; }
inline static unsigned int _0(unsigned int x) { return x & 0x7FF; }
inline static unsigned int _1(unsigned int x) { return x >> 11 & 0x7FF; }
inline static unsigned int _2(unsigned int x) { return x >> 22; }
Storage::integer clip_plane_edge(double sp0[3], double sp1[3], double p[3], double n[3], double node[3])
{
Storage::real u[3], w[3], D, N, sI;
u[0] = sp1[0] - sp0[0]; u[1] = sp1[1] - sp0[1]; u[2] = sp1[2] - sp0[2];
w[0] = sp0[0] - p[0]; w[1] = sp0[1] - p[1]; w[2] = sp0[2] - p[2];
D = (n[0] * u[0] + n[1] * u[1] + n[2] * u[2]);
N = -(n[0] * w[0] + n[1] * w[1] + n[2] * w[2]);
if (fabs(D) < 1.0e-9)
{
if (fabs(N) < 1.0e-9) return CLIP_FULL;
else return CLIP_NONE;
}
else
{
sI = N / D;
if (sI < 0 - 1.0e-9 || sI > 1 + 1.0e-9) return CLIP_NONE;
else
{
node[0] = sp0[0] + sI * u[0];
node[1] = sp0[1] + sI * u[1];
node[2] = sp0[2] + sI * u[2];
return (sI > 1.0e-9 && sI < 1.0 - 1.0e-9) ? CLIP_NODE : CLIP_ENDP;
}
}
}
void kdtree::radix_sort(int dim, struct entry * temp)
{
unsigned int i;
const unsigned int kHist = 2048;
unsigned int b0[kHist * 3];
unsigned int *b1 = b0 + kHist;
unsigned int *b2 = b1 + kHist;
memset(b0, 0, sizeof(unsigned int)*kHist * 3);
for (i = 0; i < size; i++)
{
unsigned int fi = flip((unsigned int *)&set[i].xyz[dim]);
++b0[_0(fi)]; ++b1[_1(fi)]; ++b2[_2(fi)];
}
{
unsigned int sum0 = 0, sum1 = 0, sum2 = 0;
for (i = 0; i < kHist; i++)
{
b0[kHist - 1] = b0[i] + sum0; b0[i] = sum0 - 1; sum0 = b0[kHist - 1];
b1[kHist - 1] = b1[i] + sum1; b1[i] = sum1 - 1; sum1 = b1[kHist - 1];
b2[kHist - 1] = b2[i] + sum2; b2[i] = sum2 - 1; sum2 = b2[kHist - 1];
}
}
for (i = 0; i < size; i++) temp[++b0[_0(flip((unsigned int *)&set[i].xyz[dim]))]] = set[i];
for (i = 0; i < size; i++) set[++b1[_1(flip((unsigned int *)&temp[i].xyz[dim]))]] = temp[i];
for (i = 0; i < size; i++) temp[++b2[_2(flip((unsigned int *)&set[i].xyz[dim]))]] = set[i];
for (i = 0; i < size; i++) set[i] = temp[i];
}
void kdtree::kdtree_build(int dim, int & done, int total, struct entry * temp)
{
if (size > 1)
{
if (size > 128) radix_sort(dim, temp); else
switch (dim)
{
case 0: qsort(set, size, sizeof(entry), cmpElements0); break;
case 1: qsort(set, size, sizeof(entry), cmpElements1); break;
case 2: qsort(set, size, sizeof(entry), cmpElements2); break;
}
children = static_cast<kdtree *>(malloc(sizeof(kdtree)* 2));//new kdtree[2];
assert(children != NULL);
children[0].marked = 0;
children[0].children = NULL;
children[0].set = set;
children[0].size = size / 2;
children[0].m = m;
children[1].marked = 0;
children[1].children = NULL;
children[1].set = set + size / 2;
children[1].size = size - size / 2;
children[1].m = m;
children[0].kdtree_build((dim + 1) % 3, done, total, temp);
children[1].kdtree_build((dim + 1) % 3, done, total, temp);
for (int k = 0; k < 3; k++)
{
bbox[0 + 2 * k] = std::min(children[0].bbox[0 + 2 * k], children[1].bbox[0 + 2 * k]);
bbox[1 + 2 * k] = std::max(children[0].bbox[1 + 2 * k], children[1].bbox[1 + 2 * k]);
}
}
else
{
assert(size == 1);
if (GetHandleElementType(set[0].e) == EDGE)
{
Storage::real_array n1 = Edge(m, set[0].e)->getBeg()->Coords();
Storage::real_array n2 = Edge(m, set[0].e)->getEnd()->Coords();
for (int k = 0; k < 3; k++)
{
bbox[0 + 2 * k] = std::min(n1[k], n2[k]);
bbox[1 + 2 * k] = std::max(n1[k], n2[k]);
}
done++;
if (done % 150 == 0)
{
printf("%3.1f%%\r", (done*100.0) / (total*1.0));
fflush(stdout);
}
}
else
{
ElementArray<Node> nodes = Element(m, set[0].e)->getNodes();
bbox[0] = bbox[2] = bbox[4] = 1.0e20;
bbox[1] = bbox[3] = bbox[5] = -1.0e20;
for (INMOST_DATA_ENUM_TYPE k = 0; k < nodes.size(); ++k)
{
Storage::real_array coords = nodes[k].Coords();
for (INMOST_DATA_ENUM_TYPE q = 0; q < 3; q++)
{
bbox[q * 2 + 0] = std::min<float>(bbox[q * 2 + 0], coords[q]);
bbox[q * 2 + 1] = std::max<float>(bbox[q * 2 + 1], coords[q]);
}
}
}
}
}
kdtree::kdtree() : marked(0), set(NULL), size(0), children(NULL) {}
inline int kdtree::plane_bbox(double p[3], double n[3]) const
{
Storage::real pv[3], nv[3];
for (int k = 0; k < 3; ++k)
{
if (n[k] >= 0)
{
pv[k] = bbox[1 + 2 * k]; //max
nv[k] = bbox[0 + 2 * k]; //min
}
else
{
pv[k] = bbox[0 + 2 * k]; //min
nv[k] = bbox[1 + 2 * k]; //max
}
}
Storage::real pvD, nvD;
pvD = n[0] * (pv[0] - p[0]) + n[1] * (pv[1] - p[1]) + n[2] * (pv[2] - p[2]);
nvD = n[0] * (nv[0] - p[0]) + n[1] * (nv[1] - p[1]) + n[2] * (nv[2] - p[2]);
if (nvD*pvD <= 0.0)
return 2;
else if (nvD < 0.0)
return 1;
else return 0;
}
bool kdtree::sub_intersect_plane_edge(Tag clip_point, Tag clip_state, ElementArray<Cell> & cells, MarkerType mrk, double p[3], double n[3])
{
if (size == 1)
{
assert(GetHandleElementType(set[0].e) == EDGE);
Edge ee = Edge(m, set[0].e);
Storage::real_array sp0 = ee->getBeg()->Coords();
Storage::real_array sp1 = ee->getEnd()->Coords();
Storage::integer & clip = m->IntegerDF(set[0].e, clip_state);
clip = clip_plane_edge(&sp0[0], &sp1[0], p, n, &ee->RealArrayDF(clip_point)[0]);
if (clip)
{
ElementArray<Cell> ecells = ee->getCells();
for (INMOST_DATA_ENUM_TYPE k = 0; k < ecells.size(); ++k) if (!ecells[k].GetMarker(mrk))
{
ecells[k].SetMarker(mrk);
cells.push_back(ecells[k]);
}
marked = 1;
}
}
else if (plane_bbox(p, n) == 2)
{
bool test1 = children[0].sub_intersect_plane_edge(clip_point, clip_state, cells, mrk, p, n);
bool test2 = children[1].sub_intersect_plane_edge(clip_point, clip_state, cells, mrk, p, n);
if (test1 || test2) marked = 1;
}
return marked != 0;
}
void kdtree::sub_intersect_plane_faces(Tag clip_state, double p[3], double n[3])
{
if (size == 1)
{
Storage::integer state;
Element ee(m, set[0].e);
assert(ee->GetElementDimension() == 2);
ElementArray<Node> nodes = ee->getNodes();
Storage::real_array coords = nodes[0].Coords();
Storage::real dot0 = n[0] * (coords[0] - p[0]) + n[1] * (coords[1] - p[1]) + n[2] * (coords[2] - p[2]);
if (dot0 <= 0.0) state = CLIP_FACE_INSIDE; else state = CLIP_FACE_OUTSIDE;
for (INMOST_DATA_ENUM_TYPE k = 1; k < nodes.size(); k++)
{
coords = nodes[k].Coords();
Storage::real dot = n[0] * (coords[0] - p[0]) + n[1] * (coords[1] - p[1]) + n[2] * (coords[2] - p[2]);
if (dot*dot0 <= 0.0)
{
state = CLIP_FACE_INTERSECT;
break;
}
}
m->IntegerDF(set[0].e, clip_state) = state;
}
else
{
marked = plane_bbox(p, n);
if (marked == 0)
{
for (INMOST_DATA_ENUM_TYPE k = 0; k < size; k++) m->IntegerDF(set[k].e, clip_state) = CLIP_FACE_OUTSIDE;
}
else if (marked == 1)
{
for (INMOST_DATA_ENUM_TYPE k = 0; k < size; k++) m->IntegerDF(set[k].e, clip_state) = CLIP_FACE_INSIDE;
}
else
{
children[0].sub_intersect_plane_faces(clip_state, p, n);
children[1].sub_intersect_plane_faces(clip_state, p, n);
}
}
}
void kdtree::unmark_old_edges(Tag clip_state)
{
if (size == 1)
{
assert(GetHandleElementType(set[0].e) == EDGE);
marked = 0;
if (GetHandleElementType(set[0].e) == EDGE)
m->IntegerDF(set[0].e, clip_state) = CLIP_NONE;
else if (GetHandleElementType(set[0].e) == FACE)
m->IntegerDF(set[0].e, clip_state) = CLIP_FACE_NONE;
}
else if (children)
{
if (children[0].marked) { children[0].unmark_old_edges(clip_state); marked = 0; }
if (children[1].marked) { children[1].unmark_old_edges(clip_state); marked = 0; }
}
}
void kdtree::clear_children() { if (children) { children[0].clear_children(); children[1].clear_children(); free(children); } }
kdtree::kdtree(Mesh * m) : marked(0), m(m), children(NULL)
{
double tt;
size = m->NumberOfEdges();
assert(size > 1);
set = new kdtree::entry[size];
INMOST_DATA_ENUM_TYPE k = 0;
tt = Timer();
printf("Prepearing edge set.\n");
for (Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
{
set[k].e = *it;
set[k].xyz[0] = (it->getBeg()->Coords()[0] + it->getEnd()->Coords()[0])*0.5;
set[k].xyz[1] = (it->getBeg()->Coords()[1] + it->getEnd()->Coords()[1])*0.5;
set[k].xyz[2] = (it->getBeg()->Coords()[2] + it->getEnd()->Coords()[2])*0.5;
k++;
if (k % 150 == 0)
{
printf("%3.1f%%\r", (k*100.0) / (size*1.0));
fflush(stdout);
}
}
printf("Done. Time %lg\n", Timer() - tt);
int done = 0, total = size;
printf("Building KD-tree.\n");
tt = Timer();
struct entry * temp = new entry[size];
kdtree_build(0, done, total, temp);
delete[] temp;
printf("Done. Time %lg\n", Timer() - tt);
for (int k = 0; k < 3; k++)
{
bbox[0 + 2 * k] = std::min(children[0].bbox[0 + 2 * k], children[1].bbox[0 + 2 * k]);
bbox[1 + 2 * k] = std::max(children[0].bbox[1 + 2 * k], children[1].bbox[1 + 2 * k]);
}
}
kdtree::kdtree(Mesh * m, HandleType * eset, INMOST_DATA_ENUM_TYPE size) : marked(0), m(m), size(size), children(NULL)
{
double tt;
assert(size > 1);
set = new entry[size];
tt = Timer();
printf("Prepearing elements set.\n");
for (INMOST_DATA_ENUM_TYPE k = 0; k < size; k++)
{
set[k].e = eset[k];
Storage::real cnt[3];
m->GetGeometricData(set[k].e, CENTROID, cnt);
set[k].xyz[0] = cnt[0];
set[k].xyz[1] = cnt[1];
set[k].xyz[2] = cnt[2];
if (k % 150 == 0)
{
printf("%3.1f%%\r", (k*100.0) / (size*1.0));
fflush(stdout);
}
}
printf("Done. Time %lg\n", Timer() - tt);
int done = 0, total = size;
printf("Building KD-tree.\n");
tt = Timer();
struct entry * temp = new entry[size];
kdtree_build(0, done, total, temp);
delete[] temp;
printf("Done. Time %lg\n", Timer() - tt);
for (int k = 0; k < 3; k++)
{
bbox[0 + 2 * k] = std::min(children[0].bbox[0 + 2 * k], children[1].bbox[0 + 2 * k]);
bbox[1 + 2 * k] = std::max(children[0].bbox[1 + 2 * k], children[1].bbox[1 + 2 * k]);
}
}
void kdtree::intersect_plane_edge(Tag clip_point, Tag clip_state, ElementArray<Cell> & cells, MarkerType mark_cells, double p[3], double n[3])
{
if (marked)
{
unmark_old_edges(clip_state);
cells.clear();
}
sub_intersect_plane_edge(clip_point, clip_state, cells, mark_cells, p, n);
}
void kdtree::intersect_plane_face(Tag clip_state, double p[3], double n[3])
{
sub_intersect_plane_faces(clip_state, p, n);
}
kdtree::~kdtree()
{
delete[] set;
clear_children();
}
clipper::edge_point::edge_point(Storage::real _xyz[3], Storage::integer n, float v)
{
xyz[0] = _xyz[0];
xyz[1] = _xyz[1];
xyz[2] = _xyz[2];
edge = n;
val = v;
}
clipper::edge_point::edge_point(){}
bool clipper::edge_point::operator ==(const edge_point& b) const
{
Storage::real temp = 0.0;
for (int k = 0; k < 3; k++) temp += (xyz[k] - b.xyz[k])*(xyz[k] - b.xyz[k]);
if (temp < 1.0e-8) return true; else return false;
}
bool clipper::edge_point::operator !=(const edge_point& b) const { return !(operator ==(b)); }
void clipper::edge_point::print() { printf("%g %g %g e %d\n", xyz[0], xyz[1], xyz[2], edge); }
clipper::~clipper()
{
delete tree;
for (INMOST_DATA_ENUM_TYPE k = 0; k < cells.size(); k++) cells[k]->RemMarker(marker);
mm->ReleaseMarker(marker);
mm->DeleteTag(clips);
mm->DeleteTag(clipsv);
mm->DeleteTag(clipsn);
mm->DeleteTag(clip_point);
mm->DeleteTag(clip_state);
}
clipper::clipper(Mesh * m)
{
mm = m;
cells.SetMeshLink(mm);
tree = new kdtree(m);
marker = m->CreateMarker();
clips = m->CreateTag("CLIPS", DATA_REAL, CELL, CELL);
clipsv = m->CreateTag("CLIPS_VAL", DATA_REAL, CELL, CELL);
clipsn = m->CreateTag("CLIPS_NUM", DATA_INTEGER, CELL, CELL);
clip_point = m->CreateTag("CLIP_POINT", DATA_REAL, EDGE, NONE, 3);
clip_state = m->CreateTag("CLIP_STATE", DATA_INTEGER, EDGE, NONE, 1);
}
double clipper::compute_value(Edge e, Storage::real * pnt)
{
if (isColorBarEnabled())
{
Storage::real_array c1 = e->getBeg()->Coords();
Storage::real_array c2 = e->getEnd()->Coords();
Storage::real d1, d2, t;
d1 = sqrt((pnt[0] - c1[0])*(pnt[0] - c1[0]) + (pnt[1] - c1[1])*(pnt[1] - c1[1]) + (pnt[2] - c1[2])*(pnt[2] - c1[2]));
d2 = sqrt((c2[0] - c1[0])*(c2[0] - c1[0]) + (c2[1] - c1[1])*(c2[1] - c1[1]) + (c2[2] - c1[2])*(c2[2] - c1[2]));
t = d1 / d2; //(pnt == c2, t = 1 : pnt == c1, t = 0)
return e->getBeg()->RealDF(GetVisualizationTag())*(1 - t) + e->getEnd()->RealDF(GetVisualizationTag())*t;
}
else return 0.f;
}
void clipper::clip_plane(Storage::real p[3], Storage::real n[3])
{
const bool print = false;
for (INMOST_DATA_ENUM_TYPE k = 0; k < cells.size(); ++k)
if (cells[k]->GetMarker(marker))
{
cells[k]->RealArray(clips).clear();
cells[k]->RealArray(clipsv).clear();
cells[k]->IntegerArray(clipsn).clear();
cells[k]->RemMarker(marker);
}
tree->intersect_plane_edge(clip_point, clip_state, cells, marker, p, n);
dynarray<edge_point, 128> clipcoords, loopcoords;
std::vector<bool> closed;
for (INMOST_DATA_ENUM_TYPE k = 0; k < cells.size(); ++k)
{
//assuming faces are convex we will have at most one clipping edge per polygon
//otherwise every pair of clipping nodes forming edge should appear consequently
//as long as we go through face's edges in ordered way
clipcoords.clear();
//we will gather all the pairs of nodes, then form closed loop
ElementArray<Face> faces = cells[k]->getFaces();
Face full_face;
int ntotpoints = 0, ntotedges = 0;
for (INMOST_DATA_ENUM_TYPE q = 0; q < faces.size(); ++q)
{
int last_edge_type = CLIP_NONE;
int nfulledges = 0, npoints = 0, nstartedge = ntotedges;
ElementArray<Edge> edges = faces[q].getEdges();
for (INMOST_DATA_ENUM_TYPE r = 0; r < edges.size(); ++r)
{
Storage::integer state = edges[r].IntegerDF(clip_state);
if (state == CLIP_FULL)
{
nfulledges++;
edge_point n1 = edge_point(&edges[r].getBeg()->Coords()[0], ntotedges, isColorBarEnabled() ? edges[r].getBeg()->RealDF(GetVisualizationTag()) : 0.f);
edge_point n2 = edge_point(&edges[r].getEnd()->Coords()[0], ntotedges, isColorBarEnabled() ? edges[r].getEnd()->RealDF(GetVisualizationTag()) : 0.f);
if (npoints % 2 == 0) //all privious edges are closed, just add this one
{
clipcoords.push_back(n1);
clipcoords.push_back(n2);
npoints += 2;
ntotedges++;
last_edge_type = CLIP_FULL;
}
else if (n1 == clipcoords.back()) //this may be prolongation of one point that hit one edge
{
clipcoords.push_back(n2);
npoints++;
ntotedges++;
last_edge_type = CLIP_FULL;
}
else if (n2 == clipcoords.back())
{
clipcoords.push_back(n1);
npoints++;
ntotedges++;
last_edge_type = CLIP_FULL;
}
else printf("%s:%d strange orphan node before me\n", __FILE__, __LINE__);
}
else if (state == CLIP_ENDP)
{
edge_point n = edge_point(&edges[r].RealArrayDF(clip_point)[0], ntotedges, compute_value(edges[r], &edges[r].RealArrayDF(clip_point)[0]));
bool add = true;
if (last_edge_type == CLIP_ENDP)
{
if (n == clipcoords.back())
add = false;
}
else if (last_edge_type == CLIP_FULL)