1 /*
2 * Copyright (c) 2019-2021, The OSKAR Developers.
3 * See the LICENSE file at the top-level directory of this distribution.
4 */
5
6 #include "mem/oskar_mem.h"
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <fitsio.h>
10 #include <string.h>
11
12 #ifdef __cplusplus
13 extern "C" {
14 #endif
15
16 #define MAX_AXES 10
17
oskar_mem_read_fits(oskar_Mem * data,size_t offset,size_t num_pixels,const char * file_name,int num_index_dims,const int * start_index,int * num_axes,int ** axis_size,double ** axis_inc,int * status)18 void oskar_mem_read_fits(oskar_Mem* data, size_t offset, size_t num_pixels,
19 const char* file_name, int num_index_dims, const int* start_index,
20 int* num_axes, int** axis_size, double** axis_inc, int* status)
21 {
22 int i = 0, imagetype = 0, naxis = 0, anynul = 0, type_fits = 0;
23 long firstpix[MAX_AXES], naxes[MAX_AXES], num_pixels_cube = 1;
24 double cdelt[MAX_AXES];
25 double nul = 0.0;
26 oskar_Mem *data_ptr = 0, *data_temp = 0;
27 fitsfile* fptr = 0;
28 if (*status) return;
29
30 /* Open the file and get the cube parameters. */
31 fits_open_file(&fptr, file_name, READONLY, status);
32 if (*status || !fptr)
33 {
34 *status = OSKAR_ERR_FILE_IO;
35 return;
36 }
37 fits_get_img_param(fptr, MAX_AXES, &imagetype, &naxis, naxes, status);
38 if (*status || naxis < 1 || naxis > MAX_AXES)
39 {
40 fits_close_file(fptr, status);
41 *status = OSKAR_ERR_FILE_IO;
42 return;
43 }
44
45 /* Store axis sizes and increments. */
46 if (axis_inc)
47 {
48 fits_read_keys_dbl(fptr, "CDELT", 1, naxis, cdelt, &i, status);
49 }
50 if (num_axes && *num_axes < naxis)
51 {
52 if (axis_size)
53 {
54 *axis_size = (int*) realloc(*axis_size, naxis * sizeof(int));
55 }
56 if (axis_inc)
57 {
58 *axis_inc = (double*) realloc(*axis_inc, naxis * sizeof(double));
59 }
60 }
61 if (num_axes) *num_axes = naxis;
62 for (i = 0; i < naxis; ++i)
63 {
64 firstpix[i] = 1;
65 if (axis_size) (*axis_size)[i] = naxes[i];
66 if (axis_inc) (*axis_inc)[i] = cdelt[i];
67 num_pixels_cube *= naxes[i];
68 }
69 if (!data || num_index_dims == 0 || !start_index)
70 {
71 fits_close_file(fptr, status);
72 return;
73 }
74
75 /* Get the FITS data type of the output array. */
76 switch (oskar_mem_type(data))
77 {
78 case OSKAR_INT:
79 type_fits = TINT;
80 break;
81 case OSKAR_SINGLE:
82 type_fits = TFLOAT;
83 break;
84 case OSKAR_DOUBLE:
85 type_fits = TDOUBLE;
86 break;
87 default:
88 fits_close_file(fptr, status);
89 *status = OSKAR_ERR_BAD_DATA_TYPE;
90 return;
91 }
92
93 /* Read the data. */
94 offset *= oskar_mem_element_size(oskar_mem_type(data));
95 if (num_index_dims > 0)
96 {
97 for (i = 0; i < naxis && i < num_index_dims; ++i)
98 {
99 firstpix[i] = 1 + start_index[i];
100 }
101 }
102 else if (num_pixels == 0)
103 {
104 num_pixels = (size_t) num_pixels_cube;
105 }
106 data_ptr = data;
107 if (oskar_mem_location(data) != OSKAR_CPU)
108 {
109 data_temp = oskar_mem_create(
110 oskar_mem_type(data), OSKAR_CPU, num_pixels, status);
111 data_ptr = data_temp;
112 }
113 oskar_mem_ensure(data_ptr, num_pixels, status);
114 fits_read_pix(fptr, type_fits, firstpix, (long) num_pixels,
115 &nul, oskar_mem_char(data_ptr) + offset, &anynul, status);
116 fits_close_file(fptr, status);
117 if (data_ptr != data)
118 {
119 oskar_mem_copy(data, data_ptr, status);
120 }
121 oskar_mem_free(data_temp, status);
122 }
123
124 #ifdef __cplusplus
125 }
126 #endif
127