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

Go to the documentation of this file.

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

#include "FlowGeometry.h"
#include "FlowChannel.h"
#include "FlowData.h"
#include "reverseBytes.h"

using namespace std;



FlowGeometry::FlowGeometry()
{
  geometryData     = NULL;
  normGeometryData = NULL;
}


FlowGeometry::~FlowGeometry()
{
  if ( geometryData )
    delete[] geometryData;

  if ( normGeometryData )
    delete[] normGeometryData;
}


bool FlowGeometry::readFromFile(char* header, FILE* fp, bool bigEndian)
{
  isFlipped = false;

  std::cout << "\n\n ------- Loading grid ------- " << std::endl;

  //determine the dimensions
  sscanf( header, "SN4DB %d %d %d", &dim[0], &dim[1], &dim[2] );
  std::cout << "Dimensions: " << dim[0] << " x " << dim[1] << " x " << dim[2] << std::endl;

  if ( dim[2]!=1 )
  {
    std::cerr << "Invalid Z dimension value." << std::endl;
    return false;
  }

  // create space for grid-data
  geometryData = new vec3[ dim[0] * dim[1] ];

  //read the data in one chunk
  int result = (int) fread( geometryData, sizeof(vec3), dim[0] * dim[1], fp );

  // check if everything went fine
  if ( result != dim[0]*dim[1] )
  {
    std::cerr << "+ Error reading grid file." << std::endl << std::endl;
    return false;
  }

  // we usually do not have big endian
  if ( bigEndian )
    for (int j = 0; j < getDimX()*getDimY(); j++)
      for (int k = 0; k < 3; k++)
        geometryData[j][k] = reverseBytes<float>(geometryData[j][k]);


  // determine boundaries
  boundaryMin = vec3( getPos( 0 ) );                  //first vertex
  boundaryMax = vec3( getPos( (dim[0]*dim[1]) - 1) ); //last vertex
  boundarySize = boundaryMax - boundaryMin;           // size

  // print to screen
  std::cout << "X-Boundaries : " << boundaryMin[0]  << " .. " << boundaryMax[0]  << std::endl;
  std::cout << "Y-Boundaries : " << boundaryMin[1]  << " .. " << boundaryMax[1]  << std::endl;
  std::cout << "Boundary-Size: " << boundarySize[0] << " | "  << boundarySize[1] << std::endl;


  // check if axes are flipped
  // if X and Y are flipped, going in a row the y value increases and x stays the same
  // (block dataset is flipped)
  float maxYPos = getPosY( dim[0]-1 );
  float avgYPos = ( boundaryMin[1] + boundaryMax[1] ) * 0.5f;

  if ( avgYPos < maxYPos )
  //if ( avgYPos > maxYPos )
  //if (posy > average) // gv Shouldn't it be the other way 'round?
  {
    isFlipped = true;
    std::cout << "Flipped X and Y dimensions." << std::endl;
  }
  else isFlipped = false;


  // ( gv Do not change anything -> do this with OGL rotations !!! ) -> posssibly wrong -> flip axes now

  // gv Finally get rid of the problem with inverted, switched etc.. axes!!
  if ( isFlipped )
  {
    // switch axes NOW!
    // and in the remaining programm noone has to ever think about it again!

    vec3 temp;

    for ( int i=0; i<dim[0]*dim[1]; i++ )
    {
      temp = geometryData[ i ];

      geometryData[ i ].v[0] = temp.v[1];
      geometryData[ i ].v[1] = temp.v[0];
    }

    int dimTemp = dim[X];
    dim[X] = dim[Y];
    dim[Y] = dimTemp;
  } // if flipped


  // gv Move grid to origin
  if ( boundaryMin.v[0] != 0.0 || boundaryMin.v[1] != 0.0 )
    for( int i=0; i<dim[0]*dim[1]; i++ )
      geometryData[i] = geometryData[i] - boundaryMin;


  boundaryMax = boundaryMax - boundaryMin; // alternative: getPos( (dim[0]*dim[1]) - 1)
  boundaryMin = getPos( 0 ); // should always be zero now!

  // gv
  // Create space for normalized grid-points (so that we still have access to the original data)
  // do not invert y, since we do that with OGL rotations later
  normGeometryData = new vec3[ dim[0]*dim[1] ];

  for( int i=0; i<dim[0]*dim[1]; i++ )
    normGeometryData[i] = normalizeCoords( geometryData[i] );


/* // gv Not sure if this is really such a good idea, since I think we just don't need it?!

  // nk Invert y of geometry data (convert from LHS to RHS) and normalize [0...1]
  for(int i = 0; i < getDimX() * getDimY(); i++)
  {
    vec3 pos = geometryData[i];
    pos = normalizeCoords(pos);
    pos[Y] = 1.0f - pos[Y];
    geometryData[i] = pos;
  }
*/

/*
  //first vertex
  boundaryMin = vec3( getPos(0) );

  //last vertex
  boundaryMax = vec3( getPos( (dim[0]*dim[1]) - 1) );

  // size
  boundarySize = boundaryMax - boundaryMin;
*/

  // print new boundaries
  std::cout << "New X Boundaries: " << boundaryMin[0] << " .. " << boundaryMax[0] << std::endl;
  std::cout << "New Y Boundaries: " << boundaryMin[1] << " .. " << boundaryMax[1] << std::endl;

  // gv TEST: print coords of grid-points
  //printGridPoints();

  std::cout << " ------- Done loading grid ------- " << std::endl;

  // all done (hopefully)
  return true;
}


