C:/Projekte/C++/FlowVIS_107/src/StreamLine.cpp

Go to the documentation of this file.

#include "StreamLine.h"
#include "Texture.h"
#include "Shader.h"

// needed for random-feed
#include <ctime>

// needed for writing log-file
#include <iostream>
#include <fstream>


StreamLine::StreamLine( FlowData *fd )
{
  //cout << "New streamline..." << endl;

  m_flowData = fd;
  m_flowGeometry = m_flowData->getGeometry();

  m_dSep  = 0.05f;
  m_dTest = m_dSep * 0.5f;

  reset();
}


StreamLine::~StreamLine()
{
  cout << "Delete streamline..." << endl;
}


void StreamLine::create()
{
  //cout << "Create streamlines..." << endl;

  /*
    - create sl-grid
    - find random start-point

    repeat
      - create streamline
      - create seed-points
      - take next seed-point, remove it from queue and check if (still) valid
        - if valid -> make new streamline
          else -> take next seed-point
    until no more seed-points

  */

  //cout << "\n\n -- create SL -- " << endl;

  // reset streamlines
  reset();

  // set global parameters
  m_d = m_dTest * 0.8f;
  m_maxSteps = 1000;
  m_maxLines = 5000;
  m_maxLineThickness = 0.005f;

  if ( m_dSep==0 ) m_dSep = 0.001f;

  cout << "dsep: " << m_dSep << ", dTest: " << m_dTest << ", d: " << m_d << "\n";

  // just for testing
  //setTestParameters();

  // create the initial streamline
  m_slgGrid = new StreamLineGrid( m_dSep );
  createInitialStreamLine();
  createSeedPoints();

  int nrStreamLines = 1;

  // while there are seed-points, create more streamlines
  while ( ( nrStreamLines < m_maxLines ) && (!m_seedPoints.empty()) )
  {
    vec3 startPos = m_seedPoints.at(0);
    m_seedPoints.erase( m_seedPoints.begin() );

    // check if valid seedpoint
    if ( isValidSeedPoint( startPos ) )
    {
      createStreamLine( startPos );
      createSeedPoints();
      ++nrStreamLines;
    }
  }

  // TEST
  //m_slgGrid->printGrid();

  //cout << "\n -- finished SL -- \n" << endl;

}


void StreamLine::create( float dSep, float dTest, bool linesActive, bool taperingActive, bool glyphActive )
{
  // set global parameters
  m_dSep  = dSep;
  m_dTest = dTest;

  m_isLinesActive    = linesActive;
  m_isTaperingActive = taperingActive;
  m_isGlyphActive    = glyphActive;

  // create the streamline
  create();
}


void StreamLine::reset()
{
  cout << "Reseting streamlines...\n";

  // delete all streamlines
  for ( int i=0; i<(int)m_streamLines.size(); i++ )
    m_streamLines.at(i).streamLine.clear();

  m_streamLines.clear();
  //delete m_slgGrid;

  // delete glyph-texture
  if ( m_glyphTexId > 0 )
    glDeleteTextures( 1, &m_glyphTexId );


  // set default values
  float color[] = { 1.0f, 1.0f, 1.0f, 1.0f };
  setLineColor( color );
  setGlyphColor( color );

  setGlyphSize( 100 );
  m_glyphDistance = 5;

  //m_dSep  = 0.05f;
  //m_dTest = m_dSep * 0.5f;

  m_isLinesActive    = true;
  m_isTaperingActive = false;
  m_isGlyphActive    = false;

  m_glyphTexId = 0;
  m_integrationMethod = EULER;
}


void StreamLine::draw()
{
  cout << "Draw streamlines..." << endl;

  glPushMatrix();

  // draw streamlines
  if ( m_isLinesActive )
  {
      drawStreamLines();
  }

  // draw glyphs
  if ( m_isGlyphActive )
  {
    bool haveGlyphTexture = false;

    if ( m_glyphTexId == 0 ) 
      // no texture created yet -> create it
      haveGlyphTexture = createGlyphTexture();
    else
      haveGlyphTexture = true;

    if ( haveGlyphTexture ) drawGlyphs();
  }

/*
  // TEST
  glTranslatef( -0.5f, 0.0f, 0.0f );
  glRotatef( 180.0f, 0.0f, 1.0f, 0.0f );
  glTranslatef(  0.5f, 0.0f, 0.0f );
*/
  
  glPopMatrix();

}


