streamlines.cpp

Go to the documentation of this file.

#include "streamlines.h"
#include <cmath>
#include <algorithm>

float StreamlineGenerator::GetMinDistance(const vec2 &v) const
{
    typedef std::list<vec2>::const_iterator iter_t;

    int cx, cy;
    ComputeCellPosition(v, cx, cy);

    // compute squared minimum distance to other points.
    // only consider points in the 3x3 grid surrounding the point.
    float mindist2 = HUGE;

    const iter_t beg = grid(cx, cy).begin();
    const iter_t end = grid(cx, cy).end();

    for (iter_t it = beg; it != end; ++it) {
        const vec2 diff = v - *it;
        mindist2 = std::min(mindist2, dot(diff, diff));
    }

    const int y0 = std::max(0, cy - 1);
    const int y1 = std::min(grid.GetSizeY() - 1, cy + 1);
    const int x0 = std::max(0, cx - 1);
    const int x1 = std::min(grid.GetSizeX() - 1, cx + 1);

    for (int y = y0; y <= y1; ++y) {
        for (int x = x0; x <= x1; ++x) {
            if (x == cx && y == cy) continue;

            const iter_t beg = grid(x, y).begin();
            const iter_t end = grid(x, y).end();

            for (iter_t it = beg; it != end; ++it) {
                const vec2 diff = v - *it;
                mindist2 = std::min(mindist2, dot(diff, diff));
            }
        }
    }

    return std::sqrt(mindist2);
}

void StreamlineGenerator::PutIntoGrid(const vec2 &v)
{
    int cx, cy;
    ComputeCellPosition(v, cx, cy);
    grid(cx, cy).push_back(v);
}

void StreamlineGenerator::DelayedPutIntoGrid(const vec2 &v)
{
    point_queue.push_back(v);

    if (point_queue.size() > grid_add_delay) {
        const vec2 &p = point_queue.front();
        PutIntoGrid(p);
        point_queue.pop_front();
    }
}

void StreamlineGenerator::BeginTraceLine() {}

void StreamlineGenerator::EndTraceLine()
{
    while (point_queue.size()) {
        const vec2 &p = point_queue.back();
        PutIntoGrid(p);
        point_queue.pop_back();
    }
}

bool StreamlineGenerator::InterpolateTo(vec2 &position, bool forward)
{
    vec2 steppos;
    float direction = (forward)?(1.0f):(-1.0f);

    switch(inter_type)
    {
    case ITER_RK:
        steppos = position + (v(position) * direction) * stepsize_half;
        if (IsOutsideBounds(steppos)) return false;
        else position = position + (v(steppos) *direction) * stepsize;
        break;
    case ITER_EULER:
        position = position + (v(position) * direction) * stepsize;
        break;
    }

    return true;
}

void StreamlineGenerator::TraceLine(const vec2& startpos, Streamline &sline)
{
    Streamline tmpline;

    BeginTraceLine();

    bool trace_fwd = true;
    bool trace_back = true;

    vec2 fpos = startpos;
    vec2 bpos = startpos;

    // termination criteria
    if (IsOutsideBounds(fpos)) return;

    const float mindist = GetMinDistance(fpos);
    if (mindist < dtest) return;

    if(!InterpolateTo(fpos, true)) trace_fwd = false;

    for (int i = 0; i < maxiter / 2; ++i) {
        if (trace_fwd) {
            if (IsOutsideBounds(fpos)) goto tf_break;

            const float fmindist = GetMinDistance(fpos);
            if (fmindist < dtest) goto tf_break;

            const vec2 v_fpos = v(fpos);
            tmpline.push_back(SamplePoint(fpos, v_fpos, ComputeThickness(fmindist)));

            DelayedPutIntoGrid(fpos);

            if(!InterpolateTo(fpos, true)) goto tf_break;
        }

        goto tf_end;
        tf_break: trace_fwd = false;
        tf_end: ;

        if (trace_back) {
            if (IsOutsideBounds(bpos)) goto tb_break;

            const float bmindist = GetMinDistance(bpos);
            if (bmindist < dtest) goto tb_break;

            const vec2 v_bpos = v(bpos);
            sline.push_back(SamplePoint(bpos, v_bpos, ComputeThickness(bmindist)));

            DelayedPutIntoGrid(bpos);

            if(!InterpolateTo(bpos, false)) goto tb_break;
        }

        goto tb_end;
        tb_break: trace_back = false;
        tb_end: ;
    }


    // reverse the order of the points in the line which has been traces backwards
    std::reverse(sline.begin(), sline.end());

    // Concatenate both line traces
    sline.insert(sline.end(), tmpline.begin(), tmpline.end());

    EndTraceLine();
}

void StreamlineGenerator::TraceAndStore(const vec2 &v, Streamline &sl)
{
    sl.clear();
    TraceLine(v, sl);

    // when the candidate point was invalid the line is empty
    if (sl.size()) streamlines.push_back(sl);
}

StreamlineGenerator::StreamlineGenerator(
    const SL_Rectangle &domain,
    float dsep,
    float dtest,
    float stepsize,
    int maxiter,
    VelocityFunc v,
    SLIterAlg inter_type
)
    : domain(domain)
    , dsep(dsep)
    , dtest(dtest)
    , stepsize(stepsize)
    , stepsize_half(stepsize / 2.0f)
    , maxiter(maxiter)
    , v(v)
    , grid_add_delay(std::ceil(2.0f * (0.5f * (dtest + dsep)) / stepsize))
    , grid(std::ceil(domain.GetWidth() / dsep), std::ceil(domain.GetHeight() / dsep))
    , inter_type(inter_type)
{
    //assert(stepsize < dsep);

    // variable used throughout the method
    Streamline sline;

    // trace the initial line
    const float x = (domain.left + domain.right) * 0.5f;
    const float y = (domain.bottom + domain.top) * 0.5f;
    const vec2 startpos(x, y);
    TraceAndStore(startpos, sline);


    for (size_t cli = 0; cli < streamlines.size(); ++cli) {
        for (size_t i = 0, _n = streamlines[cli].size(); i < _n; ++i) {
            const Streamline &line = streamlines[cli];

            // offset the starting point perpendicular to the direction vector
            const vec2 offset = rotate90(line[i].dir) * dsep;
            const vec2 sp1 = line[i].pos + offset;
            const vec2 sp2 = line[i].pos - offset;

            // trace both lines
            TraceAndStore(sp1, sline);
            TraceAndStore(sp2, sline);
        }
    }
}

---