FITS file processor

fits.c

#include <stdio.h>
#include <fitsio.h>
#include <dirent.h>
#include <string.h>
#include <mpi.h>

/*If status != 0, report the error and return 1.*/
#define FITS_ERROR_CHECK if(status){fits_report_error(stdout, status); return 1;}

int countNumFits(int limit, char * dir);
char ** allocateStrArr(int size);
char ** populateStrArr(int size, char * dir, char ** strArr);
long * allocateLongArr(int size);
int * allocateIntArr(int size);

int main (int argc, char ** argv) 
{
  if(argc != 4 )
  {
    printf("usage: %s <input-directory> <output-directory> <max-num-files>\n",argv[0]);
    return 0;
  }
  
  /*For fits processing*/
  int fitsLimit = atoi(argv[3]);
  int numFits, index;
  
  FILE *filePtr = NULL;      //pointers to file
  int   sumPix;          //stores sum of connected pixels
  int  valueThresh = 1150;    //threshold value of pix we want
  int   clusterThresh = 20;    //threshold number of pix we want connected  
  int *valueArrPtr = NULL;    //traverses valueArr
  int *clusterArrPtr= NULL;    //traverses clusterArr
  long rowCtr, colCtr, k, i;    //loop counters

  /*For cfitsio*/
  fitsfile *fptr = NULL;       //pointer to .fit
  fitsfile *newFptr = NULL;    //pointer to new .fit
  int status = 0;          //must be set to 0
  int hdunum;            //HDU number
  int hdutype;          //type of HDU - Image, ASCII, or Binary
  int numKeys;          //number of keywords in header
  int moreKeys;          //amount of empty space for keywords in header
  int keyLength = 81;        //80-characters for key length
  char keyword[keyLength];    //Arr for keywords in header
  int keyNum;            //Keyword number in header
  int maxDimension = 2;      //maximum number of image dimensions allowed
  int bitPix;            //bit size of each pixel in image
  int nAxes;            //number of axes in image
  long nAxis[2];          //size of each axis  in image
  long numPix;          //number of pixels in image
  long firstPixel[2];        //coor of first pixel (lower left) in image
  long lastPixel[2];        //coor of the last pixel (upper right) in image
  long increment[2];        //increment - reads every inc-th pixel of image
  long longNull;          //for fits_read_subset
    
  /*For mpi */
  int rank;
   int numProcesses;
  int fitsPerProcess;
  int fitsPerProcessRemainder;
  int localIndex;
  int localIndexStart;
  int localIndexEnd;
  double startTime, endTime, loopTime1, loopTime2;

  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &numProcesses);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  
  /*Count the number of fits files in directory provided to argv[1]*/
  numFits = countNumFits(fitsLimit, argv[1]);
  
  /*Allocate and populate array to store fits files*/
  char ** fitsNameArr = allocateStrArr(numFits);
  
  fitsNameArr = populateStrArr(numFits, argv[1], fitsNameArr);
  
  /*Determines number of .fit files each processor will handle*/
  fitsPerProcess = numFits/ numProcesses;
  //fitsPerProcessRemainder = numFits % numProcesses;
  
  /*Determines the beginning and end indices for each processor based on rank*/
  localIndexStart = fitsPerProcess * rank;
  localIndexEnd = localIndexStart + fitsPerProcess;
  
  fprintf(stderr,"%d: [%d,%d)\n",rank,localIndexStart,localIndexEnd);
  loopTime1 = 0.0;
  loopTime2 = 0.0;
  /*Begins processing .fit files*/
  for (localIndex = localIndexStart; localIndex < localIndexEnd; ++localIndex) 
  {
    startTime = MPI_Wtime();
    /*Change into directory provided in argv[1]*/
    chdir(argv[1]);
    /*open fits files argv[1] directory*/
    fits_open_file(&fptr, fitsNameArr[localIndex], READONLY, &status);
    FITS_ERROR_CHECK 
    /*Return the number of exisiting keywords*/
    fits_get_hdrspace(fptr, &numKeys, &moreKeys, &status);
    FITS_ERROR_CHECK
    /*Return values for maxDimension, bitpix, nAxes, nAxis[] */
    fits_get_img_param(fptr, maxDimension, &bitPix, &nAxes, nAxis, &status);
    FITS_ERROR_CHECK
    numPix = nAxis[0]*nAxis[1];    
    firstPixel[0] = 1;        
    firstPixel[1] =  1;
    lastPixel[0] = nAxis[0];    
    lastPixel[1] = nAxis[1];
    increment[0] = 1;  
    increment[1] = 1;
    longNull = 0;          
    /*Allocates and populates imageArray by reading a subset from fits image*/
    long * imageArr = allocateLongArr(numPix);
    fits_read_subset(fptr, TLONG, firstPixel, lastPixel, increment, 
              &longNull, imageArr, NULL, &status);
    FITS_ERROR_CHECK;
    /*Allocates and populates valueArr for pixel value filtering.
    Pixels above valueThresh are assigned 1, and pixels below valueThresh are assigned 0*/
    int * valueArr = allocateIntArr(numPix);
    for(i = 0; i < numPix; ++i)
    {
      if(imageArr[i] >= valueThresh)
      {
        valueArr[i] = 1;
      }
    }  
    /*Allocates clusterArr for pixel cluster filtering*/
    int * clusterArr = allocateIntArr(numPix);
    /*Pointers to traverse valueArr and clusterArr for populating clusterArr*/
    valueArrPtr = valueArr;          
    clusterArrPtr = clusterArr;      
    /*Staring with the bottom left pixel, traverse each row to detect
    the clusterThresh number of connected pixels by summing their values (0 or 1)
    and testing the sum. If sumPix = valueThresh, write a 1 to the corresponding 
    indices in clusterArr. Subtract the first pixel and add another. Retest and write.
    Continue until the last column.  */
    int numRows = nAxis[1];
    int numCol = nAxis[0];
    sumPix = 0;
    for(rowCtr = 0; rowCtr < numRows ; ++rowCtr) 
    {
      /*Sum the first (clusterThresh - 1) pixels in the first row.*/
      for(colCtr = 0; colCtr < (clusterThresh - 1); ++colCtr) 
      {
          sumPix += valueArrPtr[colCtr];
      }
      
      for(colCtr = (clusterThresh - 1); colCtr < numCol; ++colCtr) 
      {
        /*Add the next pixel in the clusterThresh and test sumPix */
        sumPix += valueArrPtr[colCtr];  
        if(sumPix == clusterThresh)          
        {
          /*Write a 1 corresponding indices in clusterArr*/
          for(k = colCtr; k > colCtr - clusterThresh; --k) 
          {
            clusterArrPtr[k] = 1;
          }
        }
        /*Subtract the first pixel*/
        sumPix -= valueArrPtr[colCtr - clusterThresh + 1];   
      }
      /*Add a row to the starting point for valueArrPtr and clusterArrPtr */
      valueArrPtr += numCol;
      clusterArrPtr += numCol;
    }
    endTime = MPI_Wtime();
    loopTime1 += endTime-startTime;
    startTime = endTime;
    /*Resets valueArrPtr and clusterArrPtr */
    valueArrPtr = valueArr;
    clusterArrPtr = clusterArr;
    /*Staring with the bottom left pixel, traverse each column to detect
    the clusterThresh number of connected pixels by summing their values (0 or 1)
    and testing the sum. If sumPix = valueThresh, write a 1 to the corresponding \
    indices in clusterArr. Subtract the first pixel and add another. Retest and write.
    Continue until the last row.  */
    sumPix = 0;
    for(colCtr=0; colCtr < numCol; ++colCtr) 
    {  
      /*Sum the first (clusterThresh - 1) pixels in the first column.*/
      for(rowCtr = 0; rowCtr < (clusterThresh-1); ++rowCtr) 
      {
        sumPix += valueArrPtr[rowCtr * numCol];
      }
      
      
      for(rowCtr = (clusterThresh-1); rowCtr < nAxis[1]; ++rowCtr) 
      {
        /*Add the next pixel in clusterThresh and test sumPix */
        sumPix += valueArrPtr[rowCtr* numCol]; 
        if(sumPix == clusterThresh)
        {
          /*Write a 1 to corresponding indices in clusterArr */
          for(k = rowCtr; k > rowCtr-(clusterThresh-1); --k) 
          {
            clusterArrPtr[k*numCol] = 1;
          }
        }
        /*Subtract the first pixel*/
        sumPix -= valueArrPtr[(rowCtr - clusterThresh + 1)*numCol];      
      }
       /*Add a column to the starting point valueArrPtr and clusterArrPtr*/
       valueArrPtr += 1;
       clusterArrPtr += 1;
    }
    /*Change into FITS directory passed into argv[2]*/
    chdir(argv[2]);          
    /*Create new .fit file in the newFITS directory*/
    fits_create_file(&newFptr, fitsNameArr[localIndex], &status);
    FITS_ERROR_CHECK
    /*copies .fit file to new .fit file*/
    fits_copy_file(fptr, newFptr, 0, 1, 0, &status);
    FITS_ERROR_CHECK
    /*closes new .fit file*/
    fits_close_file(newFptr, &status);
    FITS_ERROR_CHECK
    /*closes original .fit file*/
    fits_close_file(fptr, &status);
    FITS_ERROR_CHECK
    free(imageArr);
    free(valueArr);
    free(clusterArr);
    endTime = MPI_Wtime();
    loopTime2 += endTime - startTime;
  }
  fprintf(stderr,"%d: loopTime1: %f loopTime2: %f\n",rank,loopTime1,loopTime2);
  for(index = 0; index < numFits; ++index)
  {
    free(fitsNameArr[index]);
  }
  free(fitsNameArr);
  MPI_Finalize();  
  return 0;
}