void StreamLine::drawStreamLines() 
{
  // save old color
  float *oldColor = new float[4];
  glGetFloatv( GL_CURRENT_COLOR, oldColor );

  // enable blending
  glEnable( GL_BLEND );
  glBlendFunc( GL_SRC_ALPHA, GL_ONE );

  // set SL-color
  glColor4fv( m_lineColor );

  // with or without tapering
  if ( m_isTaperingActive )
    drawStreamLinesWithTapering();
  else
    drawStreamLinesWithoutTapering();

  glDisable( GL_BLEND );

  // restore old color
  glColor4fv( oldColor );
}


void StreamLine::drawStreamLinesWithTapering()
{
  // draw all streamlines with tapering
  for ( int j=0; j<(int)m_streamLines.size(); j++ )
  {
    glBegin( GL_QUAD_STRIP );
    for ( int i=0; i<(int)m_streamLines.at(j).streamLine.size(); i++ )
    {
      glVertex3f( m_streamLines.at(j).streamLine.at( i ).point[X], 
                  m_streamLines.at(j).streamLine.at( i ).point[Y],
                  0.0f );

      glVertex3f( m_streamLines.at(j).streamLine.at( i ).tPoint[X], 
                  m_streamLines.at(j).streamLine.at( i ).tPoint[Y],
                  0.0f );
    } // for all points of streamline
    glEnd();
  } // for

/*
  // TEST ---
  for ( int j=0; j<(int)m_streamLines.size(); j++ )
  {
    //TEST
    //float r =  (float)((float)j / (float)m_streamLines.size());
    //glColor4f( r, 0.0f, 0.0f, 0.0f );

    cout << " pts: " << m_streamLines.at(j).streamLine.size() << "\n";

    glBegin( GL_POINTS );
    for ( int i=0; i<(int)m_streamLines.at(j).streamLine.size(); i++ )
    {
      //TEST
      int sz = (int)m_streamLines.at(j).streamLine.size();
      float g =  (float) ( (float)i / (float)sz );
      //cout << "g: " << g << ", sz: " << sz <<"\n";
      glColor4f( 1.0f, g, 0.0f, 0.0f );

      glVertex3f( m_streamLines.at(j).streamLine.at( i ).point[X], 
                  m_streamLines.at(j).streamLine.at( i ).point[Y],
                  0.0f );
/
      glVertex3f( m_streamLines.at(j).streamLine.at( i ).tPoint[X], 
                  m_streamLines.at(j).streamLine.at( i ).tPoint[Y],
                  0.0f );*
    } // for all points of streamline
    glEnd();
  } // for
  // end TEST
*/

}


void StreamLine::drawStreamLinesWithoutTapering()
{
  // draw all streamlines
  for ( int j=0; j<(int)m_streamLines.size(); j++ )
  {
    glBegin( GL_LINE_STRIP );
    for ( int i=0; i<(int)m_streamLines.at(j).streamLine.size(); i++ )
    {
      glVertex3f( m_streamLines.at(j).streamLine.at( i ).point[X], 
                  m_streamLines.at(j).streamLine.at( i ).point[Y],
                  0.0f );
    } // for all points of streamline
    glEnd();
  } // for 
}


