path: root/src/vtkAMRStructuredPointsReader.cxx

#include "vtkAMRStructuredPointsReader.h"
#include "vtkScalars.h"
#include "vtkFloatArray.h"
#include "vtkIntArray.h"
#include "vtkShortArray.h"
#include "vtkUnsignedShortArray.h"
#include "vtkUnsignedCharArray.h"
#include <stdio.h>
#include <string.h>

vtkAMRStructuredPointsReader::vtkAMRStructuredPointsReader():hdf_file(0),file(0),index(0){
  filename[0]='\0';
  skipmap=0;
  curset=curtime=curlevel=0;
  maxset=maxtime=maxlevel=0;
  myreader=NULL;
  SetLoggedOff();
  fprintf(stderr,"before loggedrange\n");
  loggedmax=0.0;
  loggedmin=0.0;
  fprintf(stderr,"after loggedrange\n");
}


unsigned long vtkAMRStructuredPointsReader::GetMTime()
{
  unsigned long dtime = this->vtkSource::GetMTime();
  // get last modification time... need to set this myself...
  //unsigned long rtime = this->Reader.GetMTime();
  //return (dtime > rtime ? dtime : rtime);
  return dtime;
}

// Description:
// Specify file name of vtk polygonal data file to read.
// immediately constructs the file handle, but reading
// doesn't occur until Execute is called (through Update method).

void vtkAMRStructuredPointsReader::SetFileName(char *name) {
  if(file) delete file; // reopen for new filename
  file=0;
  hdf_file=0;
  if(!name){
    vtkErrorMacro(<<"Null Datafile!");
    return;
  }
  // statically allocating filename for right now
  strncpy(filename,name,sizeof(filename)); // copy the name (for printing info)
  filename[sizeof(filename)-1]='\0'; // cap it off;
  file = new IEEEIO(name,IObase::Read);
  if(!file->isValid()){
    delete file;
    file=0;
    file = hdf_file = new HDFIO(name,IObase::Read);
    if(!file->isValid()){
      vtkErrorMacro(<<"Not a valid IEEEIO or HDF-SDS file: "<<name);
      delete file;
      file=0;
      hdf_file=0;
      filename[0]='\0';
      return;
    }
  }
  // here is where we should scan to differentiate datasets from coords
  if(skipmap) delete skipmap;
  int nds=file->nDatasets();
  skipmap = new int[nds];
  for(index=0,maxindex=0;index<nds;index++){
    file->seek(index);
    if(!hdf_file || hdf_file->isCoord()) continue;
    else skipmap[maxindex++] = index;
  }
  maxindex--;
  index = skipmap[0]; // reset the file index to 0; (should scan for datasets)
  myreader=new Brads_AmrFileReader(*file);
  if (myreader==NULL){
    vtkErrorMacro(<<"Not a valid AMR file: "<<name);
    return;
  }
  myreader->setDataLoadingOff();
  maxlevel=myreader->getNumLevels()-1;
  myreader->getTimeRange(mintime,maxtime);
  int ii=0;
  for (ii;ii<=maxlevel;ii++) myreader->hideLevel(ii);
  myreader->showLevel(curlevel);
  myreader->setTime(curtime);
  maxset=myreader->getNumGrids();
  Modified(); // this has been modified
}

int vtkAMRStructuredPointsReader::SelectTimeStep(int timestep){
  if (timestep<mintime || timestep>maxtime){
	vtkErrorMacro(<<"Timestep out of range");
	return 0;
  }
  myreader->setTime(timestep);
  curtime=timestep;
  maxset=myreader->getNumGrids();
  Modified();
}

int vtkAMRStructuredPointsReader::SelectLevel(int level){
  if (level<0 || level>maxlevel){
	vtkErrorMacro(<<"Level out of range");
	return 0;
  }
  curlevel=level;
  int ii=0;
  for (ii;ii<=maxlevel;ii++) myreader->hideLevel(ii);
  myreader->showLevel(curlevel);
  SelectTimeStep(curtime);
  Modified();
}

