Commit bf10cd3e authored by Edward Vigmond's avatar Edward Vigmond
Browse files

Preliminary support for isosurfaces. Isolines for triangles are working but

nothing else.
parent 0e76ef0c
......@@ -697,7 +697,7 @@ ClipPlane::ClipPlane( TBmeshWin *tb ) {
} // Fl_Group* Clip2
{ Clip3 = new Fl_Group(5, 250, 185, 220, "Clipping Plane 4");
Clip3->box(FL_UP_BOX);
Clip3->align(Fl_Align(FL_ALIGN_BOTTOM));
Clip3->align(FL_ALIGN_BOTTOM);
{ Fl_Group* o = new Fl_Group(15, 305, 115, 145, "Plane Normal4");
o->box(FL_UP_BOX);
{ xdir3 = new Fl_Value_Input(35, 315, 80, 25, "X:");
......@@ -739,7 +739,7 @@ ClipPlane::ClipPlane( TBmeshWin *tb ) {
} // Fl_Group* Clip3
{ Clip4 = new Fl_Group(195, 250, 185, 220, "Clipping Plane 5");
Clip4->box(FL_UP_BOX);
Clip4->align(Fl_Align(FL_ALIGN_BOTTOM));
Clip4->align(FL_ALIGN_BOTTOM);
{ Fl_Group* o = new Fl_Group(205, 305, 115, 145, "Plane Normal5");
o->box(FL_UP_BOX);
{ xdir4 = new Fl_Value_Input(225, 315, 80, 25, "X:");
......@@ -781,7 +781,7 @@ ClipPlane::ClipPlane( TBmeshWin *tb ) {
} // Fl_Group* Clip4
{ Clip5 = new Fl_Group(385, 250, 185, 220, "Clipping Plane 6");
Clip5->box(FL_UP_BOX);
Clip5->align(Fl_Align(FL_ALIGN_BOTTOM));
Clip5->align(FL_ALIGN_BOTTOM);
{ Fl_Group* o = new Fl_Group(395, 305, 115, 145, "Plane Normal6");
o->box(FL_UP_BOX);
{ xdir5 = new Fl_Value_Input(415, 315, 80, 25, "X:");
......
......@@ -100,13 +100,13 @@ DataOpacity::DataOpacity( TBmeshWin *tbmw ) {
minopacval->type(3);
minopacval->value(0.5);
minopacval->callback((Fl_Callback*)cb_minopacval);
minopacval->align(Fl_Align(FL_ALIGN_LEFT));
minopacval->align(FL_ALIGN_LEFT);
} // Fl_Value_Slider* minopacval
{ maxopacval = new Fl_Value_Slider(125, 100, 120, 30, "Maximum Opacity");
maxopacval->type(3);
maxopacval->value(1);
maxopacval->callback((Fl_Callback*)cb_maxopacval);
maxopacval->align(Fl_Align(FL_ALIGN_LEFT));
maxopacval->align(FL_ALIGN_LEFT);
} // Fl_Value_Slider* maxopacval
{ onbut = new Fl_Light_Button(225, 10, 120, 30, "Data Opacity ");
onbut->labelsize(16);
......
......@@ -5,6 +5,8 @@
enum lpint_enum{ BOTH_ON_PLANE, NO_INTERSECTION };
static const int simple_index[] = { 0,1,2,3,4,5,6,7,8,9,10,11,12,13};
/** find the intersection of a plane with a line
*
* \param a first point of line
......@@ -144,7 +146,7 @@ void SurfaceElement::read_normals( int e0, int e1, const char *fnb )
* \pre \p nl must be at least _ptsPerObj*\p n long
* \post \p _n is allocated if need be
*/
void MultiPoint::define( int *nl, int n )
void MultiPoint::define( const int *nl, int n )
{
delete[] _node;
_node = new int[n*_ptsPerObj];
......@@ -308,3 +310,95 @@ VolElement::planecut( char *pd, GLfloat* cp,
return se;
}
/** linearly interpolate the position of a value along an element edge
*
* \param p0 vertex 0
* \param dat0 data on p0
* \param p1 vertex 1
* \param dat1 data on p1
* \param val data between dat0 and dat1
* \param pint location of val (computed)
*/
void
edge_interp( const GLfloat *p0, DATA_TYPE dat0, const GLfloat *p1,
DATA_TYPE dat1, DATA_TYPE val, GLfloat* pint )
{
Vector3D<GLfloat> b = p1;
Vector3D<GLfloat> e = b - p0;
e *= 1.- (val-dat0)/(dat1-dat0);
b -= e;
pint[0] = b.X();
pint[1] = b.Y();
pint[2] = b.Z();
}
/** determine the isosurface for a multipoint object
*
* A row in a table is determined from the nodes above the threshold value
* Each row is of the form \n
* \#polygon \#sides_poly0 edge0_node0 edge0_node1 edge1_node0
* edge1_ node1 ... edgeN_node1 \#sides_poly1 edge0_node0 ...
* edgeN_node2
*
* \param dat data for all the nodes
* \param val value for the isosurface
* \param npoly number of polygons
*
* \return a list of element pointers
* \post npoly is the number of elements in the list
*/
MultiPoint ** MultiPoint::isosurf( DATA_TYPE *dat, DATA_TYPE val, int &npoly )
{
// determine row index into table
unsigned int index=0;
for( int i=_ptsPerObj-1; i>=0; i-- ) {
index <<= 1;
if( dat[_node[i]]>val )
index += 1;
}
const int* poly = iso_polys(index);
npoly = poly[0];// number of polygons to create
int poly_start = 1; // first polygon defined after \#polygons
MultiPoint **isoele = new MultiPoint *[poly[0]]; //element pointer list
for( int n=0; n<npoly; n++ ) {
int npts = poly[poly_start]; // \#nodes defining polygon
GLfloat *pt = new GLfloat[npts*3]; // local point list
for( int i=0; i<npts; i++ ) {
int pindex = poly_start+1+i*2;
edge_interp( (*_pt)[_node[poly[pindex]]], dat[_node[poly[pindex]]],
(*_pt)[_node[poly[pindex+1]]], dat[_node[poly[pindex+1]]],
val, pt+3*i );
}
Point *pts = new Point;
pts->add( pt, npts );
switch(poly[poly_start]) {
case 1:
assert(0);
break;
case 2:
isoele[n] = new Connection( pts );
break;
case 3:
isoele[n] = new Triangle( pts );
break;
case 4:
isoele[n] = new Quadrilateral( pts );
break;
default:
isoele[n] = new PolyGon( pts, npts );
break;
}
isoele[n]->define( simple_index );
poly_start += npts*2+1;
}
return isoele;
}
......@@ -54,6 +54,7 @@ class Point: public DrawingObj
const GLfloat* offset() const { return _offset; }
void base1(bool b){ _base1 = b; }
void add( GLfloat *, int n=1 );
const GLfloat* operator[] (int i){ return _pts+3*i; }
private:
GLfloat* _pts; //!< point list
vector<bool>*_visible; //!< points which get drawn
......@@ -66,30 +67,36 @@ class Point: public DrawingObj
class MultiPoint : public DrawingObj
{
public:
MultiPoint( Point *p, int n ):_pt(p),_ptsPerObj(n),_node(NULL){}
MultiPoint( Point *p, int n, int e ):_pt(p),_ptsPerObj(n),_node(NULL),
_nedge(e){}
~MultiPoint(){ if ( _n ) delete[] _node; }
const int* obj( int n=0 ) { return _node+n*_ptsPerObj; }
int ptsPerObj(){ return _ptsPerObj; }
void register_vertices(int, int, vector<bool>& );
void add( int *n );
const Point* pt(){ return _pt; }
void define( int *nl, int n=1 );
const int* obj( int n=0 ) { return _node+n*_ptsPerObj; }
int ptsPerObj(){ return _ptsPerObj; }
void register_vertices(int, int, vector<bool>& );
void add( int *n );
const Point* pt(){ return _pt; }
void define( const int *nl, int n=1 );
MultiPoint **isosurf( DATA_TYPE *d, DATA_TYPE val, int & );
virtual const int*iso_polys(unsigned int)=0;
protected:
int * _node; //!< list of nodes defining objects
Point* _pt; //!< pointer to point list
int _ptsPerObj; //!< \#nodes to define one object
int _nedge; //!< \#edges
};
class Connection : public MultiPoint
{
public:
Connection(Point *p):MultiPoint(p,2) {}
Connection(Point *p):MultiPoint(p,2,1) {}
void add( int, int ); //!< add a connection
virtual void draw( int, GLfloat*, float size=1 );
virtual void draw( int, int, GLfloat*, Colourscale*, DATA_TYPE*,
int stride=1, dataOpac* dopac=NULL );
virtual DrawingObj *isosurf( DATA_TYPE *d, DATA_TYPE val ){}
virtual bool read( const char * );
const int *iso_polys(unsigned int index){return NULL;}
};
......@@ -97,7 +104,7 @@ class Connection : public MultiPoint
class ContCable : public MultiPoint
{
public:
ContCable(Point *p):MultiPoint(p,1) {}
ContCable(Point *p):MultiPoint(p,1,1) {}
const int* obj( int n=0 ) { return _node+n; }
void add( int, int );
virtual void draw( int, GLfloat*, float size=1 );
......@@ -105,13 +112,15 @@ class ContCable : public MultiPoint
int stride=1, dataOpac* dopac=NULL );
virtual bool read( const char * );
void register_vertices( int, int, vector<bool>& );
const int *iso_polys(unsigned int index){return NULL;}
};
// closed convex polygons
class SurfaceElement : public MultiPoint
{
public:
SurfaceElement(Point *p, int n):MultiPoint(p,n),_nrml(NULL) {}
SurfaceElement(Point *p, int n):MultiPoint(p,n,n),_nrml(NULL) {}
virtual void compute_normals( int, int )=0;
const void nrml( GLfloat *n ){ _nrml=n; };
inline const GLfloat* nrml( int a=0 ) {return _nrml==NULL?NULL:_nrml+3*a; }
......@@ -131,12 +140,14 @@ class PolyGon : public SurfaceElement
public:
PolyGon( Point *p, int n ):SurfaceElement(p,n) {}
virtual void compute_normals( int a, int b );
virtual void draw( int, GLfloat*, float size=1 );
virtual void draw( int, int, GLfloat*, Colourscale*, DATA_TYPE*,
virtual void draw( int, GLfloat*, float size=1 );
virtual void draw( int, int, GLfloat*, Colourscale*, DATA_TYPE*,
int stride, dataOpac* dopac, const GLfloat * );
virtual void draw( int, int, GLfloat*, Colourscale*, DATA_TYPE*,
virtual void draw( int, int, GLfloat*, Colourscale*, DATA_TYPE*,
int stride=1, dataOpac* dopac=NULL );
virtual bool read( const char * ){}};
virtual bool read( const char * ){}
const int* iso_polys(unsigned int index){return NULL;}
};
class Triangle : public SurfaceElement
......@@ -150,7 +161,9 @@ class Triangle : public SurfaceElement
int stride, dataOpac* dopac, const GLfloat *);
virtual bool read( const char * );
virtual void compute_normals( int, int );
virtual DrawingObj *isosurf( DATA_TYPE *d, DATA_TYPE val ){}
bool add( const char * );
const int* iso_polys(unsigned int);
protected:
int countInFile( const char * );
};
......@@ -168,6 +181,7 @@ class Quadrilateral : public SurfaceElement
virtual bool read( const char * );
virtual void compute_normals( int, int );
bool add( const char * );
const int* iso_polys(unsigned int);
protected:
int countTrisInFile( const char * );
};
......@@ -177,7 +191,7 @@ class Quadrilateral : public SurfaceElement
class VolElement : public MultiPoint
{
public:
VolElement( Point *p, int n ):MultiPoint( p, n ) {}
VolElement( Point *p, int n, int e ):MultiPoint( p, n, e ) {}
const int* region() const { return _region; }
int region(int a) const { return _region[a]; }
void region(int a, int b) { _region[a] = b; }
......@@ -194,50 +208,53 @@ class VolElement : public MultiPoint
class Tetrahedral : public VolElement
{
public:
Tetrahedral(Point *p ): VolElement(p,4) {}
Tetrahedral(Point *p ): VolElement(p,4,6) {}
virtual void draw( int, GLfloat*, float size=1 );
virtual void draw( int, int, GLfloat*, Colourscale*, DATA_TYPE*,
int stride=1, dataOpac* dopac=NULL );
virtual bool read( const char * );
virtual void draw_out_face( int );
virtual SurfaceElement* cut( char*, GLfloat* cp, Interpolator<DATA_TYPE>*&,int=0);
virtual SurfaceElement* cut(char*,GLfloat*,Interpolator<DATA_TYPE>*&,int=0);
const int* iso_polys( unsigned int );
};
class Prism : public VolElement
{
public:
Prism(Point *p ): VolElement(p,6) {}
Prism(Point *p ): VolElement(p,6,9) {}
virtual void draw( int, GLfloat*, float size=1 );
virtual void draw( int, int, GLfloat*, Colourscale*, DATA_TYPE*,
int stride=1, dataOpac* dopac=NULL );
virtual bool read( const char * );
virtual void draw_out_face( int );
virtual SurfaceElement* cut( char *pd, GLfloat* cp, Interpolator<DATA_TYPE>*&,int );
virtual SurfaceElement* cut(char *,GLfloat*,Interpolator<DATA_TYPE>*&,int);
const int* iso_polys(unsigned int);
};
class Hexahedron : public VolElement
{
public:
Hexahedron(Point *p ): VolElement(p,8) {}
Hexahedron(Point *p ): VolElement(p,8,12) {}
virtual void draw( int, GLfloat*, float size=1 );
virtual void draw( int, int, GLfloat*, Colourscale*, DATA_TYPE*,
int stride=1, dataOpac* dopac=NULL );
virtual bool read( const char * );
virtual void draw_out_face( int );
virtual SurfaceElement* cut( char * pd, GLfloat* cp, Interpolator<DATA_TYPE>*&,int );
virtual SurfaceElement* cut(char *,GLfloat*,Interpolator<DATA_TYPE>*&,int);
const int* iso_polys( unsigned int );
};
class Pyramid : public VolElement
{
public:
Pyramid(Point *p ): VolElement(p,5) {}
Pyramid(Point *p ): VolElement(p,5,8) {}
virtual void draw( int, GLfloat*, float size=1 );
virtual void draw( int, int, GLfloat*, Colourscale*, DATA_TYPE*,
int stride=1, dataOpac* dopac=NULL );
virtual bool read( const char * );
virtual void draw_out_face( int );
virtual SurfaceElement* cut( char* pd, GLfloat* cp, Interpolator<DATA_TYPE>*&,int );
virtual SurfaceElement* cut(char*, GLfloat*, Interpolator<DATA_TYPE>*&,int);
const int* iso_polys( unsigned int );
};
#endif
......@@ -3,9 +3,9 @@
const int hex_num_edge=12;
const int hex_edges[hex_num_edge][2] =
{
{0,1}
,{1,2},{2,3},{3,0},{0,4},{1,7},{2,6},{3,5},{4,5},{5,6},{6,7},{7,4}
{0,1} ,{1,2},{2,3},{3,0},{0,4},{1,7},{2,6},{3,5},{4,5},{5,6},{6,7},{7,4}
};
const int hex_iso_table[][3] = { {1,2,3} };
/** draw many Hexahedrons
*
......@@ -245,3 +245,14 @@ Hexahedron::cut( char*pd, GLfloat* cp,
return planecut( pd, cp, interp, hex_num_edge, hex_edges, e );
}
/* return the list of intersection polygons for isosurfaces
*
* \param index for each node, the bit is true if the isoval is exceeded
*
* \return the row in the table
*/
const int* Hexahedron::iso_polys(unsigned int index)
{
return hex_iso_table[index];
}
#include "IsoLines.h"
int
IsoLine :: process( Surfaces *s, DATA_TYPE *dat )
{
int num_lines=0;
for( int i=0; i<_nl; i++ ){
double val = _nl==1? _v0: _v0 + i*(_v1-_v0)/(float)(_nl-1.);
for( int j=0; j<s->num(); j++ ) {
int npoly;
MultiPoint **lpoly = s->ele(j)->isosurf( dat, val, npoly );
for( int j=0; j<npoly; j++ ) {
_polygon.push_back(lpoly[j]);
num_lines++;
}
}
}
return num_lines;
}
IsoLine::~IsoLine()
{
for( int i=0; i<_polygon.size(); i++ ) {
delete _polygon[i]->pt();
delete _polygon[i];
}
}
void IsoLine::draw( Colourscale *cs, GLfloat size )
{
for( int i=0; i<_polygon.size(); i++ ) {
_polygon[i]->draw(0, _color, size );
}
}
#ifndef ISOLINE_H
#define ISOLINE_H
#include "Model.h"
#include "DrawingObjects.h"
#include "Surfaces.h"
#include <vector>
class IsoLine {
public:
IsoLine( double v0, double v1, int n, int t ):_v0(v0),_v1(v1),_nl(n),
_t(t) {_color[0]=_color[1]=_color[2]=0;_color[3]=1.;}
~IsoLine();
int process( Surfaces *, DATA_TYPE * );
int nl(){ return _nl; }
int v0(){ return _v0; }
int v1(){ return _v1; }
int tm(){ return _t; }
void draw( Colourscale *, GLfloat );
private:
vector<MultiPoint *> _polygon;
GLfloat _color[4];
double _v0, _v1;
int _nl;
int _t;
};
#endif
#include "IsoSurface.h"
IsoSurface::IsoSurface(Model* m, DATA_TYPE *dat, double v, vector<bool>&member)
:_val(v)
{
for( int i=0; i<m->numVol(); i++ ) {
if( member[i] ) {
int npoly;
MultiPoint **lpoly = m->_vol[i]->isosurf( dat, _val, npoly );
for( int j=0; j<npoly; j++ )
polygon.push_back(lpoly[j]);
}
}
}
IsoSurface::~IsoSurface()
{
for( int i=0; i<polygon.size(); i++ ) {
delete polygon[i]->pt();
delete polygon[i];
}
}
void IsoSurface::draw()
{
glPolygonMode( GL_FRONT_AND_BACK, GL_FILL );
for( int i=0; i<polygon.size(); i++ ) {
polygon[i]->draw(0, colour, 1 );
}
}
#ifndef ISOSURFACE_H
#define ISOSURFACE_H
#include "Model.h"
#include "DrawingObjects.h"
#include <vector>
class IsoSurface {
public:
IsoSurface( Model *m, DATA_TYPE *dat, double v, vector<bool>& );
~IsoSurface();
void draw();
GLfloat *color(){ return colour; }
double isoval(){ return _val; }
void color(const GLfloat *c){memmove(colour,c, sizeof(colour) );}
private:
vector<MultiPoint *> polygon;
GLfloat colour[4];
double _val;
};
#endif
......@@ -36,6 +36,9 @@ OBJS = Fl_Gl_Tb_Window.o\
Hexahedron.o\
HiLiteWinInfo.o\
IGBheader.o\
isosurf.o\
IsoLines.o\
IsoSurface.o\
main.o\
Matrix4x4.o\
Model.o\
......
......@@ -186,7 +186,7 @@ int Model::add_surface_from_elem( const char *fn )
gzFile in;
string fname(fn);
if ( (in=gzopen( (fname+="elem").c_str(), "r" )) == NULL )
if ( (in=gzopen( (fname+=".gz").c_str(), "r" )) == NULL ) {
if ((in=gzopen( (fname+=".gz").c_str(), "r" )) == NULL ) {
return -1;
}
......
......@@ -3,9 +3,9 @@
const int prism_num_edge=9;
const int prism_edges[prism_num_edge][2] =
{
{0,1}
,{1,2},{2,0},{0,3},{1,5},{2,4},{3,4},{4,5},{5,3}
{0,1},{1,2},{2,0},{0,3},{1,5},{2,4},{3,4},{4,5},{5,3}
};
const int prism_iso_table[][3] = { {1,2,3},{2,3,4} };
/** draw many Prisms
*
......@@ -217,3 +217,9 @@ Prism::cut( char* pd, GLfloat* cp,
return planecut( pd, cp, interp, prism_num_edge, prism_edges, e );
}
const int*
Prism::iso_polys(unsigned int index)
{
return prism_iso_table[index];
}
......@@ -3,9 +3,9 @@
const int pyr_num_edge=8;
static int pyr_edges[pyr_num_edge][2] =
{
{0,1}
,{1,2},{2,3},{3,0},{0,4},{1,4},{2,4},{3,4}
{0,1},{1,2},{2,3},{3,0},{0,4},{1,4},{2,4},{3,4}
};
static const int pyramid_iso_table[][2] = { {2,3},{3,4} };
/** draw many Pyramids
*
......@@ -191,3 +191,15 @@ Pyramid::cut( char *pd, GLfloat* cp,
return planecut( pd, cp, interp, pyr_num_edge, pyr_edges, e );
}
/* return the list of intersection polygons for isosurfaces
*
* \param index for each node, the bit is true if the isoval is exceeded
*
* \return the row in the table
*/
const int* Pyramid::iso_polys(unsigned int index)
{
return pyramid_iso_table[index];
}
#include "DrawingObjects.h"
static int quad_iso_table[][3] = {{1,2,3}};
/** draw many Quadrilaterals
*
* \param p0 first index of point to draw
......@@ -172,3 +174,10 @@ void Quadrilateral::compute_normals(int e0, int e1)
}
}
const int*
Quadrilateral::iso_polys(unsigned int index)
{
return quad_iso_table[index];
}
......@@ -7,7 +7,7 @@ int RRegion_sort( const void *a, const void *b )
return (*(RRegion **)a)->label() - (*(RRegion **)b)->label();
}
void RRegion:: initialize( int n, int l )
void RRegion:: initialize( int n, int nvol, int l )
{
is_visible = true;
_label = l;
......@@ -21,6 +21,8 @@ void RRegion:: initialize( int n, int l )
memset( _3D, 0, maxobject*sizeof(bool) );
_member.resize(n);
_member.assign(n,false);
_elemember.resize(nvol);
_elemember.assign(nvol,false);
for( int i=0; i<maxobject; i++ ) _size[i] = 1.;
}
......@@ -36,11 +38,12 @@ RRegion::RRegion( VolElement **v, int nv, int n, int l )
{
startind[Tetrahedron] = -1;