M3_format_healpix_wtf.c

#include "M3_format_healpix.h"
#include "fitsio.h"
#include "ccFitsTools.h"


void M3_File_parseFormatString_healpix( M3_File *theFile, M3_FormatStruct_healpix *formatStruct )
{
  char *sptr;


  formatStruct->order = 'r';
  formatStruct->nside = 0;
  formatStruct->fieldIndex = 0;

  sptr = strchr( theFile->format, '_');
  if( sptr)
    sptr++;

  while( sptr )
  {
    if( strncmp( sptr, "nside=", 6 ) == 0 )
    {
      sptr += 6;
      formatStruct->nside = atoi(sptr);
    }
    else if( strncmp( sptr, "fieldIndex=", 11 ) == 0 )
    {
      sptr += 11;
      formatStruct->fieldIndex = atoi(sptr);
    }
    else if( strncmp( sptr, "order=", 6 ) == 0 )
    {
      sptr += 6;
      if( strncmp( sptr, "ring", 4 ) == 0 )
        formatStruct->order = 'r';
      else if( strncmp( sptr, "nest", 4) == 0 )
        formatStruct->order = 'n';
    }
    sptr = strchr( sptr, '_');
    if(sptr) 
      sptr++;
  }
}

void M3_File_readHeader_healpix( M3_File *theFile, void *data )
{
  switch( theFile->fileType )
  {
    case M3_COORD_FILE_TYPE:
      M3_File_readHeader_healpix_coord( theFile, (M3_CoordHeader *)data );
      break;
    case M3_MAP_FILE_TYPE:
      M3_File_readHeader_healpix_map( theFile, (M3_MapHeader *)data );
      break;
    default:
      M3_ErrorCheck(-1, "healpix file format only valid for coord and map files\n", 0, M3_EFLAG_DEFAULT);
  }

}
void M3_File_readData_healpix( M3_File *theFile, void *readArg, void *data )
{
  switch( theFile->fileType )
  {
    case M3_COORD_FILE_TYPE:
      M3_File_readData_healpix_coord( theFile, *((int32_t*)readArg), (M3_CoordEl *)data );
      break;
    case M3_MAP_FILE_TYPE:
      M3_File_readData_healpix_map( theFile, *((int32_t*)readArg), (M3_MapEl *)data );
      break;
    default:
      M3_ErrorCheck(-1, "healpix file format only valid for coord and map files\n", 0, M3_EFLAG_DEFAULT);
  }  
  
}

void M3_File_readHeader_healpix_coord( M3_File *theFile, M3_CoordHeader *header )
{
  M3_FormatStruct_healpix formatStruct;

  M3_File_parseFormatString_healpix( theFile, &formatStruct );
  header->numPixelInClass = 12*formatStruct.nside*formatStruct.nside;

}


void M3_File_readData_healpix_coord( M3_File *theFile, int32_t numReadPixel, M3_CoordEl *data )
{
  int32_t i;
  double theta, phi;
  int32_t nside;
  M3_FormatStruct_healpix formatStruct ;
  M3_File_parseFormatString_healpix(theFile, &formatStruct );
  
  if( formatStruct.order == 'n')
    for( i = 0; i < numReadPixel; i++ )
    {
      pix2ang_nest( formatStruct.nside, data[i].pixel, &theta, &phi);
      data[i].ra = M_1_PI*180.0*phi;
      data[i].dec = 90 - M_1_PI*180.0*theta;
    }
  else
    for( i = 0; i < numReadPixel; i++ )
    {
      pix2ang_ring( formatStruct.nside, data[i].pixel, &theta, &phi);
      data[i].ra = M_1_PI*180.0*phi;
      data[i].dec = 90 - M_1_PI*180.0*theta;
    }
}