int vtkAMRStructuredPointsReader::SelectDataset(int dataset){
  if (dataset<0 || dataset>=maxset){
    vtkErrorMacro(<<"Dataset out of range");
    return 0;
  }
  curset=dataset;
  index=file->seek(myreader->getActiveIndex(curset));
  Modified();
}



void vtkAMRStructuredPointsReader::Execute(){

  //fprintf(stderr,"in execute\n");
  //puts("vtkAMR:execute()");
    int numPts=0,i;
    char line[256];
    int dimsRead=0;
#ifdef VTK31
    vtkStructuredPoints *output=(vtkStructuredPoints *)this->Outputs[0];
#else
    vtkStructuredPoints *output=(vtkStructuredPoints *)this->Output;
#endif
    
    if(!file || !(file->isValid())) {
      vtkErrorMacro(<<"Invalid Data File (cannot read)!");
      return;
    }
    // seek to current 
    file->seek(index);
    // read dimensions
    int rank,dims[3]={1,1,1};
    IObase::DataType datatype;
    i=file->readInfo(datatype,rank,dims,3);
    //printf("ReadInfo: type=%u rank=%u dims[0]=%u dims[1]=%u dims[2]=%u retval=%u\n",(int)datatype,rank,dims[0],dims[1],dims[2],i);
    //for(i=0,numPts=1;i<rank;i++) numPts*=dims[i];
    numPts = IObase::nElements(rank,dims);
    output->SetDimensions(dims);
    // danger if file is >3 dimensions!!!
    
    // read x,y, and z coordinates
    // for now, lets just set them to to the dimensions of the dataset
    // ReadCoordinates(output,<dimnumber>,ncoords)
    // ReadPointData(output,npoints)
    // for now, we only support two different types
    { // push stack
      void *ptr;
      vtkScalars *scalars = vtkScalars::New();
      scalars->Allocate(numPts*8); // bogus?
      
      switch(datatype){
      case IObase::Float32:
	{
	  scalars->SetDataTypeToFloat();
	  vtkFloatArray *farray = (vtkFloatArray*)scalars->GetData();
	  ptr = (void*)farray->WritePointer(0, numPts);
	}
	break;
      case IObase::Char:
      case IObase::uChar:
	{
	  scalars->SetDataTypeToUnsignedChar();
	  vtkUnsignedCharArray *ucarray = 
	    (vtkUnsignedCharArray*)scalars->GetData();
	  ptr = (void*)ucarray->WritePointer(0, numPts);
	}
	break;
      case IObase::Int16: 
	{
	  scalars->SetDataTypeToShort();
	  vtkShortArray *sarray = (vtkShortArray*)scalars->GetData();
	  ptr = (void*)sarray->WritePointer(0, numPts);
	}
	break;
      case IObase::uInt16:
	{
	  scalars->SetDataTypeToUnsignedShort();
	  vtkUnsignedShortArray *usarray = 
	    (vtkUnsignedShortArray*)scalars->GetData();
	  ptr = (void*)usarray->WritePointer(0, numPts);
	}
	break;
      case IObase::Int32:
	{
	  scalars->SetDataTypeToInt();
	  vtkIntArray *iarray = (vtkIntArray*)scalars->GetData();
	  ptr = (void*)iarray->WritePointer(0, numPts);
	}
	break;
      }
      
      file->read(ptr); // read it directly into vtk pre-allocated buffer
      
      // in case the data is logged.  -BPM
      //fprintf(stderr,"in execute\n");
      if(logged){
	double scalar;
	int n = scalars->GetNumberOfScalars();
	float maxscalar = scalars->GetRange()[1];
	
	int i=0;
	for(i;i<n;i++){
	  scalar = scalars->GetScalar(i);
	  
	  if (scalar > 0.0) {
	    scalar = log10((double)scalars->GetScalar(i));
	    scalars->InsertScalar(i, (float) scalar);
	  }
	  else {
	    if (scalar == 0.0) {
	      scalar = log10((double) (maxscalar*1.0e-10));
	      scalars->InsertScalar(i, (float) scalar);
	    }
	    else {
	      logged=-1;
	      vtkErrorMacro(<<"Negative value!");
	      return;
	    }
	  }
	}
	
	scalars->Modified();
	float* scarange = scalars->GetRange();
	if(scarange[0]<loggedmin)
	  loggedmin=scarange[0];
	if(scarange[1]>loggedmax)
	  loggedmax=scarange[1];
      }
      //fprintf(stderr,"in execute\n");
      //    printf("read %u points\n",numPts);
      
      output->GetPointData()->SetScalars(scalars); // should set ->Modified() and Refcount()
      scalars->Delete(); // remove local copy of reference counted data
    }
    { // push the stack...
      // check for coord info (read for xcoordinate)
      float min_ext[5],max_ext[5];
      int use_extents=0;
      int nelements;
      int minindex,maxindex;
      IObase::DataType dtype;
      
      if((minindex=file->readAttributeInfo("min_ext",dtype,nelements))>0 &&
	 (maxindex=file->readAttributeInfo("max_ext",dtype,nelements))>0 &&
	 (dtype==IObase::Float32 || dtype==IObase::Float64)) {
	if(dtype==IObase::Float32){
	  file->readAttribute(maxindex,min_ext);
	  file->readAttribute(minindex,max_ext);
	}
	else if(dtype==IObase::Float64){
	  double maxext[5],minext[5];
	  file->readAttribute(maxindex,maxext);
	  file->readAttribute(minindex,minext);
	  int i=0;
	  for(i;i<rank;i++){
	    min_ext[i]=(float)(minext[i]);
	    max_ext[i]=(float)(maxext[i]);
	  }
	}
	output->SetOrigin(min_ext);
	int i=0;
	for(i;i<rank;i++){
	  max_ext[i]-=min_ext[i];
	  if(dims[i]-1 > 0)
	    max_ext[i]/=(float)(dims[i]-1);
	}
	output->SetSpacing(max_ext);
	//use_extents=1;
      } //***** OK, if we cant find extents, then try origin & delta
      else if((minindex=file->readAttributeInfo("origin",dtype,nelements))>0 &&
	      (maxindex=file->readAttributeInfo("delta",dtype,nelements))>0 &&
	      (dtype==IObase::Float32 || dtype==IObase::Float64))  {
	if(dtype==IObase::Float32){
	  file->readAttribute(maxindex,min_ext);
	  file->readAttribute(minindex,max_ext);
	  output->SetOrigin(min_ext);
	  output->SetSpacing(max_ext);
	}
	else if(dtype==IObase::Float64){
	  double maxext[5],minext[5];
	  file->readAttribute(maxindex,maxext);
	  file->readAttribute(minindex,minext);
	  int i=0;
	  for(i;i<rank;i++){
	    min_ext[i]=(float)(minext[i]);
	    max_ext[i]=(float)(maxext[i]);
	  }
	  output->SetOrigin(min_ext);
	  output->SetSpacing(max_ext);
	}
      }
      else{
	int i=0;
	for(i;i<rank;i++){
	  min_ext[i]=0;
	  max_ext[i]=(float)(dims[i]-1);
	}
	output->SetOrigin(min_ext);
	output->SetSpacing(max_ext);
      }
    }
}

void vtkAMRStructuredPointsReader::PrintSelf(ostream& os, vtkIndent indent)
{
  os << indent << "AMR File Name: " << (this->filename) << "\n";
  os << indent << "CurrentIndex:  " << (this->index) << "\n";
  //output->PrintSelf(os,indent);
  vtkStructuredPointsSource::PrintSelf(os,indent);
}