diff options
Diffstat (limited to 'src/AVSreadcpp.cc')
| -rw-r--r-- | src/AVSreadcpp.cc | 216 |
1 files changed, 216 insertions, 0 deletions
diff --git a/src/AVSreadcpp.cc b/src/AVSreadcpp.cc new file mode 100644 index 0000000..eb5f139 --- /dev/null +++ b/src/AVSreadcpp.cc @@ -0,0 +1,216 @@ +#include <stdio.h> +#include <avs/avs.h> +#include <avs/port.h> +#include <avs/field.h> +/* John's Includes */ +#include <stdio.h> +#include <stdlib.h> +#include <math.h> +#include <string.h> +/* Local Includes */ +#include <IO.hh> +#include <IEEEIO.hh> +#ifdef WITH_HDF +#include <HDFIO.hh> +#endif +// #include <Reader.hh> // an even higher level interface, but must skip for now + +extern "C" { +int readIEEE(AVSfield **out,char *Filename,int RecordNumber, + int SetExtents,char *SetCoordinates); +} + +int readIEEE(AVSfield **out,char *Filename,int RecordNumber, + int SetExtents,char *SetCoordinates){ + static IObase *infile=0; + // Handle case where Filename has changed + if(AVSparameter_changed("Filename") || !infile){ + // Filename has changed, so close any currently open data files + // (this is done by delete'ing here.) + if(infile) delete infile; + infile=0; + // Just in case the new file is bogus, hide the dial widgets so + // that they can't be messed with if the data isn't valid. + AVSparameter_visible("RecordNumber",0); +#ifdef WITH_HDF + if(strstr(strchr(Filename,','),"hdf")) + infile = new HDFIO(Filename,IObase::Read); + else +#endif + infile = new IEEEIO(Filename,IObase::Read); // open the datafile + + if(!infile->isValid()){ // if the file is not IEEEIO, exit + delete infile; + infile=0; + fprintf(stderr,"%s is not an IEEEIO file",Filename); + return 0; + } + if(infile->nDatasets()>1){// don't set widget dial for min=max=0 + AVSmodify_parameter("RecordNumber",AVS_VALUE|AVS_MAXVAL, + 0,0,infile->nDatasets()-1); + AVSparameter_visible("RecordNumber",1); + } // so dial widget remains hidden if dataset contains 1 or fewer datasets + RecordNumber=0; + } + // Failsafe sanity check (things are really hosed if this check fails) + if(!infile) { + puts("Not an IEEEIO file"); + AVSparameter_visible("RecordNumber",0); + return 0; + } + + // now do memory allocation for AVS + char AllocationDescription[64]; + char *typname; + IObase::DataType datatype; + int rank,dims[5]; + int nelements; + + infile->seek(RecordNumber); // seek to desired record + infile->readInfo(datatype,rank,dims); // get info about the dataset + // we are building the AVS memory description text string to do the allocation. + switch(datatype){ + case IObase::Float32: + puts("Data Type float"); + typname="float"; + break; + case IObase::Float64: + puts("Data Type double"); + typname="double"; + break; + case IObase::Int32: + typname="integer"; + break; + case IObase::uChar: + case IObase::Int8: + typname="byte"; + break; + } + + if(!strcmp("Set Edges",SetCoordinates)) + sprintf(AllocationDescription,"field %uD %u-space scalar rectilinear %s", + rank,rank,typname); // put it all together into the allocation string + else + sprintf(AllocationDescription,"field %uD %u-space scalar uniform %s", + rank,rank,typname); // put it all together into the allocation string + + //printf("AVS allocating %s\n",AllocationDescription); + //printf("With dims[0]=%d dims[1]=%d dims[2]=%d\n",dims[0],dims[1],dims[2]); + if(*out) AVSfield_free(*out); // free old memory if its still around + *out = (AVSfield *)AVSdata_alloc(AllocationDescription,dims); // and allocate! + if(!*out){ + fprintf(stderr,"Allocation of %s failed\n", + AllocationDescription); + return 0; // allocation failed + } + + /*--- + Note: So if this read doesn't work, you'll need to use a switch statement + with exact casts to the proper datatypes. + I skipped using the switch statement because it looks really stupid + but its necessary for some machines which align vectors + differently depending on the size of the type the vector + points to. (yes folks... this is what the "align" keyword is about) + Hopefully we won't use such poorly designed machines. + + This need not be done for the parameters though since that field + of the avsfld is always float. + ---*/ + AVSfield_float *fld=(AVSfield_float*)(*out); + puts("Preparing to read data"); + infile->read(fld->data); // read the data + puts("read data"); + int minindex,maxindex; + /*--- + In this section we are just setting the extents and edge points information + for AVS (the physical size of the dataset in floating point coordinates). + There are just two Attribute pairs that it will recognize to do this. + + 1) min_ext,max_ext + These are vectors of size <rank of dataset> which contain the + coordinates of the minimum and maximum corners of the bounding box. + 2) origin,delta + The GR people seem to like this. Its the origin (same as the min_ext) + of the dataset and then the delta (or "dx") between grid points in the + dataset. + + It is now able to handle these parameters as Float32 or Float64 data. + Integer is ignored. The dataset will simply be ranged from 0-dim[x] + if the above Attributes are not found. + ---*/ + { + float min_ext[5],max_ext[5]; + + if((minindex=infile->readAttributeInfo("min_ext",datatype,nelements))>0 && + (maxindex=infile->readAttributeInfo("max_ext",datatype,nelements))>0 && + (datatype==IObase::Float32 || datatype==IObase::Float64)) { + if(datatype==IObase::Float32){ + infile->readAttribute(maxindex,min_ext); + infile->readAttribute(minindex,max_ext); + } + else if(datatype==IObase::Float64){ + double maxext[5],minext[5]; + infile->readAttribute(maxindex,maxext); + infile->readAttribute(minindex,minext); + for(int i=0;i<rank;i++){ + min_ext[i]=(float)(minext[i]); + max_ext[i]=(float)(maxext[i]); + } + } + + } //***** OK, if we cant find extents, then try origin & delta + else if((minindex=infile->readAttributeInfo("origin",datatype,nelements))>0 && + (maxindex=infile->readAttributeInfo("delta",datatype,nelements))>0 && + (datatype==IObase::Float32 || datatype==IObase::Float64)) { + if(datatype==IObase::Float32){ + infile->readAttribute(maxindex,fld->min_extent); + infile->readAttribute(minindex,fld->max_extent); + for(int i=0;i<rank;i++){ + max_ext[i]*=(float)dims[i]; + max_ext[i]+=min_ext[i]; + } + } + else if(datatype==IObase::Float64){ + double maxext[5],minext[5]; + infile->readAttribute(maxindex,maxext); + infile->readAttribute(minindex,minext); + for(int i=0;i<rank;i++){ + maxext[i]*=(double)dims[i]; + maxext[i]+=minext[i]; + min_ext[i]=(float)(minext[i]); + max_ext[i]=(float)(maxext[i]); + } + } + } + else { // No information about edges. + for(int i=0;i<rank;i++){ + min_ext[i]=0.0; + max_ext[i]=(float)(dims[i]-1); + } + } + if(SetExtents){ + for(int i=0;i<rank;i++){ + fld->min_extent[i]=min_ext[i]; + fld->max_extent[i]=max_ext[i]; + } + } + if(!strcmp("Set Min/Max",SetCoordinates)){ + for(int i=0;i<rank;i++){ + fld->points[i*2]=min_ext[i]; + fld->points[i*2 + 1]=max_ext[i]; + } + } + else if(!strcmp("Set Edges",SetCoordinates)){ + for(int index=0,i=0;i<rank;i++){ + double dx = (max_ext[i]-min_ext[i])/((float)(fld->dimensions[i])); + printf("SetEdges[%u]: dx=%lf\n",i,dx); + for(int j=0;j<fld->dimensions[i];j++,index++){ + fld->points[index] = min_ext[i] + ((double)j)*dx; + } + } + } + } + return 1; // for now +} + + |