bool FlowGeometry::getIsFlipped()
{
  return isFlipped;
}



/* gv Obsolete?
int FlowGeometry::getXYvtx(vec3 pos)
{
  int i;
  int j;
  //search for the column left to the vertex
  for (i = 0; (i < dim[0])&&(getPosX(getVtx(i,0)) < pos[0]); i++);
  //search for the row under the vertex
  for (j = 0; (j < dim[1])&&(getPosY(getVtx(0,j)) < pos[1]); j++);

  //return the vertex ID of the found vertex
  return getVtx((i<dim[0]) ? i : dim[0]-1, (j<dim[1]) ? j : dim[1]-1);
}
*/

inline vec3 FlowGeometry::getPos(int vtxID)
{
  return geometryData[vtxID];
}

inline float FlowGeometry::getPosX(int vtxID)
{
  return geometryData[vtxID][0];
}

inline float FlowGeometry::getPosY(int vtxID)
{
  return geometryData[vtxID][1];
}


vec3 FlowGeometry::getNormPos(int vtxID)
{
  return normGeometryData[vtxID];
}


static int closest = 0;
//static bool global = true;  // gv Maybe needed for global vs. local search


int FlowGeometry::getNearestVtx(vec3 pos)
{
  //take the norm to the first vertex
  float dist = pos.dist2(getPos(0));
  //mark the vertex 0 as the closest one

  //float newd;
  //Iterate through all vertices and search for the enarest one. I know, this is so inefficient, that you'll have to rewrite it :P
  //for (int i = 1; i < dim[0]*dim[1]; i++)
  //{
  //  newd = pos.dist2(getPos(i));
  //  if (newd < dist)
  //  {
  //    dist = newd;
  //    closest = i;
  //  }
  //}
  //printf("%f, %f, id: %d\n", pos.v[0], pos.v[1], closest);
  return closest++;
}


bool FlowGeometry::getInterpolationAt_old(vec3 pos, int* vtxID, float* coef)
//bool FlowGeometry::getInterpolationAt_old(vec3 pos, int* vtxID, float* coef)
{
  std::cout << "getInterpolationAt... (old)" << std::endl;

  //if we are outside of the dataset, return false
  //please note, that this test is valid only for rectangular datasets (block, hurricane), but not for tube
  if ((pos[0]<boundaryMin[0])||(pos[1]<boundaryMin[1])||(pos[0]>boundaryMax[0])||(pos[1]>boundaryMax[1]))
    return false;

  //example of a low-quality lookup with no interpolation
  vtxID[0] = getNearestVtx(pos); //finds the nearest vertex to the given position and returns it's value as the dominanting one
  vtxID[1] = 0;
  vtxID[2] = 0;
  vtxID[3] = 0;
  coef[0] = 1.0f;
  coef[1] = 0.0f;
  coef[2] = 0.0f;
  coef[3] = 0.0f;
  return true;
}


