Commit 249408f0 authored by shadow_walker's avatar shadow_walker
Browse files

Merge branch 'mrcio' into compression

parents 31f1b6f9 1f326521
......@@ -1106,237 +1106,8 @@ int MrcIO::read_data(float *rdata, int image_index, const Region * area, bool)
return 0;
}
int MrcIO::write_data(float *data, int image_index, const Region* area,
EMUtil::EMDataType dt, bool use_host_endian)
{
ENTERFUNC;
bool append = (image_index == -1);
if (is_stack && append) {
image_index = stack_size - 1;
}
int max_images = 0;
if (! is_stack) {
max_images = 1;
image_index = 0;
}
check_write_access(rw_mode, image_index, max_images, data);
check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
int nx, ny, nz;
if (! area) {
nx = mrch.nx;
ny = mrch.ny;
nz = mrch.nz;
}
else {
nx = (int)(area->get_width());
ny = (int)(area->get_height());
nz = (int)(area->get_depth());
}
bool got_one_image = (nz > 1);
if (is_stack && ! got_one_image) {
nz = 1;
}
size_t size = (size_t)nx * ny * nz;
if (is_complex_mode()) {
nx *= 2;
if (! is_ri) {
Util::ap2ri(data, size);
is_ri = 1;
}
Util::flip_complex_phase(data, size);
Util::rotate_phase_origin(data, nx, ny, nz);
}
portable_fseek(file, sizeof(MrcHeader), SEEK_SET);
if(dt == EMUtil::EM_COMPRESSED) {
if (renderbits <= 0) mrch.mode = MRC_FLOAT;
else if (renderbits <= 8) mrch.mode = MRC_UCHAR;
else if (renderbits <= 16) mrch.mode = MRC_USHORT;
}
if ((is_big_endian != ByteOrder::is_host_big_endian()) || ! use_host_endian) {
if (mrch.mode != MRC_UCHAR && mrch.mode != MRC_CHAR) {
if (mode_size == sizeof(short)) {
ByteOrder::swap_bytes((short*) data, size);
}
else if (mode_size == sizeof(float)) {
ByteOrder::swap_bytes((float*) data, size);
}
}
}
mode_size = get_mode_size(mrch.mode);
void * ptr_data = data;
int truebits=0;
switch(mrch.mode) {
case MRC_UCHAR: truebits=8; break;
case MRC_SHORT: truebits=16; break;
case MRC_FLOAT: truebits=32; break;
case MRC_SHORT_COMPLEX: truebits=32; break;
case MRC_FLOAT_COMPLEX: truebits=25; break;
case MRC_USHORT: truebits=16; break;
case MRC_UCHAR3: truebits=8; break;
case MRC_CHAR: truebits=8; break;
case MRC_UHEX: truebits=4; break;
case MRC_UNKNOWN: truebits=0; break;
}
if (renderbits==0 || renderbits>truebits) renderbits=truebits;
float rmin = rendermin;
float rmax = rendermax;
int rbits = renderbits;
EMUtil::getRenderMinMax(data, nx, ny, rmin, rmax, rbits, nz);
signed char * scdata = NULL;
unsigned char * cdata = NULL;
short * sdata = NULL;
unsigned short * usdata = NULL;
if (is_stack && ! got_one_image) {
mrch.nz = stack_size;
}
bool dont_rescale;
if (mrch.mode == MRC_UCHAR) {
cdata = new unsigned char[size];
dont_rescale = (rmin == 0.0 && rmax == UCHAR_MAX);
for (size_t i = 0; i < size; ++i) {
if (data[i] <= rmin) {
cdata[i] = 0;
}
else if (data[i] >= rmax){
cdata[i] = UCHAR_MAX;
}
else if (dont_rescale) {
cdata[i] = (unsigned char) data[i];
}
else {
cdata[i] = (unsigned char)((data[i] - rmin) /
(rmax - rmin) *
UCHAR_MAX);
}
}
ptr_data = cdata;
update_stats(ptr_data, size);
}
else if (mrch.mode == MRC_CHAR) {
scdata = new signed char[size];
dont_rescale = (rmin == SCHAR_MIN && rmax == SCHAR_MAX);
for (size_t i = 0; i < size; ++i) {
if (data[i] <= rmin) {
scdata[i] = SCHAR_MIN;
}
else if (data[i] >= rmax){
scdata[i] = SCHAR_MAX;
}
else if (dont_rescale) {
scdata[i] = (signed char) data[i];
}
else {
scdata[i] = (signed char)((data[i] - rmin) /
(rmax - rmin) *
(SCHAR_MAX - SCHAR_MIN) + SCHAR_MIN);
}
}
ptr_data = scdata;
update_stats(ptr_data, size);
}
else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
sdata = new short[size];
dont_rescale = (rmin == SHRT_MIN && rmax == SHRT_MAX);
for (size_t i = 0; i < size; ++i) {
if (data[i] <= rmin) {
sdata[i] = SHRT_MIN;
}
else if (data[i] >= rmax) {
sdata[i] = SHRT_MAX;
}
else if (dont_rescale) {
sdata[i] = (short) data[i];
}
else {
sdata[i] = (short)(((data[i] - rmin) /
(rmax - rmin)) *
(SHRT_MAX - SHRT_MIN) + SHRT_MIN);
}
}
ptr_data = sdata;
update_stats(ptr_data, size);
}
else if (mrch.mode == MRC_USHORT) {
usdata = new unsigned short[size];
dont_rescale = (rmin == 0.0 && rmax == USHRT_MAX);
for (size_t i = 0; i < size; ++i) {
if (data[i] <= rmin) {
usdata[i] = 0;
}
else if (data[i] >= rmax) {
usdata[i] = USHRT_MAX;
}
else if (dont_rescale) {
usdata[i] = (unsigned short) data[i];
}
else {
usdata[i] = (unsigned short)((data[i] - rmin) /
(rmax - rmin) *
USHRT_MAX);
}
}
ptr_data = usdata;
update_stats(ptr_data, size);
}
if (is_stack && ! got_one_image) {
mrch.nz = 1;
}
EMUtil::process_region_io(ptr_data, file, WRITE_ONLY, image_index,
mode_size, mrch.nx, mrch.ny, mrch.nz, area);
if (cdata) {delete [] cdata; cdata = NULL;}
if (scdata) {delete [] scdata; scdata = NULL;}
if (sdata) {delete [] sdata; sdata = NULL;}
if (usdata) {delete [] usdata; usdata = NULL;}
EXITFUNC;
return 0;
}
void MrcIO::update_stats(void * data, size_t size)
template<class T>
void MrcIO::update_stats(const vector<T> &data)
{
float v; // variable to hold pixel value
double sum;
......@@ -1345,7 +1116,7 @@ void MrcIO::update_stats(void * data, size_t size)
double sigma;
double vv;
float min, max;
signed char * scdata = NULL;
unsigned char * cdata = NULL;
short * sdata = NULL;
......@@ -1355,32 +1126,33 @@ void MrcIO::update_stats(void * data, size_t size)
bool use_uchar = (mrch.mode == MRC_UCHAR);
bool use_short = (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX);
bool use_ushort = (mrch.mode == MRC_USHORT);
if (use_uchar) {
max = 0.0;
min = UCHAR_MAX;
cdata = (unsigned char *) data;
cdata = (unsigned char *) data.data();
}
else if (use_schar) {
max = SCHAR_MIN;
min = SCHAR_MAX;
scdata = (signed char *) data;
scdata = (signed char *) data.data();
}
else if (use_short) {
max = (float) SHRT_MIN;
min = (float) SHRT_MAX;
sdata = (short *) data;
sdata = (short *) data.data();
}
else if (use_ushort) {
max = 0.0f;
min = (float) USHRT_MAX;
usdata = (unsigned short *) data;
usdata = (unsigned short *) data.data();
}
else {
throw InvalidCallException("This function is used to write 8bit/16bit mrc file only.");
}
sum = 0.0;
auto size = data.size();
for (size_t i = 0; i < size; i++) {
if (use_uchar) {
......@@ -1402,12 +1174,10 @@ void MrcIO::update_stats(void * data, size_t size)
sum = sum + v;
}
if (size > 0) {
if (size > 0)
mean = sum / (double) size;
}
else {
else
mean = 0.0;
}
square_sum = 0.0;
......@@ -1430,12 +1200,10 @@ void MrcIO::update_stats(void * data, size_t size)
square_sum = square_sum + vv * vv;
}
if (size > 1) {
if (size > 1)
sigma = std::sqrt(square_sum / (double) (size-1));
}
else {
else
sigma = 0.0;
}
/* change mrch.amin / amax / amean / rms here */
......@@ -1443,16 +1211,136 @@ void MrcIO::update_stats(void * data, size_t size)
mrch.amax = max;
mrch.amean = (float) mean;
mrch.rms = (float) sigma;
portable_fseek(file, 0, SEEK_SET);
if (fwrite(& mrch, sizeof(MrcHeader), 1, file) != 1) {
throw ImageWriteException(filename, "Error writing MRC header to update statistics.");
}
portable_fseek(file, sizeof(MrcHeader), SEEK_SET);
}
template<class T>
auto MrcIO::write_compressed(float *data, size_t size, int image_index, const Region* area) {
void * ptr_data = data;
if constexpr (!std::is_same<T, float>::value) {
auto [rendered_data, count] = getRenderedDataAndRendertrunc<T>(data, size);
update_stats<T>(rendered_data);
ptr_data = rendered_data.data();
}
EMUtil::process_region_io(ptr_data, file, WRITE_ONLY, image_index,
mode_size, mrch.nx, mrch.ny, mrch.nz, area);
}
int MrcIO::write_data(float *data, int image_index, const Region* area,
EMUtil::EMDataType dt, bool use_host_endian)
{
ENTERFUNC;
if (is_stack && (image_index == -1))
image_index = stack_size - 1;
int max_images = 0;
if (! is_stack) {
max_images = 1;
image_index = 0;
}
check_write_access(rw_mode, image_index, max_images, data);
check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
int nx, ny, nz;
if (! area) {
nx = mrch.nx;
ny = mrch.ny;
nz = mrch.nz;
}
else {
nx = (int)(area->get_width());
ny = (int)(area->get_height());
nz = (int)(area->get_depth());
}
if (is_stack && ! (nz > 1)) {
nz = 1;
mrch.nz = 1;
}
size_t size = (size_t)nx * ny * nz;
if (is_complex_mode()) {
nx *= 2;
if (! is_ri) {
Util::ap2ri(data, size);
is_ri = 1;
}
Util::flip_complex_phase(data, size);
Util::rotate_phase_origin(data, nx, ny, nz);
}
portable_fseek(file, sizeof(MrcHeader), SEEK_SET);
if(dt == EMUtil::EM_COMPRESSED) {
if (renderbits <= 0) mrch.mode = MRC_FLOAT;
else if (renderbits <= 8) mrch.mode = MRC_UCHAR;
else if (renderbits <= 16) mrch.mode = MRC_USHORT;
}
if ((is_big_endian != ByteOrder::is_host_big_endian()) || ! use_host_endian) {
if (mrch.mode != MRC_UCHAR && mrch.mode != MRC_CHAR) {
if (mode_size == sizeof(short))
ByteOrder::swap_bytes((short*) data, size);
else if (mode_size == sizeof(float))
ByteOrder::swap_bytes((float*) data, size);
}
}
mode_size = get_mode_size(mrch.mode);
int truebits=0;
switch(mrch.mode) {
case MRC_UCHAR: truebits=8; break;
case MRC_SHORT: truebits=16; break;
case MRC_FLOAT: truebits=32; break;
case MRC_SHORT_COMPLEX: truebits=32; break;
case MRC_FLOAT_COMPLEX: truebits=25; break;
case MRC_USHORT: truebits=16; break;
case MRC_UCHAR3: truebits=8; break;
case MRC_CHAR: truebits=8; break;
case MRC_UHEX: truebits=4; break;
case MRC_UNKNOWN: truebits=0; break;
}
if (renderbits==0 || renderbits>truebits) renderbits=truebits;
float rmin = rendermin;
float rmax = rendermax;
int rbits = renderbits;
EMUtil::getRenderMinMax(data, nx, ny, rmin, rmax, rbits, nz);
if (mrch.mode == MRC_UCHAR)
write_compressed<unsigned char>(data, size, image_index, area);
else if (mrch.mode == MRC_CHAR)
write_compressed<char>(data, size, image_index, area);
else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX)
write_compressed<short>(data, size, image_index, area);
else if (mrch.mode == MRC_USHORT)
write_compressed<unsigned short>(data, size, image_index, area);
else
write_compressed<float>(data, size, image_index, area);
EXITFUNC;
return 0;
}
bool MrcIO::is_complex_mode()
{
init();
......
......@@ -293,7 +293,8 @@ namespace EMAN
* when write MRC file as 16 bit or 8 bit. It will write new set of
* min/max/mean/sigma to mrch.
* this function needs get the output data storage type from mrch.mode.*/
void update_stats(void* data, size_t size);
template<class T>
void update_stats(const vector<T> &data);
/** This is a utility routine to tell whether to byte swap MRC header. */
static void check_swap(const int * data, const char * filnam, bool show_errors,
......@@ -304,6 +305,9 @@ namespace EMAN
//utility funciton to tranpose x and y dimension in case the source mrc image is mapc=2,mapr=1
int transpose(float *data, int nx, int ny, int nz) const;
template<class T>
auto write_compressed(float *data, size_t size, int image_index, const Region* area);
};
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment