implemented finite differences and SG filters in stim::envi

David Mayerich
1 parent 2ce6954b
Showing 7 changed files with 521 additions and 0 deletions Show diff stats
matlab/fd_coefficients.m
matlab/sg_coefficients.m
stim/envi/bil.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/envi.h
stim/math/fd_coefficients.h
+function c = fd_coefficients(z,x,m)
+% Calculates FD weights. The parameters are:
+%  z   location where approximations are to be accurate,
+%  x   vector with x-coordinates for grid points,
+%  m   highest derivative that we want to find weights for
+%  c   array size m+1,lentgh(x) containing (as output) in 
+%      successive rows the weights for derivatives 0,1,...,m.
+
+n=length(x); c=zeros(m+1,n); c1=1; c4=x(1)-z; c(1,1)=1;
+for i=2:n
+   mn=min(i,m+1); c2=1; c5=c4; c4=x(i)-z;
+   for j=1:i-1
+      c3=x(i)-x(j);  c2=c2*c3;
+      if j==i-1 
+         c(2:mn,i)=c1*((1:mn-1)'.*c(1:mn-1,i-1)-c5*c(2:mn,i-1))/c2;
+         c(1,i)=-c1*c5*c(1,i-1)/c2;
+      end
+      c(2:mn,j)=(c4*c(2:mn,j)-(1:mn-1)'.*c(1:mn-1,j))/c3;
+      c(1,j)=c4*c(1,j)/c3;
+   end
+   c1=c2;
+end
 \ No newline at end of file
+function C = sg_coefficients(n, order)
+
+if(mod(n, 2) == 0)
+    disp('The number of coefficients must be odd');
+    return;
+end
+
+%assemble the column vector based on positions
+r = floor(n/2);         %maximum extent around zero (0)
+x = -r:1:r;
+
+%calculate J
+J = zeros(n, order + 1);
+
+%columns values are 1, x_i, x_i^2, ...
+for i = 1:order+1
+    J(:, i) = (x').^(i-1);
+end
+C = (J' * J)^(-1) * J';
 \ No newline at end of file
@@ -1174,7 +1174,106 @@ public:
 		}
 	}
  
+	/// Convolve the given band range with a kernel specified by a vector of coefficients.
  
+	/// @param outfile is an already open stream to the output file
+	/// @param C is an array of coefficients
+	/// @param start is the band to start processing (the first coefficient starts here)
+	/// @param nbands is the number of bands to process
+	/// @param center is the index for the center coefficient for the kernel (used to set the wavelengths in the output file)
+
+	void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, bool PROGRESS = false){
+		std::ofstream out(outfile.c_str(), std::ios::binary);	//open the output file for writing
+		size_t B = end - start + 1;
+		size_t Xb = X() * sizeof(T);								//calculate the size of a band (frame) in bytes
+		size_t XBb = Xb * Z();
+
+		size_t N = C.size();										//get the number of bands that the kernel spans
+		std::deque<T*> frame(N, NULL);								//create an array to store pointers to each frame
+		for(size_t f = 0; f < N; f++)
+			frame[f] = (T*)malloc(Xb);								//allocate space for the frame		
+
+		T* outline = (T*)malloc(Xb);								//allocate space for the output band
+
+		//Implementation: In order to minimize reads from secondary storage, each band is only loaded once into the 'frame' deque.
+		//					When a new band is loaded, the last band is popped, a new frame is copied to the pointer, and it is
+		//					re-inserted into the deque.
+		for(size_t y = 0; y < Y(); y++){
+			file.seekg(y * XBb + start * Xb, std::ios::beg);		//move to the beginning of the 'start' band
+			for(size_t f = 0; f < N; f++)								//for each frame
+				file.read((char*)frame[f], Xb);							//load the frame
+			for(size_t b = 0; b < B; b++){								//for each band
+				memset(outline, 0, Xb);									//set the output band to zero (0)
+				for(size_t c = 0; c < N; c++){							//for each frame (corresponding to each coefficient)
+					for(size_t x = 0; x < X(); x++){					//for each pixel
+						outline[x] += (T)(C[c] * frame[c][x]);			//calculate the contribution of the current frame (scaled by the corresponding coefficient)
+					}
+				}
+				out.write((char*)outline, Xb);							//output the band
+				file.read((char*)frame[0], Xb);							//read the next band
+				frame.push_back(frame.front());							//put the first element in the back
+				frame.pop_front();										//pop the first element
+			}
+			if(PROGRESS) progress = (double)(y+1) / (double)Y() * 100;
+		}		
+	}
+
+	/// Approximate the spectral derivative of the image
+
+	void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w = std::vector<double>(), bool PROGRESS = false){
+		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file for writing
+				
+		size_t Xb = X() * sizeof(T);								//calculate the size of a line (frame) in bytes
+		size_t XBb = Xb * Z();
+		size_t B = Z();
+
+		file.seekg(0, std::ios::beg);								//seek to the beginning of the file		
+
+		size_t N = order + d;										//approximating a derivative requires order + d samples
+		std::deque<T*> frame(N, NULL);								//create an array to store pointers to each frame
+		for(size_t f = 0; f < N; f++)								//for each frame
+			frame[f] = (T*)malloc(Xb);								//allocate space for the frame	
+
+		T* outline = (T*)malloc(Xb);								//allocate space for the output band
+
+		//Implementation: In order to minimize reads from secondary storage, each band is only loaded once into the 'frame' deque.
+		//					When a new band is loaded, the last band is popped, a new frame is copied to the pointer, and it is
+		//					re-inserted into the deque.
+		size_t mid = (size_t)(N / 2);								//calculate the mid point of the kernel
+		size_t iw;													//index to the first wavelength used to evaluate the derivative at this band
+
+		for(size_t y = 0; y < Y(); y++){
+			//file.seekg(y * XBb, std::ios::beg);							//seek to the beginning of the current Y plane
+			for(size_t f = 0; f < N; f++)								//for each frame
+				file.read((char*)frame[f], Xb);							//load a line
+
+			for(size_t b = 0; b < B; b++){								//for each band
+				if(b < mid)												//calculate the first wavelength used to evaluate the derivative at this band
+					iw = 0;
+				else if(b > B - (N - mid + 1))
+					iw = B - N;
+				else{
+					iw = b - mid;
+					file.read((char*)frame[0], Xb);											//read the next band
+					frame.push_back(frame.front());											//put the first element in the back
+					frame.pop_front();														//pop the first element
+				}
+				std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + N);			//get the wavelengths corresponding to each sample
+				std::vector<double> C = diff_coefficients(w[b], w_pts, d);				//get the optimal sample weights
+
+				memset(outline, 0, Xb);													//set the output band to zero (0)
+				for(size_t c = 0; c < N; c++){											//for each frame (corresponding to each coefficient)
+					for(size_t x = 0; x < X(); x++){									//for each pixel
+						outline[x] += (T)(C[c] * frame[c][x]);							//calculate the contribution of the current frame (scaled by the corresponding coefficient)
+					}
+				}
+				out.write((char*)outline, Xb);											//output the band
+				
+			}
+			if(PROGRESS) progress = (double)(y+1) / (double)Y() * 100;
+		}
+
+	}
  
 	/// Close the file.
 	bool close(){
@@ -1250,6 +1250,99 @@ public:
 		}
 	}
  