bool FlowGeometry::getInterpolationAt( vec3 pos, vec3* value, FlowData* data )
{
  //std::cout << "getInterpolationAt... (gv)" << std::endl;


  // lookup position in grid
  int xGridIdx = binarySearch( pos.v[0], 0, dim[X]-1, 0 );
  int yGridIdx = binarySearch( pos.v[1], 0, dim[Y]-1, 1 );

  //std::cout << "x-index found for " << pos.v[0] << ": " << xGridIdx << std::endl;
  //std::cout << "y-index found for " << pos.v[1] << ": " << yGridIdx << std::endl;

  // return false if not in grid
  if ( xGridIdx == -1 || yGridIdx == -1 ) return false;

  // interpolate cell-points
  int cell[2] = { xGridIdx, yGridIdx };  // lower-left point (=index) of cell containing the searchpoint
  vec3 interpolatedValue[1]; 
  bool retCode = interpolate( pos, cell, data, interpolatedValue );
  //vec3 interpolatedValue = interpolate( pos, cell, data );

  //std::cout << "interpolated value: " << interpolatedValue[0] << std::endl;
  value[0] = interpolatedValue[0];

  // TEST: print coords of grid-points
  //printGridPoints();

  // everything ok
  //return true;
  return retCode;
}


bool FlowGeometry::interpolate( vec3 pos, int* cell, FlowData* data, vec3 *interpolatedValue )
{
  // indices of the cellpoints (cp):
  //
  //  (0,1)    (1,1)
  //    |--------|
  //    |        |
  //    |--------|
  //  (0,0)    (1,0)

  // calculate cell-incides
  int idxCp00 = cell[X] + cell[Y] * dim[X];
  int idxCp10 = idxCp00 + 1;
  int idxCp01 = idxCp00 + dim[X];
  int idxCp11 = idxCp01 + 1;

  // lookup cell-positions
  vec3 cp00 = normGeometryData[ idxCp00 ];
  vec3 cp10 = normGeometryData[ idxCp10 ];
  vec3 cp01 = normGeometryData[ idxCp01 ];
  vec3 cp11 = normGeometryData[ idxCp11 ];
/*
  vec3 cp00 = geometryData[ idxCp00 ];
  vec3 cp10 = geometryData[ idxCp10 ];
  vec3 cp01 = geometryData[ idxCp01 ];
  vec3 cp11 = geometryData[ idxCp11 ];
*/
  // the coefficients for the interpolation
  float coeff00 = 0.0f;
  float coeff10 = 0.0f;
  float coeff01 = 0.0f;
  float coeff11 = 0.0f;

  // we need to calcualte the coefficients from the grid
  // Then interpolate the VELOCITIY VALUES from this grid-points with theri coeff

  // vec3 interpolatedVelocityValue = velocity[cp00] * coeff00 + velocity[cp10] * coeff10 
  //                                + velocity[cp01] * coeff01 + velocity[cp11] * coeff11;


  // calculate the coeff 
  // with weighted-square-method

  // cp0   v   pos      w        cp1
  //  |---------|-----------------| 

  float vx = pos[X] - cp00[X];
  float wx = cp10[X] - pos[X];

  float vy = pos[Y] - cp00[Y];
  float wy = cp01[Y] - pos[Y];

  float maxX = cp10[X] - cp00[X];
  float maxY = cp01[Y] - cp00[Y];

  coeff00 = (wx * wy) / ( maxX * maxY );
  coeff10 = (vx * wy) / ( maxX * maxY );
  coeff01 = (wx * vy) / ( maxX * maxY );
  coeff11 = (vx * vy) / ( maxX * maxY );

/*
  // TEST
  // test == pos if correct coeffcients
  vec3 test = cp00 * coeff00 + 
              cp10 * coeff10 + 
              cp01 * coeff01 + 
              cp11 * coeff11;
*/

  // calculate interpolated value
  FlowChannel* velX = data->getChannel( X );
  FlowChannel* velY = data->getChannel( Y );
  FlowChannel* velZ = data->getChannel( Z );
/*
  int ti = 152 + 203 * dim[X];
  vec3 tv( velX->getValue( ti ), velY->getValue( ti ), velZ->getValue( ti ) );
  cout << "tv: " << tv << "\n";
*/
  // lookup velocity-values
  vec3 vel00( velX->getValue( idxCp00 ), velY->getValue( idxCp00 ), velZ->getValue( idxCp00 ) );
  vec3 vel10( velX->getValue( idxCp10 ), velY->getValue( idxCp10 ), velZ->getValue( idxCp10 ) );
  vec3 vel01( velX->getValue( idxCp01 ), velY->getValue( idxCp01 ), velZ->getValue( idxCp01 ) );
  vec3 vel11( velX->getValue( idxCp11 ), velY->getValue( idxCp11 ), velZ->getValue( idxCp11 ) );
/*
  cout << "v00: " << vel00
       << "\nv10: " << vel10
       << "\nv01: " << vel01
       << "\nv11: " << vel11 << "\n";
*/

  if ( ((vel00[X] == 0) && (vel00[Y] == 0)) ||
       ((vel10[X] == 0) && (vel10[Y] == 0)) ||
       ((vel01[X] == 0) && (vel01[Y] == 0)) ||
       ((vel11[X] == 0) && (vel11[Y] == 0)) )
  {
    /*
    cout << "there is a null-vector !!! pos: " << pos << "\n";
    cout << "cp00: " << idxCp00
       << "\ncp10: " << idxCp10
       << "\ncp01: " << idxCp01
       << "\ncp11: " << idxCp11 << "\n";
    */
    return false;
  }

  /*vec3*/ interpolatedValue[0] = vel00 * coeff00 + 
                              vel10 * coeff10 + 
                              vel01 * coeff01 + 
                              vel11 * coeff11;


/*
  // print values (for debug)
  std::cout << "\n\n --- Print interpolation info --- " << std::endl;
  std::cout << "Spos: " << pos << std::endl << std::endl;
  std::cout << "TEST: " << test << std::endl << std::endl;

  std::cout << "idxCp00: " << idxCp00 << std::endl;
  std::cout << "idxCp10: " << idxCp10 << std::endl;
  std::cout << "idxCp01: " << idxCp01 << std::endl;
  std::cout << "idxCp11: " << idxCp11 << std::endl << std::endl;

  std::cout << "grid [ cp00 ] : " << geometryData[ idxCp00 ] << std::endl;
  std::cout << "grid [ cp10 ] : " << geometryData[ idxCp10 ] << std::endl;
  std::cout << "grid [ cp01 ] : " << geometryData[ idxCp01 ] << std::endl;
  std::cout << "grid [ cp11 ] : " << geometryData[ idxCp11 ] << std::endl << std::endl;

  std::cout << "vx: " << vx << ", wx: " << wx << std::endl;
  std::cout << "vy: " << vy << ", wy: " << wy << std::endl << std::endl;

  std::cout << "vel00: " << vel00 << std::endl;
  std::cout << "vel10: " << vel10 << std::endl;
  std::cout << "vel01: " << vel01 << std::endl;
  std::cout << "vel11: " << vel11 << std::endl << std::endl;

  std::cout << "coeff00: " << coeff00 << std::endl;
  std::cout << "coeff10: " << coeff10 << std::endl;
  std::cout << "coeff01: " << coeff01 << std::endl;
  std::cout << "coeff11: " << coeff11 << std::endl << std::endl;

  std::cout << "iPol: " << interpolatedValue << std::endl << std::endl;
*/

  return true;
  //return interpolatedValue;
}