void StreamLine::drawGlyphs() 
{
  // NOTE: don't expect any nice code here... very "hacky" :(

  cout << "drawGlyphs...\n";

  glPushMatrix();

  // enable blending
  glEnable( GL_BLEND );
  glBlendFunc( GL_SRC_ALPHA, GL_ONE );

  // save old color
  float *oldColor = new float[4];
  glGetFloatv( GL_CURRENT_COLOR, oldColor );

  // set color
  glColor4fv( m_glyphColor );

  // configure texture-unit
  glActiveTexture( GL_TEXTURE0 );
  glEnable( GL_TEXTURE_2D );
  glEnable( GL_POINT_SPRITE_ARB );

  // prepare texture
  glBindTexture( GL_TEXTURE_2D, m_glyphTexId );
  glTexEnvf( GL_POINT_SPRITE_ARB, GL_COORD_REPLACE_ARB, GL_TRUE );
  glTexEnvi( GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_NONE );


  // initialise some parameters
  glPointSize( m_glyphSize );

  if ( m_glyphDistance == 0 ) m_glyphDistance = 1;

  int n = m_glyphDistance; // display glyphs for every n'th point of the streamline
  int m = m_glyphDistance; // display glyphs for every m'th streamline

  vec3 pos( -1.0f, -1.0f, -1.0f );
  float angle = 0;


  // --- point-sprite-shader --- 

  Shader *shader = new Shader();
  shader->Load( "sprite.vert", "sprite.frag" );

  GLint spriteTexLoc = glGetUniformLocation( shader->GetProgram(), "SpriteTexture" );
  GLint angleLoc     = glGetUniformLocation( shader->GetProgram(), "angle" );

  int nrSls = (int)m_streamLines.size(); // number of streamlines

  // for all streamlines
  for ( int j=0; j<nrSls; j=j+m )
  {
    int nrPts = (int)m_streamLines.at(j).streamLine.size();  // number of points per streamline

    //cout << "max glyphs: " << nrPts << " / sls: " << nrSls << "\n";

    // print glyphs
    for ( int i=0; i<nrPts; i=i+n )
    {
      pos = m_streamLines.at( j ).streamLine.at( i ).point;
      angle = atan2( m_streamLines.at(j).streamLine.at(i).direction[Y], m_streamLines.at(j).streamLine.at(i).direction[X] );

      //cout << "glyph nr: " << i << ", angle: " << angle << "\n";

      shader->Bind();  

      glUniform1i( spriteTexLoc, 0 );
      glUniform1f( angleLoc, angle ); 

      glBegin( GL_POINTS );
        glVertex3fv( pos.v );
      glEnd();

      shader->Unbind();
    } // for all points in one streamline
  } // for all streamlines

  // --- ( end shader ) --- 


  // disable gl-properties
  glDisable( GL_POINT_SPRITE_ARB );
  glDisable( GL_TEXTURE_2D );
  glDisable( GL_BLEND );

  glPopMatrix();

  // restore old color
  glColor4fv( oldColor );
}


void StreamLine::createInitialStreamLine()
{
  // create random start-point
  vec3 pos( getRandomFloat( 654321 ), getRandomFloat( 123456 ), 0.0f );
  cout << "Initial streamline starting at: " << pos << endl;

  createStreamLine( pos );

  // TEST
  //createStreamLine( vec3( 0.478145f, 0.35, 0.0f ) );
  //createStreamLine( vec3( 0.478145f, 0.312874f, 0.0f ) ); // -> null-vel-vec!
}


void StreamLine::createStreamLine( vec3 pos )
{
  // create forward-integration and backward-integration
  vector<slSamplePoint> fwdInt = doIntegration( pos, FORWARD );
  vector<slSamplePoint> bwdInt = doIntegration( pos, BACKWARD );

  // merge fwdInt- and bwdInt-points to one streamline
  vector<slSamplePoint> allSamplePoints;
  allSamplePoints.clear();

  // add inverted bwdInt
  int nrBwdPts = (int)bwdInt.size();
  for ( int i=nrBwdPts-1; i>=0; i-- )
    allSamplePoints.push_back( bwdInt.at(i) );

  // add fwdInt-points
  int nrFwdPts = (int)fwdInt.size();
  for ( int i=0; i<nrFwdPts; i++ )
    // i=1 -> don't add startPoint twice!
    allSamplePoints.push_back( fwdInt.at(i) );

/*  if ( allSamplePoints.empty() ) 
  {
    cout << "don't add empty streamline\n";
  }
  else
  {*/
    // add this streamline to the streamlines-list
    slStreamLine sl;
    sl.streamLine = allSamplePoints;

    m_streamLines.push_back( sl );
  //}
}


