ViSP
 All Classes Functions Variables Enumerations Enumerator Friends Groups Pages
testSvd.cpp
1 /****************************************************************************
2  *
3  * $Id: testSvd.cpp 4056 2013-01-05 13:04:42Z fspindle $
4  *
5  * This file is part of the ViSP software.
6  * Copyright (C) 2005 - 2013 by INRIA. All rights reserved.
7  *
8  * This software is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU General Public License
10  * ("GPL") version 2 as published by the Free Software Foundation.
11  * See the file LICENSE.txt at the root directory of this source
12  * distribution for additional information about the GNU GPL.
13  *
14  * For using ViSP with software that can not be combined with the GNU
15  * GPL, please contact INRIA about acquiring a ViSP Professional
16  * Edition License.
17  *
18  * See http://www.irisa.fr/lagadic/visp/visp.html for more information.
19  *
20  * This software was developed at:
21  * INRIA Rennes - Bretagne Atlantique
22  * Campus Universitaire de Beaulieu
23  * 35042 Rennes Cedex
24  * France
25  * http://www.irisa.fr/lagadic
26  *
27  * If you have questions regarding the use of this file, please contact
28  * INRIA at visp@inria.fr
29  *
30  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
31  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
32  *
33  *
34  * Description:
35  * Test various svd decompositions.
36  *
37  * Authors:
38  * Eric Marchand
39  * Fabien Spindler
40  *
41  *****************************************************************************/
42 
43 
51 #include <visp/vpTime.h>
52 
53 #include <visp/vpMatrix.h>
54 #include <visp/vpColVector.h>
55 #include <visp/vpParseArgv.h>
56 #include <vector>
57 #include <algorithm>
58 #include <stdlib.h>
59 #include <stdio.h>
60 // List of allowed command line options
61 #define GETOPTARGS "h"
62 
68 void usage(const char *name, const char *badparam)
69 {
70  fprintf(stdout, "\n\
71 Test various svd decompositions.\n\
72 \n\
73 SYNOPSIS\n\
74  %s [-h]\n", name);
75 
76  fprintf(stdout, "\n\
77 OPTIONS: Default\n\
78  -h\n\
79  Print the help.\n");
80 
81  if (badparam)
82  fprintf(stdout, "\nERROR: Bad parameter [%s]\n", badparam);
83 }
91 bool getOptions(int argc, const char **argv)
92 {
93  const char *optarg;
94  int c;
95  while ((c = vpParseArgv::parse(argc, argv, GETOPTARGS, &optarg)) > 1) {
96 
97  switch (c) {
98  case 'h': usage(argv[0], NULL); return false; break;
99 
100  default:
101  usage(argv[0], optarg);
102  return false; break;
103  }
104  }
105 
106  if ((c == 1) || (c == -1)) {
107  // standalone param or error
108  usage(argv[0], NULL);
109  std::cerr << "ERROR: " << std::endl;
110  std::cerr << " Bad argument " << optarg << std::endl << std::endl;
111  return false;
112  }
113 
114  return true;
115 }
116 #define abs(x) ((x) < 0 ? - (x) : (x))
117 #ifdef VISP_HAVE_GSL
118 
119 
120 bool testRandom(double epsilon)
121 {
122  vpMatrix L0(6,6);
123  vpMatrix L1(6,6);
124 
125  for (unsigned int i=0 ; i < L0.getRows() ; i++)
126  for (unsigned int j=0 ; j < L0.getCols() ; j++)
127  L1[i][j] = L0[i][j] = (double)rand()/(double)RAND_MAX;
128 
129 
130 
131  vpColVector W0(L0.getCols()) ;
132  vpMatrix V0(L0.getCols(), L0.getCols()) ;
133  vpColVector W1(L1.getCols()) ;
134  vpMatrix V1(L1.getCols(), L1.getCols()) ;
135 
136  L0.svdNr(W0,V0);
137  L1.svdGsl(W1,V1);
138 
139  vpColVector _W0 = vpColVector::sort(W0);
140  vpColVector _W1 = vpColVector::sort(W1);
141 
142  vpColVector diff = _W0-_W1;
143  double error=-1.0;
144 
145  for(unsigned int i=0;i<6;i++)
146  error=std::max(abs(diff[i]),error);
147 
148  return error<epsilon;
149 
150 }
151 
152 #endif
153 
154 bool testSvdOpenCvGSLCoherence(double epsilon){
155 #if (VISP_HAVE_OPENCV_VERSION >= 0x020101) && defined (VISP_HAVE_GSL) // Require opencv >= 2.1.1
156  vpMatrix A;
157  vpMatrix vA;
158  vpColVector wA;
159  vpMatrix B;
160  vpMatrix vB;
161  vpColVector wB;
162  A.resize(6,6);
163  B.resize(6,6);
164  vA.resize(6,6);
165  vB.resize(6,6);
166  wA.resize(6);
167  wB.resize(6);
168 
169  for (unsigned int i=0 ; i < A.getRows() ; i++)
170  for (unsigned int j=0 ; j < A.getCols() ; j++)
171  B[i][j] = A[i][j] = (double)rand()/(double)RAND_MAX;
172 
173  A.svdOpenCV(wA,vA);
174  B.svdGsl(wB,vB);
175 
176  bool error=false;
177  for (unsigned int i=0 ; i < A.getRows() ; i++){
178  error = error | (abs(wA[i]-wB[i])>epsilon);
179  }
180 
181  return !error;
182 #else
183  (void)epsilon;
184  return true;
185 #endif
186 }
187 
188 int
189 main(int argc, const char ** argv)
190 {
191  // Read the command line options
192  if (getOptions(argc, argv) == false) {
193  exit (-1);
194  }
195 
196  unsigned int i,j ;
197  vpMatrix L(60000,6), Ls ;
198  for (i=0 ; i < L.getRows() ; i++)
199  for (j=0 ; j < L.getCols() ; j++)
200  L[i][j] = 2*i+j + cos((double)(i+j))+((double)(i)) ;
201  // std::cout << L << std::endl ;
202  Ls = L ;
203  std::cout << "--------------------------------------"<<std::endl ;
204 
205  vpColVector W(L.getCols()) ;
206  vpMatrix V(L.getCols(), L.getCols()) ;
207 
208  double t = vpTime::measureTimeMs() ;
209  L.svdNr(W,V) ;
210  t = vpTime::measureTimeMs() -t ;
211 
212  std::cout <<"svdNr Numerical recipes \n time " <<t << std::endl;
213  std::cout << W.t() ;
214  std::cout << "--------------------------------------"<<std::endl ;
215 
216 
217 #ifdef VISP_HAVE_GSL
218  L = Ls ;
219  t = vpTime::measureTimeMs() ;
220  L.svdGsl(W,V) ;
221  t = vpTime::measureTimeMs() -t ;
222  std::cout <<"svdGsl_mod \n time " <<t << std::endl;
223  std::cout << W.t() ;
224 
225  std::cout << "--------------------------------------"<<std::endl ;
226  std::cout << "TESTING RANDOM MATRICES:" ;
227 
228  bool ret = true;
229  for(int i=0;i<2000;i++)
230  ret = ret & testRandom(0.00001);
231  if(ret)
232  std:: cout << "Success"<< std:: endl;
233  else
234  std:: cout << "Fail"<< std:: endl;
235 
236  std::cout << "--------------------------------------"<<std::endl ;
237 #endif
238 
239  std::cout << "--------------------------------------"<<std::endl ;
240  std::cout << "TESTING OPENCV-GSL coherence:" ;
241 
242  bool ret2 = true;
243  for(int i=0;i<1;i++)
244  ret2 = ret2 & testSvdOpenCvGSLCoherence(0.00001);
245  if(ret2)
246  std:: cout << "Success"<< std:: endl;
247  else
248  std:: cout << "Fail"<< std:: endl;
249 
250  std::cout << "--------------------------------------"<<std::endl ;
251 
252  L = Ls ;
253  t = vpTime::measureTimeMs() ;
254  L.svdFlake(W,V) ;
255  t = vpTime::measureTimeMs() -t ;
256  std::cout <<"svdFlake\n time " <<t << std::endl;
257  std::cout << W.t() ;
258 }
259