int FlowGeometry::binarySearch( float searchValue, int begin, int end, int dimension )
{
  bool found    = false;
  int  leftEnd  = begin;
  int  rightEnd = end;
  int  mid      = 0;
  int  currentIndex = 0;


  while ( !found && leftEnd <= rightEnd )
  {
    float currentValue;
    mid = ( leftEnd + rightEnd ) / 2;
    currentIndex = mid;

    currentValue = getCurrentSearchValue( currentIndex, dimension );

    if ( currentValue == searchValue )
      // found value
      found = true;
    else
    {

      if ( currentValue < searchValue  )
      {
        // check if neighborhood + 1 < searchValue
        currentIndex = mid + 1;
        currentValue = getCurrentSearchValue( currentIndex, dimension );

        if ( currentValue > searchValue )
        {
          // found value
          found = true;
          --currentIndex; // we want the lower left corner of the cell
        }
      }

      else // if ( currentValue > searchValue )
      {
        // check if neighborhood - 1 < searchValue
        currentIndex = mid - 1;
        currentValue = getCurrentSearchValue( currentIndex, dimension );

        if ( currentValue < searchValue )
        {
          // found value
          found = true;
        }
      }

      if ( !found ) // not yet found
      {
        if ( searchValue < currentValue )
          // continue search in left half
          rightEnd = mid - 1;
        else
          // continue search in right half
          leftEnd = mid + 1;
      }
    }
  } // while

  if ( found )
    // return found index
    return currentIndex;
  else
    // not found
    return -1;
}


