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>
96 const std::vector<std::string>& inputSensorBands,
const std::map<std::string, std::vector<double> >& components,
97 te::rst::Raster& outputRaster,
const bool enableProgressInterface)
throw(te::rp::Exception)
101 const unsigned int nComponents = components.size();
102 const unsigned int nBands = inputRasterBands.size();
104 boost::numeric::ublas::matrix<double> transposeA = boost::numeric::ublas::matrix<double>(nComponents + 1, nBands + 1);
105 boost::numeric::ublas::matrix<double> matrixA (transposeA.size2(), transposeA.size1());
112 std::vector<double> Qmax;
113 for (
unsigned i = 0; i < inputSensorBands.size(); i++)
118 for (
unsigned j = 0; j < nBands; j++)
119 transposeA(row, j) = 1.0;
120 transposeA(row, nBands) = 0.0;
121 std::map<std::string, std::vector<double> >::const_iterator it;
122 for (row = 1, it = components.begin(); it != components.end(); row++, it++)
124 for (
unsigned j = 0; j < nBands; j++)
125 transposeA(row, j) = it->second[j] / Qmax[j];
126 transposeA(row, nBands) = 1.0;
130 matrixA = boost::numeric::ublas::trans(transposeA);
133 boost::numeric::ublas::matrix<double> productAtA = boost::numeric::ublas::matrix<double>(transposeA.size1(), matrixA.size2());
134 productAtA = boost::numeric::ublas::prod(transposeA, matrixA);
137 boost::numeric::ublas::matrix<double> invertAtA = boost::numeric::ublas::matrix<double>(productAtA.size1(), productAtA.size2());
140 boost::numeric::ublas::matrix<double> matrixR = boost::numeric::ublas::matrix<double>(matrixA.size1(), 1);
141 boost::numeric::ublas::matrix<double> productAtR = boost::numeric::ublas::matrix<double>(transposeA.size1(), matrixR.size2());
142 boost::numeric::ublas::matrix<double> matrixX = boost::numeric::ublas::matrix<double>(invertAtA.size1(), productAtR.size2());
143 boost::numeric::ublas::matrix<double> productAX = boost::numeric::ublas::matrix<double>(matrixA.size1(), matrixX.size2());
144 boost::numeric::ublas::matrix<double> matrixE = boost::numeric::ublas::matrix<double>(matrixR.size1(), matrixR.size2());
148 for (
unsigned r = 0; r < inputRaster.getNumberOfRows(); r++)
150 for (
unsigned c = 0; c < inputRaster.getNumberOfColumns(); c++)
153 for (
unsigned b = 0; b < nBands; b++)
155 inputRaster.getValue(c, r, value, inputRasterBands[b]);
156 matrixR(b, 0) = value / Qmax[b];
158 matrixR(nBands, 0) = 0.0;
160 productAtR = boost::numeric::ublas::prod(transposeA, matrixR);
162 matrixX = boost::numeric::ublas::prod(invertAtA, productAtR);
164 for (
unsigned b = 0; b < nComponents; b++)
165 outputRaster.setValue(c, r, matrixX(b + 1, 0), b);
168 if (nComponents < outputRaster.getNumberOfBands())
170 productAX = boost::numeric::ublas::prod(matrixA, matrixX);
171 matrixE = matrixR - productAX;
174 for (
unsigned b = nComponents, e = 0; b < outputRaster.getNumberOfBands(); b++, e++)
175 outputRaster.setValue(c, r, matrixE(e, 0), b);
AbstractParameters * clone() const
Create a clone copy of this instance.
bool initialize(StrategyParameters const *const strategyParams)
Initialize the segmentation strategy.
MixtureModelLinearStrategy::Parameters m_parameters
Internal execution parameters.
Raster strategy parameters base class.
MixtureModelLinearStrategy()
Raster Processing functions.
te::rp::MixtureModelStrategy * build()
Concrete factories (derived from this one) must implement this method in order to create objects...
~MixtureModelLinearStrategyFactory()
Raster mixture model strategy base class.
double GetDigitalNumberBandMax(std::string bandName)
Returns the maximum digital number of a given sensor/band.
void reset()
Clear all internal allocated resources and reset the parameters instance to its initial state...
#define TR_RP(message)
It marks a string in order to get translated. This is a special mark used in the Raster Processing mo...
~MixtureModelLinearStrategy()
MixtureModelLinearStrategyFactory()
#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...
Abstract parameters base interface.
Raster linear strategy for mixture model classification.
Raster linear mixture model strategy factory.
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.
bool getInverseMatrix(const boost::numeric::ublas::matrix< T > &inputMatrix, boost::numeric::ublas::matrix< T > &outputMatrix)
Matrix inversion.
bool m_isInitialized
True if this instance is initialized.
An abstract class for raster data strucutures.
This class can be used to inform the progress of a task.
Raster Mixture model strategy factory base class.
const Parameters & operator=(const Parameters ¶ms)