diff --git a/stim/envi/bil.h b/stim/envi/bil.h index 007356f..1cc2646 100644 --- a/stim/envi/bil.h +++ b/stim/envi/bil.h @@ -2,7 +2,7 @@ #define STIM_BIL_H #include "../envi/envi_header.h" -#include "../envi/binary.h" +#include "../envi/hsi.h" #include #include @@ -19,28 +19,29 @@ namespace stim{ template -class bil: public binary { +class bil: public hsi { protected: - std::vector w; //band wavelength + //std::vector w; //band wavelengths - unsigned long X(){ + /*unsigned long long X(){ return R[0]; } - unsigned long Y(){ + unsigned long long Y(){ return R[2]; } - unsigned long Z(){ + unsigned long long Z(){ return R[1]; - } - + }*/ + using hsi::w; //use the wavelength array in stim::hsi + using hsi::nnz; using binary::progress; /// Call the binary nnz() function for the BIL orientation - unsigned long long nnz(unsigned char* mask){ - return binary::nnz(mask, X()*Y()); - } + //unsigned long long nnz(unsigned char* mask){ + // return binary::nnz(mask, X()*Y()); + //} public: @@ -48,6 +49,8 @@ public: using binary::file; using binary::R; + bil(){ hsi::init_bil(); } + /// Open a data file for reading using the class interface. /// @param filename is the name of the binary file on disk @@ -56,11 +59,16 @@ public: /// @param B is the number of samples (bands) along dimension 3 /// @param header_offset is the number of bytes (if any) in the binary header /// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band - bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int B, unsigned int header_offset, std::vector wavelengths){ + bool open(std::string filename, + unsigned long long X, + unsigned long long Y, + unsigned long long B, + unsigned long long header_offset, + std::vector wavelengths){ w = wavelengths; - return open(filename, vec(X, B, Y), header_offset); + return open(filename, vec(X, B, Y), header_offset); } @@ -68,12 +76,12 @@ public: /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param page <= B is the integer number of the band to be copied. - bool band_index( T * p, unsigned int page, bool PROGRESS = false){ + bool band_index( T * p, unsigned long long page, bool PROGRESS = false){ //return binary::read_plane_1(p, page); if(PROGRESS) progress = 0; - unsigned int L = X() * sizeof(T); //caculate the number of bytes in a sample line - unsigned int jump = X() * (Z() - 1) * sizeof(T); + unsigned long long L = X() * sizeof(T); //caculate the number of bytes in a sample line + unsigned long long jump = X() * (Z() - 1) * sizeof(T); if (page >= Z()){ //make sure the bank number is right std::cout<<"ERROR: page out of range"<::read_line_1(p, x, y, PROGRESS); } @@ -210,18 +218,18 @@ public: /// @param p is a pointer to pre-allocated memory at least sizeof(T) in size. /// @param n is an integer index to the pixel using linear array indexing. - bool pixel(T * p, unsigned n){ + bool pixel(T * p, unsigned long long n){ //calculate the corresponding x, y - unsigned int x = n % X(); - unsigned int y = n / X(); + unsigned long long x = n % X(); + unsigned long long y = n / X(); //get the pixel return spectrum(p, x, y); } //given a Y ,return a XZ slice - bool read_plane_y(T * p, unsigned y){ + bool read_plane_y(T * p, unsigned long long y){ return binary::read_plane_2(p, y); } @@ -232,15 +240,15 @@ public: /// @param wls is the list of baseline points based on band labels. bool baseline(std::string outname, std::vector wls, unsigned char* mask = NULL, bool PROGRESS = false){ - unsigned N = wls.size(); //get the number of baseline points + unsigned long long N = wls.size(); //get the number of baseline points std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file std::string headername = outname + ".hdr"; //the header file name //simplify image resolution - unsigned int ZX = Z() * X(); //calculate the number of points in a Y slice - unsigned int L = ZX * sizeof(T); //calculate the number of bytes of a Y slice - unsigned int B = Z(); + unsigned long long ZX = Z() * X(); //calculate the number of points in a Y slice + unsigned long long L = ZX * sizeof(T); //calculate the number of bytes of a Y slice + unsigned long long B = Z(); T* c; //pointer to the current Y slice c = (T*)malloc(L); //memory allocation @@ -254,7 +262,7 @@ public: double ai, bi; //stores the two baseline points wavelength surrounding the current band double ci; //stores the current band's wavelength - unsigned control; + unsigned long long control; if (a == NULL || b == NULL || c == NULL){ std::cout<<"ERROR: error allocating memory"; @@ -262,7 +270,7 @@ public: } // loop start correct every y slice - for (unsigned k =0; k < Y(); k++) + for (unsigned long long k =0; k < Y(); k++) { //get the current y slice read_plane_y(c, k); @@ -287,7 +295,7 @@ public: //correct every YZ line - for(unsigned cii = 0; cii < B; cii++){ + for(unsigned long long cii = 0; cii < B; cii++){ //update baseline points, if necessary if( w[cii] >= bi && cii != B - 1) { @@ -317,7 +325,7 @@ public: ci = w[cii]; - unsigned jump = cii * X(); + unsigned long long jump = cii * X(); //perform the baseline correction for(unsigned i=0; i < X(); i++) { @@ -347,11 +355,11 @@ public: /// @param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers. bool normalize(std::string outname, double w, unsigned char* mask = NULL, bool PROGRESS = false) { - unsigned int B = Z(); //calculate the number of bands - unsigned int ZX = Z() * X(); - unsigned int XY = X() * Y(); //calculate the number of pixels in a band - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band - unsigned int L = ZX * sizeof(T); + unsigned long long B = Z(); //calculate the number of bands + unsigned long long ZX = Z() * X(); + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band + unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band + unsigned long long L = ZX * sizeof(T); std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file std::string headername = outname + ".hdr"; //the header file name @@ -364,12 +372,12 @@ public: band(b, w); //get the certain band into memory - for(unsigned j = 0; j < Y(); j++) + for(unsigned long long j = 0; j < Y(); j++) { read_plane_y(c, j); - for(unsigned i = 0; i < B; i++) + for(unsigned long long i = 0; i < B; i++) { - for(unsigned m = 0; m < X(); m++) + for(unsigned long long m = 0; m < X(); m++) { if( mask != NULL && !mask[m + j * X()] ) c[m + i * X()] = (T)0.0; @@ -393,7 +401,7 @@ public: /// @param outname is the name of the output BSQ file to be saved to disk. bool bsq(std::string outname, bool PROGRESS = false) { - unsigned int S = X() * Y() * sizeof(T); //calculate the number of bytes in a band + unsigned long long S = X() * Y() * sizeof(T); //calculate the number of bytes in a band std::ofstream target(outname.c_str(), std::ios::binary); std::string headername = outname + ".hdr"; @@ -401,7 +409,7 @@ public: T * p; //pointer to the current band p = (T*)malloc(S); - for ( unsigned i = 0; i < Z(); i++) + for ( unsigned long long i = 0; i < Z(); i++) { band_index(p, i); target.write(reinterpret_cast(p), S); //write a band data into target file @@ -429,13 +437,13 @@ public: T * q; //pointer to the current ZX slice for bip file q = (T*)malloc(S); - for ( unsigned i = 0; i < Y(); i++) + for ( unsigned long long i = 0; i < Y(); i++) { read_plane_y(p, i); - for ( unsigned k = 0; k < Z(); k++) + for ( unsigned long long k = 0; k < Z(); k++) { - unsigned ks = k * X(); - for ( unsigned j = 0; j < X(); j++) + unsigned long long ks = k * X(); + for ( unsigned long long j = 0; j < X(); j++) q[k + j * Z()] = p[ks + j]; if(PROGRESS) progress = (double)((i+1) * Z() + k+1) / (Z() * Y()) * 100; //store the progress for the current operation @@ -461,12 +469,12 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size. bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){ - unsigned XY = X() * Y(); + unsigned long long XY = X() * Y(); band(result, wavelength); //get band //perform the baseline correction double r = (double) (wavelength - lb) / (double) (rb - lb); - for(unsigned i=0; i < XY; i++){ + for(unsigned long long i=0; i < XY; i++){ result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] ); } return true; @@ -482,8 +490,8 @@ public: T* lp; T* rp; - unsigned XY = X() * Y(); - unsigned S = XY * sizeof(T); + unsigned long long XY = X() * Y(); + unsigned long long S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation rp = (T*) malloc(S); @@ -513,8 +521,8 @@ public: T* cur; //current band 1 T* cur2; //current band 2 - unsigned XY = X() * Y(); - unsigned S = XY * sizeof(T); + unsigned long long XY = X() * Y(); + unsigned long long S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation rp = (T*) malloc(S); @@ -524,9 +532,9 @@ public: memset(result, (char)0, S); //find the wavelenght position in the whole band - unsigned int n = w.size(); - unsigned int ai = 0; //left bound position - unsigned int bi = n - 1; //right bound position + unsigned long long n = w.size(); + unsigned long long ai = 0; //left bound position + unsigned long long bi = n - 1; //right bound position @@ -555,23 +563,23 @@ public: //calculate the beginning and the ending part baseline_band(lb, rb, lp, rp, rab, cur2); //ending part baseline_band(lb, rb, lp, rp, w[bi], cur); - for(unsigned j = 0; j < XY; j++){ - result[j] += (rab - w[bi]) * (cur[j] + cur2[j]) / 2.0; + for(unsigned long long j = 0; j < XY; j++){ + result[j] += (T)((rab - w[bi]) * ((double)cur[j] + (double)cur2[j]) / 2.0); } baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part baseline_band(lb, rb, lp, rp, w[ai], cur); - for(unsigned j = 0; j < XY; j++){ - result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0; + for(unsigned long long j = 0; j < XY; j++){ + result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0); } //calculate the area ai++; - for(unsigned i = ai; i <= bi ;i++) + for(unsigned long long i = ai; i <= bi ;i++) { baseline_band(lb, rb, lp, rp, w[ai], cur2); - for(unsigned j = 0; j < XY; j++) + for(unsigned long long j = 0; j < XY; j++) { - result[j] += (w[ai] - w[ai-1]) * (cur[j] + cur2[j]) / 2.0; + result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0); } std::swap(cur,cur2); //swap the band pointers } @@ -601,7 +609,7 @@ public: height(lb1, rb1, pos1, p1); height(lb2, rb2, pos2, p2); //calculate the ratio in result - for(unsigned i = 0; i < X() * Y(); i++){ + for(unsigned long long i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -632,7 +640,7 @@ public: area(lb1, rb1, lab1, rab1, p1); height(lb2, rb2, pos, p2); //calculate the ratio in result - for(unsigned i = 0; i < X() * Y(); i++){ + for(unsigned long long i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -665,7 +673,7 @@ public: area(lb1, rb1, lab1, rab1, p1); area(lb2, rb2, lab2, rab2, p2); //calculate the ratio in result - for(unsigned i = 0; i < X() * Y(); i++){ + for(unsigned long long i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -690,8 +698,8 @@ public: T* cur; //current band 1 T* cur2; //current band 2 - unsigned XY = X() * Y(); - unsigned S = XY * sizeof(T); + unsigned long long XY = X() * Y(); + unsigned long long S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation rp = (T*) malloc(S); @@ -701,9 +709,9 @@ public: memset(result, (char)0, S); //find the wavelenght position in the whole band - unsigned int n = w.size(); - unsigned int ai = 0; //left bound position - unsigned int bi = n - 1; //right bound position + unsigned long long n = w.size(); + unsigned long long ai = 0; //left bound position + unsigned long long bi = n - 1; //right bound position //to make sure the left and the right bound are in the bandwidth if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){ @@ -730,23 +738,23 @@ public: //calculate the beginning and the ending part baseline_band(lb, rb, lp, rp, rab, cur2); //ending part baseline_band(lb, rb, lp, rp, w[bi], cur); - for(unsigned j = 0; j < XY; j++){ - result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0; + for(unsigned long long j = 0; j < XY; j++){ + result[j] += (T)((rab - w[bi]) * (rab + w[bi]) * ((double)cur[j] + (double)cur2[j]) / 4.0); } baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part baseline_band(lb, rb, lp, rp, w[ai], cur); - for(unsigned j = 0; j < XY; j++){ - result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0; + for(unsigned long long j = 0; j < XY; j++){ + result[j] += (T)((w[ai] - lab) * (w[ai] + lab) * ((double)cur[j] + (double)cur2[j]) / 4.0); } //calculate f(x) times x ai++; - for(unsigned i = ai; i <= bi ;i++) + for(unsigned long long i = ai; i <= bi ;i++) { baseline_band(lb, rb, lp, rp, w[ai], cur2); - for(unsigned j = 0; j < XY; j++) + for(unsigned long long j = 0; j < XY; j++) { - result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0; + result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0); } std::swap(cur,cur2); //swap the band pointers } @@ -773,7 +781,7 @@ public: x_area(lb, rb, lab, rab, p1); area(lb, rb, lab, rab, p2); //calculate the ratio in result - for(unsigned i = 0; i < X() * Y(); i++){ + for(unsigned long long i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -797,7 +805,7 @@ public: T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band band(temp, mask_band, PROGRESS); - for (unsigned int i = 0; i < X() * Y(); i++) { + for (unsigned long long i = 0; i < X() * Y(); i++) { if (temp[i] < threshold) p[i] = 0; else @@ -817,17 +825,17 @@ public: std::ofstream target(outfile.c_str(), std::ios::binary); //I THINK THIS IS WRONG - unsigned XZ = X() * Z(); //calculate the number of values in a page on disk - unsigned L = XZ * sizeof(T); //calculate the size of the page (in bytes) + unsigned long long XZ = X() * Z(); //calculate the number of values in a page on disk + unsigned long long L = XZ * sizeof(T); //calculate the size of the page (in bytes) T * temp = (T*)malloc(L); //allocate memory for a temporary page - for (unsigned i = 0; i < Y(); i++) //for each value in Y() (BIP should be X) + for (unsigned long long i = 0; i < Y(); i++) //for each value in Y() (BIP should be X) { read_plane_y(temp, i); //retrieve an ZX slice, stored in temp - for ( unsigned j = 0; j < Z(); j++) //for each Z() (Y) + for ( unsigned long long j = 0; j < Z(); j++) //for each Z() (Y) { - for (unsigned k = 0; k < X(); k++) //for each band + for (unsigned long long k = 0; k < X(); k++) //for each band { if(p[i * X() + k] == 0) temp[j * X() + k] = 0; @@ -849,9 +857,9 @@ public: std::ofstream target(outfile.c_str(), std::ios::binary); //for loading pages: - unsigned long XZ = X() * Z(); //calculate the number of values in an XZ page on disk - unsigned long B = Z(); //calculate the number of bands - unsigned long L = XZ * sizeof(T); //calculate the size of the page (in bytes) + unsigned long long XZ = X() * Z(); //calculate the number of values in an XZ page on disk + unsigned long long B = Z(); //calculate the number of bands + unsigned long long L = XZ * sizeof(T); //calculate the size of the page (in bytes) //allocate temporary memory for a XZ slice T* slice = (T*) malloc(L); @@ -860,18 +868,18 @@ public: T* spec = (T*) malloc(B * sizeof(T)); //for each slice along the y axis - for (unsigned long y = 0; y < Y(); y++) //Select a page by choosing Y coordinate, Y() + for (unsigned long long y = 0; y < Y(); y++) //Select a page by choosing Y coordinate, Y() { read_plane_y(slice, y); //retrieve an ZX page, store in "slice" //for each sample along X - for (unsigned long x = 0; x < X(); x++) //Select a pixel by choosing X coordinate in the page, X() + for (unsigned long long x = 0; x < X(); x++) //Select a pixel by choosing X coordinate in the page, X() { //if the mask != 0 at that xy pixel if (p[y * X() + x] != 0) //if the mask != 0 at that XY pixel { //for each band at that pixel - for (unsigned long b = 0; b < B; b++) //Select a voxel by choosing Z coordinate at the pixel + for (unsigned long long b = 0; b < B; b++) //Select a voxel by choosing Z coordinate at the pixel { spec[b] = slice[b*X() + x]; //Pass the correct spectral value from XZ page into the spectrum to be saved. } @@ -896,20 +904,20 @@ public: T* temp = (T*)malloc(sizeof(T) * XZ); T* line = (T*)malloc(sizeof(T) * X()); - for (unsigned i = 0; i < Y(); i++){ + for (unsigned long long i = 0; i < Y(); i++){ getY(temp, i); //initialize x-line - for (unsigned j = 0; j < X(); j++){ + for (unsigned long long j = 0; j < X(); j++){ line[j] = 0; } - unsigned c = 0; - for (unsigned j = 0; j < Z(); j++){ - for (unsigned k = 0; k < X(); k++){ + unsigned long long c = 0; + for (unsigned long long j = 0; j < Z(); j++){ + for (unsigned long long k = 0; k < X(); k++){ line[k] += temp[c] / (T)Z(); c++; } } - for (unsigned j = 0; j < X(); j++){ + for (unsigned long long j = 0; j < X(); j++){ p[j + i * X()] = line[j]; } } @@ -925,7 +933,7 @@ public: unsigned long long XZ = X() * Z(); unsigned long long XY = X() * Y(); T* temp = (T*)malloc(sizeof(T) * XZ); - for (unsigned j = 0; j < Z(); j++){ + for (unsigned long long j = 0; j < Z(); j++){ p[j] = 0; } //calculate vaild number in a band @@ -1038,9 +1046,9 @@ public: //set the start position for the cropped region file.seekg((y0 * X() * Z() + b0 * X() + x0) * sizeof(T), std::ios::beg); - for (unsigned x = 0; x < lines; x++) + for (unsigned long long x = 0; x < lines; x++) { - for (unsigned z = b0; z < b1; z++) + for (unsigned long long z = b0; z < b1; z++) { file.read((char *)(temp + z * samples), sizeof(T) * samples); file.seekg(jumpb, std::ios::cur); //go to the next band diff --git a/stim/envi/binary.h b/stim/envi/binary.h index e010b50..b597745 100644 --- a/stim/envi/binary.h +++ b/stim/envi/binary.h @@ -25,7 +25,7 @@ protected: std::string name; //file name unsigned long long R[D]; //resolution - unsigned int header; //header size (in bytes) + unsigned long long header; //header size (in bytes) unsigned char* mask; //pointer to a character array: 0 = background, 1 = foreground (or valid data) double progress; //stores the progress on the current operation (accessible using a thread) @@ -33,30 +33,13 @@ protected: /// Private initialization function used to set default parameters in the data structure. void init(){ - memset(R, 0, sizeof(unsigned int) * D); //initialize the resolution to zero + memset(R, 0, sizeof(unsigned long long) * D); //initialize the resolution to zero header = 0; //initialize the header size to zero mask = NULL; progress = 0; } - //calculate the number of non-zero pixels in a mask - unsigned long long nnz(unsigned char* mask, unsigned long long N){ - - unsigned long long n = 0; //initialize the counter to 0 (zero) - if(mask == NULL) return N; //if the mask is NULL, assume all pixels are masked - - for(unsigned long long i = 0; i < N; i++){ //for each pixel - if(mask[i] != 0) n++; //increment the counter for every non-zero pixel in the mask - } - return n; //return the number of nonzero pixels - } - - /// Calculate the number of nonzero pixels in a mask over X-Y - unsigned long long nnz(unsigned char* mask){ - return nnz(mask, R[0] * R[1]); - } - /// Private helper function that returns the size of the file on disk using system functions. long long int get_file_size(){ #ifdef _WIN32 @@ -122,7 +105,7 @@ protected: public: - unsigned int get_progress(){ + double get_progress(){ return progress; } @@ -135,9 +118,9 @@ public: /// @param filename is the name of the binary file /// @param r is a STIM vector specifying the size of the binary file along each dimension /// @param h is the length (in bytes) of any header file (default zero) - bool open(std::string filename, vec r, unsigned int h = 0){ + bool open(std::string filename, vec r, unsigned long long h = 0){ - for(unsigned int i = 0; i < D; i++) //set the dimensions of the binary file object + for(unsigned long long i = 0; i < D; i++) //set the dimensions of the binary file object R[i] = r[i]; header = h; //save the header size @@ -154,17 +137,17 @@ public: /// @param filename is the name of the binary file to be created /// @param r is a STIM vector specifying the size of the file along each dimension /// @offset specifies how many bytes to offset the file (used to leave room for a header) - bool create(std::string filename, vec r, unsigned int offset = 0){ + bool create(std::string filename, vec r, unsigned long long offset = 0){ std::ofstream target(filename.c_str(), std::ios::binary); //initialize binary file T p = 0; - for(unsigned int i =0; i < r[0] * r[1] * r[2]; i++){ + for(unsigned long long i =0; i < r[0] * r[1] * r[2]; i++){ target.write((char*)(&p), sizeof(T)); } - for(unsigned int i = 0; i < D; i++) //set the dimensions of the binary file object + for(unsigned long long i = 0; i < D; i++) //set the dimensions of the binary file object R[i] = r[i]; header = offset; //save the header size @@ -178,7 +161,7 @@ public: /// @param p is a pointer to the data to be written /// @param page is the page number (index of the highest-numbered dimension) - bool write_page( T * p, unsigned int page){ + bool write_page( T * p, unsigned long long page){ if(p == NULL){ std::cout<<"ERROR: unable to write into file, empty pointer"<= R[1] || z >= R[2] ){ std::cout<<"ERROR: sample or line out of range"<= R[0] || z >= R[2] ){ @@ -271,7 +254,7 @@ public: } file.seekg((z * R[0] * R[1] + x) * sizeof(T), std::ios::beg); //seek to the start of the line - for (unsigned int i = 0; i < R[1]; i++){ //for each pixel in the line + for (unsigned long long i = 0; i < R[1]; i++){ //for each pixel in the line file.read((char *)(p + i), sizeof(T)); //read the pixel file.seekg((R[0] - 1) * sizeof(T), std::ios::cur); //seek to the next pixel in the line if(PROGRESS) progress = (double)i / (double)R[1] * 100; @@ -284,22 +267,22 @@ public: /// @param p is a pointer to pre-allocated memory of size R[1] * R[2] * sizeof(T) /// @param n is the 0-axis coordinate used to retrieve the plane - bool read_plane_0(T* p, unsigned int n, bool PROGRESS = false){ + bool read_plane_0(T* p, unsigned long long n, bool PROGRESS = false){ if(PROGRESS) progress = 0; if (n >= R[0]){ //make sure the number is within the possible range std::cout<<"ERROR read_plane_0: page out of range"<= R[1]){ //make sure the bank number is right std::cout<<"ERROR read_plane_1: page out of range"<= R[0] * R[1] * R[2]){ std::cout<<"ERROR read_pixel: n is out of range"<= R[0] || y < 0 || y >= R[1] || z < 0 || z > R[2]){ std::cout<<"ERROR read_pixel: (x,y,z) is out of range"< #include @@ -23,15 +23,15 @@ namespace stim{ */ template -class bip: public binary { +class bip: public hsi { protected: - std::vector w; //band wavelength - unsigned int offset; //header offset + //std::vector w; //band wavelength + unsigned long long offset; //header offset - unsigned long long X(){ + /*unsigned long long X(){ return R[1]; } unsigned long long Y(){ @@ -39,17 +39,18 @@ protected: } unsigned long long Z(){ return R[0]; - } - + }*/ + using hsi::w; //use the wavelength array in stim::hsi + using hsi::nnz; using binary::progress; /// Call the binary nnz() function for the BIP orientation - unsigned long long nnz(unsigned char* mask){ - return binary::nnz(mask, X()*Y()); - } + //unsigned long long nnz(unsigned char* mask){ + // return binary::nnz(mask, X()*Y()); + //} /// Linear interpolation of a spectral value given the bounding spectral samples - T lerp(double w, T low_v, double low_w, T high_v, double high_w){ + /*T lerp(double w, T low_v, double low_w, T high_v, double high_w){ if(low_w == high_w) return low_v; //if the interval is of zero length, just return one of the bounds double alpha = (w - low_w) / (high_w - low_w); //calculate the interpolation factor return (1.0 - alpha) * low_v + alpha * high_v; //interpolate @@ -104,7 +105,7 @@ protected: v[n] = interp_spectrum(s, wavelengths[n]); //interpolate the measurement } return v; - } + }*/ public: @@ -113,6 +114,8 @@ public: using binary::R; using binary::read_line_0; + bip(){ hsi::init_bip(); } + /// Open a data file for reading using the class interface. /// @param filename is the name of the binary file on disk @@ -121,14 +124,19 @@ public: /// @param B is the number of samples (bands) along dimension 3 /// @param header_offset is the number of bytes (if any) in the binary header /// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band - bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int B, unsigned int header_offset, std::vector wavelengths){ + bool open(std::string filename, + unsigned long long X, + unsigned long long Y, + unsigned long long B, + unsigned long long header_offset, + std::vector wavelengths){ //copy the wavelengths to the BSQ file structure w = wavelengths; //copy the offset to the structure offset = header_offset; - return open(filename, vec(B, X, Y), header_offset); + return open(filename, vec(B, X, Y), header_offset); } @@ -136,7 +144,7 @@ public: /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param page <= B is the integer number of the band to be copied. - bool band_index( T * p, unsigned int page, bool PROGRESS = false){ + bool band_index( T * p, unsigned long long page, bool PROGRESS = false){ return binary::read_plane_0(p, page, PROGRESS); } @@ -148,9 +156,9 @@ public: //if there are no wavelengths in the BSQ file if(w.size() == 0) - return band_index(p, (unsigned int)wavelength, PROGRESS); + return band_index(p, (unsigned long long)wavelength, PROGRESS); - unsigned int XY = X() * Y(); //calculate the number of pixels in a band + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band unsigned page=0; //bands around the wavelength @@ -180,9 +188,9 @@ public: p2=(T*)malloc( XY * sizeof(T)); band_index(p1, page - 1); band_index(p2, page, PROGRESS); - for(unsigned i=0; i < XY; i++){ + for(unsigned long long i=0; i < XY; i++){ double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]); - p[i] = (p2[i] - p1[i]) * r + p1[i]; + p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]); } free(p1); free(p2); @@ -200,7 +208,7 @@ public: /// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size. /// @param x is the x-coordinate (dimension 1) of the spectrum. /// @param y is the y-coordinate (dimension 2) of the spectrum. - bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){ + bool spectrum(T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){ return read_line_0(p, x, y, PROGRESS); //read a line in the binary YZ plane (dimension order for BIP is ZXY) } @@ -211,16 +219,16 @@ public: /// @param wavelength is the wavelength of X values to retrieve bool read_x_from_xz(T* p, T* c, double wavelength) { - unsigned int B = Z(); + unsigned long long B = Z(); - unsigned page=0; //samples around the wavelength + unsigned long long page=0; //samples around the wavelength //get the bands numbers around the wavelength //if wavelength is smaller than the first one in header file if ( w[page] > wavelength ){ - for(unsigned j = 0; j < X(); j++) + for(unsigned long long j = 0; j < X(); j++) { p[j] = c[j * B]; } @@ -232,7 +240,7 @@ public: page++; //if wavelength is larger than the last wavelength in header file if (page == B) { - for(unsigned j = 0; j < X(); j++) + for(unsigned long long j = 0; j < X(); j++) { p[j] = c[(j + 1) * B - 1]; } @@ -246,17 +254,17 @@ public: p1=(T*)malloc( X() * sizeof(T)); //memory allocation p2=(T*)malloc( X() * sizeof(T)); //band_index(p1, page - 1); - for(unsigned j = 0; j < X(); j++) + for(unsigned long long j = 0; j < X(); j++) { p1[j] = c[j * B + page - 1]; } //band_index(p2, page ); - for(unsigned j = 0; j < X(); j++) + for(unsigned long long j = 0; j < X(); j++) { p2[j] = c[j * B + page]; } - for(unsigned i=0; i < X(); i++){ + for(unsigned long long i=0; i < X(); i++){ double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]); p[i] = (p2[i] - p1[i]) * r + p1[i]; } @@ -266,7 +274,7 @@ public: else //if the wavelength is equal to a wavelength in header file { //band_index(p, page); - for(unsigned j = 0; j < X(); j++) + for(unsigned long long j = 0; j < X(); j++) { p[j] = c[j * B + page]; } @@ -276,29 +284,30 @@ public: } /// Retrieve a single pixel and store it in a pre-allocated double array. - bool pixeld(double* p, unsigned n){ - unsigned bandnum = X() * Y(); //calculate numbers in one band + bool pixeld(double* p, unsigned long long n){ + unsigned long long bandnum = X() * Y(); //calculate numbers in one band if ( n >= bandnum){ //make sure the pixel number is right std::cout<<"ERROR: sample or line out of range"<= N){ //make sure the pixel number is right std::cout<<"ERROR: sample or line out of range"<::read_plane_2(p, y); } @@ -330,7 +339,6 @@ public: std::vector base_vals; //allocate space for the values at each baseline point double aw, bw; //surrounding baseline point wavelengths T av, bv; //surrounding baseline point values - double alpha; unsigned long long ai, bi; //surrounding baseline point band indices for(unsigned long long n = 0; n < N; n++){ //for each pixel in the image if(mask != NULL && !mask[n]){ //if the pixel isn't masked @@ -393,14 +401,14 @@ public: T* s = (T*)malloc(sizeof(T) * B); //allocate memory to store a pixel T nv; //stores the value of the normalized band for(unsigned long long n = 0; n < N; n++){ //for each pixel in the image - pixel(s, n); //retrieve the spectrum s - nv = interp_spectrum(s, w); //find the value of the normalization band if(mask != NULL && !mask[n]) //if the normalization band is below threshold memset(s, 0, sizeof(T) * B); //set the output to zero else{ - for(unsigned long long b = 0; b < B; b++){ //for each band in the spectrum + pixel(s, n); //retrieve the spectrum s + nv = interp_spectrum(s, w); //find the value of the normalization band + + for(unsigned long long b = 0; b < B; b++) //for each band in the spectrum s[b] /= nv; //divide by the normalization value - } } if(PROGRESS) progress = (double)(n+1) / N * 100; //set the current progress @@ -419,23 +427,23 @@ public: /// @param outname is the name of the output BIL file to be saved to disk. bool bil(std::string outname, bool PROGRESS = false) { - unsigned int S = X() * Z() * sizeof(T); //calculate the number of bytes in a ZX slice + unsigned long long S = X() * Z() * sizeof(T); //calculate the number of bytes in a ZX slice std::ofstream target(outname.c_str(), std::ios::binary); - std::string headername = outname + ".hdr"; + //std::string headername = outname + ".hdr"; T * p; //pointer to the current ZX slice for bip file p = (T*)malloc(S); T * q; //pointer to the current XZ slice for bil file q = (T*)malloc(S); - for ( unsigned i = 0; i < Y(); i++) + for ( unsigned long long i = 0; i < Y(); i++) { read_plane_y(p, i); - for ( unsigned k = 0; k < Z(); k++) + for ( unsigned long long k = 0; k < Z(); k++) { - unsigned ks = k * X(); - for ( unsigned j = 0; j < X(); j++) + unsigned long long ks = k * X(); + for ( unsigned long long j = 0; j < X(); j++) q[ks + j] = p[k + j * Z()]; if(PROGRESS) progress = (double)(i * Z() + k+1) / (Y() * Z()) * 100; @@ -459,12 +467,12 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size. bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){ - unsigned XY = X() * Y(); + unsigned long long XY = X() * Y(); band(result, wavelength); //get band //perform the baseline correction double r = (double) (wavelength - lb) / (double) (rb - lb); - for(unsigned i=0; i < XY; i++){ + for(unsigned long long i=0; i < XY; i++){ result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] ); } return true; @@ -480,8 +488,8 @@ public: T* lp; T* rp; - unsigned XY = X() * Y(); - unsigned S = XY * sizeof(T); + unsigned long long XY = X() * Y(); + unsigned long long S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation rp = (T*) malloc(S); @@ -510,8 +518,8 @@ public: T* cur; //current band 1 T* cur2; //current band 2 - unsigned XY = X() * Y(); - unsigned S = XY * sizeof(T); + unsigned long long XY = X() * Y(); + unsigned long long S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation rp = (T*) malloc(S); @@ -521,9 +529,9 @@ public: memset(result, (char)0, S); //find the wavelenght position in the whole band - unsigned int n = w.size(); - unsigned int ai = 0; //left bound position - unsigned int bi = n - 1; //right bound position + unsigned long long n = w.size(); + unsigned long long ai = 0; //left bound position + unsigned long long bi = n - 1; //right bound position @@ -552,23 +560,23 @@ public: //calculate the beginning and the ending part baseline_band(lb, rb, lp, rp, rab, cur2); //ending part baseline_band(lb, rb, lp, rp, w[bi], cur); - for(unsigned j = 0; j < XY; j++){ - result[j] += (rab - w[bi]) * (cur[j] + cur2[j]) / 2.0; + for(unsigned long long j = 0; j < XY; j++){ + result[j] += (T)((rab - w[bi]) * ((double)cur[j] + (double)cur2[j]) / 2.0); } baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part baseline_band(lb, rb, lp, rp, w[ai], cur); - for(unsigned j = 0; j < XY; j++){ - result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0; + for(unsigned long long j = 0; j < XY; j++){ + result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0); } //calculate the area ai++; - for(unsigned i = ai; i <= bi ;i++) + for(unsigned long long i = ai; i <= bi ;i++) { baseline_band(lb, rb, lp, rp, w[ai], cur2); - for(unsigned j = 0; j < XY; j++) + for(unsigned long long j = 0; j < XY; j++) { - result[j] += (w[ai] - w[ai-1]) * (cur[j] + cur2[j]) / 2.0; + result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0); } std::swap(cur,cur2); //swap the band pointers } @@ -598,7 +606,7 @@ public: height(lb1, rb1, pos1, p1); height(lb2, rb2, pos2, p2); //calculate the ratio in result - for(unsigned i = 0; i < X() * Y(); i++){ + for(unsigned long long i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -629,7 +637,7 @@ public: area(lb1, rb1, lab1, rab1, p1); height(lb2, rb2, pos, p2); //calculate the ratio in result - for(unsigned i = 0; i < X() * Y(); i++){ + for(unsigned long long i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -662,7 +670,7 @@ public: area(lb1, rb1, lab1, rab1, p1); area(lb2, rb2, lab2, rab2, p2); //calculate the ratio in result - for(unsigned i = 0; i < X() * Y(); i++){ + for(unsigned long long i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -687,8 +695,8 @@ public: T* cur; //current band 1 T* cur2; //current band 2 - unsigned XY = X() * Y(); - unsigned S = XY * sizeof(T); + unsigned long long XY = X() * Y(); + unsigned long long S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation rp = (T*) malloc(S); @@ -698,9 +706,9 @@ public: memset(result, (char)0, S); //find the wavelenght position in the whole band - unsigned int n = w.size(); - unsigned int ai = 0; //left bound position - unsigned int bi = n - 1; //right bound position + unsigned long long n = w.size(); + unsigned long long ai = 0; //left bound position + unsigned long long bi = n - 1; //right bound position //to make sure the left and the right bound are in the bandwidth if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){ @@ -727,23 +735,23 @@ public: //calculate the beginning and the ending part baseline_band(lb, rb, lp, rp, rab, cur2); //ending part baseline_band(lb, rb, lp, rp, w[bi], cur); - for(unsigned j = 0; j < XY; j++){ - result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0; + for(unsigned long long j = 0; j < XY; j++){ + result[j] += (T)((rab - w[bi]) * (rab + w[bi]) * ((double)cur[j] + (double)cur2[j]) / 4.0); } baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part baseline_band(lb, rb, lp, rp, w[ai], cur); - for(unsigned j = 0; j < XY; j++){ - result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0; + for(unsigned long long j = 0; j < XY; j++){ + result[j] += (T)((w[ai] - lab) * (w[ai] + lab) * ((double)cur[j] + (double)cur2[j]) / 4.0); } //calculate f(x) times x ai++; - for(unsigned i = ai; i <= bi ;i++) + for(unsigned long long i = ai; i <= bi ;i++) { baseline_band(lb, rb, lp, rp, w[ai], cur2); - for(unsigned j = 0; j < XY; j++) + for(unsigned long long j = 0; j < XY; j++) { - result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0; + result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0); } std::swap(cur,cur2); //swap the band pointers } @@ -770,7 +778,7 @@ public: x_area(lb, rb, lab, rab, p1); area(lb, rb, lab, rab, p2); //calculate the ratio in result - for(unsigned i = 0; i < X() * Y(); i++){ + for(unsigned long long i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -794,7 +802,7 @@ public: T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band band(temp, mask_band, PROGRESS); - for (unsigned int i = 0; i < X() * Y();i++) { + for (unsigned long long i = 0; i < X() * Y();i++) { if (temp[i] < threshold) p[i] = 0; else @@ -814,17 +822,17 @@ public: std::ofstream target(outfile.c_str(), std::ios::binary); - unsigned ZX = Z() * X(); //calculate the number of values in a page (XZ in BIP) - unsigned L = ZX * sizeof(T); //calculate the number of bytes in a page + unsigned long long ZX = Z() * X(); //calculate the number of values in a page (XZ in BIP) + unsigned long long L = ZX * sizeof(T); //calculate the number of bytes in a page T * temp = (T*)malloc(L); //allocate space for that page - for (unsigned i = 0; i < Y(); i++) //for each page (Y in BIP) + for (unsigned long long i = 0; i < Y(); i++) //for each page (Y in BIP) { read_plane_y(temp, i); //load that page (it's pointed to by temp) - for ( unsigned j = 0; j < X(); j++) //for each X value + for ( unsigned long long j = 0; j < X(); j++) //for each X value { - for (unsigned k = 0; k < Z(); k++) //for each B value (band) + for (unsigned long long k = 0; k < Z(); k++) //for each B value (band) { if (p[i * X() + j] == 0) //if the mask value is zero temp[j * Z() + k] = 0; //set the pixel value to zero @@ -850,15 +858,15 @@ public: std::ofstream target(outfile.c_str(), std::ios::binary); //allocate space for a single spectrum - unsigned long B = Z(); + unsigned long long B = Z(); T* spectrum = (T*) malloc(B * sizeof(T)); //calculate the number of pixels in a band - unsigned long XY = X() * Y(); + unsigned long long XY = X() * Y(); //for each pixel - unsigned long skip = 0; //number of spectra to skip - for(unsigned long x = 0; x < XY; x++){ + unsigned long long skip = 0; //number of spectra to skip + for(unsigned long long x = 0; x < XY; x++){ //if the current pixel isn't masked if( mask[x] == 0){ @@ -879,9 +887,8 @@ public: //write this pixel out target.write((char *)spectrum, B * sizeof(T)); - - if(PROGRESS) progress = (double) (x+1) / XY * 100; } + if(PROGRESS) progress = (double) (x+1) / XY * 100; } @@ -892,7 +899,7 @@ public: return true; } - bool unsift(std::string outfile, unsigned char* mask, unsigned int samples, unsigned int lines){ + bool unsift(std::string outfile, unsigned char* mask, unsigned long long samples, unsigned long long lines){ // open an output stream std::ofstream target(outfile.c_str(), std::ios::binary); @@ -901,7 +908,7 @@ public: file.seekg(0, std::ios::beg); //allocate space for a single spectrum - unsigned long B = Z(); + unsigned long long B = Z(); T* spectrum = (T*) malloc(B * sizeof(T)); //allocate space for a spectrum of zeros @@ -909,11 +916,11 @@ public: memset(zeros, 0, B * sizeof(T)); //calculate the number of pixels in a band - unsigned long XY = samples * lines; + unsigned long long XY = samples * lines; //for each pixel - unsigned long skip = 0; //number of spectra to skip - for(unsigned long x = 0; x < XY; x++){ + unsigned long long skip = 0; //number of spectra to skip + for(unsigned long long x = 0; x < XY; x++){ //if the current pixel isn't masked if( mask[x] == 0){ @@ -951,24 +958,23 @@ public: /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){ - unsigned long long XY = X() * Y(); - T* temp = (T*)malloc(sizeof(T) * Z()); - //Iinitialize - for (unsigned j = 0; j < Z(); j++){ - p[j] = 0; - } - //calculate vaild number in a band - unsigned count = nnz(mask); - - //calculate average number of a band - for (unsigned i = 0; i < XY; i++){ - if (mask == NULL || mask[i] != 0){ - pixel(temp, i); - for (unsigned j = 0; j < Z(); j++){ - p[j] += (double)temp[j] / (double)count; + unsigned long long XY = X() * Y(); //calculate the total number of pixels in the HSI + T* temp = (T*)malloc(sizeof(T) * Z()); //allocate space for the current spectrum to be read + memset(p, 0, sizeof(double) * Z()); //initialize the average spectrum to zero (0) + //for (unsigned j = 0; j < Z(); j++){ + // p[j] = 0; + //} + + unsigned long long count = nnz(mask); //calculate the number of masked pixels + + for (unsigned long long i = 0; i < XY; i++){ //for each pixel in the HSI + if (mask == NULL || mask[i] != 0){ //if the pixel is masked + pixel(temp, i); //get the spectrum + for (unsigned long long j = 0; j < Z(); j++){ //for each spectral component + p[j] += (double)temp[j] / (double)count; //add the weighted value to the average } } - if(PROGRESS) progress = (double)(i+1) / XY * 100; + if(PROGRESS) progress = (double)(i+1) / XY * 100; //increment the progress } free(temp); @@ -976,6 +982,7 @@ public: } #ifdef CUDA_FOUND /// Calculate the covariance matrix for masked pixels using cuBLAS + /// Note that cuBLAS only supports integer-sized arrays, so there may be issues with large spectra bool co_matrix_cublas(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){ cudaError_t cudaStat; @@ -994,7 +1001,7 @@ public: cudaStat = cudaMalloc(&A_dev, B * B * sizeof(double)); //allocate space on the CUDA device for the covariance matrix cudaStat = cudaMemset(A_dev, 0, B * B * sizeof(double)); //initialize the covariance matrix to zero (0) cudaStat = cudaMalloc(&avg_dev, B * sizeof(double)); //allocate space on the CUDA device for the average spectrum - stat = cublasSetVector(B, sizeof(double), avg, 1, avg_dev, 1); //copy the average spectrum to the CUDA device + stat = cublasSetVector((int)B, sizeof(double), avg, 1, avg_dev, 1); //copy the average spectrum to the CUDA device double ger_alpha = 1.0/(double)XY; //scale the outer product by the inverse of the number of samples (mean outer product) double axpy_alpha = -1; //multiplication factor for the average spectrum (in order to perform a subtraction) @@ -1007,22 +1014,22 @@ public: for (unsigned long long xy = 0; xy < XY; xy++){ //for each pixel if (mask == NULL || mask[xy] != 0){ pixeld(s, xy); //retreive the spectrum at the current xy pixel location - stat = cublasSetVector(B, sizeof(double), s, 1, s_dev, 1); //copy the spectrum from the host to the device - stat = cublasDaxpy(handle, B, &axpy_alpha, avg_dev, 1, s_dev, 1); //subtract the average spectrum - stat = cublasDsyr(handle, CUBLAS_FILL_MODE_UPPER, B, &ger_alpha, s_dev, 1, A_dev, B); //calculate the covariance matrix (symmetric outer product) + stat = cublasSetVector((int)B, sizeof(double), s, 1, s_dev, 1); //copy the spectrum from the host to the device + stat = cublasDaxpy(handle, (int)B, &axpy_alpha, avg_dev, 1, s_dev, 1); //subtract the average spectrum + stat = cublasDsyr(handle, CUBLAS_FILL_MODE_UPPER, (int)B, &ger_alpha, s_dev, 1, A_dev, (int)B); //calculate the covariance matrix (symmetric outer product) } if(PROGRESS) progress = (double)(xy+1) / XY * 100; //record the current progress } - cublasGetMatrix(B, B, sizeof(double), A_dev, B, co, B); //copy the result from the GPU to the CPU + cublasGetMatrix((int)B, (int)B, sizeof(double), A_dev, (int)B, co, (int)B); //copy the result from the GPU to the CPU cudaFree(A_dev); //clean up allocated device memory cudaFree(s_dev); cudaFree(avg_dev); - for(unsigned i = 0; i < B; i++){ //copy the upper triangular portion to the lower triangular portion - for(unsigned j = i+1; j < B; j++){ + for(unsigned long long i = 0; i < B; i++){ //copy the upper triangular portion to the lower triangular portion + for(unsigned long long j = i+1; j < B; j++){ co[B * i + j] = co[B * j + i]; } } @@ -1046,14 +1053,13 @@ public: unsigned long long XY = X() * Y(); unsigned long long B = Z(); T* temp = (T*)malloc(sizeof(T) * B); - //count vaild pixels in a band - unsigned long long count = nnz(mask); + + unsigned long long count = nnz(mask); //count the number of masked pixels //initialize covariance matrix memset(co, 0, B * B * sizeof(double)); //calculate covariance matrix - T i_diff; //stores (temp[i] - avg[i]) to speed math double* co_half = (double*) malloc(B * B * sizeof(double)); //allocate space for a higher-precision intermediate matrix double* temp_precise = (double*) malloc(B * sizeof(double)); memset(co_half, 0, B * B * sizeof(double)); //initialize the high-precision matrix with zeros @@ -1083,10 +1089,10 @@ public: return true; } - bool project(std::string outfile, double* center, double* basis, unsigned long long M, bool PROGRESS = false){ + bool project(std::string outfile, double* center, double* basis, unsigned long long M, unsigned char* mask = NULL, bool PROGRESS = false){ std::ofstream target(outfile.c_str(), std::ios::binary); //open the target binary file - std::string headername = outfile + ".hdr"; //the header file name + //std::string headername = outfile + ".hdr"; //the header file name //memory allocation unsigned long long XY = X() * Y(); @@ -1096,15 +1102,17 @@ public: T* rs = (T*)malloc(sizeof(T) * M); //allocate space for the projected spectrum double* bv; //pointer to the current projection vector for(unsigned long long xy = 0; xy < XY; xy++){ //for each spectrum in the image - pixel(s, xy); //load the spectrum - memset(rs, 0, sizeof(T) * M); //initialize the rotated spectrum to zero (0) - for(unsigned long long m = 0; m < M; m++){ //for each basis vector - bv = &basis[m * B]; //assign 'bv' to the beginning of the basis vector - for(unsigned long long b = 0; b < B; b++){ //for each band - rs[m] += (s[b] - center[b]) * bv[b]; //center the spectrum and perform the projection + memset(rs, 0, sizeof(T) * M); + if(mask == NULL || mask[xy]){ + pixel(s, xy); //load the spectrum + for(unsigned long long m = 0; m < M; m++){ //for each basis vector + bv = &basis[m * B]; //assign 'bv' to the beginning of the basis vector + for(unsigned long long b = 0; b < B; b++){ //for each band + rs[m] += (T)(((double)s[b] - center[b]) * bv[b]); //center the spectrum and perform the projection + } } } - + target.write(reinterpret_cast(rs), sizeof(T) * M); //write the projected vector if(PROGRESS) progress = (double)(xy+1) / XY * 100; } @@ -1135,7 +1143,7 @@ public: for(unsigned long long c = 0; c < C; c++){ //for each basis vector coefficient bv = &basis[c * B]; //assign 'bv' to the beginning of the basis vector for(unsigned long long b = 0; b < B; b++){ //for each component of the basis vector - s[b] += coeff[c] * bv[b] + center[b]; //calculate the contribution of each element of the basis vector in the final spectrum + s[b] += (T)((double)coeff[c] * bv[b] + center[b]); //calculate the contribution of each element of the basis vector in the final spectrum } } diff --git a/stim/envi/bsq.h b/stim/envi/bsq.h index e533ce7..81f1dd3 100644 --- a/stim/envi/bsq.h +++ b/stim/envi/bsq.h @@ -2,7 +2,7 @@ #define STIM_BSQ_H #include "../envi/envi_header.h" -#include "../envi/binary.h" +#include "../envi/hsi.h" #include "../envi/bil.h" #include #include @@ -22,27 +22,27 @@ namespace stim{ */ template -class bsq: public binary { +class bsq: public hsi { protected: - std::vector w; //band wavelengths - unsigned int offset; + //std::vector w; //band wavelengths + unsigned long long offset; using binary::R; - unsigned long X(){ + /*unsigned long long X(){ return R[0]; } - unsigned long Y(){ + unsigned long long Y(){ return R[1]; } - unsigned long Z(){ + unsigned long long Z(){ return R[2]; - } - - using binary::nnz; + }*/ + using hsi::w; //use the wavelength array in stim::hsi + using hsi::nnz; using binary::progress; public: @@ -51,8 +51,8 @@ public: using binary::file; using binary::read_line_2; using binary::read_plane_2; - //using binary::getSlice; + bsq(){ hsi::init_bsq(); } /// Open a data file for reading using the class interface. @@ -62,21 +62,26 @@ public: /// @param B is the number of samples (bands) along dimension 3 /// @param header_offset is the number of bytes (if any) in the binary header /// @param wavelengths is an STL vector of size B specifying a numerical label for each band - bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int B, unsigned int header_offset, std::vector wavelengths){ + bool open(std::string filename, + unsigned long long X, + unsigned long long Y, + unsigned long long B, + unsigned long long header_offset, + std::vector wavelengths){ //copy the wavelengths to the BSQ file structure w = wavelengths; //copy the wavelengths to the structure offset = header_offset; - return open(filename, vec(X, Y, B), header_offset); + return open(filename, vec(X, Y, B), header_offset); } /// Retrieve a single band (based on index) and stores it in pre-allocated memory. /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param page <= B is the integer number of the band to be copied. - bool band_index( T * p, unsigned int page){ + bool band_index( T * p, unsigned long long page){ return read_plane_2(p, page); //call the binary read_plane function (don't let it update the progress) } @@ -88,10 +93,10 @@ public: if(PROGRESS) progress = 0; //if there are no wavelengths in the BSQ file if(w.size() == 0) - return band_index(p, (unsigned int)wavelength); + return band_index(p, (unsigned long long)wavelength); - unsigned int XY = X() * Y(); //calculate the number of pixels in a band - unsigned page = 0; + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band + unsigned long long page = 0; //get the two neighboring bands (above and below 'wavelength') @@ -122,9 +127,9 @@ public: p2=(T*)malloc( XY * sizeof(T)); band_index(p1, page - 1); band_index(p2, page ); - for(unsigned i=0; i < XY; i++){ - double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]); - p[i] = (p2[i] - p1[i]) * r + p1[i]; + for(unsigned long long i=0; i < XY; i++){ + double r = (wavelength - w[page-1]) / (w[page] - w[page-1]); + p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]); } free(p1); free(p2); @@ -142,7 +147,7 @@ public: /// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size. /// @param x is the x-coordinate (dimension 1) of the spectrum. /// @param y is the y-coordinate (dimension 2) of the spectrum. - bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){ + bool spectrum(T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){ return read_line_2(p, x, y, PROGRESS); } @@ -150,16 +155,16 @@ public: /// @param p is a pointer to pre-allocated memory at least sizeof(T) in size. /// @param n is an integer index to the pixel using linear array indexing. - bool pixel(T * p, unsigned n){ + bool pixel(T * p, unsigned long long n){ - unsigned bandnum = X() * Y(); //calculate numbers in one band + unsigned long long bandnum = X() * Y(); //calculate numbers in one band if ( n >= bandnum){ //make sure the pixel number is right std::cout<<"ERROR: sample or line out of range"<::header, std::ios::beg); //point to the certain pixel - for (unsigned i = 0; i < Z(); i++) + for (unsigned long long i = 0; i < Z(); i++) { file.read((char *)(p + i), sizeof(T)); file.seekg((bandnum - 1) * sizeof(T), std::ios::cur); //go to the next band @@ -174,20 +179,20 @@ public: /// @param wls is the list of baseline points based on band labels. bool baseline(std::string outname, std::vector wls, unsigned char* mask = NULL, bool PROGRESS = false ) { - unsigned N = wls.size(); //get the number of baseline points + size_t N = wls.size(); //get the number of baseline points std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file std::string headername = outname + ".hdr"; //the header file name //simplify image resolution - unsigned int B = Z(); //calculate the number of bands - unsigned int XY = X() * Y(); //calculate the number of pixels in a band - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band + unsigned long long B = Z(); //calculate the number of bands + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band + unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band double ai, bi; //stores the two baseline points wavelength surrounding the current band double ci; //stores the current band's wavelength - unsigned control=0; + unsigned long long control=0; T * a; //pointers to the high and low band images T * b; @@ -222,7 +227,7 @@ public: band(b, bi); //correct every band - for(unsigned cii = 0; cii < B; cii++){ + for(unsigned long long cii = 0; cii < B; cii++){ //update baseline points, if necessary if( w[cii] >= bi && cii != B - 1) { @@ -255,7 +260,7 @@ public: ci = w[cii]; //perform the baseline correction - for(unsigned i=0; i < XY; i++){ + for(unsigned long long i=0; i < XY; i++){ if(mask != NULL && !mask[i]) //if the pixel is excluded by a mask c[i] = 0; //set the value to zero else{ @@ -287,9 +292,9 @@ public: /// @param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers. bool normalize(std::string outname, double w, unsigned char* mask = NULL, bool PROGRESS = false) { - unsigned int B = Z(); //calculate the number of bands - unsigned int XY = X() * Y(); //calculate the number of pixels in a band - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band + unsigned long long B = Z(); //calculate the number of bands + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band + unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file std::string headername = outname + ".hdr"; //the header file name @@ -302,10 +307,10 @@ public: band(b, w); //get the certain band into memory - for(unsigned j = 0; j < B; j++) + for(unsigned long long j = 0; j < B; j++) { band_index(c, j); //get the current band into memory - for(unsigned i = 0; i < XY; i++) + for(unsigned long long i = 0; i < XY; i++) { if(mask != NULL && !mask[i]) c[i] = (T)0.0; @@ -333,26 +338,24 @@ public: bool bil(std::string outname, bool PROGRESS = false) { //simplify image resolution - //unsigned int L = X() * Z() * sizeof(T); //calculate the number of bytes of a ZX slice - unsigned int jump = (Y() - 1) * X() * sizeof(T); + unsigned long long jump = (Y() - 1) * X() * sizeof(T); std::ofstream target(outname.c_str(), std::ios::binary); std::string headername = outname + ".hdr"; - unsigned int L = X(); + unsigned long long L = X(); T* line = (T*)malloc(sizeof(T) * L); - for ( unsigned y = 0; y < Y(); y++) //for each y position + for ( unsigned long long y = 0; y < Y(); y++) //for each y position { file.seekg(y * X() * sizeof(T), std::ios::beg); //seek to the beginning of the xz slice - for ( unsigned z = 0; z < Z(); z++ ) //for each band + for ( unsigned long long z = 0; z < Z(); z++ ) //for each band { file.read((char *)line, sizeof(T) * X()); //read a line target.write((char*)line, sizeof(T) * X()); //write the line to the output file file.seekg(jump, std::ios::cur); //seek to the next band if(PROGRESS) progress = (double)((y+1) * Z() + z + 1) / (Z() * Y()) * 100; //update the progress counter } - //std::cout< w[n-1] || rb >w[n-1]){ @@ -639,23 +642,23 @@ public: //calculate the beginning and the ending part baseline_band(lb, rb, lp, rp, rab, cur2); //ending part baseline_band(lb, rb, lp, rp, w[bi], cur); - for(unsigned j = 0; j < XY; j++){ - result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0; + for(unsigned long long j = 0; j < XY; j++){ + result[j] += (T)((rab - w[bi]) * (rab + w[bi]) * ((double)cur[j] + (double)cur2[j]) / 4.0); } baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part baseline_band(lb, rb, lp, rp, w[ai], cur); - for(unsigned j = 0; j < XY; j++){ - result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0; + for(unsigned long long j = 0; j < XY; j++){ + result[j] += (T)((w[ai] - lab) * (w[ai] + lab) * ((double)cur[j] + (double)cur2[j]) / 4.0); } //calculate f(x) times x ai++; - for(unsigned i = ai; i <= bi ;i++) + for(unsigned long long i = ai; i <= bi ;i++) { baseline_band(lb, rb, lp, rp, w[ai], cur2); - for(unsigned j = 0; j < XY; j++) + for(unsigned long long j = 0; j < XY; j++) { - result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0; + result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0); } std::swap(cur,cur2); //swap the band pointers } @@ -682,7 +685,7 @@ public: x_area(lb, rb, lab, rab, p1); area(lb, rb, lab, rab, p2); //calculate the ratio in result - for(unsigned i = 0; i < X() * Y(); i++){ + for(unsigned long long i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -706,13 +709,13 @@ public: T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band band(temp, mask_band); - for (unsigned int i = 0; i < X() * Y(); i++) { + for (unsigned long long i = 0; i < X() * Y(); i++) { if (temp[i] < threshold) p[i] = 0; else p[i] = 255; - if(PROGRESS) progress = (double) (i+1) / (X() * Y()); + if(PROGRESS) progress = (double) (i+1) / (X() * Y()) * 100; } free(temp); @@ -728,15 +731,15 @@ public: std::ofstream target(outfile.c_str(), std::ios::binary); - unsigned XY = X() * Y(); //calculate number of a band - unsigned L = XY * sizeof(T); + unsigned long long XY = X() * Y(); //calculate number of a band + unsigned long long L = XY * sizeof(T); T * temp = (T*)malloc(L); - for (unsigned i = 0; i < Z(); i++) //for each spectral bin + for (unsigned long long i = 0; i < Z(); i++) //for each spectral bin { band_index(temp, i); //get the specified band (by index) - for ( unsigned j = 0; j < XY; j++) // for each pixel + for ( unsigned long long j = 0; j < XY; j++) // for each pixel { if(p[j] == 0){ //if the mask is 0 at that pixel temp[j] = 0; //set temp to zero @@ -787,56 +790,24 @@ public: return true; } - - /// Saves only those spectra corresponding to mask value != 0 to pre-allocated memory - /// @param matrix is the destination for the sifted pixels - /// @param p is the mask file used for sifting - /*bool sift(T* matrix, unsigned char* p){ - - // open a band (XY plane) - unsigned long XY = X() * Y(); //Number of XY pixels - unsigned long L = XY * sizeof(T); //size of XY pixels - - //count the number of masked pixels - unsigned long P = 0; //allocate space for the number of pixels - for(unsigned long xy = 0; xy < XY; xy++){ //count the number of pixels - if(p[xy] != 0) P++; - } - - T* bandvals = (T*) malloc(sizeof(T) * P); //allocate a temporary array to store the pixels for a single band - memset(matrix, 0, sizeof(T) * P * Z()); - for (unsigned long z = 0; z < Z(); z++){ //for each band - - if(!sift_band(bandvals, z, p)){ - std::cout<<"ERROR sifting band index "<()); + file = new bsq(); else if(header.data_type == envi_header::float64) - return(file = new bsq()); + file = new bsq(); } else if(header.interleave == envi_header::BIP){ if(header.data_type ==envi_header::float32) - return(file = new bip()); + file = new bip(); else if(header.data_type == envi_header::float64) - return(file = new bip()); + file = new bip(); } else if(header.interleave == envi_header::BIL){ if(header.data_type ==envi_header::float32) - return(file = new bil()); + file = new bil(); else if(header.data_type == envi_header::float64) - return(file = new bil()); + file = new bil(); } - exit(1); //if the function hasn't already returned, we don't handle this state - } /// Open a previously opened ENVI file @@ -179,11 +177,13 @@ public: /// @param filename is the name of the ENVI binary file /// @param header is an ENVI header structure bool open(std::string filename, stim::envi_header h){ - allocate(); + header = h; //store the header fname = filename; //save the filename + allocate(); + return open(); //open the ENVI file; @@ -294,7 +294,7 @@ public: } } - void project(std::string outfile, double* center, double* basis, unsigned long long M, bool PROGRESS = false){ + void project(std::string outfile, double* center, double* basis, unsigned long long M, unsigned char* mask, bool PROGRESS = false){ if(header.interleave == envi_header::BSQ){ //if the infile is bsq file std::cout<<"ERROR: BSQ projection not supported"<*)file)->project(outfile, center, basis, M, PROGRESS); + ((bip*)file)->project(outfile, center, basis, M, mask, PROGRESS); else if(header.data_type == envi_header::float64) - ((bip*)file)->project(outfile, center, basis, M, PROGRESS); + ((bip*)file)->project(outfile, center, basis, M, mask, PROGRESS); else{ std::cout<<"ERROR: unidentified data type"< 0 ) nnz++; //create a new header @@ -625,12 +627,10 @@ public: std::cout << "ERROR: unidentified file type" << std::endl; exit(1); } - - return false; } - bool unsift(std::string outfile, unsigned char* mask, unsigned int samples, unsigned int lines){ + bool unsift(std::string outfile, unsigned char* mask, unsigned long long samples, unsigned long long lines){ //create a new header envi_header new_header = header; @@ -976,7 +976,7 @@ public: /// @param ptr is a pointer to pre-allocated memory of size B*sizeof(T) /// @param x is the x-coordinate of the spectrum /// @param y is the y-coordinate of the spectrum - bool spectrum(void* ptr, unsigned int x, unsigned int y, bool PROGRESS = false){ + bool spectrum(void* ptr, unsigned long long x, unsigned long long y, bool PROGRESS = false){ if(header.interleave == envi_header::BSQ){ //if the infile is bsq file if(header.data_type ==envi_header::float32) @@ -1015,7 +1015,7 @@ public: /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param page <= B is the integer number of the band to be copied. - bool band_index(void* ptr, unsigned int b){ + bool band_index(void* ptr, unsigned long long b){ if (header.interleave == envi_header::BSQ){ //if the infile is bsq file if (header.data_type == envi_header::float32) return ((bsq*)file)->band_index((float*)ptr, b); diff --git a/stim/envi/envi_header.h b/stim/envi/envi_header.h index 3ca70e1..163ee3a 100644 --- a/stim/envi/envi_header.h +++ b/stim/envi/envi_header.h @@ -23,10 +23,10 @@ struct envi_header std::string description; - unsigned int samples; //x-axis - unsigned int lines; //y-axis - unsigned int bands; //spectral axis - unsigned int header_offset; //header offset for binary file (in bytes) + unsigned long long samples; //x-axis + unsigned long long lines; //y-axis + unsigned long long bands; //spectral axis + unsigned long long header_offset; //header offset for binary file (in bytes) std::string file_type; //should be "ENVI Standard" envi_header::dataType data_type; //data representation; common value is 4 (32-bit float) @@ -81,7 +81,7 @@ struct envi_header if(line.length() == 0) return line; //trims whitespace from the beginning and end of line - unsigned int start_i, end_i; + size_t start_i, end_i; for(start_i=0; start_i < line.length(); start_i++) if(line[start_i] != 32) { @@ -307,7 +307,7 @@ struct envi_header } //make sure the number of bands matches the number of wavelengths - unsigned int wavelengths = wavelength.size(); + size_t wavelengths = wavelength.size(); if(wavelengths && bands != wavelengths) { std::cout<<"ENVI Header Error -- Number of wavelengths doesn't match the number of bands. Bands = "< +#include + +namespace stim{ + +/** + The BIL class represents a 3-dimensional binary file stored using band interleaved by line (BIL) image encoding. The binary file is stored + such that X-Z "frames" are stored sequentially to form an image stack along the y-axis. When accessing the data sequentially on disk, + the dimensions read, from fastest to slowest, are X, Z, Y. + + This class is optimized for data streaming, and therefore supports extremely large (terabyte-scale) files. Data is loaded from disk + on request. Functions used to access data are written to support efficient reading. +*/ + +template +class hsi: public binary { + +protected: + unsigned char O[3]; //order of dimensions (orientation on disk) + //[X Y B]: [0 1 2] = bsq, [0 2 1] = bil, [1 2 0] = bip + + std::vector w; //band wavelengths + + unsigned long long X(){ return R[O[0]]; } + unsigned long long Y(){ return R[O[1]]; } + unsigned long long Z(){ return R[O[2]]; } + + /// Initialize axis orders based on common interleave formats + void init_bsq(){ O[0] = 0; O[1] = 1; O[2] = 2; } + void init_bil(){ O[0] = 0; O[1] = 2; O[2] = 1; } + void init_bip(){ O[0] = 1; O[1] = 2; O[2] = 0; } + + /// Calculate the number of masked pixels in a given mask + unsigned long long nnz(unsigned char* mask){ + unsigned long long XY = X() * Y(); //calculate the total number of pixels in the HSI + if(mask == NULL) return XY; //if the mask is null, assume that all of the pixels are masked + + unsigned long long n = 0; //initialize the number of masked pixels to zero (0) + for(unsigned long long xy = 0; xy < XY; xy++) //for each pixel in the HSI + if(mask[xy]) n++; //if the mask value is nonzero, increment the number of masked pixels + return n; //return the number of masked pixels + } + + T lerp(double w, T low_v, double low_w, T high_v, double high_w){ + if(low_w == high_w) return low_v; //if the interval is of zero length, just return one of the bounds + double alpha = (w - low_w) / (high_w - low_w); //calculate the interpolation factor + return (T)((1.0 - alpha) * low_v + alpha * high_v); //interpolate + } + + /// Gets the two band indices surrounding a given wavelength + void band_bounds(double wavelength, unsigned long long& low, unsigned long long& high){ + unsigned long long B = Z(); + for(high = 0; high < B; high++){ + if(w[high] > wavelength) break; + } + low = 0; + if(high > 0) + low = high-1; + } + + /// Get the list of band numbers that bound a list of wavelengths + void band_bounds(std::vector wavelengths, + std::vector& low_bands, + std::vector& high_bands){ + + unsigned long long W = w.size(); //get the number of wavelengths in the list + low_bands.resize(W); //pre-allocate space for the band lists + high_bands.resize(W); + + for(unsigned long long wl = 0; wl < W; wl++){ //for each wavelength + band_bounds(wavelengths[wl], low_bands[wl], high_bands[wl]); //find the low and high bands + } + } + + /// Returns the interpolated in the given spectrum based on the given wavelength + + /// @param s is the spectrum in main memory of length Z() + /// @param wavelength is the wavelength value to interpolate out + T interp_spectrum(T* s, double wavelength){ + unsigned long long low, high; //indices for the bands surrounding wavelength + band_bounds(wavelength, low, high); //get the surrounding band indices + + if(high == w.size()) return s[w.size()-1]; //if the high band is above the wavelength range, return the highest wavelength value + + return lerp(wavelength, s[low], w[low], s[high], w[high]); + } + + /// Returns the interpolated value corresponding to the given list of wavelengths + std::vector interp_spectrum(T* s, std::vector wavelengths){ + + unsigned long long N = wavelengths.size(); //get the number of wavelength measurements + + std::vector v; //allocate space for the resulting values + v.resize(wavelengths.size()); + for(unsigned long long n = 0; n < N; n++){ //for each measurement + v[n] = interp_spectrum(s, wavelengths[n]); //interpolate the measurement + } + return v; + } +}; + +} //end namespace STIM + +#endif \ No newline at end of file diff --git a/stim/image/CImg.h b/stim/image/CImg.h index 265a9b7..3889259 100644 --- a/stim/image/CImg.h +++ b/stim/image/CImg.h @@ -54,7 +54,7 @@ // Set version number of the library. #ifndef cimg_version -#define cimg_version 169 +#define cimg_version 170 /*----------------------------------------------------------- # @@ -2340,7 +2340,7 @@ namespace cimg_library_suffixed { char *_message; CImgException() { _message = new char[1]; *_message = 0; } CImgException(const char *const format, ...):_message(0) { _cimg_exception_err("CImgException",true); } - CImgException(const CImgException& e) { + CImgException(const CImgException& e):std::exception(e) { const int size = std::strlen(e._message); _message = new char[size + 1]; std::strncpy(_message,e._message,size); @@ -2364,7 +2364,7 @@ namespace cimg_library_suffixed { char *_message; CImgAbortException() { _message = new char[1]; *_message = 0; } CImgAbortException(const char *const format, ...):_message(0) { _cimg_exception_err("CImgAbortException",true); } - CImgAbortException(const CImgAbortException& e) { + CImgAbortException(const CImgAbortException& e):std::exception(e) { const int size = std::strlen(e._message); _message = new char[size + 1]; std::strncpy(_message,e._message,size); @@ -2442,7 +2442,7 @@ namespace cimg_library_suffixed { static bool is_inf(const T) { return false; } static bool is_nan(const T) { return false; } static T min() { return ~max(); } - static T max() { return (T)(1UL<<(8*sizeof(T) - 1)); } + static T max() { return (T)1<<(8*sizeof(T) - 1); } static T inf() { return max(); } static T cut(const double val) { return val<(double)min()?min():val>(double)max()?max():(T)val; } static const char* format() { return "%s"; } @@ -13726,8 +13726,9 @@ namespace cimg_library_suffixed { char *user_function; unsigned int mempos, mem_img_median, debug_indent, init_size, result_dim; + bool is_parallelizable, need_input_copy; double *result; - const char *const calling_function; + const char *const calling_function, *s_op, *ss_op; typedef double (*mp_func)(_cimg_math_parser&); #define _cimg_mp_is_constant(arg) (memtype[arg]==1) // Is constant? @@ -13737,12 +13738,14 @@ namespace cimg_library_suffixed { #define _cimg_mp_is_vector(arg) (memtype[arg]>1) // Is vector? #define _cimg_mp_vector_size(arg) (_cimg_mp_is_scalar(arg)?0U:(unsigned int)memtype[arg] - 1) // Vector size #define _cimg_mp_calling_function calling_function_s()._data -#define _cimg_mp_check_type(arg,n_arg,s_op,mode,N) check_type(arg,n_arg,s_op,mode,N,ss,se,saved_char) -#define _cimg_mp_check_constant(arg,n_arg,s_op,is_strict) check_constant(arg,n_arg,s_op,is_strict,ss,se,saved_char) -#define _cimg_mp_check_matrix_square(arg,n_arg,s_op) check_matrix_square(arg,n_arg,s_op,ss,se,saved_char) -#define _cimg_mp_check_vector0(dim,s_op) check_vector0(dim,s_op,ss,se,saved_char) +#define _cimg_mp_op(s) s_op = s; ss_op = ss +#define _cimg_mp_check_type(arg,n_arg,mode,N) check_type(arg,n_arg,mode,N,ss,se,saved_char) +#define _cimg_mp_check_constant(arg,n_arg,is_strict) check_constant(arg,n_arg,is_strict,ss,se,saved_char) +#define _cimg_mp_check_matrix_square(arg,n_arg) check_matrix_square(arg,n_arg,ss,se,saved_char) +#define _cimg_mp_check_vector0(dim) check_vector0(dim,ss,se,saved_char) +#define _cimg_mp_check_list(is_out) check_list(is_out,ss,se,saved_char) #define _cimg_mp_defunc(mp) (*(mp_func)(*(mp).opcode))(mp) -#define _cimg_mp_return(x) { *se = saved_char; return x; } +#define _cimg_mp_return(x) { *se = saved_char; s_op = previous_s_op; ss_op = previous_ss_op; return x; } #define _cimg_mp_constant(val) _cimg_mp_return(constant(val)) #define _cimg_mp_scalar0(op) _cimg_mp_return(scalar0(op)) #define _cimg_mp_scalar1(op,i1) _cimg_mp_return(scalar1(op,i1)) @@ -13763,19 +13766,22 @@ namespace cimg_library_suffixed { code(_code),imgin(img_input),listin(list_input?*list_input:CImgList::const_empty()), imgout(img_output?*img_output:CImg::empty()),listout(list_output?*list_output:CImgList::empty()), img_stats(_img_stats),list_stats(_list_stats),list_median(_list_median),user_function(0), - mem_img_median(~0U),debug_indent(0),init_size(0),result_dim(0), - calling_function(funcname?funcname:"cimg_math_parser") { + mem_img_median(~0U),debug_indent(0),init_size(0),result_dim(0),is_parallelizable(true), + need_input_copy(false),calling_function(funcname?funcname:"cimg_math_parser") { if (!expression || !*expression) throw CImgArgumentException("[_cimg_math_parser] " "CImg<%s>::%s: Empty expression.", pixel_type(),_cimg_mp_calling_function); const char *_expression = expression; - while (*_expression && *_expression<=' ') ++_expression; + while (*_expression && (*_expression<=' ' || *_expression==';')) ++_expression; CImg::string(_expression).move_to(expr); + char *ps = &expr.back() - 1; + while (ps>expr._data && (*ps==' ' || *ps==';')) --ps; + *(++ps) = 0; expr._width = ps - expr._data + 1; // Ease the retrieval of previous non-space characters afterwards. pexpr.assign(expr._width); - const char *ps; + char c, *pe = pexpr._data; for (ps = expr._data, c = ' '; *ps; ++ps) { if (*ps!=' ') c = *ps; @@ -13854,6 +13860,7 @@ namespace cimg_library_suffixed { // [11] = xm, [12] = ym, [13] = zm, [14] = cm, [15] = xM, [16] = yM, [17] = zM, [18]=cM, [19]=i0...[28]=i9, // Compile expression into a serie of opcodes. + s_op = ""; ss_op = expr._data; const unsigned int ind_result = compile(expr._data,expr._data + expr._width - 1,0,0); p_code_end = code.end(); @@ -13867,6 +13874,7 @@ namespace cimg_library_suffixed { reserved_label.assign(); expr.assign(); pexpr.assign(); + opcode.assign(); opcode._width = opcode._depth = opcode._spectrum = 1; opcode._is_shared = true; @@ -13888,7 +13896,7 @@ namespace cimg_library_suffixed { imgin(CImg::const_empty()),listin(CImgList::const_empty()), imgout(CImg::empty()),listout(CImgList::empty()), img_stats(_img_stats),list_stats(_list_stats),list_median(_list_median),debug_indent(0), - result_dim(0),calling_function(0) { + result_dim(0),is_parallelizable(true),need_input_copy(false),calling_function(0) { mem.assign(1 + _cimg_mp_c,1,1,1,0); // Allow to skip 'is_empty?' test in operator()() result = mem._data; } @@ -13896,8 +13904,9 @@ namespace cimg_library_suffixed { _cimg_math_parser(const _cimg_math_parser& mp): mem(mp.mem),code(mp.code),p_code_begin(mp.p_code_begin),p_code_end(mp.p_code_end), imgin(mp.imgin),listin(mp.listin),imgout(mp.imgout),listout(mp.listout),img_stats(mp.img_stats), - list_stats(mp.list_stats),list_median(mp.list_median),debug_indent(0), - result_dim(mp.result_dim), result(mem._data + (mp.result - mp.mem._data)),calling_function(0) { + list_stats(mp.list_stats),list_median(mp.list_median),debug_indent(0),result_dim(mp.result_dim), + is_parallelizable(mp.is_parallelizable), need_input_copy(mp.need_input_copy), + result(mem._data + (mp.result - mp.mem._data)),calling_function(0) { #ifdef cimg_use_openmp mem[17] = omp_get_thread_num(); #endif @@ -13905,30 +13914,8 @@ namespace cimg_library_suffixed { opcode._is_shared = true; } - // Return 'true' is the specified mathematical expression requires the input image to be copied. - // Set 'is_parallelizable' to 'false' if expression should be evaluated with a single thread. - static bool needs_input_copy(const char *expression, bool &is_parallelizable) { - if (!expression || *expression=='>' || *expression=='<') return is_parallelizable = false; - for (const char *s = expression; *s; ++s) - if ((*s=='i' || *s=='j' || *s=='I' || *s=='J') && (s[1]=='(' || s[1]=='[')) { - if (s[2]=='#') is_parallelizable = false; - else { - const char opening = s[1], ending = opening=='('?')':']'; - const char *ns; - int level = 0; - for (ns = s + 2; *ns; ++ns) { // Find ending ')' or ']'. - if (*ns==ending && !level) break; - if (*ns==opening) ++level; else if (*ns==ending) --level; - } - if (*ns && (ns[1]!='=' || ns[2]=='=')) return true; - } - } else if (((*s=='R' || *s=='G' || *s=='B' || *s=='A' || *s=='I' || *s=='J') && s[1]!='#') || - (*s=='i' && s[1]>='0' && s[1]<='7' && s[2]!='#')) return true; - return false; - } - // Compilation procedure. - unsigned int compile(char *ss, char *se, const unsigned int depth, unsigned int *p_ref) { + unsigned int compile(char *ss, char *se, const unsigned int depth, unsigned int *const p_ref) { if (depth>256) { cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " @@ -13954,10 +13941,15 @@ namespace cimg_library_suffixed { if (se<=ss || !*ss) { cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " - "CImg<%s>::%s: Missing item, in expression '%s'.", - pixel_type(),_cimg_mp_calling_function, - expr._data); + "CImg<%s>::%s: %s%s Missing %s, in expression '%s%s%s'.", + pixel_type(),_cimg_mp_calling_function,s_op,*s_op?":":"", + *s_op=='F'?"argument":"item", + (ss_op - 4)>expr._data?"...":"", + (ss_op - 4)>expr._data?ss_op - 4:expr._data, + ss_op + std::strlen(ss_op)<&expr.back()?"...":""); } + + const char *const previous_s_op = s_op, *const previous_ss_op = ss_op; const unsigned int depth1 = depth + 1; unsigned int pos, p1, p2, p3, arg1, arg2, arg3, arg4, arg5, arg6; char @@ -13966,7 +13958,6 @@ namespace cimg_library_suffixed { *const ss5 = ss + 5, *const ss6 = ss + 6, *const ss7 = ss + 7, *const ss8 = ss + 8, *s, *ps, *ns, *s0, *s1, *s2, *s3, sep = 0, end = 0; double val, val1, val2; - const char *s_op; mp_func op; // 'p_ref' is a 'unsigned int[7]' used to return a reference to an image or vector value @@ -14022,22 +14013,28 @@ namespace cimg_library_suffixed { if (reserved_label['i']!=~0U) _cimg_mp_return(reserved_label['i']); _cimg_mp_scalar0(mp_i); case 'I' : + _cimg_mp_op("Variable 'I'"); if (reserved_label['I']!=~0U) _cimg_mp_return(reserved_label['I']); - _cimg_mp_check_vector0(imgin._spectrum,"variable 'I'"); + _cimg_mp_check_vector0(imgin._spectrum); + need_input_copy = true; pos = vector(imgin._spectrum); CImg::vector((uptrT)mp_Joff,pos,0,0).move_to(code); _cimg_mp_return(pos); case 'R' : if (reserved_label['R']!=~0U) _cimg_mp_return(reserved_label['R']); + need_input_copy = true; _cimg_mp_scalar6(mp_ixyzc,_cimg_mp_x,_cimg_mp_y,_cimg_mp_z,0,0,0); case 'G' : if (reserved_label['G']!=~0U) _cimg_mp_return(reserved_label['G']); + need_input_copy = true; _cimg_mp_scalar6(mp_ixyzc,_cimg_mp_x,_cimg_mp_y,_cimg_mp_z,1,0,0); case 'B' : if (reserved_label['B']!=~0U) _cimg_mp_return(reserved_label['B']); + need_input_copy = true; _cimg_mp_scalar6(mp_ixyzc,_cimg_mp_x,_cimg_mp_y,_cimg_mp_z,2,0,0); case 'A' : if (reserved_label['A']!=~0U) _cimg_mp_return(reserved_label['A']); + need_input_copy = true; _cimg_mp_scalar6(mp_ixyzc,_cimg_mp_x,_cimg_mp_y,_cimg_mp_z,3,0,0); } else if (ss2==se) { // Two-chars variable @@ -14048,6 +14045,7 @@ namespace cimg_library_suffixed { if (*ss1>='0' && *ss1<='9') { // i0...i9 pos = 19 + *ss1 - '0'; if (reserved_label[pos]!=~0U) _cimg_mp_return(reserved_label[pos]); + need_input_copy = true; _cimg_mp_scalar6(mp_ixyzc,_cimg_mp_x,_cimg_mp_y,_cimg_mp_z,pos - 19,0,0); } switch (*ss1) { @@ -14109,7 +14107,7 @@ namespace cimg_library_suffixed { cimg::strpare(variable_name); const unsigned int l_variable_name = (unsigned int)std::strlen(variable_name); char *const ve1 = ss + l_variable_name - 1; - s_op = "Operator '='"; + _cimg_mp_op("Operator '='"); // Assign image value (direct). if (l_variable_name>2 && (*ss=='i' || *ss=='j' || *ss=='I' || *ss=='J') && (*ss1=='(' || *ss1=='[') && @@ -14117,9 +14115,11 @@ namespace cimg_library_suffixed { is_relative = *ss=='j' || *ss=='J'; if (*ss1=='[' && *ve1==']') { // i/j/I/J[_#ind,offset] = value + is_parallelizable = false; if (*ss2=='#') { // Index specified s0 = ss3; while (s0='i'?1:3,p2); + _cimg_mp_check_type(arg2,2,*ss>='i'?1:3,p2); if (p_ref) { *p_ref = _cimg_mp_is_vector(arg2)?4:2; @@ -14172,9 +14172,11 @@ namespace cimg_library_suffixed { } if (*ss1=='(' && *ve1==')') { // i/j/I/J(_#ind,_x,_y,_z,_c) = value + is_parallelizable = false; if (*ss2=='#') { // Index specified s0 = ss3; while (s01) { arg2 = arg1 + 1; if (p2>2) { @@ -14212,9 +14214,9 @@ namespace cimg_library_suffixed { p3 = (unsigned int)cimg::mod((int)mem[p1],listin.width()); p2 = listin[p3]._spectrum; } - _cimg_mp_check_vector0(p2,s_op); + _cimg_mp_check_vector0(p2); } else p2 = 0; - _cimg_mp_check_type(arg5,2,s_op,*ss>='i'?1:3,p2); + _cimg_mp_check_type(arg5,2,*ss>='i'?1:3,p2); if (p_ref) { *p_ref = _cimg_mp_is_vector(arg5)?5:3; @@ -14272,7 +14274,7 @@ namespace cimg_library_suffixed { arg1 = ~0U; // Vector slot arg2 = compile(++s0,ve1,depth1,0); // Index arg3 = compile(s + 1,se,depth1,0); // Value to assign - _cimg_mp_check_type(arg3,2,s_op,1,0); + _cimg_mp_check_type(arg3,2,1,0); if (variable_name[1]) { // Multi-char variable cimglist_for(variable_def,i) if (!std::strcmp(variable_name,variable_def[i])) { @@ -14307,8 +14309,8 @@ namespace cimg_library_suffixed { } } - // Assign user-defined function. - if (l_variable_name>3 && *ve1==')' && *ss!='(') { + // Assign user-defined macro. + if (l_variable_name>2 && *ve1==')' && *ss!='(') { s0 = ve1; while (s0>ss && *s0!='(') --s0; is_sth = std::strncmp(variable_name,"debug(",6) && std::strncmp(variable_name,"print(",6); // is_valid_function_name? @@ -14454,7 +14456,7 @@ namespace cimg_library_suffixed { } } else { // Variable already exists -> assign a new value - _cimg_mp_check_type(arg2,2,s_op,_cimg_mp_is_vector(arg1)?3:1,0); + _cimg_mp_check_type(arg2,2,_cimg_mp_is_vector(arg1)?3:1,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1)) { // Vector if (_cimg_mp_is_vector(arg2)) // From vector CImg::vector((uptrT)mp_vector_copy,arg1,arg2,(uptrT)_cimg_mp_vector_size(arg1)). @@ -14478,7 +14480,7 @@ namespace cimg_library_suffixed { arg2 = compile(s + 1,se,depth1,0); // Value to assign if (*ref==1) { // Vector value (scalar): V[k] = scalar - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg2,2,1,0); arg3 = ref[1]; // Vector slot arg4 = ref[2]; // Index if (p_ref) std::memcpy(p_ref,ref,ref._width*sizeof(unsigned int)); @@ -14488,7 +14490,8 @@ namespace cimg_library_suffixed { } if (*ref==2) { // Image value (scalar): i/j[_#ind,off] = scalar - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg2,2,1,0); + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // Offset @@ -14506,7 +14509,8 @@ namespace cimg_library_suffixed { } if (*ref==3) { // Image value (scalar): i/j(_#ind,_x,_y,_z,_c) = scalar - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg2,2,1,0); + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // X @@ -14527,7 +14531,8 @@ namespace cimg_library_suffixed { } if (*ref==4) { // Image value (vector): I/J[_#ind,off] = value - _cimg_mp_check_type(arg2,2,s_op,3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // Offset @@ -14553,7 +14558,8 @@ namespace cimg_library_suffixed { } if (*ref==5) { // Image value (vector): I/J(_#ind,_x,_y,_z,_c) = value - _cimg_mp_check_type(arg2,2,s_op,3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // X @@ -14581,7 +14587,7 @@ namespace cimg_library_suffixed { } if (_cimg_mp_is_vector(arg1)) { // Vector variable: V = value - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg2,2,1,0); if (_cimg_mp_is_vector(arg2)) // From vector CImg::vector((uptrT)mp_vector_copy,arg1,arg2,(uptrT)_cimg_mp_vector_size(arg1)). move_to(code); @@ -14592,7 +14598,7 @@ namespace cimg_library_suffixed { } if (_cimg_mp_is_variable(arg1)) { // Scalar variable: s = scalar - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg2,2,1,0); CImg::vector((uptrT)mp_copy,arg1,arg2).move_to(code); _cimg_mp_return(arg1); @@ -14615,28 +14621,28 @@ namespace cimg_library_suffixed { for (s = se2, ps = se3, ns = ps - 1; s>ss1; --s, --ps, --ns) // Here, ns = ps - 1 if (*s=='=' && (*ps=='*' || *ps=='/' || *ps=='^') && *ns==*ps && level[s - expr._data]==clevel) { // Self-operators for complex numbers only (**=,//=,^^=) - s_op = *ps=='*'?"Operator '**='":*ps=='/'?"Operator '//='":"Operator '^^='"; + _cimg_mp_op(*ps=='*'?"Operator '**='":*ps=='/'?"Operator '//='":"Operator '^^='"); ref.assign(7); arg1 = compile(ss,ns,depth1,ref); // Vector slot arg2 = compile(s + 1,se,depth1,0); // Right operand if (*ps!='*') { - _cimg_mp_check_type(arg1,2,s_op,2,2); - _cimg_mp_check_type(arg2,2,s_op,2,2); + _cimg_mp_check_type(arg1,2,2,2); + _cimg_mp_check_type(arg2,2,2,2); } if (_cimg_mp_is_vector(arg2)) { // Complex **= complex or Matrix **= matrix if (*ps=='*') { if (_cimg_mp_vector_size(arg1)==2 && _cimg_mp_vector_size(arg2)==2) CImg::vector((uptrT)mp_complex_mul,arg1,arg1,arg2).move_to(code); else { - _cimg_mp_check_matrix_square(arg2,2,s_op); + _cimg_mp_check_matrix_square(arg2,2); p3 = _cimg_mp_vector_size(arg1); p2 = (unsigned int)std::sqrt((float)_cimg_mp_vector_size(arg2)); p1 = p3/p2; if (p1*p2!=p3) { *se = saved_char; cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " - "CImg<%s>::%s: %s: Sizes of left-hand and right-hand operands " + "CImg<%s>::%s: %s: Types of left-hand and right-hand operands " "('%s' and '%s') do not match, in expression '%s%s%s'.", pixel_type(),_cimg_mp_calling_function,s_op, s_type(arg1)._data,s_type(arg2)._data, @@ -14651,16 +14657,14 @@ namespace cimg_library_suffixed { else CImg::vector((uptrT)mp_complex_pow_vv,arg1,arg1,arg2).move_to(code); } else { // Complex **= scalar - if (*ps=='*') - CImg::vector((uptrT)mp_self_map_vector_s,arg1,2,(uptrT)mp_self_mul,arg2).move_to(code); - else if (*ps=='/') - CImg::vector((uptrT)mp_self_map_vector_s,arg1,2,(uptrT)mp_self_div,arg2).move_to(code); - else - CImg::vector((uptrT)mp_complex_pow_vs,arg1,arg1,arg2).move_to(code); + if (*ps=='*') self_vector_s(arg1,mp_self_mul,arg2); + else if (*ps=='/') self_vector_s(arg1,mp_self_div,arg2); + else CImg::vector((uptrT)mp_complex_pow_vs,arg1,arg1,arg2).move_to(code); } // Write computed value back in image if necessary. if (*ref==4) { // Image value (vector): I/J[_#ind,off] **= value + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // Offset @@ -14676,6 +14680,7 @@ namespace cimg_library_suffixed { } } else if (*ref==5) { // Image value (vector): I/J(_#ind,_x,_y,_z,_c) **= value + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // X @@ -14702,16 +14707,16 @@ namespace cimg_library_suffixed { (*ps=='>' && *ns=='>') || (*ps=='<' && *ns=='<')) && level[s - expr._data]==clevel) { // Self-operators (+=,-=,*=,/=,%=,>>=,<<=,&=,^=,|=) switch (*ps) { - case '+' : op = mp_self_add; s_op = "Operator '+='"; break; - case '-' : op = mp_self_sub; s_op = "Operator '-='"; break; - case '*' : op = mp_self_mul; s_op = "Operator '*='"; break; - case '/' : op = mp_self_div; s_op = "Operator '/='"; break; - case '%' : op = mp_self_modulo; s_op = "Operator '%='"; break; - case '<' : op = mp_self_bitwise_left_shift; s_op = "Operator '<<='"; break; - case '>' : op = mp_self_bitwise_right_shift; s_op = "Operator '>=='"; break; - case '&' : op = mp_self_bitwise_and; s_op = "Operator '&='"; break; - case '|' : op = mp_self_bitwise_or; s_op = "Operator '|='"; break; - default : op = mp_self_pow; s_op = "Operator '^='"; break; + case '+' : op = mp_self_add; _cimg_mp_op("Operator '+='"); break; + case '-' : op = mp_self_sub; _cimg_mp_op("Operator '-='"); break; + case '*' : op = mp_self_mul; _cimg_mp_op("Operator '*='"); break; + case '/' : op = mp_self_div; _cimg_mp_op("Operator '/='"); break; + case '%' : op = mp_self_modulo; _cimg_mp_op("Operator '%='"); break; + case '<' : op = mp_self_bitwise_left_shift; _cimg_mp_op("Operator '<<='"); break; + case '>' : op = mp_self_bitwise_right_shift; _cimg_mp_op("Operator '>=='"); break; + case '&' : op = mp_self_bitwise_and; _cimg_mp_op("Operator '&='"); break; + case '|' : op = mp_self_bitwise_or; _cimg_mp_op("Operator '|='"); break; + default : op = mp_self_pow; _cimg_mp_op("Operator '^='"); break; } s1 = *ps=='>' || *ps=='<'?ns:ps; @@ -14725,7 +14730,7 @@ namespace cimg_library_suffixed { } if (*ref==1) { // Vector value (scalar): V[k] += scalar - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg2,2,1,0); arg3 = ref[1]; // Vector slot arg4 = ref[2]; // Index if (p_ref) std::memcpy(p_ref,ref,ref._width*sizeof(unsigned int)); @@ -14736,7 +14741,8 @@ namespace cimg_library_suffixed { } if (*ref==2) { // Image value (scalar): i/j[_#ind,off] += scalar - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg2,2,1,0); + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // Offset @@ -14755,7 +14761,8 @@ namespace cimg_library_suffixed { } if (*ref==3) { // Image value (scalar): i/j(_#ind,_x,_y,_z,_c) += scalar - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg2,2,1,0); + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // X @@ -14777,17 +14784,13 @@ namespace cimg_library_suffixed { } if (*ref==4) { // Image value (vector): I/J[_#ind,off] += value - _cimg_mp_check_type(arg2,2,s_op,3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // Offset if (p_ref) std::memcpy(p_ref,ref,ref._width*sizeof(unsigned int)); - if (_cimg_mp_is_scalar(arg2)) - CImg::vector((uptrT)mp_self_map_vector_s,arg1,(uptrT)_cimg_mp_vector_size(arg1),(uptrT)op,arg2). - move_to(code); - else - CImg::vector((uptrT)mp_self_map_vector_v,arg1,(uptrT)_cimg_mp_vector_size(arg1),(uptrT)op,arg2). - move_to(code); + if (_cimg_mp_is_scalar(arg2)) self_vector_s(arg1,op,arg2); else self_vector_v(arg1,op,arg2); if (p1!=~0U) { if (!listout) _cimg_mp_return(arg1); CImg::vector((uptrT)(is_relative?mp_list_set_Joff_v:mp_list_set_Ioff_v), @@ -14800,20 +14803,16 @@ namespace cimg_library_suffixed { _cimg_mp_return(arg1); } - if (*ref==5) { // Image value (vector): I/J(_#ind,_x,_y,_z,_c) = value - _cimg_mp_check_type(arg2,2,s_op,3,_cimg_mp_vector_size(arg1)); + if (*ref==5) { // Image value (vector): I/J(_#ind,_x,_y,_z,_c) += value + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // X arg4 = ref[4]; // Y arg5 = ref[5]; // Z if (p_ref) std::memcpy(p_ref,ref,ref._width*sizeof(unsigned int)); - if (_cimg_mp_is_scalar(arg2)) - CImg::vector((uptrT)mp_self_map_vector_s,arg1,(uptrT)_cimg_mp_vector_size(arg1),(uptrT)op,arg2). - move_to(code); - else - CImg::vector((uptrT)mp_self_map_vector_v,arg1,(uptrT)_cimg_mp_vector_size(arg1),(uptrT)op,arg2). - move_to(code); + if (_cimg_mp_is_scalar(arg2)) self_vector_s(arg1,op,arg2); else self_vector_v(arg1,op,arg2); if (p1!=~0U) { if (!listout) _cimg_mp_return(arg1); CImg::vector((uptrT)(is_relative?mp_list_set_Jxyz_v:mp_list_set_Ixyz_v), @@ -14827,18 +14826,14 @@ namespace cimg_library_suffixed { } if (_cimg_mp_is_vector(arg1)) { // Vector variable: V += value - _cimg_mp_check_type(arg2,2,s_op,3,_cimg_mp_vector_size(arg1)); - if (_cimg_mp_is_vector(arg2)) // Vector += vector - CImg::vector((uptrT)mp_self_map_vector_v,arg1,(uptrT)_cimg_mp_vector_size(arg1),(uptrT)op,arg2). - move_to(code); - else // Vector += scalar - CImg::vector((uptrT)mp_self_map_vector_s,arg1,(uptrT)_cimg_mp_vector_size(arg1),(uptrT)op,arg2). - move_to(code); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); + if (_cimg_mp_is_vector(arg2)) self_vector_v(arg1,op,arg2); // Vector += vector + else self_vector_s(arg1,op,arg2); // Vector += scalar _cimg_mp_return(arg1); } if (_cimg_mp_is_variable(arg1)) { // Scalar variable: s += scalar - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg2,2,1,0); CImg::vector((uptrT)op,arg1,arg2).move_to(code); _cimg_mp_return(arg1); } @@ -14857,17 +14852,19 @@ namespace cimg_library_suffixed { for (s = ss1; s::vector((uptrT)mp_if,pos,arg1,arg2,arg3, @@ -14877,12 +14874,12 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='|' && *ns=='|' && level[s - expr._data]==clevel) { // Logical or ('||') - s_op = "Operator '||'"; + _cimg_mp_op("Operator '||'"); arg1 = compile(ss,s,depth1,0); p2 = code._width; arg2 = compile(s + 2,se,depth1,0); - _cimg_mp_check_type(arg1,1,s_op,1,0); - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg1,1,1,0); + _cimg_mp_check_type(arg2,2,1,0); if (_cimg_mp_is_constant(arg1) && _cimg_mp_is_constant(arg2)) _cimg_mp_constant(mem[arg1] || mem[arg2]); pos = scalar(); @@ -14893,12 +14890,12 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='&' && *ns=='&' && level[s - expr._data]==clevel) { // Logical and ('&&') - s_op = "Operator '&&'"; + _cimg_mp_op("Operator '&&'"); arg1 = compile(ss,s,depth1,0); p2 = code._width; arg2 = compile(s + 2,se,depth1,0); - _cimg_mp_check_type(arg1,1,s_op,1,0); - _cimg_mp_check_type(arg2,2,s_op,1,0); + _cimg_mp_check_type(arg1,1,1,0); + _cimg_mp_check_type(arg2,2,1,0); if (_cimg_mp_is_constant(arg1) && _cimg_mp_is_constant(arg2)) _cimg_mp_constant(mem[arg1] && mem[arg2]); pos = scalar(); @@ -14909,9 +14906,10 @@ namespace cimg_library_suffixed { for (s = se2; s>ss; --s) if (*s=='|' && level[s - expr._data]==clevel) { // Bitwise or ('|') + _cimg_mp_op("Operator '|'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 1,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '|'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_bitwise_or,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_bitwise_or,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_bitwise_or,arg1,arg2); @@ -14922,9 +14920,10 @@ namespace cimg_library_suffixed { for (s = se2; s>ss; --s) if (*s=='&' && level[s - expr._data]==clevel) { // Bitwise and ('&') + _cimg_mp_op("Operator '&'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 1,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '&'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_bitwise_and,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_bitwise_and,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_bitwise_and,arg1,arg2); @@ -14935,9 +14934,10 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='!' && *ns=='=' && level[s - expr._data]==clevel) { // Not equal to ('!=') + _cimg_mp_op("Operator '!='"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 2,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '!='",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_neq,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_neq,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_neq,arg1,arg2); @@ -14947,9 +14947,10 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='=' && *ns=='=' && level[s - expr._data]==clevel) { // Equal to ('==') + _cimg_mp_op("Operator '=='"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 2,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '=='",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_eq,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_eq,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_eq,arg1,arg2); @@ -14959,9 +14960,10 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='<' && *ns=='=' && level[s - expr._data]==clevel) { // Less or equal than ('<=') + _cimg_mp_op("Operator '<='"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 2,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '<='",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_lte,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_lte,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_lte,arg1,arg2); @@ -14971,9 +14973,10 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='>' && *ns=='=' && level[s - expr._data]==clevel) { // Greater or equal than ('>=') + _cimg_mp_op("Operator '>='"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 2,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '>='",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_gte,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_gte,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_gte,arg1,arg2); @@ -14983,9 +14986,10 @@ namespace cimg_library_suffixed { for (s = se2, ns = se1, ps = se3; s>ss; --s, --ns, --ps) if (*s=='<' && *ns!='<' && *ps!='<' && level[s - expr._data]==clevel) { // Less than ('<') + _cimg_mp_op("Operator '<'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 1,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '<'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_lt,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_lt,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_lt,arg1,arg2); @@ -14995,9 +14999,10 @@ namespace cimg_library_suffixed { for (s = se2, ns = se1, ps = se3; s>ss; --s, --ns, --ps) if (*s=='>' && *ns!='>' && *ps!='>' && level[s - expr._data]==clevel) { // Greather than ('>') + _cimg_mp_op("Operator '>'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 1,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '>'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_gt,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_gt,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_gt,arg1,arg2); @@ -15007,9 +15012,10 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='<' && *ns=='<' && level[s - expr._data]==clevel) { // Left bit shift ('<<') + _cimg_mp_op("Operator '<<'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 2,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '<<'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_bitwise_left_shift,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) @@ -15023,9 +15029,10 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='>' && *ns=='>' && level[s - expr._data]==clevel) { // Right bit shift ('>>') + _cimg_mp_op("Operator '>>'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 2,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '>>'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_bitwise_right_shift,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) @@ -15043,9 +15050,10 @@ namespace cimg_library_suffixed { (*ps!='e' || !(ps - pexpr._data>ss - expr._data && (*(ps - 1)=='.' || (*(ps - 1)>='0' && *(ps - 1)<='9')))) && level[s - expr._data]==clevel) { // Addition ('+') + _cimg_mp_op("Operator '+'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 1,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '+'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_add,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_add,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_add,arg1,arg2); @@ -15061,9 +15069,10 @@ namespace cimg_library_suffixed { (*ps!='e' || !(ps - pexpr._data>ss - expr._data && (*(ps - 1)=='.' || (*(ps - 1)>='0' && *(ps - 1)<='9')))) && level[s - expr._data]==clevel) { // Subtraction ('-') + _cimg_mp_op("Operator '-'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 1,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '-'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_sub,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_sub,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_sub,arg1,arg2); @@ -15074,7 +15083,7 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='*' && *ns=='*' && level[s - expr._data]==clevel) { // Complex/matrix multiplication ('**') - s_op = "Operator '**'"; + _cimg_mp_op("Operator '**'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 2,se,depth1,0); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) { @@ -15089,7 +15098,7 @@ namespace cimg_library_suffixed { if (arg4*p2!=p1) { *se = saved_char; cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " - "CImg<%s>::%s: %s: Sizes of left-hand and right-hand operands " + "CImg<%s>::%s: %s: Types of left-hand and right-hand operands " "('%s' and '%s') do not match, in expression '%s%s%s'.", pixel_type(),_cimg_mp_calling_function,s_op, s_type(arg1)._data,s_type(arg2)._data, @@ -15110,11 +15119,11 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='/' && *ns=='/' && level[s - expr._data]==clevel) { // Complex division ('//') - s_op = "Operator '//'"; + _cimg_mp_op("Operator '//'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 2,se,depth1,0); - _cimg_mp_check_type(arg1,1,s_op,3,2); - _cimg_mp_check_type(arg2,2,s_op,3,2); + _cimg_mp_check_type(arg1,1,3,2); + _cimg_mp_check_type(arg2,2,3,2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) { pos = vector(2); CImg::vector((uptrT)mp_complex_div_vv,pos,arg1,arg2).move_to(code); @@ -15131,9 +15140,10 @@ namespace cimg_library_suffixed { } for (s = se2; s>ss; --s) if (*s=='*' && level[s - expr._data]==clevel) { // Multiplication ('*') + _cimg_mp_op("Operator '*'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 1,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '*'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_mul,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_mul,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_mul,arg1,arg2); @@ -15143,9 +15153,10 @@ namespace cimg_library_suffixed { for (s = se2; s>ss; --s) if (*s=='/' && level[s - expr._data]==clevel) { // Division ('/') + _cimg_mp_op("Operator '/'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 1,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '/'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_div,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_div,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_div,arg1,arg2); @@ -15155,9 +15166,10 @@ namespace cimg_library_suffixed { for (s = se2, ns = se1; s>ss; --s, --ns) if (*s=='%' && *ns!='^' && level[s - expr._data]==clevel) { // Modulo ('%') + _cimg_mp_op("Operator '%'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 1,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '%'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_modulo,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_modulo,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_modulo,arg1,arg2); @@ -15167,10 +15179,13 @@ namespace cimg_library_suffixed { } if (se1>ss) { - if (*ss=='+' && (*ss1!='+' || (ss2='0' && *ss2<='9'))) // Unary plus ('+') + if (*ss=='+' && (*ss1!='+' || (ss2='0' && *ss2<='9'))) { // Unary plus ('+') + _cimg_mp_op("Operator '+'"); _cimg_mp_return(compile(ss1,se,depth1,0)); + } if (*ss=='-' && (*ss1!='-' || (ss2='0' && *ss2<='9'))) { // Unary minus ('-') + _cimg_mp_op("Operator '-'"); arg1 = compile(ss1,se,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_minus,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(-mem[arg1]); @@ -15178,6 +15193,7 @@ namespace cimg_library_suffixed { } if (*ss=='!') { // Logical not ('!') + _cimg_mp_op("Operator '!'"); arg1 = compile(ss1,se,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_logical_not,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(!mem[arg1]); @@ -15185,6 +15201,7 @@ namespace cimg_library_suffixed { } if (*ss=='~') { // Bitwise not ('~') + _cimg_mp_op("Operator '~'"); arg1 = compile(ss1,se,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_bitwise_not,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(~(unsigned long)mem[arg1]); @@ -15194,11 +15211,11 @@ namespace cimg_library_suffixed { for (s = se3, ns = se2; s>ss; --s, --ns) if (*s=='^' && *ns=='^' && level[s - expr._data]==clevel) { // Complex power ('^^') - s_op = "Operator '^^'"; + _cimg_mp_op("Operator '^^'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 2,se,depth1,0); - _cimg_mp_check_type(arg1,1,s_op,3,2); - _cimg_mp_check_type(arg2,2,s_op,3,2); + _cimg_mp_check_type(arg1,1,3,2); + _cimg_mp_check_type(arg2,2,3,2); pos = (_cimg_mp_is_vector(arg1) || _cimg_mp_is_vector(arg2))?vector(2):0; if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) { CImg::vector((uptrT)mp_complex_pow_vv,pos,arg1,arg2).move_to(code); @@ -15226,9 +15243,10 @@ namespace cimg_library_suffixed { for (s = se2; s>ss; --s) if (*s=='^' && level[s - expr._data]==clevel) { // Power ('^') + _cimg_mp_op("Operator '^'"); arg1 = compile(ss,s,depth1,0); arg2 = compile(s + 1,se,depth1,0); - _cimg_mp_check_type(arg2,2,"operator '^'",3,_cimg_mp_vector_size(arg1)); + _cimg_mp_check_type(arg2,2,3,_cimg_mp_vector_size(arg1)); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_vv(mp_pow,arg1,arg2); if (_cimg_mp_is_vector(arg1) && _cimg_mp_is_scalar(arg2)) _cimg_mp_vector2_vs(mp_pow,arg1,arg2); if (_cimg_mp_is_scalar(arg1) && _cimg_mp_is_vector(arg2)) _cimg_mp_vector2_sv(mp_pow,arg1,arg2); @@ -15246,9 +15264,13 @@ namespace cimg_library_suffixed { is_sth = ss1ss && (*se1=='+' || *se1=='-') && *se2==*se1)) { // Pre/post-decrement and increment - if ((is_sth && *ss=='+') || (!is_sth && *se1=='+')) { op = mp_self_increment; s_op = "Operator '++'"; } - else { op = mp_self_decrement; s_op = "Operator '--'"; } - + if ((is_sth && *ss=='+') || (!is_sth && *se1=='+')) { + _cimg_mp_op("Operator '++'"); + op = mp_self_increment; + } else { + _cimg_mp_op("Operator '--'"); + op = mp_self_decrement; + } ref.assign(7); arg1 = is_sth?compile(ss2,se,depth1,ref):compile(ss,se2,depth1,ref); // Variable slot @@ -15274,6 +15296,7 @@ namespace cimg_library_suffixed { } if (*ref==2) { // Image value (scalar): i/j[_#ind,off]++ + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // Offset @@ -15292,6 +15315,7 @@ namespace cimg_library_suffixed { } if (*ref==3) { // Image value (scalar): i/j(_#ind,_x,_y,_z,_c)++ + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // X @@ -15313,13 +15337,12 @@ namespace cimg_library_suffixed { } if (*ref==4) { // Image value (vector): I/J[_#ind,off]++ + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // Offset if (p_ref) std::memcpy(p_ref,ref,ref._width*sizeof(unsigned int)); - CImg::vector((uptrT)mp_self_map_vector_s,arg1,(uptrT)_cimg_mp_vector_size(arg1), - (uptrT)(op==mp_self_increment?mp_self_add:mp_self_sub),1). - move_to(code); + self_vector_s(arg1,op==mp_self_increment?mp_self_add:mp_self_sub,1); if (p1!=~0U) { if (!listout) _cimg_mp_return(pos); CImg::vector((uptrT)(is_relative?mp_list_set_Joff_v:mp_list_set_Ioff_v), @@ -15333,15 +15356,14 @@ namespace cimg_library_suffixed { } if (*ref==5) { // Image value (vector): I/J(_#ind,_x,_y,_z,_c)++ + is_parallelizable = false; p1 = ref[1]; // Index is_relative = (bool)ref[2]; arg3 = ref[3]; // X arg4 = ref[4]; // Y arg5 = ref[5]; // Z if (p_ref) std::memcpy(p_ref,ref,ref._width*sizeof(unsigned int)); - CImg::vector((uptrT)mp_self_map_vector_s,arg1,(uptrT)_cimg_mp_vector_size(arg1), - (uptrT)(op==mp_self_increment?mp_self_add:mp_self_sub),1). - move_to(code); + self_vector_s(arg1,op==mp_self_increment?mp_self_add:mp_self_sub,1); if (p1!=~0U) { if (!listout) _cimg_mp_return(pos); CImg::vector((uptrT)(is_relative?mp_list_set_Jxyz_v:mp_list_set_Ixyz_v), @@ -15355,9 +15377,7 @@ namespace cimg_library_suffixed { } if (_cimg_mp_is_vector(arg1)) { // Vector variable: V++ - CImg::vector((uptrT)mp_self_map_vector_s,arg1,(uptrT)_cimg_mp_vector_size(arg1), - (uptrT)(op==mp_self_increment?mp_self_add:mp_self_sub),1). - move_to(code); + self_vector_s(arg1,op==mp_self_increment?mp_self_add:mp_self_sub,1); _cimg_mp_return(pos); } @@ -15382,7 +15402,7 @@ namespace cimg_library_suffixed { // Array-like access to vectors and image values 'i/j[_#ind,offset,_boundary]' and 'vector[offset]'. if (*se1==']' && *ss!='[') { - s_op = "Operator '[]'"; + _cimg_mp_op("Operator '[]'"); is_relative = *ss=='j' || *ss=='J'; if ((*ss=='I' || *ss=='J') && *ss1=='[' && @@ -15390,10 +15410,11 @@ namespace cimg_library_suffixed { if (*ss2=='#') { // Index specified s0 = ss3; while (s0::vector((uptrT)(is_relative?mp_list_Joff:mp_list_Ioff), pos,p1,arg1,arg2==~0U?reserved_label[30]:arg2).move_to(code); } else { + need_input_copy = true; CImg::vector((uptrT)(is_relative?mp_Joff:mp_Ioff), pos,arg1,arg2==~0U?reserved_label[30]:arg2).move_to(code); } @@ -15428,7 +15450,7 @@ namespace cimg_library_suffixed { } else { p1 = ~0U; s0 = ss2; } s1 = s0; while (s1 sub-vector extraction arg2 = compile(++s0,s1,depth1,0); arg3 = compile(++s1,se1,depth1,0); - _cimg_mp_check_constant(arg2,1,s_op,false); - _cimg_mp_check_constant(arg3,2,s_op,false); + _cimg_mp_check_constant(arg2,1,false); + _cimg_mp_check_constant(arg3,2,false); p1 = (unsigned int)mem[arg2]; p2 = (unsigned int)mem[arg3]; p3 = _cimg_mp_vector_size(arg1); @@ -15527,12 +15550,14 @@ namespace cimg_library_suffixed { if (*se1==')') { if (*ss=='(') _cimg_mp_return(compile(ss1,se1,depth1,p_ref)); // Simple parentheses is_relative = *ss=='j' || *ss=='J'; + _cimg_mp_op("Operator '()'"); - // I/J(_#ind,_x,_y,_z,_c,_interpolation,_boundary) + // I/J(_#ind,_x,_y,_z,_interpolation,_boundary) if ((*ss=='I' || *ss=='J') && *ss1=='(') { // Image value as scalar if (*ss2=='#') { // Index specified s0 = ss3; while (s01) { arg2 = arg1 + 1; if (p2>2) arg3 = arg2 + 1; @@ -15551,7 +15576,7 @@ namespace cimg_library_suffixed { if (s1::vector((uptrT)(is_relative?mp_list_Jxyz:mp_list_Ixyz), pos,p1,arg1,arg2,arg3, arg4==~0U?reserved_label[29]:arg4, arg5==~0U?reserved_label[30]:arg5).move_to(code); - else + else { + need_input_copy = true; CImg::vector((uptrT)(is_relative?mp_Jxyz:mp_Ixyz), pos,arg1,arg2,arg3, arg4==~0U?reserved_label[29]:arg4, arg5==~0U?reserved_label[30]:arg5).move_to(code); + } _cimg_mp_return(pos); } @@ -15614,9 +15641,9 @@ namespace cimg_library_suffixed { if (s01) { arg2 = arg1 + 1; if (p2>2) { @@ -15627,7 +15654,7 @@ namespace cimg_library_suffixed { if (s1::vector((uptrT)mp_complex_conj,pos,arg1).move_to(code); _cimg_mp_return(pos); } if (!std::strncmp(ss,"cexp(",5)) { // Complex exponential + _cimg_mp_op("Function 'cexp()'"); arg1 = compile(ss5,se1,depth1,0); - _cimg_mp_check_type(arg1,0,"Function 'cexp()'",2,2); + _cimg_mp_check_type(arg1,0,2,2); pos = vector(2); CImg::vector((uptrT)mp_complex_exp,pos,arg1).move_to(code); _cimg_mp_return(pos); } if (!std::strncmp(ss,"clog(",5)) { // Complex logarithm + _cimg_mp_op("Function 'clog()'"); arg1 = compile(ss5,se1,depth1,0); - _cimg_mp_check_type(arg1,0,"Function 'clog()'",2,2); + _cimg_mp_check_type(arg1,0,2,2); pos = vector(2); CImg::vector((uptrT)mp_complex_log,pos,arg1).move_to(code); _cimg_mp_return(pos); } + if (!std::strncmp(ss,"copy(",5)) { // Memory copy + _cimg_mp_op("Function 'copy()'"); + ref.assign(14); + s1 = ss5; while (s1(1,21).move_to(code); + code.back().get_shared_rows(0,6).fill((uptrT)mp_memcopy,p1,arg1,arg2,arg3,arg4,arg5); + code.back().get_shared_rows(7,20).fill(ref); + _cimg_mp_return(p1); + } + if (!std::strncmp(ss,"cos(",4)) { // Cosine + _cimg_mp_op("Function 'cos()'"); arg1 = compile(ss4,se1,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_cos,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(std::cos(mem[arg1])); @@ -15776,30 +15848,161 @@ namespace cimg_library_suffixed { } if (!std::strncmp(ss,"cosh(",5)) { // Hyperbolic cosine + _cimg_mp_op("Function 'cosh()'"); arg1 = compile(ss5,se1,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_cosh,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(std::cosh(mem[arg1])); _cimg_mp_scalar1(mp_cosh,arg1); } + if (!std::strncmp(ss,"crop(",5)) { // Image crop + _cimg_mp_op("Function 'crop()'"); + if (*ss5=='#') { // Index specified + s0 = ss6; while (s0::sequence((uptrT)_cimg_mp_vector_size(arg1),arg1 + 1, + arg1 + (uptrT)_cimg_mp_vector_size(arg1)); + opcode.resize(1,cimg::min(opcode._height,4U),1,1,0).move_to(_opcode); + is_sth = true; + } else { + _cimg_mp_check_type(arg1,pos + 1,1,0); + CImg::vector(arg1).move_to(_opcode); + } + s = ns; + } + (_opcode>'y').move_to(opcode); + + arg1 = 0; arg2 = p1!=~0U?1:0; + switch (opcode._height) { + case 0 : case 1 : + CImg::vector(0,0,0,0,~0U,~0U,~0U,~0U,0).move_to(opcode); + break; + case 2 : + CImg::vector(*opcode,0,0,0,opcode[1],~0U,~0U,~0U,reserved_label[30]).move_to(opcode); + arg1 = arg2?3:2; + break; + case 3 : + CImg::vector(*opcode,0,0,0,opcode[1],~0U,~0U,~0U,opcode[2]).move_to(opcode); + arg1 = arg2?3:2; + break; + case 4 : + CImg::vector(*opcode,opcode[1],0,0,opcode[2],opcode[3],~0U,~0U,reserved_label[30]). + move_to(opcode); + arg1 = (is_sth?2:1) + arg2; + break; + case 5 : + CImg::vector(*opcode,opcode[1],0,0,opcode[2],opcode[3],~0U,~0U,opcode[4]). + move_to(opcode); + arg1 = (is_sth?2:1) + arg2; + break; + case 6 : + CImg::vector(*opcode,opcode[1],opcode[2],0,opcode[3],opcode[4],opcode[5],~0U, + reserved_label[30]).move_to(opcode); + arg1 = (is_sth?2:4) + arg2; + break; + case 7 : + CImg::vector(*opcode,opcode[1],opcode[2],0,opcode[3],opcode[4],opcode[5],~0U, + opcode[6]).move_to(opcode); + arg1 = (is_sth?2:4) + arg2; + break; + case 8 : + CImg::vector(*opcode,opcode[1],opcode[2],opcode[3],opcode[4],opcode[5],opcode[6], + opcode[7],reserved_label[30]).move_to(opcode); + arg1 = (is_sth?2:5) + arg2; + break; + case 9 : + arg1 = (is_sth?2:5) + arg2; + break; + default : // Error -> too much arguments + throw CImgArgumentException("[_cimg_math_parser] " + "CImg<%s>::%s: %s: Too much arguments specified, " + "in expression '%s%s%s'.", + pixel_type(),_cimg_mp_calling_function,s_op, + (ss - 4)>expr._data?"...":"", + (ss - 4)>expr._data?ss - 4:expr._data, + se<&expr.back()?"...":""); + } + _cimg_mp_check_type(*opcode,arg2 + 1,1,0); + _cimg_mp_check_type(opcode[1],arg2 + (is_sth?0:1),1,0); + _cimg_mp_check_type(opcode[2],arg2 + (is_sth?0:2),1,0); + _cimg_mp_check_type(opcode[3],arg2 + (is_sth?0:3),1,0); + + if (opcode[4]!=(uptrT)~0U) { + _cimg_mp_check_constant(opcode[4],arg1,true); + opcode[4] = (uptrT)mem[opcode[4]]; + } + if (opcode[5]!=(uptrT)~0U) { + _cimg_mp_check_constant(opcode[5],arg1 + 1,true); + opcode[5] = (uptrT)mem[opcode[5]]; + } + if (opcode[6]!=(uptrT)~0U) { + _cimg_mp_check_constant(opcode[6],arg1 + 2,true); + opcode[6] = (uptrT)mem[opcode[6]]; + } + if (opcode[7]!=(uptrT)~0U) { + _cimg_mp_check_constant(opcode[7],arg1 + 3,true); + opcode[7] = (uptrT)mem[opcode[7]]; + } + _cimg_mp_check_type(opcode[8],arg1 + 4,1,0); + + if (opcode[4]==(uptrT)~0U || opcode[5]==(uptrT)~0U || + opcode[6]==(uptrT)~0U || opcode[7]==(uptrT)~0U) { + if (p1!=~0U) { + _cimg_mp_check_constant(p1,1,false); + p1 = (unsigned int)cimg::mod((int)mem[p1],listin.width()); + } + const CImg &img = p1!=~0U?listin[p1]:imgin; + if (!img) + throw CImgArgumentException("[_cimg_math_parser] " + "CImg<%s>::%s: %s: Cannot crop empty image when " + "some xyzc-coordinates are unspecified, in expression '%s%s%s'.", + pixel_type(),_cimg_mp_calling_function,s_op, + (ss - 4)>expr._data?"...":"", + (ss - 4)>expr._data?ss - 4:expr._data, + se<&expr.back()?"...":""); + if (opcode[4]==(uptrT)~0U) opcode[4] = (uptrT)img._width; + if (opcode[5]==(uptrT)~0U) opcode[5] = (uptrT)img._height; + if (opcode[6]==(uptrT)~0U) opcode[6] = (uptrT)img._depth; + if (opcode[7]==(uptrT)~0U) opcode[7] = (uptrT)img._spectrum; + } + + pos = vector(opcode[4]*opcode[5]*opcode[6]*opcode[7]); + CImg::vector((uptrT)mp_crop, + pos,p1, + *opcode,opcode[1],opcode[2],opcode[3], + opcode[4],opcode[5],opcode[6],opcode[7], + opcode[8]).move_to(code); + _cimg_mp_return(pos); + } + if (!std::strncmp(ss,"cross(",6)) { // Cross product - s_op = "Function 'cross()"; + _cimg_mp_op("Function 'cross()'"); s1 = ss6; while (s1::vector((uptrT)mp_cross,pos,arg1,arg2).move_to(code); _cimg_mp_return(pos); } if (!std::strncmp(ss,"cut(",4)) { // Cut + _cimg_mp_op("Function 'cut()'"); s1 = ss4; while (s1::vector((uptrT)mp_diag,pos,arg1,p1).move_to(code); + _cimg_mp_return(pos); + } + if (!std::strncmp(ss,"dot(",4)) { // Dot product - s_op = "Function 'dot()'"; + _cimg_mp_op("Function 'dot()'"); s1 = ss4; while (s1::vector((uptrT)mp_dowhile,arg1,arg2,code._width - p1).move_to(code,p1); _cimg_mp_return(arg1); } + + if (!std::strncmp(ss,"draw(",5)) { // Draw image + _cimg_mp_op("Function 'draw()'"); + is_parallelizable = false; + if (*ss5=='#') { // Index specified + s0 = ss6; while (s01) { + arg3 = arg2 + 1; + if (p2>2) { + arg4 = arg3 + 1; + if (p2>3) arg5 = arg4 + 1; + } + } + ++s0; + is_sth = true; + } else { + if (s0::vector((uptrT)mp_draw,arg1,p1,arg2,arg3,arg4,arg5,0,0,0,0,1,(uptrT)-1,0,1). + move_to(opcode); + + arg2 = arg3 = arg4 = arg5 = ~0U; + p2 = p1!=~0U?0:1; + if (s0 &img = p1!=~0U?listout[p1]:imgout; + if (arg2==~0U) arg2 = img._width; + if (arg3==~0U) arg3 = img._height; + if (arg4==~0U) arg4 = img._depth; + if (arg5==~0U) arg5 = img._spectrum; + } + if (arg2*arg3*arg4*arg5!=_cimg_mp_vector_size(arg1)) { + *se = saved_char; cimg::strellipsize(expr,64); + throw CImgArgumentException("[_cimg_math_parser] " + "CImg<%s>::%s: %s: Type of %s argument ('%s') and specified size " + "(%u,%u,%u,%u) do not match, in expression '%s%s%s'.", + pixel_type(),_cimg_mp_calling_function,s_op, + p1==~0U?"first":"second",s_type(arg1)._data, + arg2,arg3,arg4,arg5, + (ss - 4)>expr._data?"...":"", + (ss - 4)>expr._data?ss - 4:expr._data, + se<&expr.back()?"...":""); + } + opcode[7] = (uptrT)arg2; + opcode[8] = (uptrT)arg3; + opcode[9] = (uptrT)arg4; + opcode[10] = (uptrT)arg5; + + if (s0::%s: %s: Type of opacity mask ('%s') and specified size " + "(%u,%u,%u,%u) do not match, in expression '%s%s%s'.", + pixel_type(),_cimg_mp_calling_function,s_op, + s_type(p2)._data, + arg2,arg3,arg4,arg5, + (ss - 4)>expr._data?"...":"", + (ss - 4)>expr._data?ss - 4:expr._data, + se<&expr.back()?"...":""); + } + opcode[12] = p2; + opcode[13] = _cimg_mp_vector_size(p2)/(arg2*arg3*arg4); + p3 = s0::vector((uptrT)mp_eig,pos,arg1,p1).move_to(code); + _cimg_mp_return(pos); + } + if (!std::strncmp(ss,"exp(",4)) { // Exponential + _cimg_mp_op("Function 'exp()'"); arg1 = compile(ss4,se1,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_exp,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(std::exp(mem[arg1])); _cimg_mp_scalar1(mp_exp,arg1); } + + if (!std::strncmp(ss,"eye(",4)) { // Identity matrix + _cimg_mp_op("Function 'eye()'"); + arg1 = compile(ss4,se1,depth1,0); + _cimg_mp_check_constant(arg1,1,true); + p1 = (unsigned int)mem[arg1]; + pos = vector(p1*p1); + CImg::vector((uptrT)mp_eye,pos,p1).move_to(code); + _cimg_mp_return(pos); + } break; case 'f' : if (*ss1=='o' && *ss2=='r' && (*ss3=='(' || (*ss3 && *ss3<=' ' && *ss4=='('))) { // For loop + _cimg_mp_op("Function 'for()'"); if (*ss3<=' ') cimg::swap(*ss3,*ss4); // Allow space before opening brace s1 = ss4; while (s1::vector((uptrT)mp_whiledo,pos,arg1,p2 - p1,code._width - p2,arg2).move_to(code,p1); _cimg_mp_return(pos); @@ -15889,14 +16275,11 @@ namespace cimg_library_suffixed { case 'g' : if (!std::strncmp(ss,"gauss(",6)) { // Gaussian function + _cimg_mp_op("Function 'gauss()'"); s1 = ss6; while (s1=se1?0:compile(s2 + 1,se1,depth1,0); - _cimg_mp_check_type(arg1,1,s_op,1,0); - _cimg_mp_check_type(arg3,3,s_op,3,_cimg_mp_vector_size(arg2)); - if (_cimg_mp_is_constant(arg1) && _cimg_mp_is_constant(arg2) && _cimg_mp_is_constant(arg3)) - _cimg_mp_constant(mem[arg1]?mem[arg2]:mem[arg3]); + arg3 = s2::vector((uptrT)mp_if,pos,arg1,arg2,arg3, @@ -15949,12 +16334,13 @@ namespace cimg_library_suffixed { } if (!std::strncmp(ss,"init(",5)) { // Init + _cimg_mp_op("Function 'init()'"); if (ss0!=expr._data || code.width()) { // (only allowed as the first instruction) *se = saved_char; cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " - "CImg<%s>::%s: Function 'init()': Init invokation not done at the " + "CImg<%s>::%s: %s: Init invokation not done at the " "beginning of expression '%s%s%s'.", - pixel_type(),_cimg_mp_calling_function, + pixel_type(),_cimg_mp_calling_function,s_op, (ss - 4)>expr._data?"...":"", (ss - 4)>expr._data?ss - 4:expr._data, se<&expr.back()?"...":""); @@ -15965,15 +16351,27 @@ namespace cimg_library_suffixed { } if (!std::strncmp(ss,"int(",4)) { // Integer cast + _cimg_mp_op("Function 'int()'"); arg1 = compile(ss4,se1,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_int,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant((long)mem[arg1]); _cimg_mp_scalar1(mp_int,arg1); } + if (!std::strncmp(ss,"inv(",4)) { // Matrix/scalar inversion + _cimg_mp_op("Function 'inv()'"); + arg1 = compile(ss4,se1,depth1,0); + _cimg_mp_check_matrix_square(arg1,1); + p1 = (unsigned int)std::sqrt((float)_cimg_mp_vector_size(arg1)); + pos = vector(p1*p1); + CImg::vector((uptrT)mp_inv,pos,arg1,p1).move_to(code); + _cimg_mp_return(pos); + } + if (*ss1=='s') { // Family of 'is_?()' functions if (!std::strncmp(ss,"isbool(",7)) { // Is boolean? + _cimg_mp_op("Function 'isbool()'"); if (ss7==se1) _cimg_mp_return(0); arg1 = compile(ss7,se1,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_isbool,arg1); @@ -15982,6 +16380,7 @@ namespace cimg_library_suffixed { } if (!std::strncmp(ss,"isdir(",6)) { // Is directory? + _cimg_mp_op("Function 'isdir()'"); *se1 = 0; is_sth = cimg::is_directory(ss6); *se1 = ')'; @@ -15989,6 +16388,7 @@ namespace cimg_library_suffixed { } if (!std::strncmp(ss,"isfile(",7)) { // Is file? + _cimg_mp_op("Function 'isfile()'"); *se1 = 0; is_sth = cimg::is_file(ss7); *se1 = ')'; @@ -15996,6 +16396,8 @@ namespace cimg_library_suffixed { } if (!std::strncmp(ss,"isin(",5)) { // Is in sequence/vector? + if (ss5>=se1) _cimg_mp_return(0); + _cimg_mp_op("Function 'isin()'"); pos = scalar(); CImg::vector((uptrT)mp_isin,pos).move_to(_opcode); for (s = ss5; s::vector((uptrT)mp_matrix_diag,pos,arg1,p1).move_to(code); - _cimg_mp_return(pos); - } - - if (!std::strncmp(ss,"meig(",5)) { // Matrix eigenvalues/eigenvector - arg1 = compile(ss5,se1,depth1,0); - _cimg_mp_check_matrix_square(arg1,1,"Function 'meig()'"); - p1 = (unsigned int)std::sqrt((float)_cimg_mp_vector_size(arg1)); - pos = vector((p1 + 1)*p1); - CImg::vector((uptrT)mp_matrix_eig,pos,arg1,p1).move_to(code); - _cimg_mp_return(pos); - } - - if (!std::strncmp(ss,"meye(",5)) { // Matrix eigenvalues/eigenvector - arg1 = compile(ss5,se1,depth1,0); - _cimg_mp_check_constant(arg1,1,"Function 'meye()'",true); - p1 = (unsigned int)mem[arg1]; - pos = vector(p1*p1); - CImg::vector((uptrT)mp_matrix_eye,pos,p1).move_to(code); - _cimg_mp_return(pos); - } - - if (!std::strncmp(ss,"minv(",5)) { // Matrix inversion - arg1 = compile(ss5,se1,depth1,0); - _cimg_mp_check_matrix_square(arg1,1,"Function 'minv()'"); - p1 = (unsigned int)std::sqrt((float)_cimg_mp_vector_size(arg1)); - pos = vector(p1*p1); - CImg::vector((uptrT)mp_matrix_inv,pos,arg1,p1).move_to(code); - _cimg_mp_return(pos); - } - - if (!std::strncmp(ss,"mmul(",5)) { // Matrix multiplication - s_op = "Function 'mmul()'"; - s1 = ss5; while (s1::%s: %s: Sizes of first and second arguments ('%s' and '%s') " + "CImg<%s>::%s: %s: Types of first and second arguments ('%s' and '%s') " "do not match for third argument 'nb_colsB=%u', " "in expression '%s%s%s'.", pixel_type(),_cimg_mp_calling_function,s_op, @@ -16144,95 +16510,12 @@ namespace cimg_library_suffixed { CImg::vector((uptrT)mp_matrix_mul,pos,arg1,arg2,arg4,arg5,p3).move_to(code); _cimg_mp_return(pos); } - - if (!std::strncmp(ss,"mrot(",5)) { // Rotation matrix - s_op = "Function 'mrot()'"; - s1 = ss5; while (s1::vector((uptrT)mp_matrix_rot,pos,arg1,arg2,arg3,arg4).move_to(code); - _cimg_mp_return(pos); - } - - if (!std::strncmp(ss,"msolve(",7)) { // Solve linear system - s_op = "Function 'msolve()'"; - s1 = ss7; while (s1::%s: %s: Sizes of first and second arguments ('%s' and '%s') " - "do not match for third argument 'nb_colsB=%u', " - "in expression '%s%s%s'.", - pixel_type(),_cimg_mp_calling_function,s_op, - s_type(arg1)._data,s_type(arg2)._data,p3, - (ss - 4)>expr._data?"...":"", - (ss - 4)>expr._data?ss - 4:expr._data, - se<&expr.back()?"...":""); - } - pos = vector(arg5*p3); - CImg::vector((uptrT)mp_matrix_solve,pos,arg1,arg2,arg4,arg5,p3).move_to(code); - _cimg_mp_return(pos); - } - - if (!std::strncmp(ss,"mtrace(",7)) { // Matrix trace - arg1 = compile(ss7,se1,depth1,0); - _cimg_mp_check_matrix_square(arg1,1,"Function 'mtrace()'"); - p1 = (unsigned int)std::sqrt((float)_cimg_mp_vector_size(arg1)); - _cimg_mp_scalar2(mp_matrix_trace,arg1,p1); - } - - if (!std::strncmp(ss,"mtrans(",7)) { // Matrix transpose - s_op = "Function 'mtrans()'"; - s1 = ss7; while (s1::%s: %s: Size of first argument ('%s') does not match" - "for second specified argument 'nb_cols=%u', " - "in expression '%s%s%s'.", - pixel_type(),_cimg_mp_calling_function,s_op, - s_type(arg1)._data,p2, - (ss - 4)>expr._data?"...":"", - (ss - 4)>expr._data?ss - 4:expr._data, - se<&expr.back()?"...":""); - } - pos = vector(p3*p2); - CImg::vector((uptrT)mp_matrix_trans,pos,arg1,p2,p3).move_to(code); - _cimg_mp_return(pos); - } break; case 'n' : if (!std::strncmp(ss,"narg(",5)) { // Number of arguments - if (*ss5==')') _cimg_mp_return(0); + _cimg_mp_op("Function 'narg()'"); + if (ss5>=se1) _cimg_mp_return(0); arg1 = 0; for (s = ss5; s::vector((uptrT)mp_norm0,pos).move_to(_opcode); break; @@ -16272,6 +16556,7 @@ namespace cimg_library_suffixed { case 'p' : if (!std::strncmp(ss,"print(",6)) { // Print expression + _cimg_mp_op("Function 'print()'"); pos = compile(ss6,se1,depth1,p_ref); *se1 = 0; if (_cimg_mp_is_vector(pos)) // Vector @@ -16287,14 +16572,11 @@ namespace cimg_library_suffixed { case 'r' : if (!std::strncmp(ss,"rol(",4) || !std::strncmp(ss,"ror(",4)) { // Bitwise rotation + _cimg_mp_op(ss[2]=='l'?"Function 'rol()'":"Function 'ror()'"); s1 = ss4; while (s11) { + arg2 = arg1 + 1; + if (p2>2) arg3 = arg2 + 1; + } + arg4 = compile(++s1,se1,depth1,0); + } else { + s2 = s1 + 1; while (s2::vector((uptrT)mp_rot3d,pos,arg1,arg2,arg3,arg4).move_to(code); + } else { // 2d rotation + _cimg_mp_check_type(arg1,1,1,0); + pos = vector(4); + CImg::vector((uptrT)mp_rot2d,pos,arg1).move_to(code); + } + _cimg_mp_return(pos); + } + if (!std::strncmp(ss,"round(",6)) { // Value rounding - s_op = "Function 'round()'"; + _cimg_mp_op("Function 'round()'"); s1 = ss6; while (s1::vector((uptrT)mp_single,arg1,code._width - p1).move_to(code,p1); + _cimg_mp_return(arg1); + } + if (!std::strncmp(ss,"sinh(",5)) { // Hyperbolic sine + _cimg_mp_op("Function 'sinh()'"); arg1 = compile(ss5,se1,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_sinh,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(std::sinh(mem[arg1])); @@ -16352,24 +16683,76 @@ namespace cimg_library_suffixed { } if (!std::strncmp(ss,"size(",5)) { // Vector size. + _cimg_mp_op("Function 'size()'"); arg1 = compile(ss5,se1,depth1,0); _cimg_mp_constant(_cimg_mp_is_scalar(arg1)?0:_cimg_mp_vector_size(arg1)); } + if (!std::strncmp(ss,"solve(",6)) { // Solve linear system + _cimg_mp_op("Function 'solve()'"); + s1 = ss6; while (s1::%s: %s: Types of first and second arguments ('%s' and '%s') " + "do not match for third argument 'nb_colsB=%u', " + "in expression '%s%s%s'.", + pixel_type(),_cimg_mp_calling_function,s_op, + s_type(arg1)._data,s_type(arg2)._data,p3, + (ss - 4)>expr._data?"...":"", + (ss - 4)>expr._data?ss - 4:expr._data, + se<&expr.back()?"...":""); + } + pos = vector(arg4*p3); + CImg::vector((uptrT)mp_solve,pos,arg1,arg2,arg4,arg5,p3).move_to(code); + _cimg_mp_return(pos); + } + if (!std::strncmp(ss,"sort(",5)) { // Sort vector - s_op = "Function 'sort()'"; + _cimg_mp_op("Function 'sort()'"); s1 = ss6; while (s1::%s: %s: Invalid specified chunk size (%u) for first argument " + "('%s'), in expression '%s%s%s'.", + pixel_type(),_cimg_mp_calling_function,s_op, + arg3,s_type(arg1)._data, + (ss - 4)>expr._data?"...":"", + (ss - 4)>expr._data?ss - 4:expr._data, + se<&expr.back()?"...":""); + } pos = vector(p1); - CImg::vector((uptrT)mp_vector_sort,pos,arg1,p1,arg2).move_to(code); + CImg::vector((uptrT)mp_sort,pos,arg1,p1,arg2,arg3).move_to(code); _cimg_mp_return(pos); } if (!std::strncmp(ss,"sqr(",4)) { // Square + _cimg_mp_op("Function 'sqr()'"); arg1 = compile(ss4,se1,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_sqr,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(cimg::sqr(mem[arg1])); @@ -16377,6 +16760,7 @@ namespace cimg_library_suffixed { } if (!std::strncmp(ss,"sqrt(",5)) { // Square root + _cimg_mp_op("Function 'sqrt()'"); arg1 = compile(ss5,se1,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_sqrt,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(std::sqrt(mem[arg1])); @@ -16386,6 +16770,7 @@ namespace cimg_library_suffixed { case 't' : if (!std::strncmp(ss,"tan(",4)) { // Tangent + _cimg_mp_op("Function 'tan()'"); arg1 = compile(ss4,se1,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_tan,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(std::tan(mem[arg1])); @@ -16393,20 +16778,57 @@ namespace cimg_library_suffixed { } if (!std::strncmp(ss,"tanh(",5)) { // Hyperbolic tangent + _cimg_mp_op("Function 'tanh()'"); arg1 = compile(ss5,se1,depth1,0); if (_cimg_mp_is_vector(arg1)) _cimg_mp_vector1_v(mp_tanh,arg1); if (_cimg_mp_is_constant(arg1)) _cimg_mp_constant(std::tanh(mem[arg1])); _cimg_mp_scalar1(mp_tanh,arg1); } + + if (!std::strncmp(ss,"trace(",6)) { // Matrix trace + _cimg_mp_op("Function 'trace()'"); + arg1 = compile(ss6,se1,depth1,0); + _cimg_mp_check_matrix_square(arg1,1); + p1 = (unsigned int)std::sqrt((float)_cimg_mp_vector_size(arg1)); + _cimg_mp_scalar2(mp_trace,arg1,p1); + } + + if (!std::strncmp(ss,"transp(",7)) { // Matrix transpose + _cimg_mp_op("Function 'transp()'"); + s1 = ss7; while (s1::%s: %s: Size of first argument ('%s') does not match" + "for second specified argument 'nb_cols=%u', " + "in expression '%s%s%s'.", + pixel_type(),_cimg_mp_calling_function,s_op, + s_type(arg1)._data,p2, + (ss - 4)>expr._data?"...":"", + (ss - 4)>expr._data?ss - 4:expr._data, + se<&expr.back()?"...":""); + } + pos = vector(p3*p2); + CImg::vector((uptrT)mp_transp,pos,arg1,p2,p3).move_to(code); + _cimg_mp_return(pos); + } break; case 'u' : if (*ss1=='(') { // Random value with uniform distribution + _cimg_mp_op("Function 'u()'"); if (*ss2==')') _cimg_mp_scalar2(mp_u,0,1); s1 = ss2; while (s10) || !std::strncmp(ss,"vector(",7)) { // Vector + _cimg_mp_op("Function 'vector()'"); arg2 = 0; // Number of specified values. s = std::strchr(ss6,'(') + 1; - if (*s!=')' || arg1==~0U) for (; s::vector((uptrT)mp_vector_init,pos,arg1),0); (_opcode>'y').move_to(code); @@ -16441,13 +16864,14 @@ namespace cimg_library_suffixed { case 'w' : if (!std::strncmp(ss,"whiledo",7) && (*ss7=='(' || (*ss7 && *ss7<=' ' && *ss8=='('))) { // While...do + _cimg_mp_op("Function 'whiledo()'"); if (*ss7<=' ') cimg::swap(*ss7,*ss8); // Allow space before opening brace s1 = ss8; while (s1::vector((uptrT)mp_whiledo,pos,arg1,p2 - p1,code._width - p2,arg2).move_to(code,p1); _cimg_mp_return(pos); @@ -16457,15 +16881,31 @@ namespace cimg_library_suffixed { if (!std::strncmp(ss,"min(",4) || !std::strncmp(ss,"max(",4) || !std::strncmp(ss,"med(",4) || !std::strncmp(ss,"kth(",4) || - !std::strncmp(ss,"arg(",4) || + !std::strncmp(ss,"arg(",4) || !std::strncmp(ss,"sum(",4) || + !std::strncmp(ss,"std(",4) || !std::strncmp(ss,"var(",4) || + !std::strncmp(ss,"prod(",5) || !std::strncmp(ss,"mean(",5) || !std::strncmp(ss,"argmin(",7) || !std::strncmp(ss,"argmax(",7)) { // Multi-argument functions + _cimg_mp_op(*ss=='a'?(ss[3]=='('?"Function 'arg()'":ss[4]=='i'?"Function 'argmin()'": + "Function 'argmax()'"): + *ss=='s'?(ss[1]=='u'?"Function 'sum()'":"Function 'std()'"): + *ss=='k'?"Function 'kth()'": + *ss=='p'?"Function 'prod()'": + *ss=='v'?"Function 'var()'": + ss[1]=='i'?"Function 'min()'": + ss[1]=='a'?"Function 'max()'": + ss[2]=='a'?"Function 'mean()'":"Function 'med()'"); pos = scalar(); - is_sth = *ss=='a' && ss[3]!='('; CImg::vector((uptrT)(*ss=='a'?(ss[3]=='('?mp_arg:ss[4]=='i'?mp_argmin:mp_argmax): - *ss=='k'?mp_kth:ss[1]=='i'?mp_min: - ss[1]=='a'?mp_max:mp_med),pos). + *ss=='s'?(ss[1]=='u'?mp_sum:mp_std): + *ss=='k'?mp_kth: + *ss=='p'?mp_prod: + *ss=='v'?mp_var: + ss[1]=='i'?mp_min: + ss[1]=='a'?mp_max: + ss[2]=='a'?mp_mean: + mp_med),pos). move_to(_opcode); - for (s = is_sth?ss7:ss4; s Look for a user-defined function. + // No corresponding built-in function -> Look for a user-defined macro. s0 = strchr(ss,'('); if (s0) { variable_name.assign(ss,s0 - ss + 1).back() = 0; @@ -16512,7 +16952,7 @@ namespace cimg_library_suffixed { if (p1!=p2+1) { // Number of specified argument do not fit *se = saved_char; cimg::strellipsize(variable_name,64); cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " - "CImg<%s>::%s: function '%s()': Number of specified arguments does not " + "CImg<%s>::%s: Function '%s()': Number of specified arguments does not " "match function declaration (%u argument%s required), " "in expression '%s%s%s'.", pixel_type(),_cimg_mp_calling_function,variable_name._data, @@ -16550,6 +16990,7 @@ namespace cimg_library_suffixed { // Vector specification using initializer '[ ... ]'. if (*ss=='[' && *se1==']') { + _cimg_mp_op("Operator '[]'"); arg1 = 0; // Number of specified values. if (*ss1!=']') for (s = ss1; s::vector(arg2).move_to(_opcode); ++arg1; } s = ns; } - _cimg_mp_check_vector0(arg1,"operator '[]'"); + _cimg_mp_check_vector0(arg1); pos = vector(arg1); _opcode.insert(CImg::vector((uptrT)mp_vector_init,pos,arg1),0); (_opcode>'y').move_to(code); @@ -16702,7 +17143,7 @@ namespace cimg_library_suffixed { // Reached an unknown item -> error. is_sth = true; // is_valid_variable_name if (*variable_name>='0' && *variable_name<='9') is_sth = false; - else for (ns = variable_name._data + 1; *ns; ++ns) + else for (ns = variable_name._data; *ns; ++ns) if (!is_varchar(*ns)) { is_sth = false; break; } *se = saved_char; cimg::strellipsize(variable_name,64); cimg::strellipsize(expr,64); @@ -16883,7 +17324,7 @@ namespace cimg_library_suffixed { while (siz-->0) *(ptr++) = -1; } - unsigned int vector(const unsigned int siz) { // Insert new vector of specified size in memory. + unsigned int vector(const unsigned int siz) { // Insert new vector of specified size in memory if (mempos + siz>=mem._width) { mem.resize(2*mem._width + siz,1,1,1,0); memtype.resize(mem._width,1,1,1,0); @@ -16895,7 +17336,14 @@ namespace cimg_library_suffixed { return pos; } - unsigned int vector_copy(const unsigned int arg) { // Insert new copy of specified vector in memory. + unsigned int vector(const unsigned int siz, const double value) { // Insert new initialized vector + const unsigned int pos = vector(siz); + double *ptr = &mem[pos] + 1; + for (unsigned int i = 0; i24) CImg::vector((uptrT)mp_self_map_vector_s,pos,siz,(uptrT)op,arg1).move_to(code); + else { + code.insert(siz); + for (unsigned int k = 1; k<=siz; ++k) + CImg::vector((uptrT)op,pos + k,arg1).move_to(code[code._width - 1 - siz + k]); + } + } + + void self_vector_v(const unsigned int pos, const mp_func op, const unsigned int arg1) { + const unsigned int siz = _cimg_mp_vector_size(pos); + if (siz>24) CImg::vector((uptrT)mp_self_map_vector_v,pos,siz,(uptrT)op,arg1).move_to(code); + else { + code.insert(siz); + for (unsigned int k = 1; k<=siz; ++k) + CImg::vector((uptrT)op,pos + k,arg1 + k).move_to(code[code._width - 1 - siz + k]); + } + } + unsigned int vector1_v(const mp_func op, const unsigned int arg1) { const unsigned int siz = _cimg_mp_vector_size(arg1), pos = is_tmp_vector(arg1)?arg1:vector(siz); - CImg::vector((uptrT)mp_vector_map_v,pos,siz,(uptrT)op,arg1).move_to(code); + if (siz>24) CImg::vector((uptrT)mp_vector_map_v,pos,siz,(uptrT)op,arg1).move_to(code); + else { + code.insert(siz); + for (unsigned int k = 1; k<=siz; ++k) + CImg::vector((uptrT)op,pos + k,arg1 + k).move_to(code[code._width - 1 - siz + k]); + } return pos; } @@ -16915,7 +17388,12 @@ namespace cimg_library_suffixed { const unsigned int siz = _cimg_mp_vector_size(arg1), pos = is_tmp_vector(arg1)?arg1:is_tmp_vector(arg2)?arg2:vector(siz); - CImg::vector((uptrT)mp_vector_map_vv,pos,siz,(uptrT)op,arg1,arg2).move_to(code); + if (siz>24) CImg::vector((uptrT)mp_vector_map_vv,pos,siz,(uptrT)op,arg1,arg2).move_to(code); + else { + code.insert(siz); + for (unsigned int k = 1; k<=siz; ++k) + CImg::vector((uptrT)op,pos + k,arg1 + k,arg2 + k).move_to(code[code._width - 1 - siz + k]); + } return pos; } @@ -16923,7 +17401,12 @@ namespace cimg_library_suffixed { const unsigned int siz = _cimg_mp_vector_size(arg1), pos = is_tmp_vector(arg1)?arg1:vector(siz); - CImg::vector((uptrT)mp_vector_map_vs,pos,siz,(uptrT)op,arg1,arg2).move_to(code); + if (siz>24) CImg::vector((uptrT)mp_vector_map_vs,pos,siz,(uptrT)op,arg1,arg2).move_to(code); + else { + code.insert(siz); + for (unsigned int k = 1; k<=siz; ++k) + CImg::vector((uptrT)op,pos + k,arg1 + k,arg2).move_to(code[code._width - 1 - siz + k]); + } return pos; } @@ -16931,7 +17414,12 @@ namespace cimg_library_suffixed { const unsigned int siz = _cimg_mp_vector_size(arg2), pos = is_tmp_vector(arg2)?arg2:vector(siz); - CImg::vector((uptrT)mp_vector_map_sv,pos,siz,(uptrT)op,arg1,arg2).move_to(code); + if (siz>24) CImg::vector((uptrT)mp_vector_map_sv,pos,siz,(uptrT)op,arg1,arg2).move_to(code); + else { + code.insert(siz); + for (unsigned int k = 1; k<=siz; ++k) + CImg::vector((uptrT)op,pos + k,arg1,arg2 + k).move_to(code[code._width - 1 - siz + k]); + } return pos; } @@ -16940,22 +17428,29 @@ namespace cimg_library_suffixed { const unsigned int siz = _cimg_mp_vector_size(arg1), pos = is_tmp_vector(arg1)?arg1:vector(siz); - CImg::vector((uptrT)mp_vector_map_vss,pos,siz,(uptrT)op,arg1,arg2,arg3).move_to(code); + if (siz>24) CImg::vector((uptrT)mp_vector_map_vss,pos,siz,(uptrT)op,arg1,arg2,arg3).move_to(code); + else { + code.insert(siz); + for (unsigned int k = 1; k<=siz; ++k) + CImg::vector((uptrT)op,pos + k,arg1 + k,arg2,arg3).move_to(code[code._width - 1 - siz + k]); + } return pos; } // Check if a memory slot is a positive integer constant scalar value. - void check_constant(const unsigned int arg, const unsigned int n_arg, const char *const s_op, + void check_constant(const unsigned int arg, const unsigned int n_arg, const bool is_strictly_positive, const char *const ss, char *const se, const char saved_char) { - _cimg_mp_check_type(arg,n_arg,s_op,1,0); + _cimg_mp_check_type(arg,n_arg,1,0); if (!_cimg_mp_is_constant(arg) || mem[arg]<(is_strictly_positive?1:0) || (double)(int)mem[arg]!=mem[arg]) { - const char *s_arg = !n_arg?"":n_arg==1?"First ":n_arg==2?"Second ":n_arg==3?"Third ":"One "; + const char *s_arg = !n_arg?"":n_arg==1?"First ":n_arg==2?"Second ":n_arg==3?"Third ": + n_arg==4?"Fourth ":n_arg==5?"Fifth ":n_arg==6?"Sixth ":n_arg==7?"Seventh ":n_arg==8?"Eighth ": + n_arg==9?"Ninth ":"One of the "; *se = saved_char; cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " - "CImg<%s>::%s(): %s: %s%s (of type '%s') is not a %spositive integer constant, " + "CImg<%s>::%s(): %s%s %s%s (of type '%s') is not a %spositive integer constant, " "in expression '%s%s%s'.", - pixel_type(),_cimg_mp_calling_function,s_op, + pixel_type(),_cimg_mp_calling_function,s_op,*s_op?":":"", s_arg,*s_arg?"argument":"Argument",s_type(arg)._data, is_strictly_positive?"strictly ":"", (ss - 4)>expr._data?"...":"", @@ -16965,9 +17460,9 @@ namespace cimg_library_suffixed { } // Check a matrix is square. - void check_matrix_square(const unsigned int arg, const unsigned int n_arg, const char *const s_op, + void check_matrix_square(const unsigned int arg, const unsigned int n_arg, const char *const ss, char *const se, const char saved_char) { - _cimg_mp_check_type(arg,n_arg,s_op,2,0); + _cimg_mp_check_type(arg,n_arg,2,0); const unsigned int siz = _cimg_mp_vector_size(arg), n = (unsigned int)std::sqrt((float)siz); @@ -16977,9 +17472,9 @@ namespace cimg_library_suffixed { else s_arg = !n_arg?"":n_arg==1?"First ":n_arg==2?"Second ":n_arg==3?"Third ":"One "; *se = saved_char; cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " - "CImg<%s>::%s(): %s: %s%s (of type '%s') " + "CImg<%s>::%s(): %s%s %s%s (of type '%s') " "cannot be considered as a square matrix, in expression '%s%s%s'.", - pixel_type(),_cimg_mp_calling_function,s_op, + pixel_type(),_cimg_mp_calling_function,s_op,*s_op?":":"", s_arg,*s_op=='F'?(*s_arg?"argument":"Argument"):(*s_arg?"operand":"Operand"), s_type(arg)._data, (ss - 4)>expr._data?"...":"", @@ -16992,7 +17487,7 @@ namespace cimg_library_suffixed { // Bits of 'mode' tells what types are allowed: // { 1 = scalar | 2 = vectorN }. // If 'N' is not zero, it also restricts the vectors to be of size N only. - void check_type(const unsigned int arg, const unsigned int n_arg, const char *const s_op, + void check_type(const unsigned int arg, const unsigned int n_arg, const unsigned int mode, const unsigned int N, const char *const ss, char *const se, const char saved_char) { const bool @@ -17004,7 +17499,9 @@ namespace cimg_library_suffixed { if (!cond) { const char *s_arg; if (*s_op!='F') s_arg = !n_arg?"":n_arg==1?"Left-hand ":"Right-hand "; - else s_arg = !n_arg?"":n_arg==1?"First ":n_arg==2?"Second ":n_arg==3?"Third ":"One "; + else s_arg = !n_arg?"":n_arg==1?"First ":n_arg==2?"Second ":n_arg==3?"Third ": + n_arg==4?"Fourth ":n_arg==5?"Fifth ":n_arg==6?"Sixth ":n_arg==7?"Seventh ":n_arg==8?"Eighth": + n_arg==9?"Ninth":"One of the "; CImg sb_type(32); if (mode==1) cimg_snprintf(sb_type,sb_type._width,"'scalar'"); else if (mode==2) { @@ -17016,9 +17513,9 @@ namespace cimg_library_suffixed { } *se = saved_char; cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " - "CImg<%s>::%s(): %s: %s%s has invalid type '%s' (should be %s), " + "CImg<%s>::%s(): %s%s %s%s has invalid type '%s' (should be %s), " "in expression '%s%s%s'.", - pixel_type(),_cimg_mp_calling_function,s_op, + pixel_type(),_cimg_mp_calling_function,s_op,*s_op?":":"", s_arg,*s_op=='F'?(*s_arg?"argument":"Argument"):(*s_arg?"operand":"Operand"), s_type(arg)._data,sb_type._data, (ss - 4)>expr._data?"...":"", @@ -17027,24 +17524,39 @@ namespace cimg_library_suffixed { } } + // Check is listin is not empty. + void check_list(const bool is_out, + const char *const ss, char *const se, const char saved_char) { + if ((!is_out && !listin) || (is_out && !listout)) { + *se = saved_char; cimg::strellipsize(expr,64); + throw CImgArgumentException("[_cimg_math_parser] " + "CImg<%s>::%s(): %s%s Invalid call with an empty image list, " + "in expression '%s%s%s'.", + pixel_type(),_cimg_mp_calling_function,s_op,*s_op?":":"", + (ss - 4)>expr._data?"...":"", + (ss - 4)>expr._data?ss - 4:expr._data, + se<&expr.back()?"...":""); + } + } + // Check a vector is not 0-dimensional, or with unknown dimension at compile time. - void check_vector0(const unsigned int dim, const char *const s_op, + void check_vector0(const unsigned int dim, const char *const ss, char *const se, const char saved_char) { if (!dim) { *se = saved_char; cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " - "CImg<%s>::%s(): %s: Invalid construction of a 0-dimensional vector, " + "CImg<%s>::%s(): %s%s Invalid construction of a 0-dimensional vector, " "in expression '%s%s%s'.", - pixel_type(),_cimg_mp_calling_function,s_op, + pixel_type(),_cimg_mp_calling_function,s_op,*s_op?":":"", (ss - 4)>expr._data?"...":"", (ss - 4)>expr._data?ss - 4:expr._data, se<&expr.back()?"...":""); } else if (dim==~0U) { *se = saved_char; cimg::strellipsize(expr,64); throw CImgArgumentException("[_cimg_math_parser] " - "CImg<%s>::%s(): %s: Invalid construction of a vector with dynamic size, " + "CImg<%s>::%s(): %s%s Invalid construction of a vector with dynamic size, " "in expression '%s%s%s'.", - pixel_type(),_cimg_mp_calling_function,s_op, + pixel_type(),_cimg_mp_calling_function,s_op,*s_op?":":"", (ss - 4)>expr._data?"...":"", (ss - 4)>expr._data?ss - 4:expr._data, se<&expr.back()?"...":""); @@ -17253,6 +17765,26 @@ namespace cimg_library_suffixed { return std::cosh(_mp_arg(2)); } + static double mp_crop(_cimg_math_parser& mp) { + double *ptrd = &_mp_arg(1) + 1; + const int x = (int)_mp_arg(3), y = (int)_mp_arg(4), z = (int)_mp_arg(5), c = (int)_mp_arg(6); + const unsigned int + dx = (unsigned int)mp.opcode[7], + dy = (unsigned int)mp.opcode[8], + dz = (unsigned int)mp.opcode[9], + dc = (unsigned int)mp.opcode[10]; + const bool boundary_conditions = (bool)_mp_arg(11); + unsigned int ind = (unsigned int)mp.opcode[2]; + if (ind!=~0U) ind = (unsigned int)cimg::mod((int)_mp_arg(2),mp.listin.width()); + const CImg &img = ind==~0U?mp.imgin:mp.listin[ind]; + if (!img) std::memset(ptrd,0,dx*dy*dz*dc*sizeof(double)); + else CImg(ptrd,dx,dy,dz,dc,true) = img.get_crop(x,y,z,c, + x + dx - 1,y + dy - 1, + z + dz - 1,c + dc - 1, + boundary_conditions); + return cimg::type::nan(); + } + static double mp_cross(_cimg_math_parser& mp) { CImg vout(&_mp_arg(1) + 1,1,3,1,1,true), @@ -17268,52 +17800,56 @@ namespace cimg_library_suffixed { } static double mp_debug(_cimg_math_parser& mp) { -#ifdef cimg_use_openmp - const unsigned int n_thread = omp_get_thread_num(); -#else - const unsigned int n_thread = 0; -#endif CImg expr(mp.opcode._height - 3); const uptrT *ptrs = mp.opcode._data + 3; cimg_for(expr,ptrd,char) *ptrd = (char)*(ptrs++); cimg::strellipsize(expr); const uptrT g_target = mp.opcode[1]; - std::fprintf(cimg::output(), - "\n[_cimg_math_parser] %p[thread #%u]:%*c" - "Start debugging expression '%s', code length %u -> mem[%u] (memsize: %u)", - (void*)&mp,n_thread,mp.debug_indent,' ', - expr._data,(unsigned int)mp.opcode[2],(unsigned int)g_target,mp.mem._width); - std::fflush(cimg::output()); - const CImg *const p_end = (++mp.p_code) + mp.opcode[2]; - CImg _op; - mp.debug_indent+=3; - for ( ; mp.p_code &op = *mp.p_code; - mp.opcode._data = op._data; mp.opcode._height = op._height; - _op.assign(1,op._height - 1); - const uptrT *ptrs = op._data + 1; - for (uptrT *ptrd = _op._data, *const ptrde = _op._data + _op._height; ptrd mem[%u] (memsize: %u)", + (void*)&mp,n_thread,mp.debug_indent,' ', + expr._data,(unsigned int)mp.opcode[2],(unsigned int)g_target,mp.mem._width); + std::fflush(cimg::output()); + const CImg *const p_end = (++mp.p_code) + mp.opcode[2]; + CImg _op; + mp.debug_indent+=3; + for ( ; mp.p_code &op = *mp.p_code; + mp.opcode._data = op._data; mp.opcode._height = op._height; + + _op.assign(1,op._height - 1); + const uptrT *ptrs = op._data + 1; + for (uptrT *ptrd = _op._data, *const ptrde = _op._data + _op._height; ptrd mem[%u] = %g", + (void*)&mp,n_thread,mp.debug_indent,' ', + (void*)mp.opcode._data,(void*)*mp.opcode,_op.value_string().data(), + (unsigned int)target,mp.mem[target]); + std::fflush(cimg::output()); + } + mp.debug_indent-=3; std::fprintf(cimg::output(), "\n[_cimg_math_parser] %p[thread #%u]:%*c" - "Opcode %p = [ %p,%s ] -> mem[%u] = %g", + "End debugging expression '%s' -> mem[%u] = %g (memsize: %u)", (void*)&mp,n_thread,mp.debug_indent,' ', - (void*)mp.opcode._data,(void*)*mp.opcode,_op.value_string().data(), - (unsigned int)target,mp.mem[target]); + expr._data,(unsigned int)g_target,mp.mem[g_target],mp.mem._width); std::fflush(cimg::output()); + --mp.p_code; } - mp.debug_indent-=3; - std::fprintf(cimg::output(), - "\n[_cimg_math_parser] %p[thread #%u]:%*c" - "End debugging expression '%s' -> mem[%u] = %g (memsize: %u)", - (void*)&mp,n_thread,mp.debug_indent,' ', - expr._data,(unsigned int)g_target,mp.mem[g_target],mp.mem._width); - std::fflush(cimg::output()); - --mp.p_code; return mp.mem[g_target]; } @@ -17321,6 +17857,24 @@ namespace cimg_library_suffixed { return _mp_arg(2) - 1; } + static double mp_det(_cimg_math_parser& mp) { + const double *ptrs = &_mp_arg(2) + 1; + const unsigned int k = (unsigned int)mp.opcode(3); + return CImg(ptrs,k,k,1,1,true).det(); + } + + static double mp_diag(_cimg_math_parser& mp) { + double *ptrd = &_mp_arg(1) + 1; + const double *ptrs = &_mp_arg(2) + 1; + const unsigned int k = (unsigned int)mp.opcode(3); + CImg(ptrd,k,k,1,1,true) = CImg(ptrs,1,k,1,1,true).get_diagonal(); + return cimg::type::nan(); + } + + static double mp_div(_cimg_math_parser& mp) { + return _mp_arg(2)/_mp_arg(3); + } + static double mp_dot(_cimg_math_parser& mp) { const unsigned int siz = (unsigned int)mp.opcode[4]; return CImg(&_mp_arg(2) + 1,1,siz,1,1,true). @@ -17346,8 +17900,36 @@ namespace cimg_library_suffixed { return mp.mem[mem_proc]; } - static double mp_div(_cimg_math_parser& mp) { - return _mp_arg(2)/_mp_arg(3); + static double mp_draw(_cimg_math_parser& mp) { + const int x = (int)_mp_arg(3), y = (int)_mp_arg(4), z = (int)_mp_arg(5), c = (int)_mp_arg(6); + const unsigned int + dx = (unsigned int)mp.opcode[7], + dy = (unsigned int)mp.opcode[8], + dz = (unsigned int)mp.opcode[9], + dc = (unsigned int)mp.opcode[10]; + const CImg S(&_mp_arg(1) + 1,dx,dy,dz,dc,true); + const float opacity = (float)_mp_arg(11); + unsigned int ind = (unsigned int)mp.opcode[2]; + if (ind!=~0U) ind = (unsigned int)cimg::mod((int)_mp_arg(2),mp.listin.width()); + CImg &img = ind==~0U?mp.imgout:mp.listout[ind]; + if (img) { + if (mp.opcode[12]!=(uptrT)-1) { + const CImg M(&_mp_arg(12) + 1,dx,dy,dz,(unsigned int)mp.opcode[13],true); + img.draw_image(x,y,z,c,S,M,opacity,(float)_mp_arg(14)); + } else img.draw_image(x,y,z,c,S,opacity); + } + return cimg::type::nan(); + } + + static double mp_eig(_cimg_math_parser& mp) { + double *ptrd = &_mp_arg(1) + 1; + const double *ptr1 = &_mp_arg(2) + 1; + const unsigned int k = (unsigned int)mp.opcode(3); + CImg val, vec; + CImg(ptr1,k,k,1,1,true).symmetric_eigen(val,vec); + CImg(ptrd,k,1,1,1,true) = val.unroll('x'); + CImg(ptrd + k,k,k,1,1,true) = vec.get_transpose(); + return cimg::type::nan(); } static double mp_eq(_cimg_math_parser& mp) { @@ -17358,6 +17940,13 @@ namespace cimg_library_suffixed { return std::exp(_mp_arg(2)); } + static double mp_eye(_cimg_math_parser& mp) { + double *ptrd = &_mp_arg(1) + 1; + const unsigned int k = (unsigned int)mp.opcode(2); + CImg(ptrd,k,k,1,1,true).identity_matrix(); + return cimg::type::nan(); + } + static double mp_g(_cimg_math_parser& mp) { cimg::unused(mp); return cimg::grand(); @@ -17424,6 +18013,14 @@ namespace cimg_library_suffixed { return (double)(long)_mp_arg(2); } + static double mp_inv(_cimg_math_parser& mp) { + double *ptrd = &_mp_arg(1) + 1; + const double *ptr1 = &_mp_arg(2) + 1; + const unsigned int k = (unsigned int)mp.opcode(3); + CImg(ptrd,k,k,1,1,true) = CImg(ptr1,k,k,1,1,true).get_invert(); + return cimg::type::nan(); + } + static double mp_ioff(_cimg_math_parser& mp) { const unsigned int boundary_conditions = (unsigned int)_mp_arg(3); @@ -17713,9 +18310,8 @@ namespace cimg_library_suffixed { z = (int)_mp_arg(5), c = (int)_mp_arg(6); const double val = _mp_arg(1); if (x>=0 && x=0 && y=0 && z=0 && c=0 && z=0 && c=0 && x=0 && y=0 && z=0 && c=0 && z=0 && c(ptrs,k,k,1,1,true).det(); - } - - static double mp_matrix_diag(_cimg_math_parser& mp) { - double *ptrd = &_mp_arg(1) + 1; - const double *ptrs = &_mp_arg(2) + 1; - const unsigned int k = (unsigned int)mp.opcode(3); - CImg(ptrd,k,k,1,1,true) = CImg(ptrs,1,k,1,1,true).get_diagonal(); - return cimg::type::nan(); - } - - static double mp_matrix_eig(_cimg_math_parser& mp) { - double *ptrd = &_mp_arg(1) + 1; - const double *ptr1 = &_mp_arg(2) + 1; - const unsigned int k = (unsigned int)mp.opcode(3); - CImg val, vec; - CImg(ptr1,k,k,1,1,true).symmetric_eigen(val,vec); - CImg(ptrd,k,1,1,1,true) = val.unroll('x'); - CImg(ptrd + k,k,k,1,1,true) = vec.get_transpose(); - return cimg::type::nan(); - } - - static double mp_matrix_eye(_cimg_math_parser& mp) { - double *ptrd = &_mp_arg(1) + 1; - const unsigned int k = (unsigned int)mp.opcode(2); - CImg(ptrd,k,k,1,1,true).identity_matrix(); - return cimg::type::nan(); - } - - static double mp_matrix_inv(_cimg_math_parser& mp) { - double *ptrd = &_mp_arg(1) + 1; - const double *ptr1 = &_mp_arg(2) + 1; - const unsigned int k = (unsigned int)mp.opcode(3); - CImg(ptrd,k,k,1,1,true) = CImg(ptr1,k,k,1,1,true).get_invert(); - return cimg::type::nan(); - } - static double mp_matrix_mul(_cimg_math_parser& mp) { double *ptrd = &_mp_arg(1) + 1; const double @@ -18152,48 +18707,110 @@ namespace cimg_library_suffixed { return cimg::type::nan(); } - static double mp_matrix_rot(_cimg_math_parser& mp) { - double *ptrd = &_mp_arg(1) + 1; - const float x = (float)_mp_arg(2), y = (float)_mp_arg(3), z = (float)_mp_arg(4), theta = (float)_mp_arg(5); - CImg(ptrd,3,3,1,1,true) = CImg::rotation_matrix(x,y,z,theta); - return cimg::type::nan(); - } - - static double mp_matrix_solve(_cimg_math_parser& mp) { - double *ptrd = &_mp_arg(1) + 1; - const double - *ptr1 = &_mp_arg(2) + 1, - *ptr2 = &_mp_arg(3) + 1; - const unsigned int - k = (unsigned int)mp.opcode(4), - l = (unsigned int)mp.opcode(5), - m = (unsigned int)mp.opcode(6); - CImg(ptrd,m,l,1,1,true) = CImg(ptr2,m,k,1,1,true).get_solve(CImg(ptr1,l,k,1,1,true)); - return cimg::type::nan(); - } - - static double mp_matrix_trace(_cimg_math_parser& mp) { - const double *ptrs = &_mp_arg(2) + 1; - const unsigned int k = (unsigned int)mp.opcode(3); - return CImg(ptrs,k,k,1,1,true).trace(); - } - - static double mp_matrix_trans(_cimg_math_parser& mp) { - double *ptrd = &_mp_arg(1) + 1; - const double *ptr1 = &_mp_arg(2) + 1; - const unsigned int - k = (unsigned int)mp.opcode(3), - l = (unsigned int)mp.opcode(4); - CImg(ptrd,l,k,1,1,true) = CImg(ptr1,k,l,1,1,true).get_transpose(); - return cimg::type::nan(); - } - static double mp_max(_cimg_math_parser& mp) { double val = _mp_arg(2); for (unsigned int i = 3; i=mp.mem.width()) + throw CImgArgumentException("[_cimg_math_parser] CImg<%s>: 'copy()': " + "Out-of-bounds variable pointer " + "(length: %ld, increment: %ld, offset start: %ld, " + "offset end: %ld, offset max: %u).", + mp.imgin.pixel_type(),siz,inc,off,eoff,mp.mem._width - 1); + return &mp.mem[off]; + } + + static float* _mp_memcopy_float(_cimg_math_parser& mp, const uptrT *const p_ref, + const long siz, const long inc) { + const unsigned ind = p_ref[1]; + const CImg &img = ind==~0U?mp.imgin:mp.listin[cimg::mod((int)mp.mem[ind],mp.listin.width())]; + const bool is_relative = (bool)p_ref[2]; + int ox, oy, oz, oc; + long off = 0; + if (is_relative) { + ox = (int)mp.mem[_cimg_mp_x]; + oy = (int)mp.mem[_cimg_mp_y]; + oz = (int)mp.mem[_cimg_mp_z]; + oc = (int)mp.mem[_cimg_mp_c]; + off = img.offset(ox,oy,oz,oc); + } + if ((*p_ref)%2) { + const int + x = (int)mp.mem[p_ref[3]], + y = (int)mp.mem[p_ref[4]], + z = (int)mp.mem[p_ref[5]], + c = *p_ref==5?0:(int)mp.mem[p_ref[6]]; + off+=(long)img.offset(x,y,z,c); + } else off+=(long)mp.mem[p_ref[3]]; + const long eoff = off + (siz - 1)*inc; + if (off<0 || eoff>=(long)img.size()) + throw CImgArgumentException("[_cimg_math_parser] CImg<%s>: Function 'copy()': " + "Out-of-bounds image pointer " + "(length: %ld, increment: %ld, offset start: %ld, " + "offset end: %ld, offset max: %lu).", + mp.imgin.pixel_type(),siz,inc,off,eoff,img.size() - 1); + return (float*)&img[off]; + } + + static double mp_memcopy(_cimg_math_parser& mp) { + long siz = (long)_mp_arg(4); + const long inc_d = (long)_mp_arg(5), inc_s = (long)_mp_arg(6); + if (siz>0) { + const bool + is_doubled = mp.opcode[7]<=1, + is_doubles = mp.opcode[14]<=1; + if (is_doubled && is_doubles) { // (double*) <- (double*) + double *ptrd = _mp_memcopy_double(mp,mp.opcode[2],&mp.opcode[7],siz,inc_d); + const double *ptrs = _mp_memcopy_double(mp,mp.opcode[3],&mp.opcode[14],siz,inc_s); + if (inc_d==1 && inc_s==1) { + if (ptrs + siz - 1ptrd + siz - 1) std::memcpy(ptrd,ptrs,siz*sizeof(double)); + else std::memmove(ptrd,ptrs,siz*sizeof(double)); + } else { + if (ptrs + (siz - 1)*inc_sptrd + (siz - 1)*inc_d) + while (siz-->0) { *ptrd = (double)*ptrs; ptrd+=inc_d; ptrs+=inc_s; } + else { // Overlapping buffers + CImg buf(siz); + cimg_for(buf,ptr,double) { *ptr = *ptrs; ptrs+=inc_s; } + ptrs = buf; + while (siz-->0) { *ptrd = *(ptrs++); ptrd+=inc_d; } + } + } + } else if (is_doubled && !is_doubles) { // (double*) <- (float*) + double *ptrd = _mp_memcopy_double(mp,mp.opcode[2],&mp.opcode[7],siz,inc_d); + const float *ptrs = _mp_memcopy_float(mp,&mp.opcode[14],siz,inc_s); + while (siz-->0) { *ptrd = (double)*ptrs; ptrd+=inc_d; ptrs+=inc_s; } + } else if (!is_doubled && is_doubles) { // (float*) <- (double*) + float *ptrd = _mp_memcopy_float(mp,&mp.opcode[7],siz,inc_d); + const double *ptrs = _mp_memcopy_double(mp,mp.opcode[3],&mp.opcode[14],siz,inc_s); + while (siz-->0) { *ptrd = (float)*ptrs; ptrd+=inc_d; ptrs+=inc_s; } + } else { // (float*) <- (float*) + float *ptrd = _mp_memcopy_float(mp,&mp.opcode[7],siz,inc_d); + const float *ptrs = _mp_memcopy_float(mp,&mp.opcode[14],siz,inc_s); + if (inc_d==1 && inc_s==1) { + if (ptrs + siz - 1ptrd + siz - 1) std::memcpy(ptrd,ptrs,siz*sizeof(float)); + else std::memmove(ptrd,ptrs,siz*sizeof(float)); + } else { + if (ptrs + (siz - 1)*inc_sptrd + (siz - 1)*inc_d) + while (siz-->0) { *ptrd = (float)*ptrs; ptrd+=inc_d; ptrs+=inc_s; } + else { // Overlapping buffers + CImg buf(siz); + cimg_for(buf,ptr,float) { *ptr = *ptrs; ptrs+=inc_s; } + ptrs = buf; + while (siz-->0) { *ptrd = *(ptrs++); ptrd+=inc_d; } + } + } + } + } + return _mp_arg(1); + } + static double mp_min(_cimg_math_parser& mp) { double val = _mp_arg(2); for (unsigned int i = 3; i vals(mp.opcode._height - 2); double *p = vals.data(); @@ -18278,13 +18901,26 @@ namespace cimg_library_suffixed { } static double mp_print(_cimg_math_parser& mp) { + cimg::mutex(6); CImg expr(mp.opcode._height - 2); const uptrT *ptrs = mp.opcode._data + 2; cimg_for(expr,ptrd,char) *ptrd = (char)*(ptrs++); cimg::strellipsize(expr); const double val = _mp_arg(1); - std::fprintf(cimg::output(),"\n[_cimg_math_parser] %s = %g",expr._data,val); - std::fflush(cimg::output()); +#ifdef cimg_use_openmp +#pragma omp critical +#endif + { + std::fprintf(cimg::output(),"\n[_cimg_math_parser] %s = %g",expr._data,val); + std::fflush(cimg::output()); + } + cimg::mutex(6,0); + return val; + } + + static double mp_prod(_cimg_math_parser& mp) { + double val = _mp_arg(2); + for (unsigned int i = 3; i::nan(); + } + + static double mp_rot3d(_cimg_math_parser& mp) { + double *ptrd = &_mp_arg(1) + 1; + const float x = (float)_mp_arg(2), y = (float)_mp_arg(3), z = (float)_mp_arg(4), theta = (float)_mp_arg(5); + CImg(ptrd,3,3,1,1,true) = CImg::rotation_matrix(x,y,z,theta); + return cimg::type::nan(); + } + static double mp_round(_cimg_math_parser& mp) { return cimg::round(_mp_arg(2),_mp_arg(3),(int)_mp_arg(4)); } @@ -18403,9 +19059,8 @@ namespace cimg_library_suffixed { z = (int)_mp_arg(4), c = (int)_mp_arg(5); const double val = _mp_arg(1); if (x>=0 && x=0 && y=0 && z=0 && c=0 && z=0 && c=0 && x=0 && y=0 && z=0 && c=0 && z=0 && c *const p_end = ++mp.p_code + mp.opcode[2]; + mp.p_code &op = *mp.p_code; + mp.opcode._data = op._data; mp.opcode._height = op._height; + const uptrT target = mp.opcode[1]; + mp.mem[target] = _cimg_mp_defunc(mp); + } + } + --mp.p_code; + return res; + } + static double mp_sinh(_cimg_math_parser& mp) { return std::sinh(_mp_arg(2)); } + static double mp_solve(_cimg_math_parser& mp) { + double *ptrd = &_mp_arg(1) + 1; + const double + *ptr1 = &_mp_arg(2) + 1, + *ptr2 = &_mp_arg(3) + 1; + const unsigned int + k = (unsigned int)mp.opcode(4), + l = (unsigned int)mp.opcode(5), + m = (unsigned int)mp.opcode(6); + CImg(ptrd,m,k,1,1,true) = CImg(ptr2,m,l,1,1,true).get_solve(CImg(ptr1,k,l,1,1,true)); + return cimg::type::nan(); + } + + static double mp_sort(_cimg_math_parser& mp) { + double *const ptrd = &_mp_arg(1) + 1; + const double *const ptrs = &_mp_arg(2) + 1; + const unsigned int + siz = mp.opcode[3], + chunk_siz = mp.opcode[5]; + const bool is_increasing = (bool)_mp_arg(4); + CImg(ptrd,chunk_siz,siz/chunk_siz,1,1,true) = CImg(ptrs,chunk_siz,siz/chunk_siz,1,1,true). + get_sort(is_increasing,chunk_siz>1?'y':0); + return cimg::type::nan(); + } + static double mp_sqr(_cimg_math_parser& mp) { return cimg::sqr(_mp_arg(2)); } @@ -18570,10 +19267,23 @@ namespace cimg_library_suffixed { return std::sqrt(_mp_arg(2)); } + static double mp_std(_cimg_math_parser& mp) { + CImg vals(mp.opcode._height - 2); + double *p = vals.data(); + for (unsigned int i = 2; i(ptrs,k,k,1,1,true).trace(); + } + + static double mp_transp(_cimg_math_parser& mp) { + double *ptrd = &_mp_arg(1) + 1; + const double *ptr1 = &_mp_arg(2) + 1; + const unsigned int + k = (unsigned int)mp.opcode(3), + l = (unsigned int)mp.opcode(4); + CImg(ptrd,l,k,1,1,true) = CImg(ptr1,k,l,1,1,true).get_transpose(); + return cimg::type::nan(); + } + static double mp_u(_cimg_math_parser& mp) { return cimg::rand(_mp_arg(2),_mp_arg(3)); } + static double mp_var(_cimg_math_parser& mp) { + CImg vals(mp.opcode._height - 2); + double *p = vals.data(); + for (unsigned int i = 2; i::nan(); @@ -18704,15 +19437,6 @@ namespace cimg_library_suffixed { return _mp_arg(5); } - static double mp_vector_sort(_cimg_math_parser& mp) { - double *const ptrd = &_mp_arg(1) + 1; - const double *const ptrs = &_mp_arg(2) + 1; - const unsigned int siz = mp.opcode[3]; - const bool is_increasing = (bool)_mp_arg(4); - CImg(ptrd,1,siz,1,1,true) = CImg(ptrs,1,siz,1,1,true).get_sort(is_increasing); - return cimg::type::nan(); - } - static double mp_vector_print(_cimg_math_parser& mp) { CImg expr(mp.opcode._height - 3); const uptrT *ptrs = mp.opcode._data + 3; @@ -20185,7 +20909,8 @@ namespace cimg_library_suffixed { case 's' : return (double)_spectrum; case 'r' : return (double)_is_shared; } - _cimg_math_parser mp(expression + (*expression=='>' || *expression=='<' || *expression=='*'?1:0),"eval", + _cimg_math_parser mp(expression + (*expression=='>' || *expression=='<' || + *expression=='*' || *expression==':'?1:0),"eval", *this,img_output,list_inputs,list_outputs); return mp(x,y,z,c); } @@ -20229,7 +20954,8 @@ namespace cimg_library_suffixed { case 's' : output.assign(1); *output = (t)_spectrum; case 'r' : output.assign(1); *output = (t)_is_shared; } - _cimg_math_parser mp(expression + (*expression=='>' || *expression=='<' || *expression=='*'?1:0),"eval", + _cimg_math_parser mp(expression + (*expression=='>' || *expression=='<' || + *expression=='*' || *expression==':'?1:0),"eval", *this,img_output,list_inputs,list_outputs); output.assign(1,cimg::max(1U,mp.result_dim)); mp(x,y,z,c,output._data); @@ -20823,13 +21549,13 @@ namespace cimg_library_suffixed { "incompatible dimensions.", cimg_instance, A._width,A._height,A._depth,A._spectrum,A._data); - if (_width!=1) { - cimg_forX(*this,i) draw_image(i,get_column(i).solve(A)); - return *this; - } - typedef _cimg_Ttfloat Ttfloat; - if (A._width==A._height) { + if (A._width==A._height) { // Classical linear system + if (_width!=1) { + CImg res(_width,A._width); + cimg_forX(*this,i) res.draw_image(i,get_column(i).solve(A)); + return res.move_to(*this); + } #ifdef cimg_use_lapack char TRANS = 'N'; int INFO, N = _height, LWORK = 4*N, *const IPIV = new int[N]; @@ -20865,6 +21591,11 @@ namespace cimg_library_suffixed { #endif } else { // Least-square solution for non-square systems. #ifdef cimg_use_lapack + if (_width!=1) { + CImg res(_width,A._width); + cimg_forX(*this,i) res.draw_image(i,get_column(i).solve(A)); + return res.move_to(*this); + } char TRANS = 'N'; int INFO, N = A._width, M = A._height, LWORK = -1, LDA = M, LDB = M, NRHS = _width; Ttfloat WORK_QUERY; @@ -22393,24 +23124,25 @@ namespace cimg_library_suffixed { CImg& _fill(const char *const expression, const bool repeat_values, const bool allow_formula, const CImgList *const list_inputs, CImgList *const list_outputs, - const char *const calling_function, const CImg *provide_base) { + const char *const calling_function, const CImg *provides_copy) { if (is_empty() || !expression || !*expression) return *this; const unsigned int omode = cimg::exception_mode(); cimg::exception_mode(0); CImg is_error; if (allow_formula) try { // Try to fill values according to a formula - bool is_parallelizable = true; - const CImg - _base = _cimg_math_parser::needs_input_copy(expression,is_parallelizable)? - (provide_base?*provide_base:+*this):CImg(), - &base = provide_base?*provide_base:_base?_base:*this; - _cimg_math_parser mp(expression + (*expression=='>' || *expression=='<' || *expression=='*'?1:0), + CImg base = provides_copy?provides_copy->get_shared():get_shared(); + _cimg_math_parser mp(expression + (*expression=='>' || *expression=='<' || + *expression=='*' || *expression==':'?1:0), calling_function,base,this,list_inputs,list_outputs); + if (!provides_copy && expression && *expression!='>' && *expression!='<' && *expression!=':' && + mp.need_input_copy) + base.assign().assign(*this); // Needs input copy + bool do_in_parallel = false; #ifdef cimg_use_openmp - cimg_openmp_if(*expression=='*' || - (is_parallelizable && _width>=320 && _height*_depth*_spectrum>=2 && + cimg_openmp_if(*expression=='*' || *expression==':' || + (mp.is_parallelizable && _width>=320 && _height*_depth*_spectrum>=2 && std::strlen(expression)>=6)) do_in_parallel = true; #endif @@ -26443,18 +27175,20 @@ namespace cimg_library_suffixed { \param boundary Boundary conditions. Can be { 0=dirichlet | 1=neumann | 2=periodic }. \note Most of the time, size of the image is modified. **/ - CImg& rotate(const float angle, const unsigned int interpolation=1, const unsigned int boundary=0) { + CImg& rotate(const float angle, const unsigned int interpolation=1, + const unsigned int boundary_conditions=0) { const float nangle = cimg::mod(angle,360.0f); if (nangle==0.0f) return *this; - return get_rotate(angle,interpolation,boundary).move_to(*this); + return get_rotate(angle,interpolation,boundary_conditions).move_to(*this); } //! Rotate image with arbitrary angle \newinstance. - CImg get_rotate(const float angle, const unsigned int interpolation=1, const unsigned int boundary=0) const { + CImg get_rotate(const float angle, const unsigned int interpolation=1, + const unsigned int boundary_conditions=0) const { if (is_empty()) return *this; CImg res; const float nangle = cimg::mod(angle,360.0f); - if (boundary!=1 && cimg::mod(nangle,90.0f)==0) { // Optimized version for orthogonal angles. + if (boundary_conditions!=1 && cimg::mod(nangle,90.0f)==0) { // Optimized version for orthogonal angles. const int wm1 = width() - 1, hm1 = height() - 1; const int iangle = (int)nangle/90; switch (iangle) { @@ -26487,7 +27221,7 @@ namespace cimg_library_suffixed { w2 = 0.5f*_width, h2 = 0.5f*_height, dw2 = 0.5f*(ux + vx), dh2 = 0.5f*(uy + vy); res.assign((int)(ux + vx),(int)(uy + vy),_depth,_spectrum); - switch (boundary) { + switch (boundary_conditions) { case 0 : { // Dirichlet boundaries. switch (interpolation) { case 2 : { // Cubic interpolation. @@ -26577,7 +27311,7 @@ namespace cimg_library_suffixed { "rotate(): Invalid specified border conditions %d " "(should be { 0=dirichlet | 1=neumann | 2=periodic }).", cimg_instance, - boundary); + boundary_conditions); } } return res; @@ -26593,13 +27327,13 @@ namespace cimg_library_suffixed { \param interpolation_type Type of interpolation. Can be { 0=nearest | 1=linear | 2=cubic }. **/ CImg& rotate(const float angle, const float cx, const float cy, const float zoom, - const unsigned int interpolation=1, const unsigned int boundary=3) { - return get_rotate(angle,cx,cy,zoom,interpolation,boundary).move_to(*this); + const unsigned int interpolation=1, const unsigned int boundary_conditions=0) { + return get_rotate(angle,cx,cy,zoom,interpolation,boundary_conditions).move_to(*this); } //! Rotate image with arbitrary angle, around a center point \newinstance. CImg get_rotate(const float angle, const float cx, const float cy, const float zoom, - const unsigned int interpolation=1, const unsigned int boundary=3) const { + const unsigned int interpolation=1, const unsigned int boundary_conditions=0) const { if (interpolation>2) throw CImgArgumentException(_cimg_instance "rotate(): Invalid specified interpolation type %d " @@ -26613,7 +27347,7 @@ namespace cimg_library_suffixed { rad = (float)((angle*cimg::PI)/180.0), ca = (float)std::cos(rad)/zoom, sa = (float)std::sin(rad)/zoom; - switch (boundary) { + switch (boundary_conditions) { case 0 : { switch (interpolation) { case 2 : { @@ -26703,7 +27437,7 @@ namespace cimg_library_suffixed { "rotate(): Invalid specified border conditions %d " "(should be { 0=dirichlet | 1=neumann | 2=periodic }).", cimg_instance, - boundary); + boundary_conditions); } return res; } @@ -29523,12 +30257,25 @@ namespace cimg_library_suffixed { //! Compute watershed transform. /** \param priority Priority map. - \param fill_lines Tells if watershed lines must be filled or not. + \param is_high_connectivity Boolean that choose between 4(false)- or 8(true)-connectivity + in 2d case, and between 6(false)- or 26(true)-connectivity in 3d case. \note Non-zero values of the instance instance are propagated to zero-valued ones according to - specified the priority map. + specified the priority map. **/ template - CImg& watershed(const CImg& priority, const bool fill_lines=true) { + CImg& watershed(const CImg& priority, const bool is_high_connectivity=false) { +#define _cimg_watershed_init(cond,X,Y,Z) \ + if (cond && !(*this)(X,Y,Z)) Q._priority_queue_insert(labels,sizeQ,priority(X,Y,Z),X,Y,Z,nb_seeds) + +#define _cimg_watershed_propagate(cond,X,Y,Z) \ + if (cond) { \ + if ((*this)(X,Y,Z)) { \ + ns = labels(X,Y,Z) - 1; xs = seeds(ns,0); ys = seeds(ns,1); zs = seeds(ns,2); \ + d = cimg::sqr((float)x - xs) + cimg::sqr((float)y - ys) + cimg::sqr((float)z - zs); \ + if (d is_queued(_width,_height,_depth,1,0); + CImg labels(_width,_height,_depth,1,0), seeds(64,3); CImg::type> Q; unsigned int sizeQ = 0; + int px, nx, py, ny, pz, nz; + bool is_px, is_nx, is_py, is_ny, is_pz, is_nz; + const bool is_3d = _depth>1; // Find seed points and insert them in priority queue. + unsigned int nb_seeds = 0; const T *ptrs = _data; - cimg_forXYZ(*this,x,y,z) if (*(ptrs++)) { - if (x - 1>=0 && !(*this)(x - 1,y,z)) - Q._priority_queue_insert(is_queued,sizeQ,priority(x - 1,y,z),x - 1,y,z); - if (x + 1=0 && !(*this)(x,y - 1,z)) - Q._priority_queue_insert(is_queued,sizeQ,priority(x,y - 1,z),x,y - 1,z); - if (y + 1=0 && !(*this)(x,y,z - 1)) - Q._priority_queue_insert(is_queued,sizeQ,priority(x,y,z - 1),x,y,z - 1); - if (z + 1=seeds._width) seeds.resize(2*seeds._width,3,1,1,0); + seeds(nb_seeds,0) = x; seeds(nb_seeds,1) = y; seeds(nb_seeds++,2) = z; + px = x - 1; nx = x + 1; + py = y - 1; ny = y + 1; + pz = z - 1; nz = z + 1; + is_px = px>=0; is_nx = nx=0; is_ny = ny=0; is_nz = nz=0; is_nx = nx=0; is_ny = ny=0; is_nz = nz=0) { - if ((*this)(x - 1,y,z)) { - if (!label) label = (*this)(x - 1,y,z); - else if (label!=(*this)(x - 1,y,z)) is_same_label = false; - } else Q._priority_queue_insert(is_queued,sizeQ,priority(x - 1,y,z),x - 1,y,z); - } - if (x + 1=0) { - if ((*this)(x,y - 1,z)) { - if (!label) label = (*this)(x,y - 1,z); - else if (label!=(*this)(x,y - 1,z)) is_same_label = false; - } else Q._priority_queue_insert(is_queued,sizeQ,priority(x,y - 1,z),x,y - 1,z); - } - if (y + 1=0) { - if ((*this)(x,y,z - 1)) { - if (!label) label = (*this)(x,y,z - 1); - else if (label!=(*this)(x,y,z - 1)) is_same_label = false; - } else Q._priority_queue_insert(is_queued,sizeQ,priority(x,y,z - 1),x,y,z - 1); - } - if (z + 1=0 && (*this)(x - 1,y,z)) || (x + 1=0 && (*this)(x,y - 1,z)) || (y + 1=0 && (*this)(x,y,z - 1)) || (z + 1>depth() && (*this)(x,y,z + 1)))) - Q._priority_queue_insert(is_queued,sizeQ,priority(x,y,z),x,y,z); + Q._priority_queue_remove(sizeQ); - // Start line filling process. - while (sizeQ) { - const int x = (int)Q(0,1), y = (int)Q(0,2), z = (int)Q(0,3); - Q._priority_queue_remove(sizeQ); - t pmax = cimg::type::min(); - int xmax = 0, ymax = 0, zmax = 0; - if (x - 1>=0) { - if ((*this)(x - 1,y,z)) { - if (priority(x - 1,y,z)>pmax) { pmax = priority(x - 1,y,z); xmax = x - 1; ymax = y; zmax = z; } - } else Q._priority_queue_insert(is_queued,sizeQ,priority(x - 1,y,z),x - 1,y,z); - } - if (x + 1pmax) { pmax = priority(x + 1,y,z); xmax = x + 1; ymax = y; zmax = z; } - } else Q._priority_queue_insert(is_queued,sizeQ,priority(x + 1,y,z),x + 1,y,z); - } - if (y - 1>=0) { - if ((*this)(x,y - 1,z)) { - if (priority(x,y - 1,z)>pmax) { pmax = priority(x,y - 1,z); xmax = x; ymax = y - 1; zmax = z; } - } else Q._priority_queue_insert(is_queued,sizeQ,priority(x,y - 1,z),x,y - 1,z); - } - if (y + 1pmax) { pmax = priority(x,y + 1,z); xmax = x; ymax = y + 1; zmax = z; } - } else Q._priority_queue_insert(is_queued,sizeQ,priority(x,y + 1,z),x,y + 1,z); - } - if (z - 1>=0) { - if ((*this)(x,y,z - 1)) { - if (priority(x,y,z - 1)>pmax) { pmax = priority(x,y,z - 1); xmax = x; ymax = y; zmax = z - 1; } - } else Q._priority_queue_insert(is_queued,sizeQ,priority(x,y,z - 1),x,y,z - 1); - } - if (z + 1pmax) { pmax = priority(x,y,z + 1); xmax = x; ymax = y; zmax = z + 1; } - } else Q._priority_queue_insert(is_queued,sizeQ,priority(x,y,z + 1),x,y,z + 1); - } - (*this)(x,y,z) = (*this)(xmax,ymax,zmax); + unsigned int xs, ys, zs, ns, nmin = 0; + float d, dmin = cimg::type::inf(); + T label = 0; + _cimg_watershed_propagate(is_px,px,y,z); + _cimg_watershed_propagate(is_nx,nx,y,z); + _cimg_watershed_propagate(is_py,x,py,z); + _cimg_watershed_propagate(is_ny,x,ny,z); + if (is_3d) { + _cimg_watershed_propagate(is_pz,x,y,pz); + _cimg_watershed_propagate(is_nz,x,y,nz); } + if (is_high_connectivity) { + _cimg_watershed_propagate(is_px && is_py,px,py,z); + _cimg_watershed_propagate(is_nx && is_py,nx,py,z); + _cimg_watershed_propagate(is_px && is_ny,px,ny,z); + _cimg_watershed_propagate(is_nx && is_ny,nx,ny,z); + if (is_3d) { + _cimg_watershed_propagate(is_px && is_pz,px,y,pz); + _cimg_watershed_propagate(is_nx && is_pz,nx,y,pz); + _cimg_watershed_propagate(is_px && is_nz,px,y,nz); + _cimg_watershed_propagate(is_nx && is_nz,nx,y,nz); + _cimg_watershed_propagate(is_py && is_pz,x,py,pz); + _cimg_watershed_propagate(is_ny && is_pz,x,ny,pz); + _cimg_watershed_propagate(is_py && is_nz,x,py,nz); + _cimg_watershed_propagate(is_ny && is_nz,x,ny,nz); + _cimg_watershed_propagate(is_px && is_py && is_pz,px,py,pz); + _cimg_watershed_propagate(is_nx && is_py && is_pz,nx,py,pz); + _cimg_watershed_propagate(is_px && is_ny && is_pz,px,ny,pz); + _cimg_watershed_propagate(is_nx && is_ny && is_pz,nx,ny,pz); + _cimg_watershed_propagate(is_px && is_py && is_nz,px,py,nz); + _cimg_watershed_propagate(is_nx && is_py && is_nz,nx,py,nz); + _cimg_watershed_propagate(is_px && is_ny && is_nz,px,ny,nz); + _cimg_watershed_propagate(is_nx && is_ny && is_nz,nx,ny,nz); + } + } + (*this)(x,y,z) = label; + labels(x,y,z) = ++nmin; } return *this; } //! Compute watershed transform \newinstance. template - CImg get_watershed(const CImg& priority, const bool fill_lines=true) const { - return (+*this).watershed(priority,fill_lines); + CImg get_watershed(const CImg& priority, const bool is_high_connectivity=false) const { + return (+*this).watershed(priority,is_high_connectivity); } // [internal] Insert/Remove items in priority queue, for watershed/distance transforms. - template - bool _priority_queue_insert(CImg& is_queued, unsigned int& siz, const t value, - const unsigned int x, const unsigned int y, const unsigned int z) { + template + bool _priority_queue_insert(CImg& is_queued, unsigned int& siz, const tv value, + const unsigned int x, const unsigned int y, const unsigned int z, + const unsigned int n=1) { if (is_queued(x,y,z)) return false; - is_queued(x,y,z) = true; + is_queued(x,y,z) = (tq)n; if (++siz>=_width) { if (!is_empty()) resize(_width*2,4,1,1,0); else assign(64,4); } - (*this)(siz - 1,0) = (T)value; (*this)(siz - 1,1) = (T)x; (*this)(siz - 1,2) = (T)y; (*this)(siz - 1,3) = (T)z; + (*this)(siz - 1,0) = (T)value; + (*this)(siz - 1,1) = (T)x; + (*this)(siz - 1,2) = (T)y; + (*this)(siz - 1,3) = (T)z; for (unsigned int pos = siz - 1, par = 0; pos && value>(*this)(par=(pos + 1)/2 - 1,0); pos = par) { - cimg::swap((*this)(pos,0),(*this)(par,0)); cimg::swap((*this)(pos,1),(*this)(par,1)); - cimg::swap((*this)(pos,2),(*this)(par,2)); cimg::swap((*this)(pos,3),(*this)(par,3)); + cimg::swap((*this)(pos,0),(*this)(par,0)); + cimg::swap((*this)(pos,1),(*this)(par,1)); + cimg::swap((*this)(pos,2),(*this)(par,2)); + cimg::swap((*this)(pos,3),(*this)(par,3)); } return true; } CImg& _priority_queue_remove(unsigned int& siz) { - (*this)(0,0) = (*this)(--siz,0); (*this)(0,1) = (*this)(siz,1); - (*this)(0,2) = (*this)(siz,2); (*this)(0,3) = (*this)(siz,3); + (*this)(0,0) = (*this)(--siz,0); + (*this)(0,1) = (*this)(siz,1); + (*this)(0,2) = (*this)(siz,2); + (*this)(0,3) = (*this)(siz,3); const float value = (*this)(0,0); for (unsigned int pos = 0, left = 0, right = 0; ((right=2*(pos + 1),(left=right - 1))(*this)(right,0)) { - cimg::swap((*this)(pos,0),(*this)(left,0)); cimg::swap((*this)(pos,1),(*this)(left,1)); - cimg::swap((*this)(pos,2),(*this)(left,2)); cimg::swap((*this)(pos,3),(*this)(left,3)); + cimg::swap((*this)(pos,0),(*this)(left,0)); + cimg::swap((*this)(pos,1),(*this)(left,1)); + cimg::swap((*this)(pos,2),(*this)(left,2)); + cimg::swap((*this)(pos,3),(*this)(left,3)); pos = left; } else { - cimg::swap((*this)(pos,0),(*this)(right,0)); cimg::swap((*this)(pos,1),(*this)(right,1)); - cimg::swap((*this)(pos,2),(*this)(right,2)); cimg::swap((*this)(pos,3),(*this)(right,3)); + cimg::swap((*this)(pos,0),(*this)(right,0)); + cimg::swap((*this)(pos,1),(*this)(right,1)); + cimg::swap((*this)(pos,2),(*this)(right,2)); + cimg::swap((*this)(pos,3),(*this)(right,3)); pos = right; } } else { - cimg::swap((*this)(pos,0),(*this)(left,0)); cimg::swap((*this)(pos,1),(*this)(left,1)); - cimg::swap((*this)(pos,2),(*this)(left,2)); cimg::swap((*this)(pos,3),(*this)(left,3)); + cimg::swap((*this)(pos,0),(*this)(left,0)); + cimg::swap((*this)(pos,1),(*this)(left,1)); + cimg::swap((*this)(pos,2),(*this)(left,2)); + cimg::swap((*this)(pos,3),(*this)(left,3)); pos = left; } } @@ -39484,11 +40232,11 @@ namespace cimg_library_suffixed { lZ = sprite.depth() - (z0 + sprite.depth()>depth()?z0 + sprite.depth() - depth():0) + (bz?z0:0), lC = sprite.spectrum() - (c0 + sprite.spectrum()>spectrum()?c0 + sprite.spectrum() - spectrum():0) + (bc?c0:0); const t - *ptrs = sprite._data - - (bx?x0:0) - - (by?y0*(uptrT)sprite.width():0) - - (bz?z0*(uptrT)sprite.width()*sprite.height():0) - - (bc?c0*(uptrT)sprite.width()*sprite.height()*sprite.depth():0); + *ptrs = sprite._data + + (bx?-x0:0) + + (by?-y0*(uptrT)sprite.width():0) + + (bz?-z0*(uptrT)sprite.width()*sprite.height():0) + + (bc?-c0*(uptrT)sprite.width()*sprite.height()*sprite.depth():0); const uptrT offX = (unsigned long)_width - lX, soffX = (unsigned long)sprite._width - lX, @@ -39528,11 +40276,11 @@ namespace cimg_library_suffixed { lZ = sprite.depth() - (z0 + sprite.depth()>depth()?z0 + sprite.depth() - depth():0) + (bz?z0:0), lC = sprite.spectrum() - (c0 + sprite.spectrum()>spectrum()?c0 + sprite.spectrum() - spectrum():0) + (bc?c0:0); const T - *ptrs = sprite._data - - (bx?x0:0) - - (by?y0*(uptrT)sprite.width():0) - - (bz?z0*(uptrT)sprite.width()*sprite.height():0) - - (bc?c0*(uptrT)sprite.width()*sprite.height()*sprite.depth():0); + *ptrs = sprite._data + + (bx?-x0:0) + + (by?-y0*(uptrT)sprite.width():0) + + (bz?-z0*(uptrT)sprite.width()*sprite.height():0) + + (bc?-c0*(uptrT)sprite.width()*sprite.height()*sprite.depth():0); const unsigned long offX = (unsigned long)_width - lX, soffX = (unsigned long)sprite._width - lX, @@ -39623,10 +40371,10 @@ namespace cimg_library_suffixed { lZ = sprite.depth() - (z0 + sprite.depth()>depth()?z0 + sprite.depth() - depth():0) + (bz?z0:0), lC = sprite.spectrum() - (c0 + sprite.spectrum()>spectrum()?c0 + sprite.spectrum() - spectrum():0) + (bc?c0:0); const uptrT - coff = -(bx?x0:0) - - (by?y0*(uptrT)mask.width():0) - - (bz?z0*(uptrT)mask.width()*mask.height():0) - - (bc?c0*(uptrT)mask.width()*mask.height()*mask.depth():0), + coff = (bx?-x0:0) + + (by?-y0*(uptrT)mask.width():0) + + (bz?-z0*(uptrT)mask.width()*mask.height():0) + + (bc?-c0*(uptrT)mask.width()*mask.height()*mask.depth():0), ssize = (uptrT)mask.width()*mask.height()*mask.depth()*mask.spectrum(); const ti *ptrs = sprite._data + coff; const tm *ptrm = mask._data + coff; @@ -43423,7 +44171,7 @@ namespace cimg_library_suffixed { unsigned int cdx = 0, dx = 0, dy = 0; int err = 0; double val; - assign(256,256); + assign(256,256,1,1,0); while ((err = std::fscanf(nfile,"%lf%255[^0-9eEinfa.+-]",&val,delimiter._data))>0) { if (err>0) (*this)(cdx++,dy) = (T)val; if (cdx>=_width) resize(3*_width/2,_height,1,1,0); @@ -48536,7 +49284,7 @@ namespace cimg_library_suffixed { **/ const CImg& save_tiff(const char *const filename, const unsigned int compression_type=0, const float *const voxel_size=0, const char *const description=0, - const bool is_bigtiff=true) const { + const bool use_bigtiff=true) const { if (!filename) throw CImgArgumentException(_cimg_instance "save_tiff(): Specified filename is (null).", @@ -48544,7 +49292,9 @@ namespace cimg_library_suffixed { if (is_empty()) { cimg::fempty(0,filename); return *this; } #ifdef cimg_use_tiff - TIFF *tif = TIFFOpen(filename,is_bigtiff?"w8":"w4"); + const bool + _use_bigtiff = use_bigtiff && sizeof(uptrT)>=8 && size()*sizeof(T)>=1UL<<31; // No bigtiff for small images. + TIFF *tif = TIFFOpen(filename,_use_bigtiff?"w8":"w4"); if (tif) { cimg_forZ(*this,z) _save_tiff(tif,z,z,compression_type,voxel_size,description); TIFFClose(tif); @@ -48554,7 +49304,7 @@ namespace cimg_library_suffixed { filename); return *this; #else - cimg::unused(compression_type,voxel_size,description,is_bigtiff); + cimg::unused(compression_type,voxel_size,description,use_bigtiff); return save_other(filename); #endif } @@ -49078,9 +49828,9 @@ namespace cimg_library_suffixed { if (is_empty()) { cimg::fempty(file,filename); return *this; } std::FILE *const nfile = file?file:cimg::fopen(filename,"wb"); - static const unsigned char header[36] = { 'P','A','N','D','O','R','E','0','4',0,0,0, - 0,0,0,0,'C','I','m','g',0,0,0,0,0, - 'N','o',' ','d','a','t','e',0,0,0,0 }; + unsigned char header[36] = { 'P','A','N','D','O','R','E','0','4',0,0,0, + 0,0,0,0,'C','I','m','g',0,0,0,0,0, + 'N','o',' ','d','a','t','e',0,0,0,0 }; unsigned int nbdims, dims[5] = { 0 }; bool saved = false; _cimg_save_pandore_case(1,1,1,"unsigned char",2); @@ -52584,7 +53334,7 @@ namespace cimg_library_suffixed { *tmp = *str_pixeltype = *str_endian = 0; unsigned int j, N, W, H, D, C; int i, err; - j = 0; while((i=std::fgetc(nfile))!='\n' && i!=EOF && j<256) tmp[j++] = (char)i; tmp[j] = 0; + j = 0; while ((i=std::fgetc(nfile))!='\n' && i!=EOF && j<256) tmp[j++] = (char)i; tmp[j] = 0; err = cimg_sscanf(tmp,"%u%*c%255[A-Za-z_]%*c%255[sA-Za-z_ ]", &N,str_pixeltype._data,str_endian._data); if (err<2) { @@ -53773,7 +54523,7 @@ namespace cimg_library_suffixed { #define _cimg_save_cimg_case(Ts,Tss) \ if (!saved && !cimg::strcasecmp(Ts,str_pixeltype)) { \ for (unsigned int l = 0; l& save_tiff(const char *const filename, const unsigned int compression_type=0, const float *const voxel_size=0, const char *const description=0, - const bool is_bigtiff=true) const { + const bool use_bigtiff=true) const { if (!filename) throw CImgArgumentException(_cimglist_instance "save_tiff(): Specified filename is (null).", @@ -53974,14 +54724,17 @@ namespace cimg_library_suffixed { if (is_empty()) { cimg::fempty(0,filename); return *this; } #ifndef cimg_use_tiff - if (_width==1) _data[0].save_tiff(filename,compression_type,voxel_size,description,is_bigtiff); + if (_width==1) _data[0].save_tiff(filename,compression_type,voxel_size,description,use_bigtiff); else cimglist_for(*this,l) { CImg nfilename(1024); cimg::number_filename(filename,l,6,nfilename); - _data[l].save_tiff(nfilename,compression_type,voxel_size,description,is_bigtiff); + _data[l].save_tiff(nfilename,compression_type,voxel_size,description,use_bigtiff); } #else - TIFF *tif = TIFFOpen(filename,is_bigtiff?"w8":"w4"); + typename CImg::uptrT siz = 0; + cimglist_for(*this,l) siz+=_data[l].size(); + const bool _use_bigtiff = use_bigtiff && sizeof(siz)>=8 && siz*sizeof(T)>=1UL<<31; // No bigtiff for small images. + TIFF *tif = TIFFOpen(filename,_use_bigtiff?"w8":"w4"); if (tif) { for (unsigned int dir = 0, l = 0; l<_width; ++l) { const CImg& img = (*this)[l]; diff --git a/stim/image/image.h b/stim/image/image.h index e10e6b0..b92a011 100644 --- a/stim/image/image.h +++ b/stim/image/image.h @@ -35,7 +35,7 @@ public: img = cimg_library::CImg(x, y, z); }*/ - image(unsigned int x, unsigned int y = 1, unsigned int z = 1, unsigned int c = 1){ + image(size_t x, size_t y = 1, size_t z = 1, size_t c = 1){ img = cimg_library::CImg(x, y, z, c); } @@ -50,13 +50,13 @@ public: } //create an image from an interleaved buffer - void set_interleaved(T* buffer, unsigned int width, unsigned int height, unsigned int channels = 1){ + void set_interleaved(T* buffer, size_t width, size_t height, size_t channels = 1){ T* non_interleaved = (T*)malloc(width * height * 3 * sizeof(T)); - unsigned int S = width * height; + size_t S = width * height; - for(unsigned int i = 0; i < S; i++){ - for(unsigned int c = 0; c < channels; c++){ + for(size_t i = 0; i < S; i++){ + for(size_t c = 0; c < channels; c++){ non_interleaved[i + c * S] = buffer[i * channels + c]; } } @@ -71,19 +71,19 @@ public: void data_interleaved(T* data){ - unsigned int C = channels(); - unsigned int X = width() * height(); + size_t C = channels(); + size_t X = width() * height(); T* ptr = img.data(); //for each channel - for(unsigned int c = 0; c < C; c++) + for(size_t c = 0; c < C; c++) //convert each pixel - for(unsigned int x = 0; x < X; x++) + for(size_t x = 0; x < X; x++) data[x * C + c] = ptr[c * X + x]; } - image channel(unsigned int c){ + image channel(size_t c){ //create a new image image single; @@ -94,7 +94,7 @@ public: } - T& operator()(unsigned x, unsigned y, unsigned z = 0, unsigned c = 0){ + T& operator()(size_t x, size_t y, size_t z = 0, size_t c = 0){ return img(x, y, z, c); } @@ -103,9 +103,9 @@ public: /// @param v is the value used to set all values in the image image operator=(T v){ - unsigned int X = width() * height() * depth() * channels(); + size_t X = width() * height() * depth() * channels(); - for(unsigned x = 0; x < X; x++) + for(size_t x = 0; x < X; x++) img.data()[x] = v; return *this; @@ -117,10 +117,10 @@ public: /// @param c is the channel number that the data will be copied to /// @param buffer is a pointer to the image to be copied to channel c - void set_channel(unsigned int c, T* buffer){ + void set_channel(size_t c, T* buffer){ //calculate the number of pixels in a channel - unsigned int channel_size = width() * height(); + size_t channel_size = width() * height(); //retreive a pointer to the raw image data T* ptr = img.data() + channel_size * c; @@ -129,7 +129,7 @@ public: memcpy(ptr, buffer, sizeof(T) * channel_size); } - image getslice(unsigned int c){ + image getslice(size_t c){ //create a new image image slice; @@ -140,19 +140,19 @@ public: } - unsigned int channels(){ - return (unsigned int)img.spectrum(); + size_t channels(){ + return (size_t)img.spectrum(); } - unsigned int width(){ + size_t width(){ return img.width(); } - unsigned int height(){ + size_t height(){ return img.height(); } - unsigned int depth(){ + size_t depth(){ return img.depth(); } @@ -162,22 +162,22 @@ public: } //returns the size (number of values) of the image - unsigned long size(){ + size_t size(){ return img.size(); } /// Returns the number of nonzero values - unsigned int nnz(){ + size_t nnz(){ - unsigned long P = width() * height(); - unsigned long C = channels(); + size_t P = width() * height(); + size_t C = channels(); T* ptr = img.data(); - unsigned long n = 0; + size_t n = 0; - for(unsigned long p = 0; p < P; p++){ - for(unsigned long c = 0; c < C; c++){ + for(size_t p = 0; p < P; p++){ + for(size_t c = 0; c < C; c++){ if(ptr[c * P + p] > 0){ n++; @@ -191,19 +191,19 @@ public: } //this function returns indices of pixels that have nonzero values - std::vector sparse_idx(){ + std::vector sparse_idx(){ - std::vector s; //allocate an array + std::vector s; //allocate an array s.resize(nnz()); //allocate space in the array - unsigned long P = width() * height(); - unsigned long C = channels(); + size_t P = width() * height(); + size_t C = channels(); T* ptr = img.data(); //get a pointer to the image data - unsigned long i = 0; - for(unsigned long p = 0; p < P; p++){ - for(unsigned long c = 0; c < C; c++){ + size_t i = 0; + for(size_t p = 0; p < P; p++){ + for(size_t c = 0; c < C; c++){ if(ptr[c * P + p] > 0){ s[i] = p; @@ -220,9 +220,9 @@ public: /// Returns the maximum pixel value in the image T maxv(){ float max = 0; - unsigned long N = width() * height(); //get the number of pixels + size_t N = width() * height(); //get the number of pixels - for (unsigned long i=0; i max) { @@ -236,9 +236,9 @@ public: /// Returns the minimum pixel value in the image T minv(){ float min = 0; - unsigned long N = width() * height(); //get the number of pixels + size_t N = width() * height(); //get the number of pixels - for (unsigned long i=0; i result; float zoom = 1; - unsigned int interpolation = 1; - unsigned int boundary = 1; + size_t interpolation = 1; + size_t boundary = 1; result.img = img.get_rotate (angle, cx, cy, zoom, interpolation, boundary); //result.save("data_output/test_rotate_neum.bmp"); @@ -286,13 +286,13 @@ public: // leila's code for non_interleaving data in 3D //create an data set from an interleaved buffer - void set_interleaved3(T* buffer, unsigned int width, unsigned int height, unsigned int depth, unsigned int channels = 3){ + void set_interleaved3(T* buffer, size_t width, size_t height, size_t depth, size_t channels = 3){ T* non_interleaved3 = (T*)malloc(width * height * depth * 3 * sizeof(T)); - unsigned int p = width * height * depth; + size_t p = width * height * depth; - for(unsigned int i = 0; i < p; i++){ - for(unsigned int c = 0; c < channels; c++){ + for(size_t i = 0; i < p; i++){ + for(size_t c = 0; c < channels; c++){ non_interleaved3[i + c * p] = buffer[i * channels + c]; } } diff --git a/stim/parser/arguments.h b/stim/parser/arguments.h index 79c4a22..d40933c 100644 --- a/stim/parser/arguments.h +++ b/stim/parser/arguments.h @@ -76,7 +76,7 @@ namespace stim{ } - int nargs() + size_t nargs() { return vals.size(); } @@ -135,9 +135,9 @@ namespace stim{ } //get the width of the left column - int col_width() + size_t col_width() { - int n = 3; + size_t n = 3; //add the length of the option name n += name.size(); @@ -148,7 +148,7 @@ namespace stim{ n += 6; //for each default option value - for(unsigned int v=0; v args; //column width of the longest option - int col_width; + size_t col_width; //list of sections std::vector sections; @@ -327,7 +327,7 @@ namespace stim{ opt.set_ansi(ansi); opts.push_back(opt); - col_width = std::max(col_width, opt.col_width()); + col_width = std::max(col_width, opt.col_width()); } ///Specifies a section name that can be used to organize parameters in the output. @@ -373,32 +373,28 @@ namespace stim{ ///Retrieves the index for a supported argument, given that argument's name. /// @param _name is the name of the requested argument - int index(std::string _name) + size_t index(std::string _name) { - unsigned int i = find(opts.begin(), opts.end(), _name) - opts.begin(); + size_t i = find(opts.begin(), opts.end(), _name) - opts.begin(); - if(i >= opts.size()) - return -1; + if(i >= opts.size()){ + std::cout<<"ERROR stim::arglist: option name '"<<_name<<"' not found"<= opts.size()) { diff --git a/stim/parser/filename.h b/stim/parser/filename.h index 62fc8a7..cc1f178 100644 --- a/stim/parser/filename.h +++ b/stim/parser/filename.h @@ -305,7 +305,7 @@ public: return result; //return the result } - stim::filename insert(unsigned i, unsigned int n = 2){ + stim::filename insert(size_t i, size_t n = 2){ std::stringstream ss; ss << std::setw(n) << std::setfill('0') << i; diff --git a/stim/visualization/colormap.h b/stim/visualization/colormap.h index 8a43c33..6f9279c 100644 --- a/stim/visualization/colormap.h +++ b/stim/visualization/colormap.h @@ -36,7 +36,7 @@ namespace stim{ enum colormapType {cmBrewer, cmGrayscale, cmRainbow}; -static void buffer2image(unsigned char* buffer, std::string filename, unsigned int width, unsigned int height) +static void buffer2image(unsigned char* buffer, std::string filename, size_t width, size_t height) { /*unsigned char* non_interleaved = (unsigned char*)malloc(x_size * y_size * 3); unsigned int S = x_size * y_size; @@ -216,9 +216,9 @@ static void gpu2image(T* gpuSource, std::string fileDest, unsigned int x_size, u #endif template -static void cpuApplyBrewer(T* cpuSource, unsigned char* cpuDest, unsigned int N, T minVal = 0, T maxVal = 1) +static void cpuApplyBrewer(T* cpuSource, unsigned char* cpuDest, size_t N, T minVal = 0, T maxVal = 1) { - for(int i=0; i -static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T valMin, T valMax, colormapType cm = cmGrayscale) +static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, size_t nVals, T valMin, T valMax, colormapType cm = cmGrayscale) { if(cm == cmBrewer) @@ -279,15 +279,15 @@ static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T if(a < 0) a = 0; if(a > 1) a = 1; - cpuDest[i * 3 + 0] = 255 * a; - cpuDest[i * 3 + 1] = 255 * a; - cpuDest[i * 3 + 2] = 255 * a; + cpuDest[i * 3 + 0] = (unsigned char)(255 * a); + cpuDest[i * 3 + 1] = (unsigned char)(255 * a); + cpuDest[i * 3 + 2] = (unsigned char)(255 * a); } } } template -static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale) +static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned long long nVals, colormapType cm = cmGrayscale) { //computes the max and min range automatically @@ -313,7 +313,7 @@ static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, co template -static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, T valMin, T valMax, colormapType cm = cmGrayscale) +static void cpu2image(T* cpuSource, std::string fileDest, size_t x_size, size_t y_size, T valMin, T valMax, colormapType cm = cmGrayscale) { //allocate a color buffer unsigned char* cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size); @@ -329,7 +329,7 @@ static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, u } template -static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, colormapType cm = cmGrayscale) +static void cpu2image(T* cpuSource, std::string fileDest, size_t x_size, size_t y_size, colormapType cm = cmGrayscale) { //allocate a color buffer unsigned char* cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size); -- libgit2 0.21.4