1 /*
2 * Copyright (c) 2017-2021, The OSKAR Developers.
3 * See the LICENSE file at the top-level directory of this distribution.
4 */
5
6 #include "imager/private_imager.h"
7 #include "imager/private_imager_read_coords.h"
8 #include "imager/oskar_imager.h"
9 #include "binary/oskar_binary.h"
10 #include "convert/oskar_convert_station_uvw_to_baseline_uvw.h"
11 #include "math/oskar_cmath.h"
12 #include "mem/oskar_binary_read_mem.h"
13 #include "ms/oskar_measurement_set.h"
14 #include "utility/oskar_timer.h"
15 #include "vis/oskar_vis_block.h"
16 #include "vis/oskar_vis_header.h"
17
18 #include <float.h>
19 #include <stdlib.h>
20
21 #ifdef __cplusplus
22 extern "C" {
23 #endif
24
oskar_imager_read_coords_ms(oskar_Imager * h,const char * filename,int i_file,int num_files,int * percent_done,int * percent_next,int * status)25 void oskar_imager_read_coords_ms(oskar_Imager* h, const char* filename,
26 int i_file, int num_files, int* percent_done, int* percent_next,
27 int* status)
28 {
29 #ifndef OSKAR_NO_MS
30 oskar_MeasurementSet* ms = 0;
31 oskar_Mem *uvw = 0, *u = 0, *v = 0, *w = 0, *weight = 0, *time_centroid = 0;
32 size_t start_row = 0;
33 double *uvw_ = 0, *u_ = 0, *v_ = 0, *w_ = 0;
34 if (*status) return;
35
36 /* Read the header. */
37 oskar_log_message(h->log, 'M', 0, "Opening Measurement Set '%s'", filename);
38 ms = oskar_ms_open_readonly(filename);
39 if (!ms)
40 {
41 *status = OSKAR_ERR_FILE_IO;
42 return;
43 }
44 const size_t num_rows = (size_t) oskar_ms_num_rows(ms);
45 const size_t num_stations = (size_t) oskar_ms_num_stations(ms);
46 const size_t num_baselines = num_stations * (num_stations - 1) / 2;
47 const int num_pols = (int) oskar_ms_num_pols(ms);
48 const int num_channels = (int) oskar_ms_num_channels(ms);
49
50 /* Set visibility meta-data. */
51 oskar_imager_set_vis_frequency(h,
52 oskar_ms_freq_start_hz(ms),
53 oskar_ms_freq_inc_hz(ms), num_channels);
54 oskar_imager_set_vis_phase_centre(h,
55 oskar_ms_phase_centre_ra_rad(ms) * 180/M_PI,
56 oskar_ms_phase_centre_dec_rad(ms) * 180/M_PI);
57
58 /* Create arrays. */
59 uvw = oskar_mem_create(OSKAR_DOUBLE, OSKAR_CPU, 3 * num_baselines, status);
60 u = oskar_mem_create(OSKAR_DOUBLE, OSKAR_CPU, num_baselines, status);
61 v = oskar_mem_create(OSKAR_DOUBLE, OSKAR_CPU, num_baselines, status);
62 w = oskar_mem_create(OSKAR_DOUBLE, OSKAR_CPU, num_baselines, status);
63 weight = oskar_mem_create(OSKAR_SINGLE, OSKAR_CPU,
64 num_baselines * num_pols, status);
65 time_centroid = oskar_mem_create(OSKAR_DOUBLE, OSKAR_CPU, num_baselines,
66 status);
67 uvw_ = oskar_mem_double(uvw, status);
68 u_ = oskar_mem_double(u, status);
69 v_ = oskar_mem_double(v, status);
70 w_ = oskar_mem_double(w, status);
71
72 /* Loop over visibility blocks. */
73 for (start_row = 0; start_row < num_rows; start_row += num_baselines)
74 {
75 size_t allocated = 0, required = 0, block_size = 0, i = 0;
76 if (*status) break;
77
78 /* Read coordinates and weights from Measurement Set. */
79 oskar_timer_resume(h->tmr_read);
80 block_size = num_rows - start_row;
81 if (block_size > num_baselines) block_size = num_baselines;
82 allocated = oskar_mem_length(uvw) *
83 oskar_mem_element_size(oskar_mem_type(uvw));
84 oskar_ms_read_column(ms, "UVW", start_row, block_size,
85 allocated, oskar_mem_void(uvw), &required, status);
86 allocated = oskar_mem_length(weight) *
87 oskar_mem_element_size(oskar_mem_type(weight));
88 oskar_ms_read_column(ms, "WEIGHT", start_row, block_size,
89 allocated, oskar_mem_void(weight), &required, status);
90 allocated = oskar_mem_length(time_centroid) *
91 oskar_mem_element_size(oskar_mem_type(time_centroid));
92 oskar_ms_read_column(ms, "TIME_CENTROID", start_row, block_size,
93 allocated, oskar_mem_void(time_centroid), &required, status);
94
95 /* Split up baseline coordinates. */
96 for (i = 0; i < block_size; ++i)
97 {
98 u_[i] = uvw_[3*i + 0];
99 v_[i] = uvw_[3*i + 1];
100 w_[i] = uvw_[3*i + 2];
101 }
102
103 /* Update the imager with the data. */
104 oskar_timer_pause(h->tmr_read);
105 oskar_imager_update(h, block_size, 0, num_channels - 1,
106 num_pols, u, v, w, 0, weight, time_centroid, status);
107 *percent_done = (int) round(100.0 * (
108 (start_row + block_size) / (double)(num_rows * num_files) +
109 i_file / (double)num_files));
110 if (percent_next && *percent_done >= *percent_next)
111 {
112 oskar_log_message(h->log, 'S', -2, "%3d%% ...", *percent_done);
113 *percent_next = 10 + 10 * (*percent_done / 10);
114 }
115 }
116 oskar_mem_free(uvw, status);
117 oskar_mem_free(u, status);
118 oskar_mem_free(v, status);
119 oskar_mem_free(w, status);
120 oskar_mem_free(weight, status);
121 oskar_mem_free(time_centroid, status);
122 oskar_ms_close(ms);
123 #else
124 (void) filename;
125 (void) i_file;
126 (void) num_files;
127 (void) percent_done;
128 (void) percent_next;
129 oskar_log_error(h->log,
130 "OSKAR was compiled without Measurement Set support.");
131 *status = OSKAR_ERR_FUNCTION_NOT_AVAILABLE;
132 #endif
133 }
134
135
oskar_imager_read_coords_vis(oskar_Imager * h,const char * filename,int i_file,int num_files,int * percent_done,int * percent_next,int * status)136 void oskar_imager_read_coords_vis(oskar_Imager* h, const char* filename,
137 int i_file, int num_files, int* percent_done, int* percent_next,
138 int* status)
139 {
140 oskar_Binary* vis_file = 0;
141 oskar_VisHeader* hdr = 0;
142 oskar_Mem *u = 0, *v = 0, *w = 0, *uu = 0, *vv = 0, *ww = 0;
143 oskar_Mem *weight = 0, *time_centroid = 0;
144 int i_block = 0;
145 double time_start_mjd = 0.0, time_inc_sec = 0.0;
146 if (*status) return;
147
148 /* Read the header. */
149 oskar_log_message(h->log, 'M', 0, "Opening '%s'", filename);
150 vis_file = oskar_binary_create(filename, 'r', status);
151 hdr = oskar_vis_header_read(vis_file, status);
152 if (*status)
153 {
154 oskar_vis_header_free(hdr, status);
155 oskar_binary_free(vis_file);
156 return;
157 }
158 const int coord_prec = oskar_vis_header_coord_precision(hdr);
159 const int max_times_per_block = oskar_vis_header_max_times_per_block(hdr);
160 const int tags_per_block = oskar_vis_header_num_tags_per_block(hdr);
161 const int num_stations = oskar_vis_header_num_stations(hdr);
162 const int num_baselines = num_stations * (num_stations - 1) / 2;
163 const int num_pols =
164 oskar_type_is_matrix(oskar_vis_header_amp_type(hdr)) ? 4 : 1;
165 const int num_weights = num_baselines * num_pols * max_times_per_block;
166 const int num_blocks = oskar_vis_header_num_blocks(hdr);
167 const double freq_inc_hz = oskar_vis_header_freq_inc_hz(hdr);
168 const double freq_start_hz = oskar_vis_header_freq_start_hz(hdr);
169 time_start_mjd = oskar_vis_header_time_start_mjd_utc(hdr) * 86400.0;
170 time_inc_sec = oskar_vis_header_time_inc_sec(hdr);
171
172 /* Set visibility meta-data. */
173 oskar_imager_set_vis_frequency(h, freq_start_hz, freq_inc_hz,
174 oskar_vis_header_num_channels_total(hdr));
175 oskar_imager_set_vis_phase_centre(h,
176 oskar_vis_header_phase_centre_ra_deg(hdr),
177 oskar_vis_header_phase_centre_dec_deg(hdr));
178
179 /* Create scratch arrays. Weights are all 1. */
180 u = oskar_mem_create(coord_prec, OSKAR_CPU, 0, status);
181 v = oskar_mem_create(coord_prec, OSKAR_CPU, 0, status);
182 w = oskar_mem_create(coord_prec, OSKAR_CPU, 0, status);
183 uu = oskar_mem_create(coord_prec, OSKAR_CPU, 0, status);
184 vv = oskar_mem_create(coord_prec, OSKAR_CPU, 0, status);
185 ww = oskar_mem_create(coord_prec, OSKAR_CPU, 0, status);
186 time_centroid = oskar_mem_create(OSKAR_DOUBLE,
187 OSKAR_CPU, num_baselines * max_times_per_block, status);
188 weight = oskar_mem_create(h->imager_prec, OSKAR_CPU, num_weights, status);
189 oskar_mem_set_value_real(weight, 1.0, 0, num_weights, status);
190
191 /* Loop over visibility blocks. */
192 for (i_block = 0; i_block < num_blocks; ++i_block)
193 {
194 int c = 0, t = 0, dim_start_and_size[6], tag_error = 0;
195 if (*status) break;
196
197 /* Read block metadata. */
198 oskar_timer_resume(h->tmr_read);
199 oskar_binary_set_query_search_start(vis_file,
200 i_block * tags_per_block, status);
201 oskar_binary_read(vis_file, OSKAR_INT,
202 OSKAR_TAG_GROUP_VIS_BLOCK,
203 OSKAR_VIS_BLOCK_TAG_DIM_START_AND_SIZE, i_block,
204 sizeof(dim_start_and_size), dim_start_and_size, status);
205 const int start_time = dim_start_and_size[0];
206 const int start_chan = dim_start_and_size[1];
207 const int num_times = dim_start_and_size[2];
208 const int num_channels = dim_start_and_size[3];
209 const size_t num_rows = num_times * num_baselines;
210
211 /* Fill in the time centroid values. */
212 for (t = 0; t < num_times; ++t)
213 {
214 oskar_mem_set_value_real(time_centroid,
215 time_start_mjd + (start_time + t + 0.5) * time_inc_sec,
216 t * num_baselines, num_baselines, status);
217 }
218
219 /* Try to read station coordinates in the block. */
220 oskar_binary_read_mem(vis_file, u, OSKAR_TAG_GROUP_VIS_BLOCK,
221 OSKAR_VIS_BLOCK_TAG_STATION_U, i_block, &tag_error);
222 if (!tag_error)
223 {
224 oskar_binary_read_mem(vis_file, v, OSKAR_TAG_GROUP_VIS_BLOCK,
225 OSKAR_VIS_BLOCK_TAG_STATION_V, i_block, status);
226 oskar_binary_read_mem(vis_file, w, OSKAR_TAG_GROUP_VIS_BLOCK,
227 OSKAR_VIS_BLOCK_TAG_STATION_W, i_block, status);
228
229 /* Convert from station to baseline coordinates. */
230 oskar_mem_ensure(uu, num_times * num_baselines, status);
231 oskar_mem_ensure(vv, num_times * num_baselines, status);
232 oskar_mem_ensure(ww, num_times * num_baselines, status);
233 for (t = 0; t < num_times; ++t)
234 {
235 oskar_convert_station_uvw_to_baseline_uvw(num_stations,
236 num_stations * t, u, v, w,
237 num_baselines * t, uu, vv, ww, status);
238 }
239 }
240 else
241 {
242 /* Station coordinates not present,
243 * so read the baseline coordinates directly. */
244 oskar_binary_read_mem(vis_file, uu, OSKAR_TAG_GROUP_VIS_BLOCK,
245 OSKAR_VIS_BLOCK_TAG_BASELINE_UU, i_block, status);
246 oskar_binary_read_mem(vis_file, vv, OSKAR_TAG_GROUP_VIS_BLOCK,
247 OSKAR_VIS_BLOCK_TAG_BASELINE_VV, i_block, status);
248 oskar_binary_read_mem(vis_file, ww, OSKAR_TAG_GROUP_VIS_BLOCK,
249 OSKAR_VIS_BLOCK_TAG_BASELINE_WW, i_block, status);
250 }
251
252 /* Update the imager with the data. */
253 oskar_timer_pause(h->tmr_read);
254 for (c = 0; c < num_channels; ++c)
255 {
256 /* Update per channel. */
257 const double freq_hz =
258 freq_start_hz + (start_chan + c) * freq_inc_hz;
259 if (freq_hz >= h->freq_min_hz &&
260 (freq_hz <= h->freq_max_hz || h->freq_max_hz == 0.0))
261 {
262 oskar_imager_update(h, num_rows,
263 start_chan + c, start_chan + c, num_pols,
264 uu, vv, ww, 0, weight, time_centroid, status);
265 }
266 }
267 *percent_done = (int) round(100.0 * (
268 (i_block + 1) / (double)(num_blocks * num_files) +
269 i_file / (double)num_files));
270 if (percent_next && *percent_done >= *percent_next)
271 {
272 oskar_log_message(h->log, 'S', -2, "%3d%% ...", *percent_done);
273 *percent_next = 10 + 10 * (*percent_done / 10);
274 }
275 }
276 oskar_mem_free(u, status);
277 oskar_mem_free(v, status);
278 oskar_mem_free(w, status);
279 oskar_mem_free(uu, status);
280 oskar_mem_free(vv, status);
281 oskar_mem_free(ww, status);
282 oskar_mem_free(weight, status);
283 oskar_mem_free(time_centroid, status);
284 oskar_vis_header_free(hdr, status);
285 oskar_binary_free(vis_file);
286 }
287
288 #ifdef __cplusplus
289 }
290 #endif
291