From 740a489821dc89e4b6baf379146194b7c4b5507e Mon Sep 17 00:00:00 2001 From: Marcin Sokolowski Date: Wed, 4 Feb 2026 15:01:55 +0800 Subject: [PATCH] initial revision of version where swapped order of STOKES and FREQUENCY axis is discovered (3rd axis STOKES and 4th axis FREQ) and program works appropriately to convert FITS to HDF5 file without going through interemediate step of swapaxis in python. --- src/Converter.cc | 78 ++++++++++++++++++++++++++++++++++++++++++-- src/Converter.h | 4 +++ src/FastConverter.cc | 2 +- src/SlowConverter.cc | 4 +-- src/Util.cc | 7 +++- src/Util.h | 2 +- 6 files changed, 90 insertions(+), 7 deletions(-) diff --git a/src/Converter.cc b/src/Converter.cc index b011fc0..131ebc2 100644 --- a/src/Converter.cc +++ b/src/Converter.cc @@ -5,7 +5,7 @@ #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); @@ -13,14 +13,23 @@ Converter::Converter(std::string inputFileName, std::string outputFileName, bool 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 wcsStokesKeywords = {"CTYPE3","CRVAL3","CDELT3","CRPIX3","CUNIT3"}; + std::vector 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"); diff --git a/src/Converter.h b/src/Converter.h index 78a45ad..796aada 100644 --- a/src/Converter.h +++ b/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; diff --git a/src/FastConverter.cc b/src/FastConverter.cc index 325eeeb..5730ca4 100644 --- a/src/FastConverter.cc +++ b/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) { diff --git a/src/SlowConverter.cc b/src/SlowConverter.cc index 5d2506a..07a8d47 100644 --- a/src/SlowConverter.cc +++ b/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");); diff --git a/src/Util.cc b/src/Util.cc index e338b00..968e968 100644 --- a/src/Util.cc +++ b/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); diff --git a/src/Util.h b/src/Util.h index be9519e..9f8c2a1 100644 --- a/src/Util.h +++ b/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);