40 #include <visp3/core/vpConfig.h>
42 #if defined(VISP_HAVE_CATCH2)
44 #include <catch_amalgamated.hpp>
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>
51 #ifdef ENABLE_VISP_NAMESPACE
55 vpMatrix orthogonalBasis(
unsigned n,
unsigned seed)
62 for (
unsigned int row = 0; row < n; ++row) {
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];
68 norm = 1.0 / sqrt(norm);
69 for (
unsigned int col = 0; col < n; ++col) {
70 basis[row][col] *= norm;
75 norms[0] = basis.getRow(0).sumSquare();
76 for (
unsigned i = 1; i < n; ++i) {
79 for (
unsigned j = 0; j < i; ++j) {
82 for (
unsigned k = 0; k < n; ++k) {
83 basis[i][k] -= res[k];
88 for (
unsigned int row = 0; row < n; ++row) {
89 double norm = sqrt(norms[row]);
91 for (
unsigned int col = 0; col < n; ++col) {
92 basis[row][col] /= norm;
99 SCENARIO(
"Using PCA features",
"[visual_features]")
101 GIVEN(
"A matrix containing simple data")
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;
109 const vpMatrix orthoFull = (orthogonalBasis(dataDim, 42) +
vpMatrix(dataDim, dataDim, 1.0)) * 127.5;
110 const vpMatrix ortho(orthoFull, 0, 0, trueComponents, dataDim);
111 const vpMatrix coefficients(numDataPoints, trueComponents);
113 for (
unsigned int i = 0; i < coefficients.getRows(); ++i) {
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];
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;
125 vpMatrix data = coefficients * ortho;
127 WHEN(
"Learning PCA basis with too many components")
129 unsigned int k = data.
getCols() + 1;
130 THEN(
"An exception is thrown")
135 WHEN(
"Learning with more images than pixels")
138 THEN(
"An exception is thrown")
143 WHEN(
"Learning PCA basis")
145 for (
unsigned int k = 2; k <= trueComponents; ++k) {
150 THEN(
"Basis has correct dimensions")
153 REQUIRE(basis.
getCols() == dataDim);
155 THEN(
"The basis is orthonormal")
157 const vpMatrix Iapprox = basis * basis.
t();
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) {
171 THEN(
"Mean vector has correct dimensions")
173 REQUIRE(pca.
getMean()->getRows() == dataDim);
174 REQUIRE(pca.
getMean()->getCols() == 1);
176 THEN(
"Modifying the basis size (number of inputs) by hand and saving")
183 REQUIRE_THROWS(pca.
save(basisFile, meanFile, varFile));
185 THEN(
"Modifying the mean Columns by hand")
192 std::shared_ptr<vpColVector> mean = pca.
getMean();
193 mean->resize(mean->getRows() + 1,
false);
194 REQUIRE_THROWS(pca.
save(basisFile, meanFile, varFile));
197 THEN(
"Saving and loading pca leads to same basis and mean")
204 pca.
save(basisFile, meanFile, varFile);
210 bool basisSame =
true;
211 bool meanSame =
true;
212 bool explainedVarSame =
true;
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) {
224 for (
unsigned int i = 0; i < meanDiff.
getRows(); ++i) {
225 if (fabs(meanDiff[i]) > 1e-10) {
226 std::cout << meanDiff << std::endl;
233 for (
unsigned int i = 0; i < explainedVarDiff.
getRows(); ++i) {
234 if (fabs(explainedVarDiff[i]) > 1e-10) {
235 explainedVarSame =
false;
239 REQUIRE(explainedVarSame);
242 THEN(
"Explained variance is below 1 and sorted in descending order")
245 REQUIRE(var.
sum() < 1.0);
246 for (
int i = 1; i < (int)var.
getRows() - 1; ++i) {
247 REQUIRE(var[i] >= var[i + 1]);
250 if (k == trueComponents) {
251 WHEN(
"K is the true manifold dimensionality")
253 THEN(
"explained variance is close to 1")
257 THEN(
"Inverse mapping leads back to the same data")
259 for (
unsigned int i = 0; i < numDataPoints; ++i) {
261 for (
unsigned int j = 0; j < data.
getCols(); ++j) {
262 I.
bitmap[j] =
static_cast<unsigned char>(data[i][j]);
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);
277 THEN(
"Projecting data is correct")
288 const unsigned border = 3;
299 WHEN(
"Saving unintialized PCA")
306 THEN(
"an exception is thrown")
308 REQUIRE_THROWS(pca.
save(basisFile, meanFile, varFile));
313 #if (VISP_CXX_STANDARD > VISP_CXX_STANDARD_11)
314 SCENARIO(
"Using DCT features",
"[visual_features]")
324 { 0.0, 1.0, 3.0, 6.0, 4.0, 2.0, 5.0, 7.0, 8.0 }
331 std::tuple<vpMatrix, vpColVector, vpMatrix> data(M1, v, M2);
333 WHEN(
"Building the associated zigzag indexing matrix")
337 const vpMatrix mAfterWriterVec = std::get<2>(data);
341 THEN(
"Calling getValues with wrong matrix rows throws")
344 REQUIRE_THROWS(zigzag.
getValues(wrongM, 0, 2, s));
346 THEN(
"Calling getValues with wrong matrix cols throws")
349 REQUIRE_THROWS(zigzag.
getValues(wrongM, 0, 2, s));
351 THEN(
"Calling getValues with wrong start and end arguments throws")
353 REQUIRE_THROWS(zigzag.
getValues(m, 2, 1, s));
355 THEN(
"Calling getValues and querying all values returns correct result")
358 REQUIRE(s == contentAsZigzag);
360 THEN(
"Calling getValues and querying a subset of the values is correct")
363 REQUIRE(s == contentAsZigzag.
extract(0, m.
size() / 2));
367 THEN(
"Calling setValues with wrong matrix rows throws")
370 REQUIRE_THROWS(zigzag.
setValues(contentAsZigzag, 0, wrongM));
372 THEN(
"Calling setValues with wrong matrix cols throws")
375 REQUIRE_THROWS(zigzag.
setValues(contentAsZigzag, 0, wrongM));
378 THEN(
"Calling setValues with wrong start and vector size arguments throws")
380 REQUIRE_THROWS(zigzag.
setValues(contentAsZigzag, m.
size() - contentAsZigzag.
size() + 1, m));
383 THEN(
"Calling setValues leads to expected result")
387 for (
unsigned i = 0; i < powered.
size(); ++i) {
388 powered[i] *= powered[i];
391 REQUIRE_NOTHROW(zigzag.
setValues(powered, 0, mWrite));
392 REQUIRE_NOTHROW(zigzag.
getValues(mWrite, 0, mWrite.size(), poweredRead));
393 REQUIRE(powered == poweredRead);
396 for (
unsigned i = 0; i < powered.
size(); ++i) {
397 indices[i] =
static_cast<double>(i);
400 REQUIRE_NOTHROW(zigzag.
setValues(indices, 0, mWrite));
401 REQUIRE(mWrite == mAfterWriterVec);
408 REQUIRE(s2 == contentAsZigzag.
extract(0, 3));
413 GIVEN(
"A constant image")
416 WHEN(
"Computing DCT")
422 THEN(
"resulting feature vector has correct size")
424 REQUIRE(s.
size() == 32);
426 THEN(
"The only non zero component is the first")
428 REQUIRE(s.
sum() == Catch::Approx(s[0]).margin(1e-5));
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]));
438 INFO(
"i = " + std::to_string(i) +
", j = " + std::to_string(j));
447 int main(
int argc,
char *argv[])
449 Catch::Session session;
450 session.applyCommandLine(argc, argv);
452 int numFailed = session.run();
unsigned int getCols() const
unsigned int size() const
Return the number of elements of the 2D array.
unsigned int getRows() const
Implementation of column vector and the associated operations.
vpColVector extract(unsigned int r, unsigned int colsize) const
unsigned int getCols() const
Type * bitmap
points toward the bitmap
unsigned int getRows() const
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.
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...
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.
vpMatrix transpose() const
Implementation of row vector and the associated operations.
Class for generating random numbers with uniform probability density.