62 fail(
"invalid end of string");
63 return std::string(
reinterpret_cast<const char*
>(
header_word(w)), len);
81 std::memcpy(
header_word(w), str.c_str(), str.size());
87 std::array<int, 3> pos{{-1, -1, -1}};
88 for (
int i = 0;
i != 3; ++
i) {
91 fail(
"Incorrect MAPC/MAPR/MAPS records");
108 for (
int i = 0;
i < 3; ++
i) {
109 double scale = 1. /
sampl[
i];
111 box.minimum.at(
i) = scale * start[
p] - 1
e-9;
112 box.maximum.at(
i) = scale * (start[
p] + size[
p] - 1) + 1
e-9;
142template<
typename T=
float>
189 fail(
"Only modes 0, 1, 2 and 6 are supported.");
191 fail(
"update_ccp4_header(): set the grid first (it has size 0)");
193 fail(
"update_ccp4_header(): run setup() first");
202 fail(
"update_ccp4_header: specify map mode explicitly (usually 2)");
213 if (std::is_same<T, std::int8_t>::value)
215 if (std::is_same<T, std::int16_t>::value)
217 if (std::is_same<T, float>::value)
219 if (std::is_same<T, std::uint16_t>::value)
234 template<
typename Stream>
236 const size_t hsize = 256;
239 fail(
"Failed to read map header: " + path);
241 fail(
"Not a CCP4 map: " + path);
244 fail(
"Unsupported machine stamp (endianness) in the file?");
251 fail(
"Unexpectedly long extended header: " + path);
254 fail(
"Failed to read extended header: " + path);
259 for (
int i = 0;
i < 3; ++
i) {
262 fail(
"Unexpected axis value in word " + std::to_string(17 +
i)
263 +
": " + std::to_string(
axis));
272 if (pos[0] == 0 && pos[1] == 1 && pos[2] == 2 &&
full_cell()) {
281 template<
typename Stream>
290 const std::string& name) {
294 template<
typename Input>
296 if (
input.is_stdin())
298 else if (
input.is_compressed())
310template<
typename From,
typename To>
311To translate_map_point(From f) {
return static_cast<To
>(f); }
314std::int8_t translate_map_point<float,std::int8_t>(
float f) {
return f != 0; }
316template<
typename Stream,
typename TFile,
typename TMem>
317void read_data(Stream& f, std::vector<TMem>& content) {
318 if (std::is_same<TFile, TMem>::value) {
319 size_t len = content.size();
320 if (!f.read(content.data(),
sizeof(TMem) * len))
321 fail(
"Failed to read all the data from the map file.");
323 constexpr size_t chunk_size = 64 * 1024;
324 std::vector<TFile> work(chunk_size);
325 for (
size_t i = 0; i < content.size(); i += chunk_size) {
326 size_t len = std::min(chunk_size, content.size() - i);
327 if (!f.read(work.data(),
sizeof(TFile) * len))
328 fail(
"Failed to read all the data from the map file.");
329 for (
size_t j = 0; j < len; ++j)
330 content[i+j] = translate_map_point<TFile,TMem>(work[j]);
335template<
typename TFile,
typename TMem>
336void write_data(
const std::vector<TMem>& content, FILE* f) {
337 if (std::is_same<TMem, TFile>::value) {
338 size_t len = content.size();
339 if (std::fwrite(content.data(),
sizeof(TFile), len, f) != len)
340 sys_fail(
"Failed to write data to the map file");
342 constexpr size_t chunk_size = 64 * 1024;
343 std::vector<TFile> work(chunk_size);
344 for (
size_t i = 0; i < content.size(); i += chunk_size) {
345 size_t len = std::min(chunk_size, content.size() - i);
346 for (
size_t j = 0; j < len; ++j)
347 work[j] =
static_cast<TFile
>(content[i+j]);
348 if (std::fwrite(work.data(),
sizeof(TFile), len, f) != len)
349 sys_fail(
"Failed to write data to the map file");
358template<
typename T>
template<
typename Stream>
360 read_ccp4_header(f, path);
361 grid.data.resize(grid.point_count());
362 int mode = header_i32(4);
364 impl::read_data<Stream, std::int8_t>(f, grid.data);
366 impl::read_data<Stream, std::int16_t>(f, grid.data);
368 impl::read_data<Stream, float>(f, grid.data);
370 impl::read_data<Stream, std::uint16_t>(f, grid.data);
372 fail(
"Mode " + std::to_string(
mode) +
" is not supported "
373 "(only 0, 1, 2 and 6 are supported).");
377 if (!same_byte_order) {
379 for (
T&
value : grid.data)
381 else if (
sizeof(
T) == 4)
382 for (
T&
value : grid.data)
392 const std::array<int, 3>
sampl = header_3i32(8);
394 const std::array<int, 3> pos = axis_positions();
395 std::array<int, 3> start = header_3i32(5);
396 int end[3] = { start[0] + grid.nu, start[1] + grid.nv, start[2] + grid.nw };
399 set_header_3i32(5, start[pos[0]], start[pos[1]], start[pos[2]]);
400 for (
int i = 0;
i < 3; ++
i) {
404 int crs[3] = { grid.nu, grid.nv, grid.nw };
405 grid.nu =
crs[pos[0]];
406 grid.nv =
crs[pos[1]];
407 grid.nw =
crs[pos[2]];
412 set_header_3i32(5, 0, 0, 0);
414 set_header_3i32(1, grid.nu, grid.nv, grid.nw);
415 set_header_3i32(17, 1, 2, 3);
418 grid.calculate_spacing();
425 for (
it[2] = start[2];
it[2] < end[2];
it[2]++)
426 for (
it[1] = start[1];
it[1] < end[1];
it[1]++)
427 for (
it[0] = start[0];
it[0] < end[0];
it[0]++) {
428 T val = grid.data[idx++];
432 grid.data = std::move(
full);
437 (end[pos[0]] - start[pos[0]] <
sampl[0] ||
438 end[pos[1]] - start[pos[1]] <
sampl[1] ||
439 end[pos[2]] - start[pos[2]] <
sampl[2]))
445 if (ccp4_header.empty())
446 fail(
"set_extent(): no header in the map. Call update_ccp4_header() first");
448 fail(
"Ccp4::set_extent() works only after setup()");
450 fail(
"Ccp4::set_extent() works only with XYZ order");
451 int u0 = (
int)std::ceil(
box.minimum.x * grid.nu);
452 int v0 = (
int)std::ceil(
box.minimum.y * grid.nv);
453 int w0 = (
int)std::ceil(
box.minimum.z * grid.nw);
454 int nu = (
int)std::floor(
box.maximum.x * grid.nu) -
u0 + 1;
455 int nv = (
int)std::floor(
box.maximum.y * grid.nv) -
v0 + 1;
456 int nw = (
int)std::floor(
box.maximum.z * grid.nw) -
w0 + 1;
458 std::vector<T>
new_data((
size_t)nu * nv * nw);
459 grid.get_subarray(
new_data.data(), {u0, v0, w0}, {nu, nv, nw});
465 set_header_3i32(1, grid.nu, grid.nv, grid.nw);
466 set_header_3i32(5,
u0,
v0,
w0);
473 assert(ccp4_header.size() >= 256);
475 std::fwrite(ccp4_header.data(), 4, ccp4_header.size(), f.get());
476 int mode = header_i32(4);
478 impl::write_data<std::int8_t>(grid.data, f.get());
480 impl::write_data<std::int16_t>(grid.data, f.get());
482 impl::write_data<float>(grid.data, f.get());
484 impl::write_data<std::uint16_t>(grid.data, f.get());
void swap_four_bytes(void *start)
std::unique_ptr< std::FILE, decltype(&std::fclose)> fileptr_t
const SpaceGroup * find_spacegroup_by_number(int ccp4) noexcept
fileptr_t file_open(const char *path, const char *mode)
void swap_two_bytes(void *start)
GEMMI_COLD void sys_fail(const std::string &msg)
DataStats calculate_data_statistics(const std::vector< T > &data)
void fail(const std::string &msg)
void set_header_str(int w, const std::string &str)
std::string header_str(int w, size_t len=80) const
std::array< int, 3 > axis_positions() const
void set_header_i32(int w, int32_t value)
void set_header_3i32(int w, int32_t x, int32_t y, int32_t z)
Position get_origin() const
std::array< int, 3 > header_3i32(int w) const
Box< Fractional > get_extent() const
int32_t header_i32(int w) const
float header_float(int w) const
void * header_word(int w)
double header_rfloat(int w) const
void set_header_float(int w, float value)
const void * header_word(int w) const
Transform get_skew_transformation() const
bool has_skew_transformation() const
std::vector< int32_t > ccp4_header
void read_ccp4_from_memory(const char *data, size_t size, const std::string &name)
void update_ccp4_header(int mode=-1, bool update_stats=true)
If the header is empty, prepare it; otherwise, update only MODE and, if update_stats==true,...
void setup(T default_value, MapSetup mode=MapSetup::Full)
static int mode_for_data()
void read_ccp4_file(const std::string &path)
void write_ccp4_map(const std::string &path) const
void prepare_ccp4_header_except_mode_and_stats()
void set_extent(const Box< Fractional > &box)
void read_ccp4_header(Stream &f, const std::string &path)
void read_ccp4(Input &&input)
void read_ccp4_stream(Stream f, const std::string &path)
void calculate_spacing()
set spacing and orth_n
Coordinates in Angstroms - orthogonal (Cartesian) coordinates.
void set(double a_, double b_, double c_, double alpha_, double beta_, double gamma_)