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

Go to the documentation of this file.

#include "FlowData.h"
#include "reverseBytes.h"

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


FlowData::FlowData()
{
  geometry = new FlowGeometry();

  //mark all the channel slots as free
  for(int i = 0; i < max_channels; i++)
  {
    channels[i] = NULL;
    freeChannel[i] = true;
    scalarChannel[i] = NONE;
  }
}

FlowData::~FlowData()
{
  delete geometry;

  //delete all the channels
  for(int i = 0; i < max_channels; i++)
  {
    if (!freeChannel[i])
    {
      deleteChannel(i);
    }
  }
}

bool FlowData::loadDataset( string filename, bool bigEndian )
{
  FILE* griFile = NULL;
  FILE* datFile = NULL; 

  char header[40];  // storeage for grid-file-header


  std::cout << "\n\n ------- READ FILE ------- " << std::endl;

  //localize the last dot in the filename 
  int lastdot = (int) filename.find_last_of(".",filename.length()-1);
  if (lastdot != string::npos)
    //if there is a dot, remove everything behind it
    filename = filename.substr(0,lastdot);      


  // GRID FILE
  string griName = filename + ".gri";
  std::cout << "- Loading grid file '" << griName << "' ... " << std::endl;

  //open the grid file
  griFile = fopen(griName.c_str(),"rb");

  if (!griFile)
  {
    std::cerr << "+ Error loading grid file:" << griName << std::endl << std::endl;
    return false;
  }

  //read the header from the grid-file
  fread( header, 40, 1, griFile );

  //pass the grid file to the geometry class to process it
  if ( !geometry->readFromFile( header, griFile, bigEndian ) )
    return false; // exit on error

  //close the file
  fclose( griFile );


  // --- done with grid-file


  // DAT FILE

  std::cout << "\n\n ------- Loading data file ------- " << std::endl;

  float DT; // gv Actually not used?

  //read some necessary data from the grid-file-header
  sscanf( header, "SN4DB %d %d %d %d %d %f", &dimX, &dimY, &dimZ, &numChannels, &timesteps, &DT );
  std::cout << "Dimensions: (" << dimX << ", " << dimY << ", " << dimZ << ")" << std::endl;
  std::cout << "Channels: " << numChannels << "\nTimesteps: " << timesteps << std::endl;         

  //since this framework handles only one timestep, we use only a single dat file
  char suffix[16];
  sprintf(suffix,".%.5u.dat",0); //the second dot and the following 5 specify that a minimum of 5 numbers will be written
  string datName = filename.append(suffix);
  std::cout << "- Loading grid file '" << datName << "' ... " << std::endl;    


  //open the dat-file
  datFile = fopen( datName.c_str(), "rb" );
  if ( !datFile )
  {
    std::cerr << "+ Error loading dat file:" << datName << std::endl << std::endl;
    return false;
  }

  numScalarChannels = numChannels;
  numChannels += 3; //add the 3 components of the velocity vector to the number of additional chanenls

  //because reading big chunks of data is much faster than single values, 
  //we read the data into a temporary array and then copy it to the channels

  //create temporary storage
  float* tmpArray = new float[ numChannels * geometry->getDimX() * geometry->getDimY() ];

  //read the data
  int result = (int) fread( tmpArray, sizeof(float), numChannels * geometry->getDimX() * geometry->getDimY(), datFile ); 

  //have we read the whole data file?
  if (result != numChannels * geometry->getDimX() * geometry->getDimY() )
  {
    std::cerr << "+ Error reading dat file:" << datName << std::endl << std::endl;
    return false;
  }

  // close the file
  fclose( datFile );

  // swap the byte order, if the file is encoded big-endian
  // (usually not necessary)
  if ( bigEndian )
    for(int j = 0; j < numChannels*geometry->getDimX()*geometry->getDimY(); j++)
      tmpArray[j] = reverseBytes<float>(tmpArray[j]);


  // --- begin processing data ---

  int currentChannel = -1;
  string channelName = "unkonwn channel data";

  int channelsToLoad = numChannels;  // gv Number of channels to load from file
  numChannels = 0;  // gv Total number of channels created in channels[]

  // assign the data to the appropriate channels
  for ( int j=0; j<channelsToLoad; j++ )
  {
   
    // just give the channels names for better debugging
    switch ( j )
    {
      case 0:
        channelName = "Velocity-X";
        break;

      case 1:
        channelName = "Velocity-Y";
        break;

      case 2:
        channelName = "Velocity-Z";
        break;

      case 3:
        channelName = "Pressure";
        scalarChannel[j] = ORIGINAL;
        break;

      case 4:
        channelName = "Vorticity";
        scalarChannel[j] = ORIGINAL;
        break;

      default:
        channelName = "unknown data";
        break;
    } // switch
    
    //create the new channel
    currentChannel = createChannel( channelName );

    //copy the values of the jth channel from tmpArray, which carries numChannels
    channels[ currentChannel ]->copyValues( tmpArray, channelsToLoad, j );

    //std::cout << "min: " << channels[ currentChannel ]->getMin() << std::endl;


/* gv We do not manipulate data here -> rotations with OGL !!!
    // gv Not sure, but do we really need to do this??
    // please read comments below...!

    if ( j==Y ) // velocity channel y - invert y
    {
      for( int i=0; i<dimX * dimY; i++ )
      {
        float val = channels[ j ]->getValue( i );
        channels[ j ]->setValue( i, -val );
      }
    } // if j==Y
*/


/*  
    // gv Wrong normalisazion. Correct, see below
    // If minimum > 1.0 => normalisation returns false results!

    if ( j>Z ) // scalar data in channel -> normalize to [0..1] 
    {
      float val;

      // normalize values
      for( int i=0; i<dimX * dimY; i++ )
      {
        val = channels[j]->getValue( i );
        val = channels[j]->normalizeValue( val );
        channels[j]->setValue( i, val );
      }
      // scalar data for background
      scalarChannel[j] = true;
    }
*/

  /* --- READ !! ---
    (gv)

    Not sure if this is correct now, since this time my brain seems to be overloaded! :(

    German now: ;)

    Wenn ich beim Grid de Achsen x<->y vertausche, bleiben de Indizes der Gitter-Punkte gleich. 
    Es ändert sich nur die Auslegung des Gitters.

    Wenn ich nun Skalar-Daten habe, haben die auch immer den gleichen Index im Array. Es ändert sich 
    für diese Daten ja nicht die Position auf dem Gitter (also der Index), sondern es ändert sich nur de Gitter-Position 
    auf denen diese Skalar-Werte liegen. 
    D.h., laut meiner Überlegung müßte es für skalare Daten egal sein, ob die Achsen geflippt werden oder nicht. 
    => Für Skalar-Daten keine Behandlung wenn flipped notwendig.

    Anders bei Vektor-Daten, da sich mit den gedrehten Achse auch die Komponenten der Vektoren vertauschen, 
    sonst zeigen sie in die falsche Richtung. D.h., man muß bei Vektor-Daten ebenfalls die Achsen tauschen, damit diese 
    in die richtige Richtung zeigen.
    
    => Dies ist bei uns nur bei den Velocity-Daten mit X u. Y der Fall. Wenn im Grid de Achsen vertauscht werden, müssen 
    wir de Velocity-Channels X und Y vertauschen. 

    Die Y-Achse invertiere ich zur Zeit nicht, da sie im Grid auch nicht mehr invertiert wird, weil somit 
    de Achsen von links-unten weggehen:

      y
      ^
      |
      |-------> x
    (0,0)

    Deshalb invertiere ich die Daten auch nicht mehr in Y-Richtung, weil sonst würd sich das Bild meiner Meinung
    nach umdrehen.

    Das sind meine Überlegungen dazu.
    Ich bin mir aber nun echt nicht sicher, ob das stimmt!!

    Bitte überleg dir das mal durch und sag mir Bescheid!

    Danke!!

    LG,
      Gü.

   --- */



  // gv If axes are flipped, flip vector-data
  if ( geometry->getIsFlipped() )
  {
    std::cout << "Flipping data-channel..." << std::endl;

/* gv Possibily wrong idea
    // switching velocity x and y data

    // switch x and y axes
    // for vector-data only!
    // scalar-data not affected by flipped axes, since the indices for the data stay all the same!
    // (hope this is true!)
    
    float tempX, tempY;

    for ( int i=0; i<dimX*dimY; i++ ) 
    {
      tempX = channels[ 0 ]->getValue( i );
      tempY = channels[ 1 ]->getValue( i );

      channels[ 0 ]->setValue( i, tempY );
      channels[ 1 ]->setValue( i, tempX );
    }
*/

    // Generate new array and switch indices
    float *newArray = new float[ dimX*dimY ];
    int newIndex = -1;
    int oldIndex = -1;

    for ( int y=0; y<dimY; y++ )
      for ( int x=0; x<dimX; x++ )
      {
        newIndex = y + x * dimY;
        oldIndex = x + y * dimX;

        newArray[ newIndex ] = channels[ currentChannel ]->getValue( oldIndex );

        // gv TEST if indices are written correct
        //newArray[ newIndex ] = oldIndex;
      }

  /* gv Old version - possibly wrong
    for ( int x=0; x<dimX; x++ )
      for ( int y=0; y<dimY; y++ )
      {
        newIndex = x + y * dimX;
        oldIndex = y + x * dimY;

        //newArray[ newIndex ] = channels[ currentChannel ]->getValue( oldIndex );
        newArray[ newIndex ] = newIndex;
      }
    */

    // copy the values from the new array to the channel: channels[currentChannel] <= newArray
    channels[ currentChannel ]->copyValues( newArray );

    delete [] newArray;
  } // if flipped


  } // for channelsToLoad


  // gv After we processed all the channels, we need to switch the dimensions
  // if the axes are flipped
  if ( geometry->getIsFlipped() )
  {
    int temp = dimX;
    dimX = dimY;
    dimY = temp;
  }


  // Normalize data and store in new channel (so that we still have the original ones)
  for ( int i=3; i<(numScalarChannels+3); i++ ) // for all scalar-channels (usually there are two of them)
  {
    std::cout << "Creating a new channel for normalized data of " << channelNames[i] << "..." << std::endl;

    //channels[ i ]->moveMinimumToZero();  // gv Actually, not neccessary

    float *normArray = new float[ dimX*dimY ];

    for ( int j=0; j<dimX*dimY; j++ )
    {
      // normalize data
      float data = channels[ i ]->getValue( j );
/*
      if(channels[i]->getMin() < 0 && data == 0.0)
      {
        normArray[ j ] = data;
      }
      else
      {
        normArray[ j ] = channels[ i ]->normalizeValue( data );
      }
*/
      // gv the code above is not correct -> min-values and zero-values are normalised alltogether to zero !!
      normArray[ j ] = channels[ i ]->normalizeValue( data );
    }

    // copy the values from the array to a new channel
    int newChannel = createChannel( "Normalized " + channelNames[i] );
    channels[ newChannel ]->copyValues( normArray );

    // nk - normalized scalar data for background
    scalarChannel[ newChannel ] = NORMALIZED;

    delete [] normArray;
  } // for normalize


  // create channel with magnitude of velocity-data ( gv Seems to be ok - not tested )
  velocityMagnitudeId = createChannelVectorLength( X, Y, Z );

  // normalize velocity-vectors (gv Seems to be ok - not tested )
  normalizeVelocityChannels( X, Y, Z );

  delete[] tmpArray;
  
  // gv TEST: write log-file
  //printData();

/*
  // gv TEST
  // check, if velZ contains any data != 0
  // => only hurricane has velZ-data! => we need to take care of that!
  for ( int i=0; i<dimX*dimY; i++ )
  {
    if ( channels[2]->getValue(i) != 0 ) std::cout << " ! VelZ: " << channels[2]->getValue(i) << std::endl;
  }
*/

  std::cout << " ------- Done loading data file ------- \n\n" << std::endl;

  // hopefully everything done well
  return true;
}



