File: purify_fitsio.cc

package info (click to toggle)
purify 5.0.1-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 186,836 kB
  • sloc: cpp: 17,731; python: 510; xml: 182; makefile: 7; sh: 6
file content (139 lines) | stat: -rw-r--r-- 6,484 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
#include "catch2/catch_all.hpp"

#include "purify/config.h"
#include "purify/types.h"
#include "purify/logging.h"

#include "purify/directories.h"
#include "purify/pfitsio.h"

using namespace purify;

/*TEST_CASE("Purify fitsio", "[readwrite]") {
  Image<t_complex> input = pfitsio::read2d(image_filename("M31.fits"));
  pfitsio::write2d(input.real(), output_filename("fits_output.fits"));
  Image<t_complex> input2 = pfitsio::read2d(output_filename("fits_output.fits"));
  CAPTURE(input2);
  CHECK(input.isApprox(input2, 1e-12));
}

TEST_CASE("readwrite2dheader", "purify fitsio") {
  Image<t_complex> input = pfitsio::read2d(image_filename("M31.fits"));
  pfitsio::header_params header_example;
  header_example.fits_name = output_filename("fits_output.fits");
  pfitsio::write2d(input.real(), header_example);
  Image<t_complex> input2 = pfitsio::read2d(output_filename("fits_output.fits"));
  CHECK(input.isApprox(input2, 1e-12));
}

TEST_CASE("readwrite3d", "purify fitsio") {
  purify::logging::set_level("debug");
  std::vector<Image<t_complex>> input = pfitsio::read3d(image_filename("cube_example.fits"));
  CHECK(input.size() == 4);
  CHECK(input[0].size() == 200 * 200);
  std::vector<Image<t_real>> input_real;
  for (int i = 0; i < input.size(); i++) {
    input_real.push_back(input[i].real());
  }
  pfitsio::write3d(input_real, output_filename("cube_output.fits"));
  std::vector<Image<t_complex>> input2 = pfitsio::read3d(output_filename("cube_output.fits"));
  for (int i = 0; i < input.size(); i++) {
    CHECK(input[i].isApprox(input2[i], 1e-12));
  }
}
TEST_CASE("readwrite3dheader", "purify fitsio") {
  std::vector<Image<t_complex>> input = pfitsio::read3d(image_filename("cube_example.fits"));
  CHECK(input.size() == 4);
  CHECK(input.at(0).size() == 200 * 200);
  std::vector<Image<t_real>> input_real;
  for (int i = 0; i < input.size(); i++) {
    CAPTURE(input.size());
    CAPTURE(i);
    input_real.push_back(input.at(i).real());
  }
  pfitsio::header_params header_example;
  header_example.fits_name = output_filename("cube_output.fits");
  pfitsio::write3d(input_real, header_example);
  std::vector<Image<t_complex>> input2 = pfitsio::read3d(output_filename("cube_output.fits"));
  CAPTURE(input.size());
  for (int i = 0; i < input.size(); i++) {
    CAPTURE(i);
    CHECK(input.at(i).isApprox(input2.at(i), 1e-12));
  }
}
TEST_CASE("readwrite3dwith2d", "purify fitsio") {
  std::vector<Image<t_complex>> input = pfitsio::read3d(image_filename("M31.fits"));
  CHECK(input.size() == 1);
  CHECK(input[0].size() == 256 * 256);
  std::vector<Image<t_real>> input_real;
  for (int i = 0; i < input.size(); i++) {
    input_real.push_back(input[i].real());
  }
  pfitsio::write3d(input_real, output_filename("2dcube_output.fits"));
  std::vector<Image<t_complex>> input2 = pfitsio::read3d(output_filename("2dcube_output.fits"));
  for (int i = 0; i < input.size(); i++) CHECK(input[i].isApprox(input2[i], 1e-12));
}
TEST_CASE("readwrite3dheaderwith2d", "purify fitsio") {
  std::vector<Image<t_complex>> input = pfitsio::read3d(image_filename("M31.fits"));
  CHECK(input.size() == 1);
  CHECK(input[0].size() == 256 * 256);
  std::vector<Image<t_real>> input_real;
  for (int i = 0; i < input.size(); i++) input_real.push_back(input[i].real());

  pfitsio::header_params header_example;
  header_example.fits_name = output_filename("2dcube_output.fits");
  pfitsio::write3d(input_real, header_example);
  std::vector<Image<t_complex>> input2 = pfitsio::read3d(output_filename("2dcube_output.fits"));
  for (int i = 0; i < input.size(); i++) CHECK(input[i].isApprox(input2[i], 1e-12));
}*/

TEST_CASE("header") {
  const std::string fits_name = "test_image.fits";
  const t_real mean_frequency = 123;  // in MHz
  const t_real cell_x = 5;            // in arcseconds
  const t_real cell_y = 2;            // in arcseconds
  const t_real ra = 12;               // in radians, converted to decimal degrees before write
  const t_real dec = 98;              // in radians, converted to decimal degrees before write
  const std::string pix_units = "Jy/PIXEL";
  const t_real channels_total = 2;
  const t_real channel_width = 11;  // in MHz
  const t_real polarisation = stokes_int.at(stokes::I);
  const int niters = 10;           // number of iterations
  const bool hasconverged = true;  // stating if model has converged
  const t_real relative_variation = 1e-3;
  const t_real residual_convergence = 1e-4;
  const t_real epsilon = 10;
  pfitsio::header_params header;
  header.fits_name = fits_name;
  header.mean_frequency = mean_frequency;
  header.cell_x = cell_x;
  header.cell_y = cell_y;
  header.ra = ra;
  header.dec = dec;
  header.pix_units = pix_units;
  header.channels_total = channels_total;
  header.channel_width = channel_width;
  header.polarisation = polarisation;
  header.niters = niters;
  header.hasconverged = hasconverged;
  header.relative_variation = relative_variation;
  header.residual_convergence = residual_convergence;
  header.epsilon = epsilon;
  const auto header_test = pfitsio::header_params(
      fits_name, pix_units, channels_total, ra, dec, stokes::I, cell_x, cell_y, mean_frequency,
      channel_width, niters, hasconverged, relative_variation, residual_convergence, epsilon);
  CHECK(header_test == header);
  CHECK(header_test == pfitsio::header_params(fits_name, pix_units, channels_total, ra, dec,
                                              stokes_string.at("p"), cell_x, cell_y, mean_frequency,
                                              channel_width, niters, hasconverged,
                                              relative_variation, residual_convergence, epsilon));
  CHECK_THROWS(header_test == pfitsio::header_params(fits_name, pix_units, channels_total, ra, dec,
                                                     stokes_string.at("z"), cell_x, cell_y,
                                                     mean_frequency, channel_width, niters,
                                                     hasconverged, relative_variation,
                                                     residual_convergence, epsilon));
  CHECK(header_test != pfitsio::header_params(fits_name, pix_units, channels_total, ra, dec,
                                              stokes::Q, cell_x, cell_y, mean_frequency,
                                              channel_width, niters, hasconverged,
                                              relative_variation, residual_convergence, epsilon));
}