27 #include "../common/MatrixUtils.h"
28 #include "../common/progress/TaskProgress.h"
29 #include "../raster/Grid.h"
30 #include "../raster/Utils.h"
37 #include <boost/numeric/ublas/io.hpp>
38 #include <boost/numeric/ublas/matrix.hpp>
40 bool gaussElimination(boost::numeric::ublas::matrix<double>& Z, std::vector<int>& row,
unsigned int nComponents);
41 bool fowardBackSubstitution(boost::numeric::ublas::matrix<double>& Z, std::vector<double>& y,
unsigned int ys, std::vector<int>& row,
unsigned int rows, std::vector<double>& x,
unsigned int xs);
101 const std::vector<std::string>& inputSensorBands,
const std::map<std::string, std::vector<double> >& components,
102 te::rst::Raster& outputRaster,
const bool enableProgressInterface)
throw(te::rp::Exception)
107 const unsigned int nComponents = components.size();
108 const unsigned int nBands = inputSensorBands.size();
110 bool computeErrorRasters =
false;
111 if (nComponents < outputRaster.getNumberOfBands())
112 computeErrorRasters =
true;
115 unsigned int i, j, k, m;
117 std::vector<double> Qmax;
118 for (i = 0; i < inputSensorBands.size(); i++)
122 boost::numeric::ublas::matrix<double> SpectralBandsComponents(nBands, nComponents);
123 std::map<std::string, std::vector<double> >::const_iterator it;
124 for(i = 0; i < nBands; i++)
125 for (j = 0, it = components.begin(); it != components.end(); j++, it++)
126 SpectralBandsComponents(i, j) = it->second[i] / Qmax[i];
131 boost::numeric::ublas::matrix<double> SpectralBandsComponentsTransposed(nComponents, nBands);
134 SpectralBandsComponentsTransposed = boost::numeric::ublas::trans(SpectralBandsComponents);
138 std::vector<double> componentsMean(nBands, 0.0);
139 std::vector<double> componentsMeanAfter(nBands, 0.0);
141 for (j = 0; j < nBands; j++)
144 componentsMean[j] = 0;
145 for(i = 0; i < nComponents; i++)
146 componentsMean[j] = componentsMean[j] + SpectralBandsComponentsTransposed(i, j);
147 componentsMean[j] = componentsMean[j] / (double) nComponents;
150 componentsMeanAfter[j] = 0;
151 for(i = 0; i < nComponents; i++)
153 SpectralBandsComponentsTransposed(i, j) = SpectralBandsComponentsTransposed(i, j) - componentsMean[j];
154 componentsMeanAfter[j] = componentsMeanAfter[j] + SpectralBandsComponentsTransposed(i, j);
156 componentsMeanAfter[j] = componentsMeanAfter[j] / (double) nComponents;
160 std::vector<double> covarianceVector(nBands * nBands, 0.0);
161 for(k = 0; k < nBands * nBands; k++)
162 covarianceVector[k] = 0;
165 for(j = 0; j < nBands; j++)
166 for(m = 0; m < j + 1; m++)
168 for(i = 0; i < nComponents; i++)
169 covarianceVector[k] = covarianceVector[k] + (SpectralBandsComponentsTransposed(i, j) - componentsMeanAfter[j]) *
170 (SpectralBandsComponentsTransposed(i, m) - componentsMeanAfter[m]);
171 covarianceVector[k] = covarianceVector[k] / (double) nComponents;
177 boost::numeric::ublas::matrix<double> covarianceMatrix(nBands, nBands);
178 for (i = 0; i < nBands; i++)
179 for (j = 0; j <= i; j++)
180 covarianceMatrix(i, j) = covarianceVector[k++];
183 boost::numeric::ublas::matrix<double> auxMatrix(covarianceMatrix);
184 boost::numeric::ublas::matrix<double> tmpMatrix;
190 boost::numeric::ublas::matrix<double> eigenReduced(nBands, nComponents - 1);
191 for (j = 0; j < nBands; j++)
192 for (i = 0; i < nComponents - 1; i++)
193 eigenReduced(j, i) = auxMatrix(j, i);
196 boost::numeric::ublas::matrix<double> SpectralBandsComponentsPCA(SpectralBandsComponentsTransposed.size1(), eigenReduced.size2());
197 SpectralBandsComponentsPCA = boost::numeric::ublas::prod(SpectralBandsComponentsTransposed, eigenReduced);
200 auxMatrix.resize(SpectralBandsComponentsPCA.size1(), SpectralBandsComponentsPCA.size2() + 1);
203 for (j = 0; j < SpectralBandsComponentsPCA.size1(); j++)
205 for (i = 0; i < SpectralBandsComponentsPCA.size2(); i++)
206 auxMatrix(j, i) = SpectralBandsComponentsPCA(j, i);
208 auxMatrix(j, SpectralBandsComponentsPCA.size2()) = 1.0;
212 SpectralBandsComponentsPCA.clear();
213 SpectralBandsComponentsPCA = boost::numeric::ublas::trans(auxMatrix);
216 std::vector<int> rows(nComponents, 0);
224 auxMatrix.resize(1, nBands);
226 boost::numeric::ublas::matrix<double> ymat(nComponents - 1, nComponents - 1);
229 std::vector<double> x(nComponents, 0.0);
230 std::vector<double> y(nComponents, 0.0);
233 std::vector<unsigned int> rowsIn(nBands, 0);
234 std::vector<unsigned int> columnIn(nBands, 0);
235 std::vector<unsigned int> columnsIn(nBands, 0);
236 std::vector<double> spectralBandsError(nBands, 0.0);
238 boost::numeric::ublas::matrix<double> rasterRowsIn(nBands, inputRaster.getNumberOfColumns());
239 boost::numeric::ublas::matrix<double> rasterRowsOut(nComponents, inputRaster.getNumberOfColumns());
243 for (
unsigned int rowout = 0; rowout < inputRaster.getNumberOfRows(); rowout++)
246 for (i = 0; i < nBands; i++)
248 for (j = 0; j < inputRaster.getNumberOfColumns(); j++)
250 inputRaster.getValue(j, rowsIn[i], value, inputRasterBands[i]);
251 rasterRowsIn(i, j) = value;
258 columnIn[i] = columnsIn[i];
262 for (
unsigned int colout = 0; colout < inputRaster.getNumberOfColumns(); colout++)
265 for (i = 0; i < nBands; i++)
266 auxMatrix(0, i) = (double) (rasterRowsIn(i, columnIn[i])) / 255.0 - componentsMean[i];
268 ymat = boost::numeric::ublas::prod(auxMatrix, eigenReduced);
270 for (i = 0; i < nComponents - 1; i++)
274 if (!
fowardBackSubstitution(SpectralBandsComponentsPCA, y, nComponents, rows, nComponents, x, nComponents))
278 for (i = 0; i < nComponents; i++)
280 prop = (short) (x[i] * 100 + 100);
281 rasterRowsOut(i, colout) = (
unsigned char) prop;
288 if (computeErrorRasters)
290 for (i = 0; i < nBands; i++)
294 for (j = 0; j < nComponents; j++)
295 aux += x[j] * SpectralBandsComponents(i, j);
298 error = (long) (rasterRowsIn(i, columnIn[i])) / 255.0 - aux;
302 if (computeErrorRasters)
303 rasterRowsIn(i, colout) = error * 255.0;
308 for (i = 0; i < nBands; i++)
314 for (i = 0; i < nComponents; i++)
315 for (j = 0; j < inputRaster.getNumberOfColumns(); j++)
316 outputRaster.setValue(j, rowout, rasterRowsOut(i, j), i);
319 if (computeErrorRasters)
321 for (i = 0; i < nComponents; i++)
322 for (j = 0; j < inputRaster.getNumberOfColumns(); j++)
323 outputRaster.setValue(j, rowout, rasterRowsIn(i, j), i + nComponents);
345 bool gaussElimination(boost::numeric::ublas::matrix<double>& Z, std::vector<int>& row,
unsigned int nComponents)
348 assert(Z.size1() >= Z.size2());
349 unsigned int nrows = Z.size1();
350 assert(nComponents >= nrows);
360 for (i = 0; i < nrows; i++)
364 for (k = 0; k < nrows; k++)
367 while (Z(i, k) == 0.0)
383 for (j=0; j < nrows; j ++)
392 for (i = k + 1; i < nrows; i++)
396 for (j = k + 1; j < nrows; j++)
397 Z (i, j) = Z (i, j) - m * Z(k, j);
402 if (Z(nrows - 1, nrows - 1) == 0.0)
408 bool fowardBackSubstitution(boost::numeric::ublas::matrix<double>& Z, std::vector<double>& y,
unsigned int ys, std::vector<int>& row,
unsigned int rows, std::vector<double>& x,
unsigned int xs)
411 assert(Z.size1() >= Z.size2());
412 unsigned int nrows = Z.size1();
414 if (((rows < nrows) || (ys < nrows) || (xs < nrows)))
419 std::vector<double> F(nrows, 0.0);
422 for (k = 0; k < (int) nrows; k++)
425 for (j = 0; j <= k - 1; j++)
426 aux = aux + Z(k, j) * F[j];
428 F[k] = y[row[k]] - aux;
432 for (k = ((
int) nrows) - 1; k >= 0; k--)
435 for (j = k + 1; j < (int) nrows; j++)
436 aux = aux + Z(k, j) * x[j];
438 x[k] = (F[k] - aux) / Z(k, k);
double GetDigitalNumberBandMax(std::string bandName)
Returns the maximum digital number of a given sensor/band.
Raster Mixture model strategy factory base class.
MixtureModelPCAStrategyFactory()
bool fowardBackSubstitution(boost::numeric::ublas::matrix< double > &Z, std::vector< double > &y, unsigned int ys, std::vector< int > &row, unsigned int rows, std::vector< double > &x, unsigned int xs)
bool initialize(StrategyParameters const *const strategyParams)
Initialize the segmentation strategy.
This class can be used to inform the progress of a task.
MixtureModelPCAStrategy()
void reset()
Clear all internal allocated resources and reset the parameters instance to its initial state...
#define TE_TR(message)
It marks a string in order to get translated.
~MixtureModelPCAStrategyFactory()
MixtureModelPCAStrategy::Parameters m_parameters
Internal execution parameters.
Raster Processing functions.
#define TERP_TRUE_OR_RETURN_FALSE(value, message)
Checks if value is true. For false values a warning message will be logged and a return of context wi...
Raster mixture model strategy base class.
An abstract class for raster data strucutures.
const Parameters & operator=(const Parameters ¶ms)
Raster strategy parameters base class.
Raster PCA mixture model strategy factory.
~MixtureModelPCAStrategy()
bool EigenVectors(const boost::numeric::ublas::matrix< T > &inputMatrix, boost::numeric::ublas::matrix< T > &eigenVectorsMatrix, boost::numeric::ublas::matrix< T > &eigenValuesMatrix)
Computes the eigenvectors of a given matrix.
Abstract parameters base interface.
bool gaussElimination(boost::numeric::ublas::matrix< double > &Z, std::vector< int > &row, unsigned int nComponents)
AbstractParameters * clone() const
Create a clone copy of this instance.
bool m_isInitialized
True if this instance is initialized.
bool execute(const te::rst::Raster &inputRaster, const std::vector< unsigned int > &inputRasterBands, const std::vector< std::string > &inputSensorBands, const std::map< std::string, std::vector< double > > &components, te::rst::Raster &outputRaster, const bool enableProgressInterface)
Executes the segmentation strategy.
PCA (Principal Component Analysis) strategy for mixture model.
te::rp::MixtureModelStrategy * build()
Concrete factories (derived from this one) must implement this method in order to create objects...