+	/// Convolve the given band range with a kernel specified by a vector of coefficients.
+
+	/// @param outfile is an already open stream to the output file
+	/// @param C is an array of coefficients
+	/// @param start is the band to start processing (the first coefficient starts here)
+	/// @param nbands is the number of bands to process
+	/// @param center is the index for the center coefficient for the kernel (used to set the wavelengths in the output file)
+
+	void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, bool PROGRESS = false){
+		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file for writing
+		
+		size_t N = end - start + 1;									//number of bands in the output spectrum
+		size_t Nb = N * sizeof(T);									//size of the output spectrum in bytes
+		size_t B = Z();												//calculate the number of values in a spectrum
+		size_t Bb = B * sizeof(T);									//calculate the size of a spectrum in bytes
+
+		file.seekg(0, std::ios::beg);								//move to the beginning of the input file
+
+		size_t nC = C.size();										//get the number of bands that the kernel spans
+		T* inspec = (T*)malloc(Bb);									//allocate space for the input spectrum
+		T* outspec = (T*)malloc(Nb);								//allocate space for the output spectrum
+
+		size_t XY = X() * Y();										//number of pixels in the image
+		for(size_t xy = 0; xy < XY; xy++){							//for each pixel
+			file.read((char*)inspec, Bb);							//read an input spectrum
+			memset(outspec, 0, Nb);									//set the output spectrum to zero (0)
+			for(size_t b = 0; b < N; b++){							//for each component of the spectrum
+				for(size_t c = 0; c < nC; c++){						//for each coefficient in the kernel
+					outspec[b] += (T)(inspec[b + start + c] * C[c]);		//perform the sum/multiply part of the convolution
+				}
+			}
+			out.write((char*)outspec, Nb);							//output the band
+			if(PROGRESS) progress = (double)(xy+1) / (double)XY * 100;
+		}
+	}
+
+	void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w, bool PROGRESS = false){
+		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file for writing
+		
+		
+		size_t B = Z();												//calculate the number of values in a spectrum
+		size_t Bb = B * sizeof(T);									//calculate the size of a spectrum in bytes
+
+		bool UNIFORM = true;
+		double ds = w[1] - w[0];									//initialize ds
+		for(size_t b = 1; b < B; b++)								//test to see if the spectral spacing is uniform (if it is, convolution is much faster)
+			if(w[b] - w[b-1] != ds) UNIFORM = false;
+
+		size_t nC = order + d;									//approximating a derivative requires order + d samples
+
+		file.seekg(0, std::ios::beg);								//move to the beginning of the input file
+
+		T* inspec = (T*)malloc(Bb);									//allocate space for the input spectrum
+		T* outspec = (T*)malloc(Bb);								//allocate space for the output spectrum
+
+		size_t XY = X() * Y();										//number of pixels in the image
+		size_t mid = (size_t)(nC / 2);							//calculate the mid point of the kernel
+		size_t iw;													//index to the first wavelength used to evaluate the derivative at this band
+		for(size_t xy = 0; xy < XY; xy++){							//for each pixel
+			file.read((char*)inspec, Bb);							//read an input spectrum
+			memset(outspec, 0, Bb);									//set the output spectrum to zero (0)
+
+			iw = 0;
+			for(size_t b = 0; b < mid; b++){							//for each component of the spectrum
+				std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + nC);			//get the wavelengths corresponding to each sample
+				std::vector<double> C = diff_coefficients(w[b], w_pts, d);					//get the optimal sample weights
+				for(size_t c = 0; c < nC; c++)						//for each coefficient in the kernel
+					outspec[b] += (T)(inspec[iw + c] * C[c]);		//perform the sum/multiply part of the convolution
+			}
+			std::vector<double> w_pts(w.begin(), w.begin() + nC);			//get the wavelengths corresponding to each sample
+			std::vector<double> C = diff_coefficients(w[0], w_pts, d);					//get the optimal sample weights
+			for(size_t b = mid; b <= B - (nC - mid); b++){
+				iw = b - mid;
+				if(!UNIFORM){																//if the spacing is non-uniform, we have to re-calculate these points every iteration
+					std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + nC);			//get the wavelengths corresponding to each sample
+					std::vector<double> C = diff_coefficients(w[b], w_pts, d);					//get the optimal sample weights
+				}
+				for(size_t c = 0; c < nC; c++)						//for each coefficient in the kernel
+					outspec[b] += (T)(inspec[iw + c] * C[c]);		//perform the sum/multiply part of the convolution
+			}
+			iw = B - nC;
+			for(size_t b = B - (nC - mid) + 1; b < B; b++){
+				std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + nC);			//get the wavelengths corresponding to each sample
+				std::vector<double> C = diff_coefficients(w[b], w_pts, d);					//get the optimal sample weights
+				for(size_t c = 0; c < nC; c++)						//for each coefficient in the kernel
+					outspec[b] += (T)(inspec[iw + c] * C[c]);		//perform the sum/multiply part of the convolution
+			}
+
+			out.write((char*)outspec, Bb);							//output the band
+			if(PROGRESS) progress = (double)(xy+1) / (double)XY * 100;
+		}
+	}
+
  
  
 	/// Close the file.
@@ -7,6 +7,7 @@
 #include <cstring>
 #include <utility>
 #include <vector>
+#include <deque>
  
  
  
@@ -1058,6 +1059,128 @@ public:
 			out.write((char*)dst, band_dst_bytes);														//write the combined image to an output file
 			if(PROGRESS) progress = (double)(b + 1) / (double) Z() * 100;
 		}
+		out.close();
+	}
+
+	/// Convolve the given band range with a kernel specified by a vector of coefficients.
+
+	/// @param outfile is an already open stream to the output file
+	/// @param C is an array of coefficients
+	/// @param start is the band to start processing (the first coefficient starts here)
+	/// @param nbands is the number of bands to process
+	/// @param center is the index for the center coefficient for the kernel (used to set the wavelengths in the output file)
+
+	void convolve(std::ofstream& out, std::vector<double> C, size_t start, size_t end, bool PROGRESS = false){
+		size_t nbands = end - start + 1;
+		size_t XY = X() * Y();										//calculate the number of values in a band
+		size_t XYb = XY * sizeof(T);								//calculate the size of a band (frame) in bytes
+
+		file.seekg(XYb * start, std::ios::beg);						//move to the beginning of the 'start' band
+
+		size_t nframes = C.size();									//get the number of bands that the kernel spans
+		std::deque<T*> frame(nframes, NULL);						//create an array to store pointers to each frame
+		for(size_t f = 0; f < nframes; f++){						//for each frame
+			frame[f] = (T*)malloc(XYb);								//allocate space for the frame
+			file.read((char*)frame[f], XYb);						//load the frame
+		}
+
+		T* outband = (T*)malloc(XYb);								//allocate space for the output band
+
+		//Implementation: In order to minimize reads from secondary storage, each band is only loaded once into the 'frame' deque.
+		//					When a new band is loaded, the last band is popped, a new frame is copied to the pointer, and it is
+		//					re-inserted into the deque.
+		for(size_t b = 0; b < nbands; b++){											//for each band
+			memset(outband, 0, XYb);												//set the output band to zero (0)
+			for(size_t c = 0; c < nframes; c++){									//for each frame (corresponding to each coefficient)
+				for(size_t xy = 0; xy < XY; xy++){									//for each pixel
+					outband[xy] += (T)(C[c] * frame[c][xy]);				//calculate the contribution of the current frame (scaled by the corresponding coefficient)
+				}
+			}
+			out.write((char*)outband, XYb);							//output the band
+			file.read((char*)frame[0], XYb);						//read the next band
+			frame.push_back(frame.front());							//put the first element in the back
+			frame.pop_front();										//pop the first element
+			if(PROGRESS) progress = (double)(b+1) / (double)nbands * 100;
+		}
+	}
+
+	/// Performs a single convolution and saves it to an output file
+
+	/// @param outfile is the convolved file to be output
+	/// @param C is an array of coefficients
+	/// @param start is the band to start processing (the first coefficient starts here)
+	/// @param nbands is the number of bands to process
+	void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, bool PROGRESS = false){
+		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file for writing
+		convolve(out, C, start, end, PROGRESS);						//start the convolution
+		out.close();
+	}
+
+	/// Performs a set of convolutions and chains the results together in a single file
+
+	/// @param outfile is the convolved file to be output
+	/// @param C is an array containing an array of coefficients for each kernel
+	/// @param start is the list of start bands for each kernel
+	/// @param end is the list of end bands for each kernel
+	void convolve(std::string outfile, std::vector< std::vector<double> > C, std::vector<size_t> start, std::vector<size_t> end, bool PROGRESS = false){
+		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file for writing
+
+		size_t K = C.size();										//get the number of kernels
+		for(size_t k = 0; k < K; k++){
+			size_t b0 = start[k];									//calculate the range of the convolution
+			size_t b1 = end[k];
+			convolve(out, C[k], b0, b1, PROGRESS);					//perform the convolution with the current kernel in the given range
+		}
+		out.close();
+	}
+
+	/// Approximate the spectral derivative of the image
+
+	void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w = std::vector<double>(), bool PROGRESS = false){
+		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file for writing
+				
+		size_t XY = X() * Y();										//calculate the number of values in a band
+		size_t XYb = XY * sizeof(T);								//calculate the size of a band (frame) in bytes
+		size_t B = Z();
+		file.seekg(0, std::ios::beg);								//move to the beginning of the file
+
+		size_t N = order + d;									//approximating a derivative requires order + d samples
+		std::deque<T*> frame(N, NULL);						//create an array to store pointers to each frame
+		for(size_t f = 0; f < N; f++){						//for each frame
+			frame[f] = (T*)malloc(XYb);								//allocate space for the frame
+			file.read((char*)frame[f], XYb);						//load the frame
+		}
+
+		T* outband = (T*)malloc(XYb);								//allocate space for the output band
+
+		//Implementation: In order to minimize reads from secondary storage, each band is only loaded once into the 'frame' deque.
+		//					When a new band is loaded, the last band is popped, a new frame is copied to the pointer, and it is
+		//					re-inserted into the deque.
+		size_t mid = (size_t)(N / 2);							//calculate the mid point of the kernel
+		size_t iw;													//index to the first wavelength used to evaluate the derivative at this band
+		for(size_t b = 0; b < B; b++){								//for each band
+			if(b < mid)												//calculate the first wavelength used to evaluate the derivative at this band
+				iw = 0;
+			else if(b > B - (N - mid + 1))
+				iw = B - N;
+			else{
+				iw = b - mid;
+				file.read((char*)frame[0], XYb);						//read the next band
+				frame.push_back(frame.front());							//put the first element in the back
+				frame.pop_front();										//pop the first element
+			}
+			std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + N);			//get the wavelengths corresponding to each sample
+			std::vector<double> C = diff_coefficients(w[b], w_pts, d);					//get the optimal sample weights
+
+			memset(outband, 0, XYb);												//set the output band to zero (0)
+			for(size_t c = 0; c < N; c++){									//for each frame (corresponding to each coefficient)
+				for(size_t xy = 0; xy < XY; xy++){									//for each pixel
+					outband[xy] += (T)(C[c] * frame[c][xy]);				//calculate the contribution of the current frame (scaled by the corresponding coefficient)
+				}
+			}
+			out.write((char*)outband, XYb);							//output the band			
+			if(PROGRESS) progress = (double)(b+1) / (double)B * 100;
+		}
  
 	}
  
