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(
TR_RP(
"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);
Raster strategy parameters base class.
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.
std::vector< te::gm::Point * > GetRandomPointsInRaster(const te::rst::Raster &inputRaster, unsigned int numberOfPoints)
Creates a vector of random positions (points) inside the raster.
Raster Processing functions.
std::vector< std::vector< double > > m_clustersMeans
The previously estimated means of the clusters (optional).
void setTotalSteps(int value)
Set the task total stepes.
This class implements and iterator to "navigate" over a raster, with a predefined number of bands...
~ClassifierEMStrategyFactory()
#define TR_RP(message)
It marks a string in order to get translated. This is a special mark used in the Raster Processing mo...
void setCurrentStep(int value)
Set the task current step.
Raster EM Classifier strategy factory.
ClassifierEMStrategyFactory()
te::rp::ClassifierStrategy * build()
Concrete factories (derived from this one) must implement this method in order to create objects...
unsigned int getRow() const
Returns the current row in iterator.
#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...
unsigned int m_maxIterations
The maximum of iterations (E/M steps) to perform if convergence is not achieved.
Abstract parameters base interface.
double m_epsilon
The stop criteria. When the clusters change in a value smaller then epsilon, the convergence is achie...
void pulse()
Calls setCurrentStep() function using getCurrentStep() + 1.
void reset()
Clear all internal allocated resources and reset the parameters instance to its initial state...
Raster classifier strategy factory base class.
bool initialize(StrategyParameters const *const strategyParams)
Initialize the classification strategy.
EM strategy for pixel-based classification. This is an unsupervised and pixel-based classification al...
bool getInverseMatrix(const boost::numeric::ublas::matrix< T > &inputMatrix, boost::numeric::ublas::matrix< T > &outputMatrix)
Matrix inversion.
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.
bool m_isInitialized
True if this instance is initialized.
This class implements the strategy to iterate with spatial restriction, the iteration occurs inside a...
ClassifierEMStrategy::Parameters m_parameters
Internal execution parameters.
bool GetDeterminant(const boost::numeric::ublas::matrix< T > &inputMatrix, double &determinant)
Get the Matrix determinant value.
unsigned int m_maxInputPoints
The maximum number of points used to estimate the clusters (default = 1000).
static RasterIterator end(Raster *r, const std::vector< unsigned int > &bands)
Returns an iterator referring to after the end of the iterator.
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.
AbstractParameters * clone() const
Create a clone copy of this instance.
const Parameters & operator=(const Parameters ¶ms)
unsigned int m_numberOfClusters
The number of clusters (classes) to estimate in the image.
void setMessage(const std::string &message)
Set the task message.
An abstract class for raster data strucutures.
This class can be used to inform the progress of a task.
unsigned int getColumn() const
Returns the current column in iterator.
Raster classifier strategy base class.
static RasterIterator begin(Raster *r, const std::vector< unsigned int > &bands)
Returns an iterator referring to the first value.
EM (Expectation-Maximization) strategy for pixel-based classification.