vector<slSamplePoint> StreamLine::doIntegration( vec3 pos, bool integration )
{
  //cout << "New Streamline starting at: " << pos << endl;

  // do some initialisation
  int i=0;
  bool stopLine = false;
  vec3 lastSamplePoint(-1, -1, -1);
  vector<slSamplePoint> samplePoints;


  // create samplepoints
  while ( i<m_maxSteps && !stopLine )
  {
    //cout << " ---\n";
    vec3 velocity[1];

    // get interpolated velocity-vector at given position
    if ( !m_flowGeometry->getInterpolationAt( pos, velocity, m_flowData ) )
    {
      // we seem to be out of the image or velocity-vector == (0, 0, 0) -> stop streamline
      //cout << "Stop SL - interpolation out of range!" << endl;
      //cout << "stop at i: " << i <<  "\n";
      stopLine = true;
    }
    else
    {
      if ( velocity[0][X] == 0.0f && velocity[0][Y] == 0.0f ) 
      {
        //cout << "Null-Vector!!\n";
        stopLine = true;
      }
      else {

        // if backward-integration, change direction
        if ( integration == BACKWARD )
          velocity[0] *= -1.0f;

        // normalize vector
        velocity[0].normalize();

        // calculate new position
        if ( m_integrationMethod == EULER )
        {
          pos = calculateEulerPosition( pos, velocity[0] );
        }
        else // RUNGE_KUTTA2
        {
          pos = calculateRungeKutta2Position( pos, velocity[0], integration );
        }
        //cout << "newPos: " << pos << " at iteration " << i << endl;


        // check if still in range
        if ( pos[X] < 0.0f || pos[X] > 1.0f ||
             pos[Y] < 0.0f || pos[Y] > 1.0f )
        {
          //cout << "Stop SL - new position out of range at sample-point " << i << "!" << endl;
          stopLine = true;

          // TODO: Move point to border, so that line doesn't stop before image-border
          // Primitive variant: just clamp outside-value to border value (we make a little error than, but quick to implement)

        }
        else 
        {
          // still in range...

          // check distance
          float minDist = m_slgGrid->getMinimalDistance( pos, lastSamplePoint );
          //cout << "minDist: " << minDist << "\n";

          if ( (minDist != -1) && (minDist < m_dTest) ) 
          {
            // there are points too close to this one -> stop line
            //cout << "Stop SL - other points too close to sample-point " << i << "!" << endl;
            stopLine = true;
          }
          else
          {
            // valid sample-point...

            slSamplePoint sp;
            sp.point = pos;

            // calculate tapering
            if ( true /*m_isTaperingActive*/ ) // -> always do this, so we don't have to recalculate everytime on (un-)checking
            {
              // calculate a 2nd point for tapering

              // rotate velocity-vector 90°
              vec3 tDist = rotateVector( velocity[0], PI/2 );

              // calculate line-thickness
              float currentThickness = m_maxLineThickness * calculateThicknessCoeff( minDist );
            
              // add offset and save point
              sp.tPoint = pos + ( tDist * currentThickness );
            }

            // store direction for glyphs
            sp.direction = velocity[0];

            // if backward-integration, change direction
            if ( integration == BACKWARD )
              sp.direction *= -1.0f;

            // finally add points
            samplePoints.push_back( sp );
            m_slgGrid->addPoint( pos );

            lastSamplePoint = pos;
            ++i;
          }
        } 
      }
    }
  } // while

  //cout << "finished a line with i: " << i << "\n";
  return samplePoints;
}