@@ -5,6 +5,7 @@
 #include "../envi/bsq.h"
 #include "../envi/bip.h"
 #include "../envi/bil.h"
+#include "../math/fd_coefficients.h"
 #include <iostream>
 //#include "../image/image.h"
  
@@ -1341,6 +1342,103 @@ public:
 			}
 		}
 	}
+
+	/// Convolve the given band range with a kernel specified by a vector of coefficients.
+
+	/// @param outfile is the combined file to be output
+	/// @param c is an array of coefficients
+	/// @param start is the band to start processing (the first coefficient starts here)
+	/// @param nbands is the number of bands to process
+	/// @param center is the index for the center coefficient for the kernel (used to set the wavelengths in the output file)
+	void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, size_t center = 0, bool PROGRESS = false){
+		size_t nbands = end - start + 1;
+		envi_header h = header;												//copy the current header
+		h.bands = nbands;													//set the number of new bands
+		if(header.wavelength.size() != 0){
+			h.wavelength.resize(nbands);									//set the number of wavelengths to the number of bands
+			for(size_t b = 0; b < nbands; b++)
+				h.wavelength[b] = header.wavelength[b+center];
+		}
+		if(header.band_names.size() != 0){
+			h.band_names.resize(nbands);
+			for(size_t b = 0; b < nbands; b++)
+				h.band_names[b] = header.band_names[b+center];
+		}
+		h.save(outfile + ".hdr");											//save the new header
+
+		if (header.interleave == envi_header::BSQ){
+			if (header.data_type == envi_header::float32)
+				((bsq<float>*)file)->convolve(outfile, C, start, end, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				((bsq<double>*)file)->convolve(outfile, C, start, end, PROGRESS);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		else if (header.interleave == envi_header::BIL){
+			if (header.data_type == envi_header::float32)
+				((bil<float>*)file)->convolve(outfile, C, start, end, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				((bil<double>*)file)->convolve(outfile, C, start, end, PROGRESS);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		else if (header.interleave == envi_header::BIP){
+			if (header.data_type == envi_header::float32)
+				((bip<float>*)file)->convolve(outfile, C, start, end, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				((bip<double>*)file)->convolve(outfile, C, start, end, PROGRESS);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+	}
+
+	/// Approximates the nth derivative of the spectra to the specified order
+
+	/// @param outfile is the file where the derivative approximation will be saved
+	/// @n is the derivative to be calculated
+	/// @order is the order of the error (must be even)
+	void deriv(std::string outfile, size_t d, size_t order, bool PROGRESS = false){
+		header.save(outfile + ".hdr");
+		if (header.interleave == envi_header::BSQ){
+			if (header.data_type == envi_header::float32)
+				((bsq<float>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				((bsq<double>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+
+		else if (header.interleave == envi_header::BIL){
+			if (header.data_type == envi_header::float32)
+				((bil<float>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				((bil<double>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+
+		else if (header.interleave == envi_header::BIP){
+			if (header.data_type == envi_header::float32)
+				((bip<float>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				((bip<double>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		exit(1);
+	}
 };
  
 }	//end namespace rts
+#ifndef STIM_FINITE_DIFFERENCE_COEFFICIENTS
+#define STIM_FINITE_DIFFERENCE_COEFFICIENTS
+
+#include <vector>
+#include <algorithm>
+
+/// Calculate the coefficients used for finite differences on arbitrary grids. This code is based on:
+///	Fornberg, "Generation of Finite Difference Formulas on Arbitrarily Spaced Grids", Mathematics of Computation, 51(184) 699-706, 1988
+
+/// @param Z is the location where the approximation is to be accurate
+/// @param X is a vector of coordinates for grid points
+/// @param M is the highest derivative to be computed
+/// This function returns a matrix C[a][b] where a is the derivative and b is the coefficient for x = b
+
+template <typename T>
+std::vector< std::vector<T> > diff_coefficient_matrix(T Z, std::vector<T> X, size_t M){
+	size_t n = X.size();								//calculate the number of data points available
+	std::vector< std::vector<T> > C(M+1);							//create an array of m arrays
+	for(size_t m = 0; m < M+1; m++)						//for each derivative
+		C[m].resize(n);									//allocate space for n coefficients
+	
+	T c1, c2, c3, c4, c5;
+	c1 = 1;												//initialize the matrix
+	c4 = X[0] - Z;
+	C[0][0] = 1;
+	for(size_t i = 2; i <= n; i++){
+		size_t mn = (std::min)(i, M+1);
+		c2 = 1;
+		c5 = c4;
+		c4 = X[i-1] - Z;
+		for(size_t j = 1; j <= i-1; j++){
+			c3 = X[i-1] - X[j-1];
+			c2 = c2 * c3;
+			if(j == i-1){
+				for(size_t k = 2; k <= mn; k++){
+					C[k-1][i-1] = c1 * ((k-1) * C[k-2][i-2] - c5 * C[k-1][i-2]) / c2;
+				}
+				C[0][i-1] = -c1 * c5 * C[0][i-2]/c2;
+			}
+			T c_p0 = C[0][j-1];
+			T c_p1 = C[1][j-1];
+			for(size_t k = 2; k <= mn; k++){
+				c_p1 = C[k-1][j-1];
+				C[k-1][j-1] = (c4 * c_p1 - (k-1) * c_p0) / c3;
+				c_p0 = c_p1;
+			}
+			C[0][j-1] = c4 * C[0][j-1] / c3;
+		}
+		c1 = c2;
+	}
+	return C;
+}
+
+/// Returns the coefficients used to estimate the nth derivative at x0 given values at points in x
+
+/// @param x0 is the point where the derivative f'(x0) is to be estimated
+/// @param x is the set of points where the values of f(x) are known
+/// @param n specifies the derivative, ex. 2 estimates the second derivative f''(x0)
+template <typename T>
+std::vector<T> diff_coefficients(T x0, std::vector<T> x, size_t n){
+	std::vector< std::vector<T> > C = diff_coefficient_matrix(x0, x, n);		//calculate all derivative coefficients
+	return C.back();
+}
+
+
+
+#endif
 \ No newline at end of file