diff --git a/matlab/fd_coefficients.m b/matlab/fd_coefficients.m new file mode 100644 index 0000000..863c5d8 --- /dev/null +++ b/matlab/fd_coefficients.m @@ -0,0 +1,22 @@ +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 diff --git a/matlab/sg_coefficients.m b/matlab/sg_coefficients.m new file mode 100644 index 0000000..86cfd71 --- /dev/null +++ b/matlab/sg_coefficients.m @@ -0,0 +1,19 @@ +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 diff --git a/stim/envi/bil.h b/stim/envi/bil.h index 0d7412a..d491310 100644 --- a/stim/envi/bil.h +++ b/stim/envi/bil.h @@ -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 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 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 w = std::vector(), 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 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 w_pts(w.begin() + iw, w.begin() + iw + N); //get the wavelengths corresponding to each sample + std::vector 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(){ diff --git a/stim/envi/bip.h b/stim/envi/bip.h index dbf01a4..8314052 100644 --- a/stim/envi/bip.h +++ b/stim/envi/bip.h @@ -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 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 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 w_pts(w.begin() + iw, w.begin() + iw + nC); //get the wavelengths corresponding to each sample + std::vector 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 w_pts(w.begin(), w.begin() + nC); //get the wavelengths corresponding to each sample + std::vector 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 w_pts(w.begin() + iw, w.begin() + iw + nC); //get the wavelengths corresponding to each sample + std::vector 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 w_pts(w.begin() + iw, w.begin() + iw + nC); //get the wavelengths corresponding to each sample + std::vector 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. diff --git a/stim/envi/bsq.h b/stim/envi/bsq.h index 82a6d6c..c131299 100644 --- a/stim/envi/bsq.h +++ b/stim/envi/bsq.h @@ -7,6 +7,7 @@ #include #include #include +#include @@ -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 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 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 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 > C, std::vector start, std::vector 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 w = std::vector(), 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 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 w_pts(w.begin() + iw, w.begin() + iw + N); //get the wavelengths corresponding to each sample + std::vector 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; + } } diff --git a/stim/envi/envi.h b/stim/envi/envi.h index f2add36..9706f6e 100644 --- a/stim/envi/envi.h +++ b/stim/envi/envi.h @@ -5,6 +5,7 @@ #include "../envi/bsq.h" #include "../envi/bip.h" #include "../envi/bil.h" +#include "../math/fd_coefficients.h" #include //#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 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*)file)->convolve(outfile, C, start, end, PROGRESS); + else if (header.data_type == envi_header::float64) + ((bsq*)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*)file)->convolve(outfile, C, start, end, PROGRESS); + else if (header.data_type == envi_header::float64) + ((bil*)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*)file)->convolve(outfile, C, start, end, PROGRESS); + else if (header.data_type == envi_header::float64) + ((bip*)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*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS); + else if (header.data_type == envi_header::float64) + ((bsq*)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*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS); + else if (header.data_type == envi_header::float64) + ((bil*)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*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS); + else if (header.data_type == envi_header::float64) + ((bip*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS); + else{ + std::cout << "ERROR: unidentified data type" << std::endl; + exit(1); + } + } + exit(1); + } }; } //end namespace rts diff --git a/stim/math/fd_coefficients.h b/stim/math/fd_coefficients.h new file mode 100644 index 0000000..764d0ce --- /dev/null +++ b/stim/math/fd_coefficients.h @@ -0,0 +1,67 @@ +#ifndef STIM_FINITE_DIFFERENCE_COEFFICIENTS +#define STIM_FINITE_DIFFERENCE_COEFFICIENTS + +#include +#include + +/// 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 +std::vector< std::vector > diff_coefficient_matrix(T Z, std::vector X, size_t M){ + size_t n = X.size(); //calculate the number of data points available + std::vector< std::vector > 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 +std::vector diff_coefficients(T x0, std::vector x, size_t n){ + std::vector< std::vector > C = diff_coefficient_matrix(x0, x, n); //calculate all derivative coefficients + return C.back(); +} + + + +#endif \ No newline at end of file -- libgit2 0.21.4