c++ - how to save a double matrix quickly and compactly if some loss of precision is acceptable? -
i have matrix of double size 1024x1024. want write file. here solution this. accepted answer 50-60% more time efficient second answer in question expected (according simple test write file in both methods).
there solution write csv file(accepted answer in question), slower(3-4 times slower)
while i'm writing file number of floating point of each element in matrix 16 , output following:
-1.6819883882999420e-001 -3.5269531607627869e-001 2.4137189984321594e-001 -3.9325976371765137e-001 -2.2069962322711945e-001 -5.9525445103645325e-002
when write file, file size becomes 24 mb in first way(first link, accepted answer) , 37 mb in third way(second link, accepted answer) both of unacceptable.
i need set precision of matrix in fast way , output become -1.6819e-01 -3.5269e-01. appreciated.
what i'm doing read image of 1024x1024, process it, write output mat(which double) file. consider have, thousands of image, images less 1 mb, running time each image less 1 sec without saving.
edit: when save same matrix in matlab, becomes 6.75 mb
edit
- now has option store in 16bit floats.
- you lose precision(depends on data range,vary rnd distribution see different error rates).
- slower(20-30ms)(see link below faster tricks if matters)
- 2mb
- if know nature of data(range,distribution surface(maybe variable length encoding)) more can done
- see here(32-bit 16-bit floating point conversion).
code
//using code lifted http://www.mathworks.com/matlabcentral/fileexchange/23173 #include "opencv2/opencv.hpp" #include <string.h> #include <math.h> using namespace cv; #define int16_type short #define uint16_type unsigned short #define int32_type long #define uint32_type unsigned long int doubles2halfp(void *target, void *source, int numel); int halfp2doubles(void *target, void *source, int numel); void writemat(char* fpath,mat& data,bool isf16) { file* fp = fopen(fpath,"wb"); if (!fp)perror("fopen"); double dbuf[1024]; for(int i=0;i<1024;++i) { for(int j=0;j<1024;++j) dbuf[j]=data.at<double>(i,j); if(isf16) { uint16_type hbuf[1024]; doubles2halfp(&hbuf,&dbuf,1024); fwrite(&hbuf,sizeof(uint16_type),1024,fp); }else { fwrite(&dbuf,sizeof(double),1024,fp); } } fclose(fp); } void readmat(char* fpath,mat& data,bool isf16) { file* fp = fopen(fpath,"rb"); if (!fp)perror("fopen"); double dbuf[1024]; for(int i=0;i<1024;++i) { if(isf16) { uint16_type hbuf[1024]; fread(&hbuf,sizeof(uint16_type),1024,fp); halfp2doubles(&dbuf,&hbuf,1024); }else{ fread(&dbuf,sizeof(double),1024,fp); } for(int j=0;j<1024;++j) { data.at<double>(i,j)=dbuf[j]; } } fclose(fp); } int main() { rng rng= therng(); mat data = mat::zeros(size(1024,1024),cv_64fc1); for(int i=0;i<1024;++i) for(int j=0;j<1024;++j) data.at<double>(i,j)=rng.uniform(-1.,1.); writemat("img.bin",data,true); mat res = mat::zeros(size(1024,1024),cv_64fc1); readmat("img.bin",res,true); double error=0; for(int i=0;i<1024;++i) for(int j=0;j<1024;++j) { //printf("%f %f\n",data.at<double>(i,j),res.at<double>(i,j)); error+=abs(data.at<double>(i,j)-res.at<double>(i,j)); } printf("err=%f avgerr=%f\n",error,error/1024/1024); getchar(); return 0; } //////////////////////////////////////// //////////////////////////////////////// //////////////////////////////////////// //////////////////////////////////////// int doubles2halfp(void *target, void *source, int numel) { uint16_type *hp = (uint16_type *) target; // type pun output unsigned 16-bit int uint32_type *xp = (uint32_type *) source; // type pun input unsigned 32-bit int uint16_type hs, he, hm; uint32_type x, xs, xe, xm; int hes; static int next; // little endian adjustment static int checkieee = 1; // flag check ieee754, endian, , word size double 1 = 1.0; // used checking ieee754 floating point format uint32_type *ip; // used checking ieee754 floating point format if( checkieee ) { // 1st call, check ieee754, endian, , word size ip = (uint32_type *) &one; if( *ip ) { // if big endian, no adjustment next = 0; } else { // if little endian, adjustment necessary next = 1; ip++; } if( *ip != 0x3ff00000u ) { // check exact ieee 754 bit pattern of 1.0 return 1; // floating point bit pattern not ieee 754 } if( sizeof(int16_type) != 2 || sizeof(int32_type) != 4 ) { return 1; // short not 16-bits, or long not 32-bits. } checkieee = 0; // checks out ok } xp += next; // little endian adjustment if necessary if( source == null || target == null ) { // nothing convert (e.g., imag part of pure real) return 0; } while( numel-- ) { x = *xp++; xp++; // xp++ skip on remaining 32 bits of mantissa if( (x & 0x7fffffffu) == 0 ) { // signed 0 *hp++ = (uint16_type) (x >> 16); // return signed 0 } else { // not 0 xs = x & 0x80000000u; // pick off sign bit xe = x & 0x7ff00000u; // pick off exponent bits xm = x & 0x000fffffu; // pick off mantissa bits if( xe == 0 ) { // denormal underflow, return signed 0 *hp++ = (uint16_type) (xs >> 16); } else if( xe == 0x7ff00000u ) { // inf or nan (all exponent bits set) if( xm == 0 ) { // if mantissa 0 ... *hp++ = (uint16_type) ((xs >> 16) | 0x7c00u); // signed inf } else { *hp++ = (uint16_type) 0xfe00u; // nan, 1st mantissa bit set } } else { // normalized number hs = (uint16_type) (xs >> 16); // sign bit hes = ((int)(xe >> 20)) - 1023 + 15; // exponent unbias double, bias halfp if( hes >= 0x1f ) { // overflow *hp++ = (uint16_type) ((xs >> 16) | 0x7c00u); // signed inf } else if( hes <= 0 ) { // underflow if( (10 - hes) > 21 ) { // mantissa shifted way off & no rounding possibility hm = (uint16_type) 0u; // set mantissa 0 } else { xm |= 0x00100000u; // add hidden leading bit hm = (uint16_type) (xm >> (11 - hes)); // mantissa if( (xm >> (10 - hes)) & 0x00000001u ) // check rounding hm += (uint16_type) 1u; // round, might overflow exp bit, ok } *hp++ = (hs | hm); // combine sign bit , mantissa bits, biased exponent 0 } else { = (uint16_type) (hes << 10); // exponent hm = (uint16_type) (xm >> 10); // mantissa if( xm & 0x00000200u ) // check rounding *hp++ = (hs | | hm) + (uint16_type) 1u; // round, might overflow inf, ok else *hp++ = (hs | | hm); // no rounding } } } } return 0; } int halfp2doubles(void *target, void *source, int numel) { uint16_type *hp = (uint16_type *) source; // type pun input unsigned 16-bit int uint32_type *xp = (uint32_type *) target; // type pun output unsigned 32-bit int uint16_type h, hs, he, hm; uint32_type xs, xe, xm; int32_type xes; int e; static int next; // little endian adjustment static int checkieee = 1; // flag check ieee754, endian, , word size double 1 = 1.0; // used checking ieee754 floating point format uint32_type *ip; // used checking ieee754 floating point format if( checkieee ) { // 1st call, check ieee754, endian, , word size ip = (uint32_type *) &one; if( *ip ) { // if big endian, no adjustment next = 0; } else { // if little endian, adjustment necessary next = 1; ip++; } if( *ip != 0x3ff00000u ) { // check exact ieee 754 bit pattern of 1.0 return 1; // floating point bit pattern not ieee 754 } if( sizeof(int16_type) != 2 || sizeof(int32_type) != 4 ) { return 1; // short not 16-bits, or long not 32-bits. } checkieee = 0; // checks out ok } xp += next; // little endian adjustment if necessary if( source == null || target == null ) // nothing convert (e.g., imag part of pure real) return 0; while( numel-- ) { h = *hp++; if( (h & 0x7fffu) == 0 ) { // signed 0 *xp++ = ((uint32_type) h) << 16; // return signed 0 } else { // not 0 hs = h & 0x8000u; // pick off sign bit = h & 0x7c00u; // pick off exponent bits hm = h & 0x03ffu; // pick off mantissa bits if( == 0 ) { // denormal convert normalized e = -1; // following loop figures out how adjust exponent { e++; hm <<= 1; } while( (hm & 0x0400u) == 0 ); // shift until leading bit overflows exponent bit xs = ((uint32_type) hs) << 16; // sign bit xes = ((int32_type) (he >> 10)) - 15 + 1023 - e; // exponent unbias halfp, bias double xe = (uint32_type) (xes << 20); // exponent xm = ((uint32_type) (hm & 0x03ffu)) << 10; // mantissa *xp++ = (xs | xe | xm); // combine sign bit, exponent bits, , mantissa bits } else if( == 0x7c00u ) { // inf or nan (all exponent bits set) if( hm == 0 ) { // if mantissa 0 ... *xp++ = (((uint32_type) hs) << 16) | ((uint32_type) 0x7ff00000u); // signed inf } else { *xp++ = (uint32_type) 0xfff80000u; // nan, 1st mantissa bit set } } else { xs = ((uint32_type) hs) << 16; // sign bit xes = ((int32_type) (he >> 10)) - 15 + 1023; // exponent unbias halfp, bias double xe = (uint32_type) (xes << 20); // exponent xm = ((uint32_type) hm) << 10; // mantissa *xp++ = (xs | xe | xm); // combine sign bit, exponent bits, , mantissa bits } } xp++; // skip on remaining 32 bits of mantissa } return 0; }
Comments
Post a Comment