SimpleStreamline.cpp

This example calculates stream lines of a vector field. A startpoint or a fan of startpoints can be specified. The update function calcultes the streamlines and stored them as a std::vector of Lines, which are cached in respect to the vectorfield and stream line parameters.

Start e.g.:

../../../bin/vish SimpleStreamline.vis

See also:

Required tutorials: - InterpolatedVector.cpp

Further tutorials: - TexturedStreamLines.cpp

#include <ocean/GLvish/VRenderObject.hpp>
#include <fish/lakeview/bone/FishSlice.hpp>
#include <fish/lakeview/bone/FishField.hpp>

#include <ocean/GLvish/Colormap.hpp> // use of color_t

#include <eagle/FixedArray.hpp>
#include <field/DirectProductArray.hpp>
#include <ocean/GLvish/BoundingBox.hpp>

#include <fish/fiber/vector/Interpolate.hpp>
#include <fish/fiber/vector/LinearIpol.hpp>
#include <fish/fiber/vector/CubicIpol.hpp>

#include <ocean/plankton/VOperatorCache.hpp>



using namespace Wizt;
using namespace Fiber;

class SimpleStreamline : public virtual VRenderObject,  public virtual Fish<Slice>, public virtual Fish<Field>
{
public:
        typedef MemArray<3, tvector>           VectorArray_t; 
        typedef MemArray<3, point>             CoordsArray_t; 
        typedef DirectProductMemArray<point>   ProceduralArray_t; 

        typedef FixedArray<float,3> rgb; 

        // define data type for storing one line
        typedef std::vector<point> Line_t;

        // prepare storing and caching of computed lines in an std::vector of lines
        typedef VOperatorCache<std::vector<Line_t> > OpCache_t;

        TypedSlot<int> LineLength;
        TypedSlot<int> LineWidth;
        TypedSlot<int> PointSize;
        TypedSlot<int> LineNumber;
        TypedSlot<double> StepSize;
        TypedSlot<double> FanDelta;

        TypedSlot<double> X,
                          Y,
                          Z;

        struct  MyState : State
        {};

        override RefPtr<State> newState() const
        {
                return new MyState();
        }

        SimpleStreamline(const string&name, int p, const RefPtr<VCreationPreferences>&VP)
        : VRenderObject(name, p, VP)
        , Fish<VObject>(this)
        , Fish<Field>("datafield")
        , LineLength(this, "LineLength", 40)
        , LineWidth(this, "LineWidth", 2)
        , PointSize(this, "PointSize", 2)
        , LineNumber(this, "LineNumber", 10)
        , StepSize(this, "StepSize", 0.1)
        , FanDelta(this, "FanDelta", 0.05)
        , X(this, "X", 0.271)
        , Y(this, "Y", 0.294)
        , Z(this, "Z", 0.531)
        {
                LineLength.setProperty("max", 100); 
                LineLength.setProperty("min", 1); 
                LineWidth.setProperty("max", 20); 
                PointSize.setProperty("max", 20); 
                LineNumber.setProperty("min", 1);
                X.setProperty("max", 1.0);
                Y.setProperty("max", 1.0);
                Z.setProperty("max", 1.0); 

                acceptType<Eagle::tvector3>();
        }


        override bool update( VRequest&R, double precision );
        override void render( VRenderContext&R ) const;

        void getCoordinates(const point& float_index, point& position_cordinates);
};


bool SimpleStreamline::update( VRequest&R, double precision )
{
        printf("SimpleStreamline::update(): updating\n"); 

FieldSelector FS = getFieldSelector(R);

// get the grid of the fiber bundle 
RefPtr<Grid>  G = FS.GridSource();
        if (!G)
                return true;

// get the coordinates of the grid
RefPtr<Field> Coords  = G->CartesianPositions();
        if (!Coords)  { puts("SimpleStreamline::update(): No coord field!"); return true; } 


// get the vector data of the grid
RefPtr<Field> Vectors = FS.getField();
        if (!Vectors) { puts("SimpleStreamline::update(): No vector field!"); return true; } 


// associate the caching object to the field. As long as the field is valid the 
// cached lines can by used. If a new field is load (e.g. from another time step) 
// OC is a reference pointer to a std::vector of lines, where they are going to be stored.
RefPtr<OpCache_t>&OC = OpCache_t::retrieve( *Vectors ) [ self() ]; 
        if (!OC) OC = new OpCache_t(); 


RefPtr<VectorArray_t> Vctrs = Vectors->getData(); 
        if( !Vctrs)
        {
                printf("SimpleStreamline::update(): Could get vector field data yet!\n"); 
                if (OC)
                {
                        OC->resize( 0 );
                        GLCache::clear( *OC );
                }
                return true;
        } 
        else
                printf("SimpleStreamline::update(): Gotvector field data!\n"); 

        RefPtr<ProceduralArray_t>   PCrds = Coords->getData(); 
        RefPtr<CoordsArray_t>       Crds  = Coords->getData(); 
        if (!(Crds || PCrds))
                {
                        printf("SimpleStreamline::update(): Did not find any coordinates :  %s\n",
                               Typename( Coords->getType() ).c_str() ); 
                        return true;
                } 



RefPtr<ValueSet> Changeables = new ValueSet(); 
        // Register parameters that shall invoke an update call, if changed.
        Changeables & R += X, Y, Z, LineLength, PointSize, StepSize, FanDelta, LineNumber; 

        // if parameters didn't change, don't do a new calculation
        if (OC->unchanged( Changeables ) ) 
        {
                printf("SimpleStreamline::update(): Line Calc Skipped, Lines in Cache...\n");
                return true;
        } 


int line_number = 1; 
        LineNumber << R >> line_number; 

// resize std::vector to number of lines
        OC->resize( line_number ); 
        GLCache::clear( *OC );

// Now do the Interpolation of the Coordinates 
// Get index Range of Data (and Coordinates)
        MultiIndex<3> FieldSize = Vctrs->Size(); 


double x, y, z; 
point float_index; 
point max_float_index; 
point vertex; 
int line_length; 
double hh; 
double fan_delta;

        LineLength << R >> line_length; 
        StepSize << R >> hh; 
        FanDelta<< R >> fan_delta;

// Get user specified position ranging from 0 to 1.0 in each coordinate covering the whole dataset 
        X << R >> x; 
        Y << R >> y; 
        Z << R >> z; 

        max_float_index = FieldSize[0] - 1, FieldSize[1] - 1, FieldSize[2] - 1;
        float_index = x * FieldSize[0] - 1, y * FieldSize[1] -1 , z * FieldSize[2] - 1; 

// Clamp wrong inputs
        for (int i = 0; i < 3; i++)
        {
                if( float_index[i] < 0 )
                        float_index[i] = 0.0; 

                if( float_index[i] >= FieldSize[i] )
                        float_index[i] = FieldSize[i] - 1;
        } 

// assume uniform PCrds here (should be extended to other types of coordinates (curvlinear, ...) as well ) 
// get maximum and minimum coordinates
point DomainStart = PCrds->first(); 
point DomaintEnd = PCrds->last(); 

tvector DomainDiagonal = DomaintEnd - DomainStart; 


// Calc Positions of the StartPoints of the fan of streamlines
std::vector<point> StartVertex; 
        StartVertex.resize(line_number); 

        // First calculate the start position of the first of the start points 
        // This is done using DomainDiagonal as a basis and componentwise array operations, so 
        // it can be optimized using the vector SMID unit of the CPU.
        StartVertex[0] = ( DomainStart + tvector(component_wise( component_wise(DomainDiagonal, max_float_index, Eagle::Operator<'/'>() ), float_index, Eagle::Operator<'*'>() ) ) ); 


tvector CurrentDir; 
point CurrentVertex; 
point CurrentFloatIndex; 
tvector FanOffset; 
        FanOffset = 1, 0, 0;

        resetBBox( R ); 
        embrace( R, StartVertex[0]); 

        // Calculate stream lines
        for(int l = 0; l < line_number; l++)
        {
                // Calucalte startpoints of the fan  by adding the offest to the x-coordinate
                StartVertex[l] = StartVertex[0] + FanOffset * l * fan_delta;
                
        Line_t&myLine = (*OC)[l]; 
                myLine.clear(); 
                
                for(int i = 0; i < line_length; i++)
                {
                        if (i == 0)
                        {
                                CurrentVertex = StartVertex[l]; 

                                // We have to use the float_index space for the interpolator. We already know 
                                // it for the startpoint.
                                Interpolate<3, tvector, CubicIpol <tvector> > VecField(*Vctrs, float_index); 
                                // Get direction at current point
                                CurrentDir = VecField.eval(); 
                        }
                        else
                        {
                                // Calculate the floatindex of the current position from data coordinates
                                CurrentFloatIndex = point( tvector( component_wise(component_wise(CurrentVertex - DomainStart, DomainDiagonal, Eagle::Operator<'/'>() ), max_float_index, Eagle::Operator<'*'>()) )); 

                                // Stop line calculation, when going out of the data range
                                if(CurrentFloatIndex [0] < 0.0 ||CurrentFloatIndex [0] > max_float_index[0])
                                        break; 
                                if(CurrentFloatIndex [1] < 0.0 ||CurrentFloatIndex [1] > max_float_index[1])
                                        break; 
                                if(CurrentFloatIndex [2] < 0.0 ||CurrentFloatIndex [2] > max_float_index[2])
                                        break; 

                                // Get direction at current point.
                                Interpolate<3, tvector, CubicIpol <tvector> > VecField(*Vctrs, CurrentFloatIndex); 
                                CurrentDir = VecField.eval(); 
                        } 

                        // calc new point
                        CurrentVertex += CurrentDir * hh; 
                        myLine.push_back(CurrentVertex); 
                        embrace( R, CurrentVertex); 
                }
        } 

        printf("SimpleStreamline::update(): Integrated into Cache\n");
        closeBBox( R ); 
        
        return true;
}