void StreamLine::createSeedPoints()
{
  int currentStreamLine = (int)m_streamLines.size()-1; 
  int nrPts             = (int)m_streamLines.at( currentStreamLine ).streamLine.size();


  // for every sample-point, generate two seed-points
  for ( int i=0; i<nrPts; i++ )
  {
    vec3 samplePoint = m_streamLines.at( currentStreamLine ).streamLine.at( i ).point;
    vec3 direction[1];

    // get the direction, which is the velocity-vector, to calculate the normals
    if ( !m_flowGeometry->getInterpolationAt( samplePoint, direction, m_flowData ) )
      // continue with next sample-point on error (should never happen anyway, but u never know)... 
      // dirty, but works :)
      continue;

    //cout << "samplepoint: " << samplePoint << ", direction: " << direction[0] << "\n";

    // rotate vel vec 90° and 270°
    vec3 normal = rotateVector( direction[0], PI/2 );
    vec3 seedPoint = samplePoint + (normal * m_dSep);
    
    //cout << "new seedpoint: " << seedPoint << "\n";
    m_seedPoints.push_back( seedPoint );

    normal = rotateVector( direction[0], 3*PI/2 );
    seedPoint = samplePoint + (normal * m_dSep);

    //cout << "new seedpoint: " << seedPoint << "\n";
    m_seedPoints.push_back( seedPoint );
  }

  //cout << "done with creating seed-points.\n";
}


bool StreamLine::isValidSeedPoint( vec3 point )
{
  float minDist = m_slgGrid->getMinimalDistance( point, vec3(-1, -1, -1) );

  if ( minDist == -1 ) return true;
  if ( minDist < m_dSep ) return false;
  return true;
}


vec3 StreamLine::calculateEulerPosition( vec3 pos, vec3 vel )
{
  return pos + ( vel * m_d );
}


vec3 StreamLine::calculateRungeKutta2Position( vec3 pos1, vec3 vel1, bool integration )
{
  vec3 pos = pos1;
  vec3 pos2 = calculateEulerPosition( pos1, vel1 );
  vec3 vel2[1];

  // get interpolated velocity-vector at position 2
  if ( !m_flowGeometry->getInterpolationAt( pos2, vel2, m_flowData ) )
  {
    // we seem to be out of the image -> just use euler at the border
    pos = calculateEulerPosition( pos1, vel1 );
  }
  else
  {
    // if backward-integration, change direction
    if ( integration == BACKWARD )
      vel2[0] *= -1.0f;

    pos = pos1 + ( vel1 + vel2[0] ) * ( m_d * 0.5f );
  }
  
  return pos;
}


float StreamLine::calculateThicknessCoeff( float distance )
{
  float thicknessCoeff = 1.0f;

  if ( (distance != -1) && (distance < m_dSep) )
  {
    thicknessCoeff = ( distance - m_dTest ) / ( m_dSep - m_dTest );
  }

  return thicknessCoeff;
}


bool StreamLine::createGlyphTexture()
{
  cout << "Create glyph-texture...\n";

  Texture tex;
  m_glyphTexId = tex.CreateGlyph();

  if ( m_glyphTexId != 0 ) return true;

  return false;
}


float StreamLine::getRandomFloat( int rnd )
{
  srand( (unsigned)time(0) * rnd );
  return rand() / ( float(RAND_MAX) + 1 );
}


float StreamLine::toRad( float grad )
{
  return (float)(2*PI* grad /360 );
}


vec3 StreamLine::rotateVector( vec3 vector, float phi )
{
  vec3 rotVec( 0.0f, 0.0f, 0.0f );

  rotVec[X] = vector[X] * cos( phi ) - vector[Y] * sin( phi );
  rotVec[Y] = vector[X] * sin( phi ) + vector[Y] * cos( phi );

  rotVec.normalize();
  return rotVec;
}


void StreamLine::setTestParameters()
{
  // TEST

  m_dSep  = 0.01f;
  m_dTest = m_dSep  * 0.5f;
  m_d     = m_dTest * 1.0f;

  //m_isTaperingActive = true;
  //m_isGlyphActive = true;
  //m_integrationMethod = EULER; // RUNGE_KUTTA2; // EULER; 

/*
  // set global parameters
  m_maxSteps = 10000;
  m_maxLines = 10000;
  m_maxLineThickness = 0.001f;
*/
  /*
  float c[] = { 1.0f, 0.0f, 0.0f, 0.0f };
  setLineColor( c );
  setGlyphColor( c );

*/
}

---