void FlowData::printData() 
{
  std::cout << "Writing log-file (this may take a while) ... please be patient!";

  ofstream dataLogFile;

  dataLogFile.open ("data.log");
  dataLogFile << "Channel-Data with their corresponding indices in channels:\n";

  // write all data in all channels 
  for ( int j=0; j<numChannels; j++ )
  {
    dataLogFile << "\nChannel " << j << " - " <<  channelNames[j] << ":\n";

    // print x/y 
    for ( int y=0; y<dimY; y++ )
    {
      dataLogFile << "\n\nY: " << y << ":\n";

      for ( int x=0; x<dimX; x++ )
      {
        dataLogFile << x << ": " << channels[j]->getValue( y*dimX + x ) << "; ";
        if ( (x % 10) == 9 )
          dataLogFile << "\n";
      }
    }


    /*
    // print array
    for ( int i=0; i<dimX * dimY; i++ )
    {
      dataLogFile << i << ": " << channels[j]->getValue( i ) << "\n";
    }
*/

  }

  dataLogFile.close();

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

ScalarDataT FlowData::getIsScalar( int i )
{
  if( i >= max_channels )
  {
    return NONE;
  }
  return scalarChannel[i];   
}


int FlowData::getChannelCount()
{
  return numChannels;
}


//VelocityT* FlowData::getVelocity() const
//{
//  return velocity;
//}


int FlowData::createChannel( string channelName )
{
  //find the first unused channel slot
  int i = 0;

  while ( (i<max_channels) && (!freeChannel[i]) ) i++;

  
  if ( i<max_channels ) //if there is a free slot
  {
    std::cout << "Creating channel " << channelName << " at " << i << std::endl;

    //create a new channel
    channels[i] = new FlowChannel(geometry);

    //remember the slot
    freeChannel[i] = false;

    // name the channel
    channelNames[i] = channelName;

    // increase numChannels
    ++numChannels;

    //return the adress of the new channel
    return i;
  }
  else //there is no free channel
  {
    std::cerr << "There is no free channel slot!" << std::endl;
    return -1;
  }
}


void FlowData::deleteChannel( int i )
{
  if ( !freeChannel[i] ) //if the address is really occupied
  {
    std::cout << "Deleting channel at " << i << " ... ";

    //delete the channel instance
    delete channels[i];

    channels[i] = NULL;

    //free the slot
    freeChannel[i] = true;
  }
  else 
    std::cout << "Tried to delete a non-existing channel at " << i << "." << std::endl;
}


FlowChannel* FlowData::getChannel( int i )
{
  return channels[i];
}

/* gv Obsolete?
int FlowData::createChannelGeometry( int dimension )
{
  int result = createChannel();

  //just take the dimension as if it was an offset to the geometryData array
  channels[result]->copyValues( (float*) geometry->geometryData, 3, dimension );

  return result;
}
*/

int FlowData::createChannelVectorLength( FlowChannel* chX, FlowChannel* chY, FlowChannel* chZ )
{
  int newChannel = createChannel( "Velocity-Length" );

  //check whether we deal with 2D or 3D vectors
  if (chZ)
  {
    for (int i = 0; i < geometry->getDimX()*geometry->getDimY(); i++)
    {
      //save the vector length
      channels[newChannel]->setValue(i,sqrt(chX->getValue(i)*chX->getValue(i) + chY->getValue(i)*chY->getValue(i) + chZ->getValue(i)*chZ->getValue(i)));
    }
  }
  else
  {
    for (int i = 0; i < geometry->getDimX()*geometry->getDimY(); i++)
    {
      channels[newChannel]->setValue(i,sqrt(chX->getValue(i)*chX->getValue(i) + chY->getValue(i)*chY->getValue(i)));
    }
  }
  return newChannel;
}


int FlowData::createChannelVectorLength( int chX, int chY, int chZ )
{
        //just a wrapper for the method above
        if (chZ >= 0)
  {
                return createChannelVectorLength(getChannel(chX),getChannel(chY),getChannel(chZ));
  }
        else 
  {
    return createChannelVectorLength(getChannel(chX),getChannel(chY));
  }
}


void FlowData::normalizeVelocityChannels(FlowChannel* chX, FlowChannel* chY, FlowChannel* chZ)
{
  if (chZ)
  {
    for (int i = 0; i < geometry->getDimX()*geometry->getDimY(); i++)
    {
      // create vec3 from data, normalize and save back
      vec3 velocity = vec3( chX->getValue(i), chY->getValue(i), chZ->getValue(i) );

      velocity.normalize();

      chX->setValue( i, velocity[X] );
      chY->setValue( i, velocity[Y] );
      chZ->setValue( i, velocity[Z] );
    }
  }
  else
  {
    for (int i = 0; i < geometry->getDimX()*geometry->getDimY(); i++)
    {
      // create vec3 from data, normalize and save back
      vec3 velocity = vec3( chX->getValue(i), chY->getValue(i) );
      velocity.normalize();

      chX->setValue( i, velocity[X] );
      chY->setValue( i, velocity[Y] );
    }
  }
}


void FlowData::normalizeVelocityChannels(int chX, int chY, int chZ)
{
  std::cout << "Normalizing (velocity-)channels " 
    << chX << ", " 
    << chY << ", " 
    << chZ << "." << std::endl;

        if (chZ >= 0)
  {
                normalizeVelocityChannels(getChannel(chX),getChannel(chY),getChannel(chZ));
  }
        else 
  {
    normalizeVelocityChannels(getChannel(chX),getChannel(chY));
  }
}

#if 0
int FlowData::createChannelVectorLength(FlowChannel* chX, FlowChannel* chY, FlowChannel* chZ)
{
  int result = createChannel();
  float mag;
  vec3 vel;
  //check whether we deal with 2D or 3D vectors
  if (chZ)
  {
    for (int i = 0; i < geometry->getDimX() * geometry->getDimY(); i++)
    {
      // get velocity vectors -> INVERT Y-Value due to LHS->RHS
      vel.v[X] = chX->getValue(i);
      vel.v[Y] = -chY->getValue(i);
      vel.v[Z] = chZ->getValue(i);

      // get magnitude of velocity
      mag = vel.length();

      //save the vector length to velocity struct & channel
      velocity[i].magnitude = mag;
      channels[result]->setValue(i, mag);

      // normalize velocity vector & save to velocity struct
      vel.normalize();
      velocity[i].normVelocity = vel;
    }
  }
  else
  {
    for (int i = 0; i < geometry->getDimX() * geometry->getDimY(); i++)
    {
      // get velocity vectors
      vel.v[X] = chX->getValue(i);
      vel.v[Y] = chY->getValue(i);

      // get magnitude of velocity
      mag = vel.length();

      //save the vector length to velocity struct & channel
      velocity[i].magnitude = mag;
      channels[result]->setValue(i, mag);

      // normalize velocity vector & save to velocity struct
      velocity[i].normVelocity = vel.normalize();
    }
  }

  return result;
}

int FlowData::createChannelVectorLength(int chX, int chY, int chZ)
{
  //just a wrapper for the method above
  if (chZ >= 0)
  {
    return createChannelVectorLength(getChannel(chX),getChannel(chY),getChannel(chZ));
  }
  else 
  {
    return createChannelVectorLength(getChannel(chX),getChannel(chY));
  }
}
#endif 


int FlowData::getScalarChannelCount()
{
  return numScalarChannels;
}


int FlowData::getNumTimesteps()
{
  return timesteps;
}


FlowGeometry* FlowData::getGeometry() const
{
  return geometry; 
}


//int FlowData::getGeomID (int dim)
//{
//  switch(dim)
//  {
//  case X: return geomX;
//          break;
//
//  case Y: return geomY;
//          break;
//
//  default: printf("Wrong geometry dimension specified\n");
//          return -1;
//          break;
//  }
//}


int FlowData::getVelMagID ()
{
  return velocityMagnitudeId;
}


void FlowData::rescaleScalarData()
{
}

---