fixed asynchronous read bugs

David Mayerich
1 parent 8cfba32f
Showing 4 changed files with 142 additions and 1 deletions Show diff stats
stim/biomodels/cellset.h
stim/envi/binary.h
stim/envi/bsq.h
stim/envi/hsi.h
@@ -15,6 +15,7 @@ protected:
 	std::vector<double*> cells;							//vector storing field data for each cell
 	std::unordered_map<std::string, size_t> fields;		//unordered map storing field->index information for each field
 	size_t ip[3];										//hard code to position indices (for speed)
+
 	void init(){
  
 	}
@@ -80,6 +81,13 @@ public:
 		return cells[c][idx];
 	}
  
+	/// returns an ID used to look up a field
+	bool exists(std::string f){
+		std::unordered_map<std::string, size_t>::iterator iter = fields.find(f);
+		if(iter == fields.end()) return false;
+		else return true;
+	}
+
 	/// Return the position of cell [i]
 	stim::vec3<double> p(size_t i){
 		stim::vec3<double> pos(cells[i][ip[0]], cells[i][ip[1]], cells[i][ip[2]]);
@@ -403,6 +403,48 @@ public:
 		return read_pixel(p, i);
 	}
  
+	/// Reads a block specified by an (x, y, z) position and size using the largest possible contiguous reads
+	bool read_block(T* dest, size_t x, size_t y, size_t z, size_t sx, size_t sy, size_t sz){
+
+		size_t size_bytes = sx * sy * sz * sizeof(T);					//size of the block to read in bytes
+
+		size_t start = z * R[0] * R[1] + y * R[0] + x;						//calculate the start postion
+		size_t start_bytes = start * sizeof(T);							//start position in bytes
+		file.seekg(start * sizeof(T), std::ios::beg);					//seek to the start position
+
+		
+		if(sx == R[0] && sy == R[1]){				//if sx and sy result in a contiguous volume along z
+			file.read((char*)dest, size_bytes);			//read the block in one pass
+			return true;
+		}
+
+		if(sx == R[0]){												//if sx is contiguous, read each z-axis slice can be read in one pass
+			size_t jump_bytes = (R[1] - sy) * R[0] * sizeof(T);		//jump between each slice
+			size_t slice_bytes = sx * sy * sizeof(T);				//size of the slice to be read
+			for(size_t zi = 0; zi < sz; zi++){						//for each z-axis slice
+				file.read((char*)dest, slice_bytes);						//read the slice
+				dest += sx * sy;									//move the destination pointer to the next slice
+				file.seekg(jump_bytes, std::ios::cur);				//skip to the next slice in the file
+			}
+			return true;
+		}
+
+		//in this case, x is not contiguous so the volume must be read line-by-line
+		size_t jump_x_bytes = (R[0] - sx) * sizeof(T);				//number of bytes skipped in the x direction
+		size_t jump_y_bytes = (R[1] - sy) * R[0] * sizeof(T) + jump_x_bytes;	//number of bytes skipped between slices
+		size_t line_bytes = sx * sizeof(T);							//size of the line to be read
+		size_t zi, yi;
+		for(zi = 0; zi < sz; zi++){									//for each slice
+			file.read((char*)dest, line_bytes);							//read the first line
+			for(yi = 1; yi < sy; yi++){								//read each additional line
+				dest += sx;											//move the pointer in the destination block to the next line
+				file.seekg(jump_x_bytes, std::ios::cur);			//skip to the next line in the file
+				file.read((char*)dest, line_bytes);						//read the line to the destination block
+			}
+			file.seekg(jump_y_bytes, std::ios::cur);				//skip to the beginning of the next slice
+		}
+	}
+
 };
  
 }
@@ -401,6 +401,79 @@ public:
 		size_t batches = (size_t)ceil((double)(Y()) / (double)batch_slices);			//calculate the number of batches
 		T* ptrDst;
 		T* ptrSrc;
