Visual Servoing Platform  version 3.6.1 under development (2025-02-18)
catchLuminanceMapping.cpp
1 
2 /*
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2024 by Inria. All rights reserved.
5  *
6  * This software is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  * See the file LICENSE.txt at the root directory of this source
11  * distribution for additional information about the GNU GPL.
12  *
13  * For using ViSP with software that can not be combined with the GNU
14  * GPL, please contact Inria about acquiring a ViSP Professional
15  * Edition License.
16  *
17  * See https://visp.inria.fr for more information.
18  *
19  * This software was developed at:
20  * Inria Rennes - Bretagne Atlantique
21  * Campus Universitaire de Beaulieu
22  * 35042 Rennes Cedex
23  * France
24  *
25  * If you have questions regarding the use of this file, please contact
26  * Inria at visp@inria.fr
27  *
28  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30  *
31  * Description:
32  * Performs various tests on the point class.
33  */
34 
40 #include <visp3/core/vpConfig.h>
41 
42 #if defined(VISP_HAVE_CATCH2)
43 
44 #include <catch_amalgamated.hpp>
45 
46 #include <visp3/core/vpSubMatrix.h>
47 #include <visp3/core/vpUniRand.h>
48 #include <visp3/core/vpIoTools.h>
49 #include <visp3/visual_features/vpFeatureLuminanceMapping.h>
50 
51 #ifdef ENABLE_VISP_NAMESPACE
52 using namespace VISP_NAMESPACE_NAME;
53 #endif
54 
55 vpMatrix orthogonalBasis(unsigned n, unsigned seed)
56 {
57  vpUniRand rand(seed);
58  vpMatrix basis(n, n);
59  vpColVector norms(n);
60 
61  //start with random basis
62  for (unsigned int row = 0; row < n; ++row) {
63  double norm = 0.0;
64  for (unsigned int col = 0; col < n; ++col) {
65  basis[row][col] = rand.uniform(-1.0, 1.0);
66  norm += basis[row][col] * basis[row][col];
67  }
68  norm = 1.0 / sqrt(norm);
69  for (unsigned int col = 0; col < n; ++col) {
70  basis[row][col] *= norm;
71  }
72 
73  }
74  // Apply gram schmidt process
75  norms[0] = basis.getRow(0).sumSquare();
76  for (unsigned i = 1; i < n; ++i) {
77  vpColVector uit = basis.getRow(i).t();
78 
79  for (unsigned j = 0; j < i; ++j) {
80  vpRowVector vj = basis.getRow(j);
81  vpRowVector res = vj * ((vj * uit) / (norms[j]));
82  for (unsigned k = 0; k < n; ++k) {
83  basis[i][k] -= res[k];
84  }
85  }
86  norms[i] = basis.getRow(i).sumSquare();
87  }
88  for (unsigned int row = 0; row < n; ++row) {
89  double norm = sqrt(norms[row]);
90 
91  for (unsigned int col = 0; col < n; ++col) {
92  basis[row][col] /= norm;
93  }
94  }
95  return basis;
96 
97 }
98 
99 SCENARIO("Using PCA features", "[visual_features]")
100 {
101  GIVEN("A matrix containing simple data")
102  {
103  const unsigned h = 16, w = 16;
104  const unsigned numDataPoints = 4;
105  const unsigned int dataDim = h * w;
106  const unsigned int trueComponents = 3;
107  // Generate numDataPoints vectors in a "dataDim"-dimensional space.
108  // The data is generated from "trueComponents" vectors, that are orthogonal
109  const vpMatrix orthoFull = (orthogonalBasis(dataDim, 42) + vpMatrix(dataDim, dataDim, 1.0)) * 127.5; // dataDim x dataDim
110  const vpMatrix ortho(orthoFull, 0, 0, trueComponents, dataDim); // trueComponents X dataDim
111  const vpMatrix coefficients(numDataPoints, trueComponents);
112  vpUniRand rand(17);
113  for (unsigned int i = 0; i < coefficients.getRows(); ++i) {
114  double sum = 0.0;
115  for (unsigned int j = 0; j < coefficients.getCols(); ++j) {
116  coefficients[i][j] = rand.uniform(0.0, 1.0);
117  sum += coefficients[i][j] * coefficients[i][j];
118  }
119  const double inv_norm = 1.0 / sqrt(sum);
120  for (unsigned int j = 0; j < coefficients.getCols(); ++j) {
121  coefficients[i][j] *= inv_norm;
122  }
123  }
124 
125  vpMatrix data = coefficients * ortho;
126 
127  WHEN("Learning PCA basis with too many components")
128  {
129  unsigned int k = data.getCols() + 1;
130  THEN("An exception is thrown")
131  {
132  REQUIRE_THROWS(vpLuminancePCA::learn(data.transpose(), k));
133  }
134  }
135  WHEN("Learning with more images than pixels")
136  {
137  vpMatrix wrongData(20, 50);
138  THEN("An exception is thrown")
139  {
140  REQUIRE_THROWS(vpLuminancePCA::learn(wrongData.transpose(), 32));
141  }
142  }
143  WHEN("Learning PCA basis")
144  {
145  for (unsigned int k = 2; k <= trueComponents; ++k) {
146 
148  const vpMatrix &basis = *pca.getBasis();
149 
150  THEN("Basis has correct dimensions")
151  {
152  REQUIRE(basis.getRows() == k);
153  REQUIRE(basis.getCols() == dataDim);
154  }
155  THEN("The basis is orthonormal")
156  {
157  const vpMatrix Iapprox = basis * basis.t();
158  vpMatrix I;
159  I.eye(basis.getRows());
160  bool matrixSame = true;
161  for (unsigned int row = 0; row < I.getRows(); ++row) {
162  for (unsigned int col = 0; col < I.getCols(); ++col) {
163  if (fabs(I[row][col] - Iapprox[row][col]) > 1e-6) {
164  matrixSame = false;
165  break;
166  }
167  }
168  }
169  REQUIRE(matrixSame);
170  }
171  THEN("Mean vector has correct dimensions")
172  {
173  REQUIRE(pca.getMean()->getRows() == dataDim);
174  REQUIRE(pca.getMean()->getCols() == 1);
175  }
176  THEN("Modifying the basis size (number of inputs) by hand and saving")
177  {
178  const std::string tempDir = vpIoTools::makeTempDirectory("visp_test_pca_wrong");
179  const std::string basisFile = vpIoTools::createFilePath(tempDir, "basis.txt");
180  const std::string meanFile = vpIoTools::createFilePath(tempDir, "mean.txt");
181  const std::string varFile = vpIoTools::createFilePath(tempDir, "var.txt");
182  pca.getBasis()->resize(pca.getBasis()->getRows(), pca.getBasis()->getCols() - 1);
183  REQUIRE_THROWS(pca.save(basisFile, meanFile, varFile));
184  }
185  THEN("Modifying the mean Columns by hand")
186  {
187  const std::string tempDir = vpIoTools::makeTempDirectory("visp_test_pca_wrong");
188  const std::string basisFile = vpIoTools::createFilePath(tempDir, "basis.txt");
189  const std::string meanFile = vpIoTools::createFilePath(tempDir, "mean.txt");
190  const std::string varFile = vpIoTools::createFilePath(tempDir, "var.txt");
191 
192  std::shared_ptr<vpColVector> mean = pca.getMean();
193  mean->resize(mean->getRows() + 1, false);
194  REQUIRE_THROWS(pca.save(basisFile, meanFile, varFile));
195 
196  }
197  THEN("Saving and loading pca leads to same basis and mean")
198  {
199  const std::string tempDir = vpIoTools::makeTempDirectory("visp_test_pca");
200  const std::string basisFile = vpIoTools::createFilePath(tempDir, "basis.txt");
201  const std::string meanFile = vpIoTools::createFilePath(tempDir, "mean.txt");
202  const std::string varFile = vpIoTools::createFilePath(tempDir, "var.txt");
203 
204  pca.save(basisFile, meanFile, varFile);
205 
206  const vpLuminancePCA pca2 = vpLuminancePCA::load(basisFile, meanFile, varFile);
207  const vpMatrix basisDiff = *pca.getBasis() - *pca2.getBasis();
208  const vpColVector meanDiff = *pca.getMean() - *pca2.getMean();
209  const vpColVector explainedVarDiff = pca.getExplainedVariance() - pca2.getExplainedVariance();
210  bool basisSame = true;
211  bool meanSame = true;
212  bool explainedVarSame = true;
213 
214  for (unsigned int i = 0; i < basisDiff.getRows(); ++i) {
215  for (unsigned int j = 0; j < basisDiff.getCols(); ++j) {
216  if (fabs(basisDiff[i][j]) > 1e-10) {
217  basisSame = false;
218  break;
219  }
220  }
221  }
222  REQUIRE(basisSame);
223 
224  for (unsigned int i = 0; i < meanDiff.getRows(); ++i) {
225  if (fabs(meanDiff[i]) > 1e-10) {
226  std::cout << meanDiff << std::endl;
227  meanSame = false;
228  break;
229  }
230  }
231  REQUIRE(meanSame);
232 
233  for (unsigned int i = 0; i < explainedVarDiff.getRows(); ++i) {
234  if (fabs(explainedVarDiff[i]) > 1e-10) {
235  explainedVarSame = false;
236  break;
237  }
238  }
239  REQUIRE(explainedVarSame);
240  }
241 
242  THEN("Explained variance is below 1 and sorted in descending order")
243  {
244  const vpColVector var = pca.getExplainedVariance();
245  REQUIRE(var.sum() < 1.0);
246  for (int i = 1; i < (int)var.getRows() - 1; ++i) {
247  REQUIRE(var[i] >= var[i + 1]);
248  }
249  }
250  if (k == trueComponents) {
251  WHEN("K is the true manifold dimensionality")
252  {
253  THEN("explained variance is close to 1")
254  {
255  REQUIRE(pca.getExplainedVariance().sum() > 0.99);
256  }
257  THEN("Inverse mapping leads back to the same data")
258  {
259  for (unsigned int i = 0; i < numDataPoints; ++i) {
260  vpImage<unsigned char> I(h, w);
261  for (unsigned int j = 0; j < data.getCols(); ++j) {
262  I.bitmap[j] = static_cast<unsigned char>(data[i][j]);
263  }
264  vpColVector s;
265  pca.setBorder(0);
266  pca.map(I, s);
268  pca.inverse(s, Irec);
269  for (unsigned int j = 0; j < data.getCols(); ++j) {
270  REQUIRE(abs(static_cast<int>(I.bitmap[j]) - static_cast<int>(Irec.bitmap[j])) < 2);
271  }
272  }
273  }
274  }
275  }
276 
277  THEN("Projecting data is correct")
278  {
279  {
280  vpColVector s;
281  pca.setBorder(0);
282  vpImage<unsigned char> I(h, w);
283  pca.map(I, s);
284  REQUIRE(s.size() == pca.getProjectionSize());
285  }
286  {
287  vpColVector s;
288  const unsigned border = 3;
289  pca.setBorder(border);
290  REQUIRE(pca.getBorder() == border);
291  vpImage<unsigned char> I(h + 2 * border, w + 2 * border);
292  pca.map(I, s);
293  REQUIRE(s.size() == pca.getProjectionSize());
294  }
295  }
296  }
297  }
298  }
299  WHEN("Saving unintialized PCA")
300  {
301  vpLuminancePCA pca;
302  const std::string tempDir = vpIoTools::makeTempDirectory("visp_test_pca");
303  const std::string basisFile = vpIoTools::createFilePath(tempDir, "basis.txt");
304  const std::string meanFile = vpIoTools::createFilePath(tempDir, "mean.txt");
305  const std::string varFile = vpIoTools::createFilePath(tempDir, "var.txt");
306  THEN("an exception is thrown")
307  {
308  REQUIRE_THROWS(pca.save(basisFile, meanFile, varFile));
309  }
310  }
311 }
312 
313 #if (VISP_CXX_STANDARD > VISP_CXX_STANDARD_11)
314 SCENARIO("Using DCT features", "[visual_features]")
315 {
316  GIVEN("A matrix")
317  {
318  vpMatrix M1({
319  {0.0, 1.0, 2.0},
320  {3.0, 4.0, 5.0},
321  {6.0, 7.0, 8.0}
322  });
323  vpColVector v(
324  { 0.0, 1.0, 3.0, 6.0, 4.0, 2.0, 5.0, 7.0, 8.0 }
325  );
326  vpMatrix M2({
327  {0.0, 1.0, 5.0},
328  {2.0, 4.0, 6.0},
329  {3.0, 7.0, 8.0}
330  });
331  std::tuple<vpMatrix, vpColVector, vpMatrix> data(M1, v, M2);
332 
333  WHEN("Building the associated zigzag indexing matrix")
334  {
335  vpMatrix m = std::get<0>(data);
336  vpColVector contentAsZigzag = std::get<1>(data);
337  const vpMatrix mAfterWriterVec = std::get<2>(data);
339  zigzag.init(m.getRows(), m.getCols());
340  vpColVector s;
341  THEN("Calling getValues with wrong matrix rows throws")
342  {
343  vpMatrix wrongM(m.getRows() + 1, m.getCols());
344  REQUIRE_THROWS(zigzag.getValues(wrongM, 0, 2, s));
345  }
346  THEN("Calling getValues with wrong matrix cols throws")
347  {
348  vpMatrix wrongM(m.getRows(), m.getCols() + 1);
349  REQUIRE_THROWS(zigzag.getValues(wrongM, 0, 2, s));
350  }
351  THEN("Calling getValues with wrong start and end arguments throws")
352  {
353  REQUIRE_THROWS(zigzag.getValues(m, 2, 1, s));
354  }
355  THEN("Calling getValues and querying all values returns correct result")
356  {
357  REQUIRE_NOTHROW(zigzag.getValues(m, 0, m.size(), s));
358  REQUIRE(s == contentAsZigzag);
359  }
360  THEN("Calling getValues and querying a subset of the values is correct")
361  {
362  REQUIRE_NOTHROW(zigzag.getValues(m, 0, m.size() / 2, s));
363  REQUIRE(s == contentAsZigzag.extract(0, m.size() / 2));
364  REQUIRE_NOTHROW(zigzag.getValues(m, m.size() / 2, m.size(), s));
365  REQUIRE(s == contentAsZigzag.extract(m.size() / 2, m.size() - m.size() / 2));
366  }
367  THEN("Calling setValues with wrong matrix rows throws")
368  {
369  vpMatrix wrongM(m.getRows() + 1, m.getCols());
370  REQUIRE_THROWS(zigzag.setValues(contentAsZigzag, 0, wrongM));
371  }
372  THEN("Calling setValues with wrong matrix cols throws")
373  {
374  vpMatrix wrongM(m.getRows(), m.getCols() + 1);
375  REQUIRE_THROWS(zigzag.setValues(contentAsZigzag, 0, wrongM));
376  }
377 
378  THEN("Calling setValues with wrong start and vector size arguments throws")
379  {
380  REQUIRE_THROWS(zigzag.setValues(contentAsZigzag, m.size() - contentAsZigzag.size() + 1, m));
381  }
382 
383  THEN("Calling setValues leads to expected result")
384  {
385  vpMatrix mWrite(m.getRows(), m.getCols());
386  vpColVector powered = contentAsZigzag;
387  for (unsigned i = 0; i < powered.size(); ++i) {
388  powered[i] *= powered[i];
389  }
390  vpColVector poweredRead;
391  REQUIRE_NOTHROW(zigzag.setValues(powered, 0, mWrite));
392  REQUIRE_NOTHROW(zigzag.getValues(mWrite, 0, mWrite.size(), poweredRead));
393  REQUIRE(powered == poweredRead);
394 
395  vpColVector indices = contentAsZigzag;
396  for (unsigned i = 0; i < powered.size(); ++i) {
397  indices[i] = static_cast<double>(i);
398  }
399  vpColVector indicesRead;
400  REQUIRE_NOTHROW(zigzag.setValues(indices, 0, mWrite));
401  REQUIRE(mWrite == mAfterWriterVec);
402 
403  vpMatrix m2(m.getRows(), m.getCols(), 0.0);
404  zigzag.setValues(contentAsZigzag.extract(0, 3), 0, m2);
405 
406  vpColVector s2;
407  zigzag.getValues(m2, 0, 3, s2);
408  REQUIRE(s2 == contentAsZigzag.extract(0, 3));
409  }
410  }
411 
412 
413  GIVEN("A constant image")
414  {
415  vpImage<unsigned char> I(32, 64, 20);
416  WHEN("Computing DCT")
417  {
418  vpLuminanceDCT dct(32);
419  dct.setBorder(0);
420  vpColVector s;
421  dct.map(I, s);
422  THEN("resulting feature vector has correct size")
423  {
424  REQUIRE(s.size() == 32);
425  }
426  THEN("The only non zero component is the first")
427  {
428  REQUIRE(s.sum() == Catch::Approx(s[0]).margin(1e-5));
429  }
430 
432  dct.inverse(s, Ir);
433  REQUIRE((Ir.getRows() == I.getRows() && Ir.getCols() == I.getCols()));
434  for (unsigned i = 0; i < I.getRows(); ++i) {
435  for (unsigned j = 0; j < I.getCols(); ++j) {
436  const int diff = abs(static_cast<int>(I[i][j]) - static_cast<int>(Ir[i][j]));
437  REQUIRE(diff < 2);
438  INFO("i = " + std::to_string(i) + ", j = " + std::to_string(j));
439  }
440  }
441  }
442  }
443  }
444 }
445 #endif
446 
447 int main(int argc, char *argv[])
448 {
449  Catch::Session session; // There must be exactly one instance
450  session.applyCommandLine(argc, argv);
451 
452  int numFailed = session.run();
453  return numFailed;
454 }
455 
456 #else
457 
458 int main()
459 {
460  return EXIT_SUCCESS;
461 }
462 #endif
unsigned int getCols() const
Definition: vpArray2D.h:417
unsigned int size() const
Return the number of elements of the 2D array.
Definition: vpArray2D.h:429
unsigned int getRows() const
Definition: vpArray2D.h:427
Implementation of column vector and the associated operations.
Definition: vpColVector.h:191
vpColVector extract(unsigned int r, unsigned int colsize) const
Definition: vpColVector.h:410
vpRowVector t() const
double sum() const
unsigned int getCols() const
Definition: vpImage.h:171
Type * bitmap
points toward the bitmap
Definition: vpImage.h:135
unsigned int getRows() const
Definition: vpImage.h:212
static std::string createFilePath(const std::string &parent, const std::string &child)
Definition: vpIoTools.cpp:1427
static std::string makeTempDirectory(const std::string &dirname)
Definition: vpIoTools.cpp:708
Helper class to iterate and get/set the values from a matrix, following a zigzag pattern.
void init(unsigned rows, unsigned cols)
Initialize the ZigZag object. Computes and stores the zigzag indexing for a given matrix size.
void setValues(const vpColVector &s, unsigned int start, vpMatrix &m) const
set the values in the matrix, according to the values stored in the vector s and the zigzag indexing ...
void getValues(const vpMatrix &m, unsigned int start, unsigned int end, vpColVector &s) const
Fill the vector s with (end - start) values, according to the zigzag matrix indexing strategy.
Implementation of .
unsigned int getProjectionSize() const
Returns the size of the space to which an image is mapped to.
unsigned int getBorder() const
Returns the number of pixels that are removed by the photometric VS computation.
void setBorder(unsigned border)
Set the number of pixels that are removed by the photometric VS computation This function should be c...
Implementation of .
static vpLuminancePCA learn(const std::vector< vpImage< unsigned char >> &images, const unsigned int projectionSize, const unsigned int imageBorder=0)
Compute a new Principal Component Analysis on set of images.
void inverse(const vpColVector &s, vpImage< unsigned char > &I) VP_OVERRIDE
Reconstruct I from a representation s.
void save(const std::string &basisFilename, const std::string &meanFileName, const std::string &explainedVarianceFile) const
Save the PCA basis to multiple text files, for later use via the load function.
std::shared_ptr< vpColVector > getMean() const
Get , the mean image computed from the dataset.
void map(const vpImage< unsigned char > &I, vpColVector &s) VP_OVERRIDE
Map an image I to a representation s. This representation s has getProjectionSize() rows.
std::shared_ptr< vpMatrix > getBasis() const
Get , the subspace projection matrix ( )
static vpLuminancePCA load(const std::string &basisFilename, const std::string &meanFileName, const std::string &explainedVarianceFile)
Save the PCA basis to multiple text files, for later use via the load function.
vpColVector getExplainedVariance() const
Get the values of explained variance by each of the eigen vectors.
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:169
vpMatrix transpose() const
vpMatrix t() const
Implementation of row vector and the associated operations.
Definition: vpRowVector.h:124
double sumSquare() const
Class for generating random numbers with uniform probability density.
Definition: vpUniRand.h:127