int FlowGeometry::binarySearch2( float searchValue, int begin, int end, int dimension )
{
  bool found    = false;
  int  leftEnd  = begin;
  int  rightEnd = end;
  int  mid      = 0;
  int  currentIndex = 0;

  // Problem: What, if value is not in array? -> returns not found
  // But: We want that the next lower value is returned!
  // => adapt this to be correct!
  //
  // Idea: keep closest match + it's distance (really neccessary?)
  // Search as long as leftEnd < rightEnd
  // If leftEnd > rightEnd, we have gone too far
  // Then, do sequential search to maxDim, stop if currentLower < searchValue && currentBigger > searchValue
  // (should work, shouldn't it?)
  //
  // We than have the index of the lower-left corner of the grid-cell containing the given point
  // Return -1 immediately if searchVAlue out of boundaries!

/*
  if ( dimension == X )
  {
    float minVal = geometryData[ leftEnd  ].v[0];
    float maxVal = geometryData[ rightEnd ].v[0];

    float dummy = 0.0f;

  }
*/

  while ( !found && leftEnd <= rightEnd )
  {
    float currentValue;
    mid = ( leftEnd + rightEnd ) / 2;
    currentIndex = mid;

    currentValue = getCurrentSearchValue( currentIndex, dimension );

    if ( currentValue == searchValue )
      // found value
      found = true;
    else
    {
      // check neighborhood + 1
      currentIndex = mid + 1;
      currentValue = getCurrentSearchValue( currentIndex, dimension );

      if ( currentValue == searchValue )
        // found value
        found = true;
      else
      {
        // check neighborhood - 1
        currentIndex = mid - 1;
        currentValue = getCurrentSearchValue( currentIndex, dimension );

        if ( currentValue == searchValue )
          // found value
          found = true;

        else {
          // not yet found

          if ( searchValue < currentValue )
            // continue search in left half
            rightEnd = mid - 1;
          else
            // continue search in right half
            leftEnd = mid + 1;
        }
      }
    }
  } // while

  if ( found )
    // return found index
    return currentIndex;
  else
    // not found
    return -1;
}


float FlowGeometry::getCurrentSearchValue( int index, int dimension )
{
  float value;

  if ( index == -1 )
  {
    //std::cout << "! Error: index -1 in FlowGeometry.getCurrentSearchValue!" << std::endl;
    return 0.0f;
  }

  if ( dimension == X )
    value = normGeometryData[ index ].v[0];
  else // dimension == Y
    value = normGeometryData[ dim[X] * index ].v[1];
/*
  if ( dimension == X )
    value = geometryData[ index ].v[0];
  else // dimension == Y
    value = geometryData[ dim[X] * index ].v[1];
*/
  return value;
}


