Skip to content
GitLab
Commit cf66ab4f authored by shadow_walker's avatar shadow_walker
Browse files

Merge branch 'compression'

parents caa86b8a 2fb17a8a
Hide whitespace changes
Inline Side-by-side
libEM/io/df3io.cpp
View file @ cf66ab4f
......@@ -173,9 +173,6 @@ int Df3IO::write_data(float *data, int, const Region*,
ENTERFUNC;
size_t img_size = (size_t)nx*ny*nz;
unsigned int * uidata = 0;
unsigned short * usdata = 0;
unsigned char * ucdata = 0;
if(dt == EMUtil::EM_COMPRESSED) {
if (renderbits <= 8) dt = EMUtil::EM_UCHAR;
......@@ -186,66 +183,25 @@ int Df3IO::write_data(float *data, int, const Region*,
if (renderbits==0 || renderbits>truebits) renderbits=truebits;
EMUtil::getRenderMinMax(data, nx, ny, rendermin, rendermax, renderbits, nz);
switch(dt) {
case EMUtil::EM_UINT:
uidata = new unsigned int[img_size];
for (size_t i = 0; i < img_size; ++i) {
if(data[i] <= rendermin) {
uidata[i] = 0;
}
else if(data[i] >= rendermax) {
uidata[i] = UINT_MAX;
}
else {
uidata[i]=(unsigned int)((data[i]-rendermin)/(rendermax-rendermin)*UINT_MAX);
}
}
ByteOrder::become_big_endian(uidata, img_size);
if(fwrite(uidata, sizeof(unsigned int), img_size, file) != img_size) {
if(dt == EMUtil::EM_UINT) {
auto [rendered_data, count] = getRenderedDataAndRendertrunc<unsigned int>(data, img_size);
ByteOrder::become_big_endian(rendered_data.data(), img_size);
if(fwrite(rendered_data.data(), sizeof(unsigned int), img_size, file) != img_size)
throw ImageWriteException(filename, "DF3 unsigned int data");
}
if(uidata) {delete [] uidata; uidata=0;}
break;
case EMUtil::EM_USHORT:
usdata = new unsigned short[img_size];
for (size_t i = 0; i < img_size; ++i) {
if(data[i] <= rendermin) {
usdata[i] = 0;
}
else if(data[i] >= rendermax) {
usdata[i] = USHRT_MAX;
}
else {
usdata[i]=(unsigned short)((data[i]-rendermin)/(rendermax-rendermin)*USHRT_MAX);
}
}
ByteOrder::become_big_endian(usdata, img_size);
if(fwrite(usdata, sizeof(unsigned short), img_size, file) != img_size) {
}
else if(dt == EMUtil::EM_USHORT) {
auto [rendered_data, count] = getRenderedDataAndRendertrunc<unsigned short>(data, img_size);
ByteOrder::become_big_endian(rendered_data.data(), img_size);
if(fwrite(rendered_data.data(), sizeof(unsigned short), img_size, file) != img_size)
throw ImageWriteException(filename, "DF3 unsigned short data");
}
if(usdata) {delete [] usdata; usdata=0;}
break;
case EMUtil::EM_UCHAR:
ucdata = new unsigned char[img_size];
for (size_t i = 0; i < img_size; ++i) {
if(data[i] <= rendermin) {
ucdata[i] = 0;
}
else if(data[i] >= rendermax){
ucdata[i] = UCHAR_MAX;
}
else {
ucdata[i]=(unsigned char)((data[i]-rendermin)/(rendermax-rendermin)*UCHAR_MAX);
}
}
if(fwrite(ucdata, sizeof(unsigned char), img_size, file) != img_size) {
}
else if(dt == EMUtil::EM_UCHAR) {
auto [rendered_data, count] = getRenderedDataAndRendertrunc<unsigned char>(data, img_size);
if(fwrite(rendered_data.data(), sizeof(unsigned char), img_size, file) != img_size)
throw ImageWriteException(filename, "DF3 unsigned char data");
}
if(ucdata) {delete [] ucdata; ucdata=0;}
break;
default:
throw ImageWriteException(filename,"DF3 does not support this data format");
}
else
throw ImageWriteException(filename,"DF3 does not support this data format");
EXITFUNC;
return 0;
......
libEM/io/jpegio.cpp
View file @ cf66ab4f
......@@ -169,7 +169,7 @@ int JpegIO::read_data(float *, int, const Region *, bool)
}
int JpegIO::write_data(float *data, int image_index, const Region* area,
EMUtil::EMDataType, bool)
EMUtil::EMDataType dt, bool)
{
ENTERFUNC;
......@@ -188,6 +188,7 @@ int JpegIO::write_data(float *data, int image_index, const Region* area,
if (renderbits==0 || renderbits>8) renderbits=8;
EMUtil::getRenderMinMax(data, nx, ny, rendermin, rendermax, renderbits);
auto [rendered_data, count] = getRenderedDataAndRendertrunc<unsigned char>(data, nx*ny);
unsigned char *cdata = (unsigned char *)malloc(nx+1);
/* Flip the image vertically, since EMAN use top-left corner as image origin
......@@ -201,10 +202,7 @@ int JpegIO::write_data(float *data, int image_index, const Region* area,
for (int i=ny-1; i>=0; i--) {
for (int j=0; j<nx; j++) {
if (data[i*nx+j] <= rendermin) cdata[j] = 0;
else if (data[i*nx+j] >= rendermax) cdata[j] = 255;
else cdata[j]= (int)((data[i*nx+j]-rendermin)/
(rendermax-rendermin)*256.0);
cdata[j]= (int)rendered_data[i*nx+j];
}
jpeg_write_scanlines(&cinfo, rp, 1);
......
libEM/io/mrcio.cpp
View file @ cf66ab4f
......@@ -1106,309 +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);
// int xlen = 0, ylen = 0, zlen = 0;
// EMUtil::get_region_dims(area, nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
// int size = xlen * ylen * zlen;
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;
}
// New way to write data which includes region writing.
// If it is tested to be OK, remove the old code in the
// #if 0 ... #endif block.
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;}
#if 0
int row_size = nx * get_mode_size(mrch.mode);
int sec_size = nx * ny;
unsigned char * cbuf = new unsigned char[row_size];
unsigned short * sbuf = (unsigned short *) cbuf;
for (int i = 0; i < nz; i++) {
int i2 = i * sec_size;
for (int j = 0; j < ny; j++) {
int k = i2 + j * nx;
void * pbuf = 0;
switch (mrch.mode) {
case MRC_UCHAR:
for (int l = 0; l < nx; l++) {
cbuf[l] = static_cast < unsigned char >(data[k + l]);
}
pbuf = cbuf;
fwrite(cbuf, row_size, 1, file);
break;
case MRC_SHORT:
case MRC_SHORT_COMPLEX:
for (int l = 0; l < nx; l++) {
sbuf[l] = static_cast < short >(data[k + l]);
}
pbuf = sbuf;
fwrite(sbuf, row_size, 1, file);
break;
case MRC_USHORT:
for (int l = 0; l < nx; l++) {
sbuf[l] = static_cast < unsigned short >(data[k + l]);
}
pbuf = sbuf;
fwrite(sbuf, row_size, 1, file);
break;
case MRC_FLOAT:
case MRC_FLOAT_COMPLEX:
pbuf = &data[k];
break;
}
if (pbuf) {
fwrite(pbuf, row_size, 1, file);
}
}
}
if (cbuf)
{
delete [] cbuf;
cbuf = NULL;
}
#endif
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;
......@@ -1417,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;
......@@ -1427,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) {
......@@ -1474,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;
......@@ -1502,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 */
......@@ -1515,31 +1211,134 @@ void MrcIO::update_stats(void * data, size_t size)
mrch.amax = max;
mrch.amean = (float) mean;
mrch.rms = (float) sigma;
// MrcHeader mrch2 = mrch;
//
// endian issue, can't get use_host_endian argument
// bool opposite_endian = false;
//
// if (!is_new_file) {
// if (is_big_endian != ByteOrder::is_host_big_endian()) {
// opposite_endian = true;
// }
//
// portable_fseek(mrcfile, 0, SEEK_SET);
// }
//
// if (opposite_endian || !use_host_endian) {
// swap_header(mrch2);
// }
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;
}