void SimpleStreamline::render( VRenderContext& Context ) const
{

// Get handle to vectorfield 
 RefPtr<Field> Vectors = getField( Context ); 
        if (!Vectors) 
        {
                puts("void SimpleStreamline::render(VRenderContext&RenderContext): NO VECTORS");
                return;
        } 

// Get cached lines
const RefPtr<OpCache_t>&OC = OpCache_t::retrieve( *Vectors ) ( self() );
        if (!OC)
        {
                puts("void SimpleStreamline::render(VRenderContext&RenderContext): NO LINE CACHE"); 
                return;
        } 
        if (OC->size()<1)
        {
                puts("void SimpleStreamline::render(VRenderContext&RenderContext): NO LINES in CACHE"); 
                return;
        } 

int line_width, point_size; 
int line_length  = 1; 

        LineWidth  << Context >> line_width; 
        PointSize  << Context >> point_size; 
        LineLength << Context >> line_length;

        if(line_width < 1)
                line_width = 1; 

        glEnable(GL_DEPTH_TEST); 
        glDisable(GL_LIGHTING); 

        glLineWidth( line_width ); 

        if(point_size > 0)
                glPointSize( point_size ); 
        else
                glPointSize( 1 ); 

        RefPtr<MyState>  state = myState(Context); 

        RefPtr<ValueSet> RenderParameterSpace = new ValueSet(); 
        RefPtr<DisplayList> DL;
                try{
                        DL = Context(*OC)(this)( RenderParameterSpace ); 
                } 
                catch(...){}
        RefPtr<ValueSet> Changeables = new ValueSet();
Changeables & Context += LineLength, LineNumber, StepSize; 

        if ( DL && OC->unchanged( Changeables  ))   //call display list if availlable, check if changables have changed         
        {
                printf("SimpleStreamline::render(): Calling List\n"); 
                DL->call(); 
        } 
        else                                        //prepare display list
        {
                printf("SimpleStreamline::render(): Rendering List\n");
 
                DL = Context[*OC][this][ RenderParameterSpace ]; 
                if (!DL) throw "SimpleStreamline::render() No Display List!?"; 

                glCompile( *DL ) 
                {
                int l = 0; 
                int p = 0; 

                        // iterate over all lines in cache and draw a line_strip
                        for( OpCache_t::const_iterator lbit = OC->begin(); lbit != OC->end(); lbit++ )
                        {
                                p = 0; 
                                //printf("SimpleStreamline::render()> :Line %d\n", l);
                        Line_t::const_iterator lit = lbit->begin(); 

                                if (lit == lbit->end() )
                                        break; 

                                glBegin( GL_LINE_STRIP );    // GLCHECK must not be used in between a glBegin and glEnd statement
                        bivector dir;

                                for(; lit != lbit->end();)
                                {
//                                      printf("SimpleStreamline::render()> :point %d\n", p);
                                point   P = *lit; 
                                        lit++; 
                                        if (lit != lbit->end()) 
                                                dir = (*lit - P).unit(); 

                                        // use length of the line as color information
                                        glColor3f((double)p/lbit->size() , 1-(double)p/lbit->size(), 0.0);
                                        glVertex( P ); 

                                        p++;
                                } 

                                l++; 
                                glEnd(); 
                        } 


                        // draw points 
                        if(point_size > 0 )
                        {
                                for( OpCache_t::const_iterator lbit = OC->begin(); lbit != OC->end(); lbit++ ) 
                                {
                                        p = 0; 
                                        //printf("SimpleStreamline::render()> :Line %d\n", l); 
                                        Line_t::const_iterator lit = lbit->begin(); 

                                        if (lit == lbit->end() ) 
                                                break; 

                                tvector dir; 
                                        glBegin(GL_POINTS); 

                                        for(; lit != lbit->end();)
                                        {
                                                //printf("SimpleStreamline::render()> :point %d\n", p);
                                        point   P = *lit; 
                                                lit++; 
                                                if (lit != lbit->end())
                                                        dir = (*lit - P).unit(); 

                                                glColor3f((double)p/lbit->size(), 1-(double)p/lbit->size(), 0.7); 
                                                glVertex( P ); 

                                                p++;
                                        } 
                                        glEnd(); 
                                } 
                        }
                }
        } 

        state->touch(); // synchinozes state, prevents unneccessary calls of the update function.
}


static Ref<VCreator<SimpleStreamline, AcceptList<Fiber::Field> >  > myCreator("Tutorial/SimpleStreamline");

VISH_DEFAULT_INIT

---