40 #include <visp3/visual_features/vpFeatureLuminanceMapping.h>
42 #include <visp3/core/vpSubMatrix.h>
44 #include <visp3/core/vpUniRand.h>
45 #include <visp3/core/vpIoTools.h>
46 #if defined(VISP_HAVE_CATCH2)
48 #define CATCH_CONFIG_RUNNER
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 SCENARIO(
"Using DCT features",
"[visual_features]")
318 std::vector<std::tuple<vpMatrix, vpColVector, vpMatrix>> data = {
326 { 0.0, 1.0, 3.0, 6.0, 4.0, 2.0, 5.0, 7.0, 8.0 }
335 for (
unsigned int i = 0; i < data.size(); ++i) {
336 WHEN(
"Building the associated zigzag indexing matrix")
339 vpColVector contentAsZigzag = std::get<1>(data[i]);
340 const vpMatrix mAfterWriterVec = std::get<2>(data[i]);
344 THEN(
"Calling getValues with wrong matrix rows throws")
347 REQUIRE_THROWS(zigzag.
getValues(wrongM, 0, 2, s));
349 THEN(
"Calling getValues with wrong matrix cols throws")
352 REQUIRE_THROWS(zigzag.
getValues(wrongM, 0, 2, s));
354 THEN(
"Calling getValues with wrong start and end arguments throws")
356 REQUIRE_THROWS(zigzag.
getValues(m, 2, 1, s));
358 THEN(
"Calling getValues and querying all values returns correct result")
361 REQUIRE(s == contentAsZigzag);
363 THEN(
"Calling getValues and querying a subset of the values is correct")
366 REQUIRE(s == contentAsZigzag.
extract(0, m.
size() / 2));
370 THEN(
"Calling setValues with wrong matrix rows throws")
373 REQUIRE_THROWS(zigzag.
setValues(contentAsZigzag, 0, wrongM));
375 THEN(
"Calling setValues with wrong matrix cols throws")
378 REQUIRE_THROWS(zigzag.
setValues(contentAsZigzag, 0, wrongM));
381 THEN(
"Calling setValues with wrong start and vector size arguments throws")
383 REQUIRE_THROWS(zigzag.
setValues(contentAsZigzag, m.
size() - contentAsZigzag.
size() + 1, m));
386 THEN(
"Calling setValues leads to expected result")
390 for (
unsigned i = 0; i < powered.
size(); ++i) {
391 powered[i] *= powered[i];
394 REQUIRE_NOTHROW(zigzag.
setValues(powered, 0, mWrite));
395 REQUIRE_NOTHROW(zigzag.
getValues(mWrite, 0, mWrite.size(), poweredRead));
396 REQUIRE(powered == poweredRead);
399 for (
unsigned i = 0; i < powered.
size(); ++i) {
400 indices[i] =
static_cast<double>(i);
403 REQUIRE_NOTHROW(zigzag.
setValues(indices, 0, mWrite));
404 REQUIRE(mWrite == mAfterWriterVec);
411 REQUIRE(s2 == contentAsZigzag.
extract(0, 3));
419 GIVEN(
"A constant image")
422 WHEN(
"Computing DCT")
428 THEN(
"resulting feature vector has correct size")
430 REQUIRE(s.
size() == 32);
432 THEN(
"The only non zero component is the first")
434 REQUIRE(s.
sum() == Approx(s[0]).margin(1e-5));
440 for (
unsigned i = 0; i < I.
getRows(); ++i) {
441 for (
unsigned j = 0; j < I.
getCols(); ++j) {
442 const int diff = abs(
static_cast<int>(I[i][j]) -
static_cast<int>(Ir[i][j]));
444 INFO(
"i = " + std::to_string(i) +
", j = " + std::to_string(j));
452 int main(
int argc,
char *argv[])
454 Catch::Session session;
455 session.applyCommandLine(argc, argv);
457 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)
Initalize 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< std::string > &imageFiles, const unsigned int projectionSize, const unsigned int imageBorder=0)
Compute a new Principal Component Analysis on set of images, stored on disk.
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.