From 8b899c24248604acddbaefeebd1a6a8e0bdd9dcf Mon Sep 17 00:00:00 2001 From: David Date: Wed, 7 Sep 2016 23:29:43 -0500 Subject: [PATCH] fixed asynchronous read bugs --- stim/biomodels/cellset.h | 8 ++++++++ stim/envi/binary.h | 42 ++++++++++++++++++++++++++++++++++++++++++ stim/envi/bsq.h | 75 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- stim/envi/hsi.h | 18 ++++++++++++++++++ 4 files changed, 142 insertions(+), 1 deletion(-) diff --git a/stim/biomodels/cellset.h b/stim/biomodels/cellset.h index 1dccc1c..3d02028 100644 --- a/stim/biomodels/cellset.h +++ b/stim/biomodels/cellset.h @@ -15,6 +15,7 @@ protected: std::vector cells; //vector storing field data for each cell std::unordered_map 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::iterator iter = fields.find(f); + if(iter == fields.end()) return false; + else return true; + } + /// Return the position of cell [i] stim::vec3 p(size_t i){ stim::vec3 pos(cells[i][ip[0]], cells[i][ip[1]], cells[i][ip[2]]); diff --git a/stim/envi/binary.h b/stim/envi/binary.h index 2af925f..ff8a725 100644 --- a/stim/envi/binary.h +++ b/stim/envi/binary.h @@ -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 + } + } + }; } diff --git a/stim/envi/bsq.h b/stim/envi/bsq.h index 2d3fb67..f8eb890 100644 --- a/stim/envi/bsq.h +++ b/stim/envi/bsq.h @@ -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::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(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(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(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"<BIL calculations: "<<(double)calc_time / (1000 * 60)<<" min"<BIL writes: "<<(double)out_time / (1000 * 60)<<" min"<::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"<::read_block(dest, d[0], d[1], d[2], sd[0], sd[1], sd[2])){ + std::cout<<"error reading block in stim::hsi: ("<