CARTAvis / carta-backend

Source code repository for the backend component of CARTA, a new visualization tool designed for the ALMA, the VLA and the SKA pathfinders.
https://cartavis.github.io/
GNU General Public License v3.0
22 stars 11 forks source link

Optimise reading of 2D slices from FITS files #965

Open veggiesaurus opened 3 years ago

veggiesaurus commented 3 years ago

A quick implementation of the simplest way to read a 2D slice from a FITS image with cfitsio is shown below:


#include <vector>
#include <chrono>

#include <fmt/format.h>
#include <fitsio.h>

int main(int argc, char** argv) {
    int status = 0;
    fitsfile* fptr;
    int dims = 4;
    std::vector<int64_t> imageSize(dims);
    std::vector<float> data_vector;
    float* data_ptr = nullptr;

    std::vector<int64_t> startPix(dims);
    for (int i = 0; i < dims; i++) {
        startPix[i] = 1;
    }

    // Open image and get data size. File is assumed to be 2D
    fits_open_file(&fptr, argv[1], READONLY, &status);
    fits_get_img_size(fptr, dims, imageSize.data(), &status);
    int64_t num_pixels = imageSize[0] * imageSize[1];
    float test_val = 0;

    // Frame::FillImageCache timing is just the memory allocation and reading part
    auto t_start = std::chrono::high_resolution_clock::now();
    data_ptr = new float[num_pixels];

    fits_read_pix(fptr, TFLOAT, startPix.data(), num_pixels, nullptr, data_ptr, nullptr, &status);
    test_val = data_ptr[num_pixels / 2];
    auto t_end = std::chrono::high_resolution_clock::now();

    auto dt = std::chrono::duration_cast<std::chrono::microseconds>(t_end - t_start).count();
    auto rate = num_pixels / (double) dt;
    fmt::print("{:.1f} MB file read in {:.1f} ms at {:.1f} MPix/s. Test val={}\n", num_pixels * sizeof(float) * 1.0e-6, dt * 1.0e-3, rate, test_val);

    fits_close_file(fptr, &status);
    delete[] data_ptr;

    return status;
}

The approach above is quite a bit (50%) faster than Frame::FillImageCache, and it can be optimised further.

veggiesaurus commented 3 years ago

Estimating as 2 story points because this will require a little refactoring.