+		size_t y = 0;			//initialize the current y-slice position
+		for(size_t c = 0; c < batches; c++){
+			if(c == (batches - 1)){
+				batch_slices = Y() - (batches - 1) * batch_slices;				//if this is the last batch, calculate the remaining # of bands
+				jump = (Y() - batch_slices) * X() * sizeof(T);
+				batchN = XB * batch_slices;
+				batch_bytes = batchN * sizeof(T);
+			}
+
+			auto in_begin = std::chrono::high_resolution_clock::now();			
+			hsi<T>::read_block(ptrIn, 0, y, 0, X(), batch_slices, Z());			//read the input block
+			y += batch_slices;
+			auto in_end = std::chrono::high_resolution_clock::now();
+			in_time += std::chrono::duration_cast<std::chrono::milliseconds>(in_end-in_begin).count();
+
+			auto calc_begin = std::chrono::high_resolution_clock::now();
+			
+			for(size_t b = 0; b < Z(); b++){									//for each line, store an XB slice in ptrDest
+				ptrSrc = ptrIn + (b * X() * batch_slices);
+				ptrDst = ptrOut + (b * X());									//initialize ptrDst to the start of the XB output slice
+				
+				for(size_t y = 0; y < batch_slices; y++){						//for each band in the current line
+					memcpy(ptrDst, ptrSrc, X() * sizeof(T));					//copy the band line from the source to the destination
+					ptrSrc += X();												//increment the pointer within the current buffer array (batch)	
+					ptrDst += X() * Z();										//increment the pointer within the XB slice (to be output)									
+				}
+			}
+			auto calc_end = std::chrono::high_resolution_clock::now();
+			calc_time += std::chrono::duration_cast<std::chrono::milliseconds>(calc_end-calc_begin).count();
+
+			auto out_begin = std::chrono::high_resolution_clock::now();
+			target.write((char*)ptrOut, batch_bytes);							//write the batch to disk
+			auto out_end = std::chrono::high_resolution_clock::now();
+			out_time += std::chrono::duration_cast<std::chrono::milliseconds>(out_end-out_begin).count();
+			if(PROGRESS) progress = (double)( c + 1 ) / (batches) * 100;
+		}
+
+		free(ptrIn);
+		free(ptrOut);
+		target.close();
+
+		std::cout<<"BSQ->BIL reads: "<<(double)in_time / (1000 * 60)<<" min"<<std::endl;
+		std::cout<<"BSQ->BIL calculations: "<<(double)calc_time / (1000 * 60)<<" min"<<std::endl;
+		std::cout<<"BSQ->BIL writes: "<<(double)out_time / (1000 * 60)<<" min"<<std::endl;
+
+		return true;
+	}
+
+	/*bool bil(std::string outname, bool PROGRESS = false){
+
+		size_t in_time, out_time, calc_time;									//initialize the timing variables
+		in_time = out_time = calc_time = 0;
+
+		size_t XY = X() * Y();													//number of elements in an input slice
+		size_t XB = X() * Z();													//number of elements in an output slice
+		size_t XYbytes = XY * sizeof(T);										//number of bytes in an input slice
+		size_t XBbytes = XB * sizeof(T);										//number of bytes in an output slice
+		size_t batch_slices = binary<T>::buffer_size / (2*XBbytes);				//calculate the number of slices that can fit in memory
+		if(Y() < batch_slices) batch_slices = Y();								//if the entire data set will fit in memory, do it
+		size_t batchN = XB * batch_slices;										//number of elements in a batch
+		size_t batch_bytes = batchN * sizeof(T);								//calculate the number of bytes in a batch
+
+		T* ptrIn = (T*) malloc(batch_bytes);									//allocate a large buffer storing the read data
+		T* ptrOut = (T*) malloc(batch_bytes);									//allocate space for storing an output buffer
+
+		size_t jump = (Y() - batch_slices) * X() * sizeof(T);					//jump between reads in the input file
+
+		std::ofstream target(outname.c_str(), std::ios::binary);
+		std::string headername = outname + ".hdr";
+
+		size_t batches = (size_t)ceil((double)(Y()) / (double)batch_slices);			//calculate the number of batches
+		T* ptrDst;
+		T* ptrSrc;
 		for(size_t c = 0; c < batches; c++){
 			file.seekg(c * X() * batch_slices * sizeof(T), std::ios::beg);
  
@@ -450,7 +523,7 @@ public:
 		std::cout<<"BSQ->BIL writes: "<<(double)out_time / (1000 * 60)<<" min"<<std::endl;
  
 		return true;
-	}
+	}*/
  
 	/*/// Convert the current BSQ file to a BIL file with the specified file name.
 	bool bil(std::string outname, bool PROGRESS = false){
@@ -202,6 +202,24 @@ public:
 		}
 	}
  
+	void read_block(T* dest, size_t x, size_t y, size_t z, size_t sx, size_t sy, size_t sz){
+		size_t d[3];					//position in the binary coordinate system
+		size_t sd[3];					//size in the binary coordinate system
+
+		d[O[0]] = x;					//set the point in the binary coordinate system
+		d[O[1]] = y;
+		d[O[2]] = z;
+
+		sd[O[0]] = sx;					//set the size in the binary coordinate system
+		sd[O[1]] = sy;
+		sd[O[2]] = sz;
+
+		if(!binary<T>::read_block(dest, d[0], d[1], d[2], sd[0], sd[1], sd[2])){
+			std::cout<<"error reading block in stim::hsi: ("<<d[0]<<", "<<d[1]<<", "<<d[2]<<") - ["<<sd[0]<<", "<<sd[1]<<", "<<sd[2]<<"]"<<std::endl;
+			exit(1);
+		}
+	}
+
 };
  
 }		//end namespace STIM