void M3_File_readHeader_healpix_map( M3_File *theFile, M3_MapHeader *header )
{
  M3_FormatStruct_healpix formatStruct;
  fitsfile *thisFitsFile;
  int32_t nside;
  char errorString[256];
  char key[256];
  int status = 0;

  fits_open_file( &thisFitsFile, theFile->name, READONLY, &status );
  M3_ErrorCheck(-1, theFile->name, status == 0, M3_EFLAG_FOPEN_READ );

  fits_movabs_hdu( thisFitsFile, 2, NULL, &status);
  sprintf( errorString, "Could not go to first fits header in file %s", theFile->name );
  M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

  fits_read_key_str( thisFitsFile, "NSIDE", key, NULL, &status );
  sprintf( errorString, "Could not read NSIDE flag from healpix map named %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT);
  nside = atoi(key);

  fits_close_file( thisFitsFile, &status );
  sprintf( errorString, "Error closing fits file named %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
  
  header->numPixelInMap = header->numPixelInClass = 12*nside*nside;
}

typedef struct 
{
  char fileName[1024];
  char fileFormat[1024];
  int32_t lastPixel;
  fitsfile *fitsFilePtr;
} M3_staticMapFileRecord;

void M3_File_readData_healpix_map( M3_File *theFile, int32_t numReadPixel, M3_MapEl *data )
{
  fitsfile *thisFitsFile;
  int status = 0;
  int byteSwap, blocks;
  char errorString[256] = {0};
  char key[256] = {0};
  M3_FormatStruct_healpix formatStruct ;
  static int64_t diskBufferSize = -1;
  char *tempPtr;
  int32_t nside, tfields, pixelIndex, firstPixel, lastPixel;
  int64_t blockSize, numBlock, bufferLength, numReads, thisBufferLength, i, j, k, l, m;
  FILE *fid;
  float *buffer, *bufferPtr;
  M3_Interval readInterval, bufferInterval;
  double *dptr;
  static M3_staticMapFileRecord staticMapFileRecord = {0};

  if( diskBufferSize == -1 )
  {
    tempPtr = getenv("M3_DISK_BUFFER_SIZE");
    if( tempPtr )
      diskBufferSize = strtoll( tempPtr, NULL, 10);
    else
      diskBufferSize = M3_DISK_BUFFER_SIZE;
  }
  
  if( strcmp( theFile->name, staticMapFileRecord->fileName) == 0 )
  {
    thisFitsFile = staticMapFileRecord->fitsFilePtr;
    /* FIX ME */

  }  


  else
  {
    fits_open_file( &thisFitsFile, theFile->name, READONLY, &status );
    M3_ErrorCheck(-1, theFile->name, status == 0, M3_EFLAG_FOPEN_READ );

    fits_movabs_hdu( thisFitsFile, 2, NULL, &status);
    sprintf( errorString, "Could not go to first fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );
  }

  M3_File_parseFormatString_healpix(theFile, &formatStruct );

  firstPixel = data[0].pixel; 
  lastPixel = data[numReadPixel-1].pixel;
  
  if( lastPixel - firstPixel + 1 == numReadPixel  )
  {
    /* Pixels in request are consecutive so use read_fits_tableblock() */
    dptr = (double*)(data + numReadPixel) - numReadPixel;
    read_fits_tableblock( thisFitsFile, firstPixel, lastPixel, NULL, 1, 0, NULL, &formatStruct.fieldIndex, dptr);

    for( i = 0, j = firstPixel; i < numReadPixel; i++, j++ )
    {
      data[i].value = dptr[i];
      data[i].pixel = j;
    }

    fits_close_file( thisFitsFile, &status );
    sprintf( errorString, "Error closing fits file named %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );

  }
  else
  {
    fits_read_key_str( thisFitsFile, "NSIDE", key, NULL, &status );
    sprintf( errorString, "Could not read NSIDE flag from healpix map named %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT);
    nside = atoi(key);

    fits_read_key_str( thisFitsFile, "TFIELDS", key, NULL, &status );
    sprintf( errorString, "Could not read TFIELDS flag from healpix map named %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT);
    tfields = atoi(key);

    fits_read_key_str( thisFitsFile, "TFORM1", key, NULL, &status );
    sprintf( errorString, "Could not read TFORM1 flag from healpix map named %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT);
    if(strncmp( key, "1024E", 5) ==  0 )
      blocks = 1;
    else
      blocks = 0;

    fits_close_file( thisFitsFile, &status );
    sprintf( errorString, "Error closing fits file named %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );

    buffer = (float *)malloc(diskBufferSize);
    M3_ErrorCheck(-1, "M3_File_readData_healpix_map()", (buffer?1:0), M3_EFLAG_MALLOC_FSTRING );

    byteSwap = (testEndian()?0:1);
    fid = fopen( theFile->name, "r");
    skipFitsHeaders( fid, 2 );

    if( blocks == 1 )
    {
      blockSize = tfields*1024*sizeof(float);
      numBlock = 12*nside*nside/1024;
      
      readInterval.firstSample = 0; 
      readInterval.lastSample = numBlock - 1;
      
      M3_bufferIntervalStart( readInterval, diskBufferSize, blockSize, &bufferLength, &numReads );

      m = 0;
      pixelIndex = 0;
    
      for( j = 0; j < numReads; j++ )
      {
        M3_bufferIntervalStep( readInterval, bufferLength, j, &bufferInterval );
        thisBufferLength = bufferInterval.lastSample - bufferInterval.firstSample + 1;
        fread(buffer, blockSize, thisBufferLength, fid );
      
        if( byteSwap )
          flipByteOrder( buffer, thisBufferLength*tfields*1024, sizeof(float));

        bufferPtr = buffer + 1024*formatStruct.fieldIndex;

        for( k = 0; k < thisBufferLength && m < numReadPixel; k++ )
        {
          for( l = 0; l < 1024 && m < numReadPixel; l++ )
          {
            if( data[m].pixel == pixelIndex )
            {
              data[m].value = bufferPtr[l];
              m++;
            }
            pixelIndex++;
          }
          bufferPtr += tfields*1024;
        }
      }
    }
    else 
    {
      readInterval.firstSample = 0;
      readInterval.lastSample = 12*nside*nside - 1;

      M3_bufferIntervalStart( readInterval, diskBufferSize, tfields*sizeof(float), &bufferLength, &numReads );
      m = 0;
      pixelIndex = 0;
      for( j = 0; j < numReads && m < numReadPixel; j++ )
      {
        M3_bufferIntervalStep( readInterval, bufferLength, j, &bufferInterval );
        thisBufferLength = bufferInterval.lastSample - bufferInterval.firstSample + 1;
        fread(buffer, tfields*sizeof(float), thisBufferLength, fid );
      
        if( byteSwap )
          flipByteOrder( buffer, thisBufferLength*tfields, sizeof(float) );

        for( k = 0; k < thisBufferLength && m < numReadPixel; k++ )
        {
          if( data[m].pixel == pixelIndex )
          {
            data[m].value = buffer[k*tfields + formatStruct.fieldIndex];
            m++;
          }
          pixelIndex++;
        }
      }

    }
    fclose(fid);
    free(buffer);
  }
}

---