Handle files with stokes and channel axes switched (addresses issue #26)

src/Converter.cc

@@ -5,22 +5,31 @@
 
 #include "Converter.h"
 
-Converter::Converter(std::string inputFileName, std::string outputFileName, bool progress, bool zMips) : timer(), progress(progress), zMips(zMips) {
+Converter::Converter(std::string inputFileName, std::string outputFileName, bool progress, bool zMips) : timer(), progress(progress), zMips(zMips), swapStokesFreqAxis(false) {
     TIMER(timer.start("Setup"););
 
     openFitsFile(&inputFilePtr, inputFileName);
 
     long dims[4];
 
     getFitsDims(inputFilePtr, N, dims);
+
+    // check relevant FITS keywords if the axis order is STOKES,FREQ -> swap is required:
+    checkIfSwapAxisRequired();
 
     stokes = N == 4 ? dims[3] : 1;
     depth = N >= 3 ? dims[2] : 1;
+    if ( swapStokesFreqAxis ) {
+       stokes = N >= 4 ? dims[2] : 1;
+       depth = N == 4 ? dims[3] : 1;
+    }
     height = dims[1];
     width = dims[0];
 
     swizzledName = N == 3 ? "ZYX" : "ZYXW";
 
+    DEBUG(std::cout << "Stokes = " << stokes << ", depth = " << depth << ", image dimensions " << height << " x " << width << " , swizzledName = " << swizzledName.c_str(););
+
     standardDims = trimAxes({stokes, depth, height, width}, N);
     tileDims = trimAxes({1, 1, TILE_SIZE, TILE_SIZE}, N);
 
@@ -72,6 +81,39 @@ void Converter::reportMemoryUsage() {
     std::cout << "TOTAL:\t" << m.total * 1e-9 << "GB" << m.note << std::endl;
 }
 
+bool Converter::checkIfSwapAxisRequired() {
+    int numAttributes;
+    readFitsHeader(inputFilePtr, numAttributes);
+
+    bool bFreqAxisFound = false;
+
+    swapStokesFreqAxis = false;
+    // Check if both STOKES and FREQ axis are present and the order is STOKES,FREQ:
+    for (int i = 1; i <= numAttributes; i++) {
+       std::string attributeName;
+       std::string attributeValue;
+       readFitsAttribute(inputFilePtr, i, attributeName, attributeValue);
+       DEBUG(std::cout << "|" << attributeName.c_str() << "| = |" << attributeValue.c_str() << "|" << std::endl;);
+
+       if (attributeName == "CTYPE3" && attributeValue.find("STOKES") != std::string::npos) {
+          std::cout << "INFO : detected 3rd axis CTYPE3 = STOKES" << std::endl;
+          swapStokesFreqAxis = true;
+       }
+       if (attributeName == "CTYPE4" && attributeValue.find("FREQ") != std::string::npos) {
+          std::cout << "INFO : detected 4th axis CTYPE4 = FREQ" << std::endl;
+          bFreqAxisFound = true;
+       }
+    }
+    if (swapStokesFreqAxis) {
+       if (!bFreqAxisFound) {
+          swapStokesFreqAxis = false; // If Stokes is 3rd axis, but 4th axis is not frequency -> no swap
+       }
+    }
+    std::cout << "INFO : swapStokesFreqAxis = " << swapStokesFreqAxis << std::endl;
+
+    return swapStokesFreqAxis;
+}
+
 void Converter::convert() {
     // CREATE OUTPUT FILE
 
@@ -114,11 +156,36 @@ void Converter::convert() {
     int numAttributes;
     readFitsHeader(inputFilePtr, numAttributes);
 
+    std::vector<std::string> wcsStokesKeywords = {"CTYPE3","CRVAL3","CDELT3","CRPIX3","CUNIT3"};
+    std::vector<std::string> wcsFreqKeywords = {"CTYPE4","CRVAL4","CDELT4","CRPIX4","CUNIT4"};
+
     // IMPORTANT: This is 1-indexed!
-    for (int i = 1; i <= numAttributes; i++) {        
+    for (int i = 1; i <= numAttributes; i++) {
+        bool bSwappingStokesFreqAxisNow = false;
         std::string attributeName;
         std::string attributeValue;
         readFitsAttribute(inputFilePtr, i, attributeName, attributeValue);
+
+        if (swapStokesFreqAxis) { // check if swap of STOKES <-> FREQ axis is required:
+           // If Stokes is 3rd axis and Freq 4th axis -> SWAP Stokes and Frequency axis in FITS header
+           if (std::find(wcsStokesKeywords.begin(), wcsStokesKeywords.end(), attributeName) != wcsStokesKeywords.end()) {
+              if(attributeName.back() == '3') {
+                 // if swap of STOKES and FREQ axis is required change 3 -> 4:
+                 std::cout << "INFO : swapping attribute " << attributeName << " (3 -> 4) " << std::endl;
+                 attributeName.back() = '4';
+                 bSwappingStokesFreqAxisNow = true;
+              }
+           }else{
+              if ( std::find(wcsFreqKeywords.begin(), wcsFreqKeywords.end(), attributeName) != wcsFreqKeywords.end() ){
+                 if(attributeName.back() == '4') {
+                    // if swap of STOKES and FREQ axis is required change 4 -> 3:
+                    std::cout << "INFO : swapping attribute " << attributeName << " (4 -> 3) " << std::endl;
+                    attributeName.back() = '3';
+                    bSwappingStokesFreqAxisNow = true;
+                 }
+              }
+           }
+        }
 
         if (attributeName.empty() || attributeName.find("COMMENT") == 0 || attributeName.find("HISTORY") == 0) {
             // TODO we should actually do something about these
@@ -133,7 +200,14 @@ void Converter::convert() {
                     // STRING
                     std::string attributeValueStr;
                     readFitsStringAttribute(inputFilePtr, attributeName, attributeValueStr);
+                    // this may not be required 
+                    if( bSwappingStokesFreqAxisNow ) {
+                       // if swapping axis now, we have to use attributeValue as the attributes's name has been swapped 
+                       // for example CTYPE3 -> CTYPE4 which means value of CTYPE4 would be read (FREQ) and no swap would happen!
+                       attributeValueStr = attributeValue;
+                    }
                     writeHdf5Attribute(outputGroup, attributeName, attributeValueStr);
+                    DEBUG(std::cout << "Written attribute " << attributeName.c_str() << " = " << attributeValueStr.c_str() << std::endl;);
                 } else if (attributeValue == "T" || attributeValue == "F") {
                     // BOOLEAN
                     bool attributeValueBool = (attributeValue == "T");

src/Converter.h

@@ -30,6 +30,9 @@ class Converter {
     void convert();
     void reportMemoryUsage();
     virtual MemoryUsage calculateMemoryUsage() = 0;
+
+    // checks the order of STOKES and FREQUENCY axis an if it is STOKES,FREQ sets flat swapStokesFreqAxis to true:
+    bool checkIfSwapAxisRequired();
 
 protected:
     virtual void copyAndCalculate() = 0;
@@ -41,6 +44,7 @@ class Converter {
     std::string tempOutputFileName;
     std::string outputFileName;
     fitsfile* inputFilePtr;
+    bool      swapStokesFreqAxis;
 
     // Main HDF5 objects
     H5::H5File outputFile;

src/FastConverter.cc

@@ -61,7 +61,7 @@ void FastConverter::copyAndCalculate() {
         // Read data into memory space
         TIMER(timer.start("Read"););
         DEBUG(std::cout << "+ Reading main dataset..." << std::flush;);
-        readFitsData(inputFilePtr, 0, currentStokes, cubeSize, standardCube);
+        readFitsData(inputFilePtr, 0, currentStokes, cubeSize, standardCube, swapStokesFreqAxis);
 
         // We have to allocate the swizzled cube for each stokes because we free it to make room for mipmaps
         if (depth > 1) {

src/SlowConverter.cc

@@ -71,7 +71,7 @@ void SlowConverter::copyAndCalculate() {
             DEBUG(std::cout << "+ Processing channel " << c << "... " << std::flush;);
             DEBUG(std::cout << " Reading main dataset..." << std::flush;);
             TIMER(timer.start("Read"););
-            readFitsData(inputFilePtr, c, s, cubeSize, standardCube);
+            readFitsData(inputFilePtr, c, s, cubeSize, standardCube, swapStokesFreqAxis);
 
             // Write the standard dataset
 
@@ -223,7 +223,7 @@ void SlowConverter::copyAndCalculate() {
             DEBUG(std::cout << " Reading main dataset..." << std::flush;);
             TIMER(timer.start("Read"););
 
-            readFitsData(inputFilePtr, c, s, cubeSize, standardCube);
+            readFitsData(inputFilePtr, c, s, cubeSize, standardCube, swapStokesFreqAxis);
 
             DEBUG(std::cout << " Calculating histogram(s)..." << std::endl;);
             TIMER(timer.start("Histograms"););

src/Util.cc

@@ -150,8 +150,13 @@ void readFitsStringAttribute(fitsfile* filePtr, const std::string& name, std::st
     value = strValueTmp;
 }
 
-void readFitsData(fitsfile* filePtr, hsize_t channel, unsigned int stokes, hsize_t size, float* destination) {
+void readFitsData(fitsfile* filePtr, hsize_t channel, unsigned int stokes, hsize_t size, float* destination, bool bSwapStokesFreqAxis) {
     long fpixel[] = {1, 1, (long)channel + 1, stokes + 1};
+    if (bSwapStokesFreqAxis) {
+       // basically should instead be : long fpixel[] = {1, 1, stokes + 1, (long)channel + 1};
+       fpixel[2] = stokes + 1;
+       fpixel[3] = (long)channel + 1;
+    }
     int status(0);
 
     fits_read_pix(filePtr, TFLOAT, fpixel, size, NULL, destination, NULL, &status);

src/Util.h

@@ -32,7 +32,7 @@ void getFitsDims(fitsfile* filePtr, int& N, long* dims);
 void readFitsHeader(fitsfile* filePtr, int& numAttributes);
 void readFitsAttribute(fitsfile* filePtr, int i, std::string& name, std::string& value);
 void readFitsStringAttribute(fitsfile* filePtr, const std::string& name, std::string& value);
-void readFitsData(fitsfile* filePtr, hsize_t channel, unsigned int stokes, hsize_t size, float* destination);
+void readFitsData(fitsfile* filePtr, hsize_t channel, unsigned int stokes, hsize_t size, float* destination, bool bSwapStokesFreqAxis);
 
 // Only available in C++ API from 1.10.1
 bool hdf5Exists(H5::H5Location& location, const std::string& name);