27 #include "../common/MatrixUtils.h"
28 #include "../common/progress/TaskProgress.h"
29 #include "../raster/Grid.h"
30 #include "../raster/PositionIterator.h"
31 #include "../raster/RasterIterator.h"
32 #include "../raster/Utils.h"
43 #include <boost/numeric/ublas/io.hpp>
44 #include <boost/numeric/ublas/lu.hpp>
45 #include <boost/numeric/ublas/matrix.hpp>
76 m_numberOfClusters = 0;
80 m_clustersMeans.clear();
121 const std::vector<te::gm::Polygon*>& inputPolygons,
te::rst::Raster& outputRaster,
122 const unsigned int outputRasterBand,
const bool enableProgressInterface)
throw(te::rp::Exception)
132 const unsigned int N = randomPoints.size();
134 const unsigned int S = inputRasterBands.size();
137 boost::numeric::ublas::matrix<double> Xk = boost::numeric::ublas::matrix<double>(N, S);
142 double max_pixel_value = 0.0;
143 while (pit != pitend)
145 for (
unsigned int l = 0; l < S; l++)
147 inputRaster.getValue(pit.getColumn(), pit.getRow(), Xk(k, l), inputRasterBands[l]);
148 if (Xk(k, l) > max_pixel_value)
149 max_pixel_value = Xk(k, l);
156 srand((
unsigned) time(0));
158 boost::numeric::ublas::matrix<double> MUj = boost::numeric::ublas::matrix<double>(M, S);
159 boost::numeric::ublas::matrix<double> previous_MUj = boost::numeric::ublas::matrix<double>(M, S);
162 for (
unsigned int j = 0; j < M; j++)
163 for (
unsigned int l = 0; l < S; l++)
169 for (
unsigned int j = 0; j < M; j++)
170 for (
unsigned int l = 0; l < S; l++)
171 MUj(j, l) = rand() % (int) ceil(max_pixel_value);
176 std::vector<boost::numeric::ublas::matrix<double> > SIGMAj;
177 for (
unsigned int j = 0; j < M; j++)
179 boost::numeric::ublas::matrix<double> tmp_sigma(S, S);
180 for (
unsigned int l = 0; l < S; l++)
182 for (
unsigned int m = 0; m < S; m++)
183 tmp_sigma(l, m) = 0.0;
184 tmp_sigma(l, l) = 10.0;
186 SIGMAj.push_back(tmp_sigma);
190 boost::numeric::ublas::matrix<double> Pj = boost::numeric::ublas::matrix<double>(M, 1);
191 boost::numeric::ublas::matrix<double> PCj_Xk = boost::numeric::ublas::matrix<double>(M, N);
192 for (
unsigned int j = 0; j < M; j++)
194 Pj(j, 0) = 1 / (double) M;
195 for (
unsigned int k = 0; k < N; k++)
201 double numerator_PCj_Xk;
202 double denominator_PCj_Xk;
207 boost::numeric::ublas::matrix<double> Xk_minus_MUj(S, 1);
208 boost::numeric::ublas::matrix<double> Xk_minus_MUj_T(1, S);
209 boost::numeric::ublas::matrix<double> product_NETAj(1, 1);
210 boost::numeric::ublas::matrix<double> product_NETAj2(1, 1);
211 boost::numeric::ublas::matrix<double> inverse_SIGMAj(S, S);
212 boost::numeric::ublas::matrix<double> inverse_SIGMAj2(S, S);
213 boost::numeric::ublas::matrix<double> product_Xk_minusMUj(S, S);
214 boost::numeric::ublas::matrix<double> sum_product_Xk_minusMUj(S, S);
220 for (
unsigned int j = 0; j < M; j++)
223 if (det_SIGMAj >= 0.0)
224 det_SIGMAj = pow(det_SIGMAj, -0.5);
228 for (
unsigned int k = 0; k < N; k++)
231 numerator_PCj_Xk = 0.0;
232 for (
unsigned int l = 0; l < S; l++)
233 Xk_minus_MUj(l, 0) = Xk(k, l) - MUj(j, l);
234 Xk_minus_MUj_T = boost::numeric::ublas::trans(Xk_minus_MUj);
238 product_NETAj = prod(Xk_minus_MUj_T, inverse_SIGMAj);
239 product_NETAj = prod(product_NETAj, Xk_minus_MUj);
240 product_NETAj *= -0.5;
242 numerator_PCj_Xk = det_SIGMAj * exp(product_NETAj(0, 0)) * Pj(j, 0);
245 denominator_PCj_Xk = 0.0;
246 for (
unsigned int j2 = 0; j2 < M; j2++)
249 if (det_SIGMAj2 >= 0.0)
250 det_SIGMAj2 = pow(det_SIGMAj2, -0.5);
254 for (
unsigned int l = 0; l < S; l++)
255 Xk_minus_MUj(l, 0) = Xk(k, l) - MUj(j2, l);
256 Xk_minus_MUj_T = boost::numeric::ublas::trans(Xk_minus_MUj);
260 product_NETAj2 = prod(Xk_minus_MUj_T, inverse_SIGMAj2);
261 product_NETAj2 = prod(product_NETAj2, Xk_minus_MUj);
262 product_NETAj2 *= -0.5;
264 denominator_PCj_Xk += det_SIGMAj2 * exp(product_NETAj2(0, 0)) * Pj(j2, 0);
266 if (denominator_PCj_Xk == 0.0)
267 denominator_PCj_Xk = 0.0000000001;
270 PCj_Xk(j, k) = numerator_PCj_Xk / denominator_PCj_Xk;
275 for (
unsigned int j = 0; j < M; j++)
278 for (
unsigned int l = 0; l < S; l++)
279 for (
unsigned int l2 = 0; l2 < S; l2++)
280 sum_product_Xk_minusMUj(l, l2) = 0.0;
281 for (
unsigned int k = 0; k < N; k++)
283 sum_PCj_Xk += PCj_Xk(j, k);
285 for (
unsigned int l = 0; l < S; l++)
286 Xk_minus_MUj(l, 0) = Xk(k, l) - MUj(j, l);
287 Xk_minus_MUj_T = boost::numeric::ublas::trans(Xk_minus_MUj);
289 product_Xk_minusMUj = prod(Xk_minus_MUj, Xk_minus_MUj_T);
290 product_Xk_minusMUj *= PCj_Xk(j, k);
292 sum_product_Xk_minusMUj += product_Xk_minusMUj;
294 if (sum_PCj_Xk == 0.0)
295 sum_PCj_Xk = 0.0000000001;
296 SIGMAj[j] = sum_product_Xk_minusMUj / sum_PCj_Xk;
300 for (
unsigned int j = 0; j < M; j++)
303 for (
unsigned int l = 0; l < S; l++)
305 for (
unsigned int k = 0; k < N; k++)
307 for (
unsigned int l = 0; l < S; l++)
308 MUj(j, l) += PCj_Xk(j, k) * Xk(k, l);
309 sum_PCj_Xk += PCj_Xk(j, k);
311 if (sum_PCj_Xk == 0.0)
312 sum_PCj_Xk = 0.0000000001;
313 for (
unsigned int l = 0; l < S; l++)
314 MUj(j, l) /= sum_PCj_Xk;
318 for (
unsigned int j = 0; j < M; j++)
321 for (
unsigned int k = 0; k < N; k++)
322 Pj(j, 0) += PCj_Xk(j, k);
329 for (
unsigned int j = 0; j < M; j++)
330 for (
unsigned int l = 0; l < S; l++)
332 a_minus_b = MUj(j, l) - previous_MUj(j, l);
333 distance_MUj += a_minus_b * a_minus_b;
335 distance_MUj = sqrt(distance_MUj);
347 boost::numeric::ublas::matrix<double> X = boost::numeric::ublas::matrix<double>(1, S);
348 boost::numeric::ublas::matrix<double> PCj_X = boost::numeric::ublas::matrix<double>(M, 1);
351 unsigned int cluster;
353 double numerator_PCj_X;
354 double denominator_PCj_X;
356 boost::numeric::ublas::matrix<double> X_minus_MUj(S, 1);
357 boost::numeric::ublas::matrix<double> X_minus_MUj_T(1, S);
359 task.
setTotalSteps(inputRaster.getNumberOfColumns() * inputRaster.getNumberOfRows());
360 task.
setMessage(
TE_TR(
"Expectation Maximization algorithm - classifying image"));
364 for (
unsigned int l = 0; l < S; l++)
368 for (
unsigned int j = 0; j < M; j++)
372 if (det_SIGMAj >= 0.0)
373 det_SIGMAj = pow(det_SIGMAj, -0.5);
377 numerator_PCj_X = 0.0;
378 for (
unsigned int l = 0; l < S; l++)
379 X_minus_MUj(l, 0) = X(0, l) - MUj(j, l);
380 X_minus_MUj_T = boost::numeric::ublas::trans(X_minus_MUj);
384 product_NETAj = prod(X_minus_MUj_T, inverse_SIGMAj);
385 product_NETAj = prod(product_NETAj, X_minus_MUj);
386 product_NETAj *= -0.5;
388 numerator_PCj_X = det_SIGMAj * exp(product_NETAj(0, 0)) * Pj(j, 0);
391 denominator_PCj_X = 0.0;
392 for (
unsigned int j2 = 0; j2 < M; j2++)
395 if (det_SIGMAj2 >= 0.0)
396 det_SIGMAj2 = pow(det_SIGMAj2, -0.5);
400 for (
unsigned int l = 0; l < S; l++)
401 X_minus_MUj(l, 0) = X(0, l) - MUj(j2, l);
402 X_minus_MUj_T = boost::numeric::ublas::trans(X_minus_MUj);
406 product_NETAj2 = prod(X_minus_MUj_T, inverse_SIGMAj2);
407 product_NETAj2 = prod(product_NETAj2, X_minus_MUj);
408 product_NETAj2 *= -0.5;
410 denominator_PCj_X += det_SIGMAj2 * exp(product_NETAj2(0, 0)) * Pj(j2, 0);
412 if (denominator_PCj_X == 0.0)
413 denominator_PCj_X = 0.0000000001;
415 PCj_X(j, 0) = numerator_PCj_X / denominator_PCj_X;
420 for (
unsigned int j = 0; j < M; j++)
421 if (PCj_X(j, 0) > max_PCj_X)
423 max_PCj_X = PCj_X(j, 0);
428 outputRaster.setValue(it.
getColumn(), it.
getRow(), cluster, outputRasterBand);
void setMessage(const std::string &message)
Set the task message.
unsigned int m_numberOfClusters
The number of clusters (classes) to estimate in the image.
unsigned int getRow() const
Returns the current row in iterator.
ClassifierEMStrategyFactory()
void reset()
Clear all internal allocated resources and reset the parameters instance to its initial state...
This class implements and iterator to "navigate" over a raster, with a predefined number of bands...
bool initialize(StrategyParameters const *const strategyParams)
Initialize the classification strategy.
This class can be used to inform the progress of a task.
static PointSetIterator end(const te::rst::Raster *r, const std::vector< te::gm::Point * > p)
Returns an iterator referring to after the end of the iterator.
EM (Expectation-Maximization) strategy for pixel-based classification.
const Parameters & operator=(const Parameters ¶ms)
#define TE_TR(message)
It marks a string in order to get translated.
This class implements the strategy to iterate with spatial restriction, the iteration occurs inside a...
Raster Processing functions.
void setTotalSteps(int value)
Set the task total stepes.
~ClassifierEMStrategyFactory()
#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 classifier strategy factory base class.
bool GetDeterminant(const boost::numeric::ublas::matrix< T > &inputMatrix, double &determinant)
Get the Matrix determinant value.
Raster EM Classifier strategy factory.
bool execute(const te::rst::Raster &inputRaster, const std::vector< unsigned int > &inputRasterBands, const std::vector< te::gm::Polygon * > &inputPolygons, te::rst::Raster &outputRaster, const unsigned int outputRasterBand, const bool enableProgressInterface)
Executes the classification strategy.
EM strategy for pixel-based classification. This is an unsupervised and pixel-based classification al...
An abstract class for raster data strucutures.
static RasterIterator begin(Raster *r, const std::vector< unsigned int > &bands)
Returns an iterator referring to the first value.
ClassifierEMStrategy::Parameters m_parameters
Internal execution parameters.
Raster strategy parameters base class.
bool m_isInitialized
True if this instance is initialized.
AbstractParameters * clone() const
Create a clone copy of this instance.
void pulse()
Calls setCurrentStep() function using getCurrentStep() + 1.
unsigned int m_maxInputPoints
The maximum number of points used to estimate the clusters (default = 1000).
Abstract parameters base interface.
unsigned int m_maxIterations
The maximum of iterations (E/M steps) to perform if convergence is not achieved.
void setCurrentStep(int value)
Set the task current step.
Raster classifier strategy base class.
bool GetInverseMatrix(const boost::numeric::ublas::matrix< T > &inputMatrix, boost::numeric::ublas::matrix< T > &outputMatrix)
Matrix inversion.
TERASTEREXPORT std::vector< te::gm::Point * > GetRandomPointsInRaster(const te::rst::Raster &inputRaster, unsigned int numberOfPoints=1000)
Creates a vector of random positions (points) inside the raster.
static RasterIterator end(Raster *r, const std::vector< unsigned int > &bands)
Returns an iterator referring to after the end of the iterator.
unsigned int getColumn() const
Returns the current column in iterator.
std::vector< std::vector< double > > m_clustersMeans
The previously estimated means of the clusters (optional).
te::rp::ClassifierStrategy * build()
Concrete factories (derived from this one) must implement this method in order to create objects...
static PointSetIterator begin(const te::rst::Raster *r, const std::vector< te::gm::Point * > p)
Returns an iterator referring to the first value of the band.
double m_epsilon
The stop criteria. When the clusters change in a value smaller then epsilon, the convergence is achie...