diff --git a/stim/envi/bsq.h b/stim/envi/bsq.h
index 8be33fb..5ebf65a 100644
--- a/stim/envi/bsq.h
+++ b/stim/envi/bsq.h
@@ -1237,6 +1237,53 @@ public:
 		return true;
 	}
 
+	///Crop out several subimages and assemble a new image from these concatenated subimages
+
+	/// @param outfile is the file name for the output image
+	/// @param sx is the width of each subimage
+	/// @param sy is the height of each subimage
+	/// @mask is the mask used to define subimage positions extracted from the input file
+	void subimages(std::string outfile, size_t sx, size_t sy, unsigned char* mask, bool PROGRESS = false){
+
+		size_t N = nnz(mask);									//get the number of subimages
+		T* dst = (T*) malloc(N * sx * sy * sizeof(T));			//allocate space for a single band of the output image
+		memset(dst, 0, N*sx*sy*sizeof(T));						//initialize the band image to zero
+
+		std::ofstream out(outfile, std::ios::binary);			//open a file for writing
+
+		T* src = (T*) malloc(X() * Y() * sizeof(T));
+
+		for(size_t b = 0; b < Z(); b++){						//for each band
+			band_index(src, b);									//load the band image
+			size_t i = 0;										//create an image index and initialize it to zero
+			size_t n = 0;
+			while(n < N){										//for each subimage
+				if(mask[i]){									//if the pixel is masked, copy the surrounding pixels into the destination band
+					size_t yi = i / X();						//determine the y position of the current pixel
+					size_t xi = i - yi * X();					//determine the x position of the current pixel
+					if( xi > sx/2 && xi < X() - sx/2 &&			//if the subimage is completely within the bounds of the original image
+						yi > sy/2 && yi < Y() - sy/2){
+						size_t cx = xi - sx/2;					//calculate the corner position for the subimage
+						size_t cy = yi - sy/2;
+						size_t cxi, cyi;
+						for(size_t syi = 0; syi < sy; syi++){					//for each line in the subimage
+							size_t src_i = (cy + syi) * X() + cx;
+							//size_t dst_i = syi * (N * sx) + n * sx;
+							size_t dst_i = (n * sy + syi) * sx;
+							memcpy(&dst[dst_i],  &src[src_i], sx * sizeof(T));	//copy one line from the subimage to the destination image
+						}
+						n++;
+					}
+				}
+				i++;
+				if(PROGRESS) progress = (double)( (n+1) * (b+1) ) / (N * Z()) * 100;
+			}//end while n
+			out.write((const char*)dst, N * sx * sy * sizeof(T));			//write the band to memory
+		}
+		free(dst);												//free memory
+		free(src);
+	}
+
 	/// Remove a list of bands from the ENVI file
 
 	/// @param outfile is the file name for the output hyperspectral image (with trimmed bands)
diff --git a/stim/envi/envi.h b/stim/envi/envi.h
index 5ebba33..66b2e60 100644
--- a/stim/envi/envi.h
+++ b/stim/envi/envi.h
@@ -1558,6 +1558,41 @@ public:
 		return false;
 	}
 
+	void subimages(std::string outfile, size_t nx, size_t ny, unsigned char* mask, bool PROGRESS = false){
+		
+		size_t nnz = 0;													//initialize the number of subimages to zero
+		for(size_t i = 0; i < header.lines * header.samples; i++)		//for each pixel in the mask
+			if(mask[i]) nnz++;											//if the pixel is valid, add a subimage
+
+
+		//save the header for the cropped file
+		stim::envi_header new_header = header;
+		new_header.samples = nx;									//calculate the width of the output image (concatenated subimages)
+		new_header.lines = nnz * ny;											//calculate the height of the output image (height of subimages)
+		
+
+		if (header.interleave == envi_header::BSQ){
+			if (header.data_type == envi_header::float32)
+				((bsq<float>*)file)->subimages(outfile, nx, ny, mask, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				((bsq<double>*)file)->subimages(outfile, nx, ny, mask, PROGRESS);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		else if (header.interleave == envi_header::BIL){
+			std::cout << "ERROR: unidentified data type" << std::endl;
+			exit(1);
+		}
+		else if (header.interleave == envi_header::BIP){
+			std::cout << "ERROR: unidentified data type" << std::endl;
+			exit(1);
+		}
+
+		new_header.save(outfile + ".hdr");									//save the header for the output file
+	}
+
 	/// Remove a list of bands from the ENVI file
 
 	/// @param outfile is the file name for the output hyperspectral image (with trimmed bands)
--
libgit2 0.21.4