streamlines.h

Go to the documentation of this file.

#ifndef STREAMLINES_H
#define STREAMLINES_H

#include "flowvis_data.h"

#include <vector>
#include <list>
#include <cassert>

struct vec2 {
    float x; 
    float y; 

    vec2() {}
    vec2(float x, float y) : x(x), y(y) {}
};

inline vec2 rotate90(const vec2 &v)
{ return vec2(-v.y, v.x); }

inline const float dot(const vec2 &lhs, const vec2 &rhs)
{ return lhs.x * rhs.x + lhs.y * rhs.y; }

inline const vec2 operator*(const vec2 &v, float f)
{ return vec2(v.x * f, v.y * f); }

inline const vec2 operator+(const vec2 &lhs, const vec2 &rhs)
{ return vec2(lhs.x + rhs.x, lhs.y + rhs.y); }

inline const vec2 operator-(const vec2 &lhs, const vec2 &rhs)
{ return vec2(lhs.x - rhs.x, lhs.y - rhs.y); }

inline const float length(const vec2 &vec)
{ return sqrt(vec.x*vec.x + vec.y*vec.y); }
struct SamplePoint {
    vec2 pos; 
    vec2 dir; 
    float thickness; 
    float rotation; 

    SamplePoint()
        : thickness(0.0f)
    {
    }

    SamplePoint(const vec2 &pos, const vec2 &dir, const float &thickness = 1.0f)
        : pos(pos), dir(dir), thickness(thickness)
    {
        float pi = 3.14159265f;

        vec2 basic(0, 1);

        float radians = acos(dot(basic, dir) / (length(basic) * length(dir)));
        
        if(dir.x < 0)
        {
            radians = (2 * pi) - radians;
        }

        radians += (pi / 4);

        if(radians >= (2 * pi))
            radians -= (2 * pi);

        rotation = radians;
    }
};

struct SL_Rectangle {
    float left, bottom;
    float right, top;

    SL_Rectangle() {}

    SL_Rectangle(float left, float bottom, float right, float top)
        : left(left)
        , bottom(bottom)
        , right(right)
        , top(top)
    {}

    float GetWidth() const
    { return right - left; }

    float GetHeight() const
    { return top - bottom; }
};

namespace detail
{
    template <typename T>
    class Array2d {
    private:

        int sizex;
        int sizey;

        std::vector<T> data;

    public:
        Array2d(int dimx, int dimy)
            : sizex(dimx)
            , sizey(dimy)
            , data(dimx * dimy)
        {}

        const T& operator()(int x, int y) const
        {
            assert(0 <= x && x < sizex);
            assert(0 <= y && y < sizey);
            return data[x + y * sizex];
        }

        T& operator()(int x, int y)
        {
            assert(0 <= x && x < sizex);
            assert(0 <= y && y < sizey);
            return data[x + y * sizex];
        }

        const int GetSizeX() const
        { return sizex; }

        const int GetSizeY() const
        { return sizey; }

        const std::vector<T>& GetData() const
        { return data; }
    };
}

typedef std::vector<SamplePoint> Streamline;

class VelocityFunctor
{
public:
    VelocityFunctor(FlowVisData *data)
        : m_data(data)
    {
    }

    ~VelocityFunctor()
    {
    }

    vec2 operator() (const vec2 &pos)
    {
        vec3 erg = !(m_data->GetVelocityAt(pos.x, pos.y));

        return rotate90(vec2(erg[1], erg[0]));
    }

private:
    FlowVisData *m_data;
};

enum SLIterAlg{
    ITER_RK,

    ITER_EULER
};

class StreamlineGenerator {
public:
    //typedef vec2 (*VelocityFunc)(const vec2 &);
    typedef VelocityFunctor VelocityFunc;

private:
    SL_Rectangle domain;

    float dsep;
    float dtest;

    float stepsize;
    float stepsize_half;

    int maxiter;
    
    SLIterAlg inter_type;

    VelocityFunc v;

    std::vector<Streamline> streamlines;

    size_t grid_add_delay;

    detail::Array2d<std::list<vec2> > grid;

    std::list<vec2> point_queue;

private:

    bool IsOutsideBounds(const vec2 &v) const
    {
        return !(domain.left <= v.x && v.x <= domain.right
            && domain.bottom <= v.y && v.y <= domain.top);
    }

    void ComputeCellPosition(const vec2 &v, int &cx, int &cy) const
    {
        cx = (v.x - domain.left) / dsep;
        cy = (v.y - domain.bottom) / dsep;
    }

    void BeginTraceLine();

    void EndTraceLine();

    float GetMinDistance(const vec2 &v) const;

    void PutIntoGrid(const vec2 &v);

    void DelayedPutIntoGrid(const vec2 &v);

    void TraceLine(const vec2 &startpos, Streamline &sline);

    void TraceAndStore(const vec2 &v, Streamline &sl);

    float ComputeThickness(float mindist)
    { return mindist < dsep ? (mindist - dtest) / (dsep - dtest) : 1.0f; }

    bool InterpolateTo(vec2 &position, bool forward);

public:
    StreamlineGenerator(
        const SL_Rectangle &domain,
        float dsep,
        float dtest,
        float stepsize,
        int maxiter,
        VelocityFunc v,
        SLIterAlg inter_type
    );

    const std::vector<Streamline>& GetStreamlines() const
    { return streamlines; }
};


#endif

---