void FlowGeometry::printGridPoints()
{
  std::cout << "Writing grid-point-coords to log-file (this may take a while - please be patient!) ... ";

  ofstream gridLogFile;

  gridLogFile.open ("grid_points.log");
  gridLogFile << "Grid-Point-Positions with their corresponding indices in geometryData:\n\n";

/*
  // all
  for ( int i=0; i<dim[X]*dim[Y]; i++ )
    gridLogFile << i << ": " << geometryData[ i ] << "\n";
*/


  // X-dim
  gridLogFile << "X:\n";
  for ( int i=0; i<dim[X]; i++ )
  {
    gridLogFile << i << ": " << normGeometryData[ i ].v[0] << "\n";
  }

  // Y-dim
  gridLogFile << "\n\nY:\n";
  for ( int i=0; i<dim[Y]; i++ )
  {
    gridLogFile << i << ": " << normGeometryData[ i * dim[X] ].v[1] << "\n";
  }


  gridLogFile.close();

  std::cout << " done. " << std::endl;
}



inline float FlowGeometry::getMinX()
{
  return boundaryMin[0];
}

inline float FlowGeometry::getMaxX()
{
  return boundaryMax[0];
}

inline float FlowGeometry::getMinY()
{
  return boundaryMin[1];
}

inline float FlowGeometry::getMaxY()
{
  return boundaryMax[1];
}

inline int FlowGeometry::getDimX()
{
  return dim[0];
  //return (isFlipped) ? dim[1] : dim[0];
}

inline int FlowGeometry::getDimY()
{
  return dim[1];
  //return (isFlipped) ? dim[0] : dim[1];
}

inline int FlowGeometry::getDimZ()
{
  return dim[2];
}

/* gv Obsolete?
int FlowGeometry::getVtx(int x, int y)
{
  return ( y*dim[0] ) + x;

  //if we need to flip the rows and columns, we do it here
  //return (isFlipped)? (x*dim[1]) + y : (y*dim[0]) + x;
}

int FlowGeometry::getVtxX(int vtxID)
{
  return vtxID % dim[0];

  //if we need to flip the rows and columns, we do it here
  //return (isFlipped)? vtxID / dim[1] : vtxID % dim[0];
}

int FlowGeometry::getVtxY(int vtxID)
{
  return vtxID / dim[0];

  //if we need to flip the rows and columns, we do it here
  //return (isFlipped)? vtxID % dim[1] : vtxID / dim[0];
}

int FlowGeometry::getRightNeigh(int vtxID)
{
  int x = getVtxX(vtxID);
  return (x+1 < dim[0]) ? getVtx(x+1,getVtxY(vtxID)) : -1;
}

int FlowGeometry::getTopNeigh(int vtxID)
{
  //remember that the data is structured with (0,0) in the upper-left corner and (1,1) in the lower-right
  //that's why we are subtracting 1 to find the top neighbour
  int y = getVtxY(vtxID);
  return (y+1 < dim[1]) ? getVtx(getVtxX(vtxID), y-1) : -1;
}

int FlowGeometry::getLeftNeigh(int vtxID)
{
  int x = getVtxX(vtxID);
  return (x > 1) ? getVtx(x-1,getVtxY(vtxID)) : -1;
}

int FlowGeometry::getBottomNeigh(int vtxID)
{
  //remember that the data is structured with (0,0) in the upper-left corner and (1,1) in the lower-right
  int y = getVtxY(vtxID);
  return (y < 1) ? getVtx(getVtxX(vtxID), y+1) : -1;
}
*/
vec3 FlowGeometry::normalizeCoords(vec3 pos)
{
  vec3 u(pos - boundaryMin);
  //scale each component according to the side length
  u[0] /= boundarySize[0];
  u[1] /= boundarySize[1];

  return u;
}

vec3 FlowGeometry::unNormalizeCoords(vec3 pos)
{
  vec3 u(pos);
  //multiply each component according to the side length
  u[0] *= boundarySize[0];
  u[1] *= boundarySize[1];

  return u += boundaryMin;
}

---