int countNumFits(int limit, char * dir)  
{
  DIR *tempDirPtr = NULL;          //directory pointer
  struct dirent *tempEntryPtr = NULL;    //entry pointer
  int tempCount;
  
  tempDirPtr = opendir (dir);
  if(tempDirPtr != NULL)
  {
    tempCount = 0;
    while((tempEntryPtr = readdir(tempDirPtr)) && (tempCount < limit))
    {
      if(tempEntryPtr->d_type != DT_REG)
        continue;
      tempCount++;
    }
    (void) closedir (tempDirPtr);
  }
  else
  {
    printf("Directory %s is empty \n", dir);
    exit(EXIT_FAILURE);
    }  
    return tempCount;
}

char ** allocateStrArr(int size)
{
  int tempIndex;
  char ** tempCharArr = (char ** )malloc(sizeof(char * ) *size);
  for (tempIndex = 0; tempIndex < size; ++tempIndex)
    tempCharArr[tempIndex] = NULL;
  return tempCharArr;
}

char ** populateStrArr(int size, char * dir, char ** strArr)
{
  DIR *tempDirPtr = NULL;          //directory pointer
  struct dirent *tempEntryPtr = NULL;    //entry pointer
  int tempSize, tempIndex;
  
  tempDirPtr = opendir (dir);  
   if (tempDirPtr != NULL)      
  {
    tempIndex = 0;
    while((tempEntryPtr = readdir(tempDirPtr)) && (tempIndex < size))
    {
      if(tempEntryPtr->d_type != DT_REG)
        continue;
      strArr[tempIndex] = (char*)malloc(tempEntryPtr->d_reclen + 1);
      strcpy(strArr[tempIndex], tempEntryPtr->d_name);
      tempIndex++;
    }
    (void) closedir (tempDirPtr);  
  }
  else
  {
    printf("Directory %s is empty \n", dir);
    exit(EXIT_FAILURE);
  }  
return strArr;
}

long * allocateLongArr(int size)
{
  int tempIndex;
  long * tempLongArr = (long *)malloc(sizeof(long)*size);
  for(tempIndex = 0; tempIndex < size; ++tempIndex)
    tempLongArr[tempIndex] = 0;
  return tempLongArr;
}  

int * allocateIntArr(int size)
{
  int tempIndex;
  int * tempIntArr = (int *)malloc(sizeof(int)*size);
  for(tempIndex = 0; tempIndex < size; ++tempIndex)
    tempIntArr[tempIndex] = 0;
  return tempIntArr;
}