Visual Servoing Platform  version 3.2.0 under development (2019-01-22)
vpColVector.cpp
1 /****************************************************************************
2  *
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2019 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 http://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  * Provide some simple operation on column vectors.
33  *
34  * Authors:
35  * Eric Marchand
36  *
37  *****************************************************************************/
38 
45 #include <assert.h>
46 #include <cmath> // std::fabs
47 #include <limits> // numeric_limits
48 #include <math.h>
49 #include <sstream>
50 #include <stdio.h>
51 #include <stdlib.h>
52 #include <string.h>
53 
54 #include <visp3/core/vpCPUFeatures.h>
55 #include <visp3/core/vpColVector.h>
56 #include <visp3/core/vpDebug.h>
57 #include <visp3/core/vpException.h>
58 #include <visp3/core/vpMath.h>
59 #include <visp3/core/vpRotationVector.h>
60 
61 #if defined __SSE2__ || defined _M_X64 || (defined _M_IX86_FP && _M_IX86_FP >= 2)
62 #include <emmintrin.h>
63 #define VISP_HAVE_SSE2 1
64 #endif
65 
68 {
69  if (getRows() != v.getRows()) {
70  throw(vpException(vpException::dimensionError, "Cannot add (%dx1) column vector to (%dx1) column vector", getRows(),
71  v.getRows()));
72  }
74 
75  for (unsigned int i = 0; i < rowNum; i++)
76  r[i] = (*this)[i] + v[i];
77  return r;
78 }
101 {
102  if (getRows() != 3) {
103  throw(vpException(vpException::dimensionError, "Cannot add %d-dimension column vector to a translation vector",
104  getRows()));
105  }
107 
108  for (unsigned int i = 0; i < 3; i++)
109  s[i] = (*this)[i] + t[i];
110 
111  return s;
112 }
113 
116 {
117  if (getRows() != v.getRows()) {
118  throw(vpException(vpException::dimensionError, "Cannot add (%dx1) column vector to (%dx1) column vector", getRows(),
119  v.getRows()));
120  }
121 
122  for (unsigned int i = 0; i < rowNum; i++)
123  (*this)[i] += v[i];
124  return (*this);
125 }
128 {
129  if (getRows() != v.getRows()) {
130  throw(vpException(vpException::dimensionError, "Cannot substract (%dx1) column vector to (%dx1) column vector",
131  getRows(), v.getRows()));
132  }
133 
134  for (unsigned int i = 0; i < rowNum; i++)
135  (*this)[i] -= v[i];
136  return (*this);
137 }
138 
146 double vpColVector::operator*(const vpColVector &v) const
147 {
148  if (size() != v.size()) {
150  "Cannot compute the dot product between column vectors "
151  "with different dimensions (%d) and (%d)",
152  size(), v.size()));
153  }
154  double r = 0;
155 
156  for (unsigned int i = 0; i < rowNum; i++)
157  r += (*this)[i] * v[i];
158  return r;
159 }
160 
171 {
172  vpMatrix M(rowNum, v.getCols());
173  for (unsigned int i = 0; i < rowNum; i++) {
174  for (unsigned int j = 0; j < v.getCols(); j++) {
175  M[i][j] = (*this)[i] * v[j];
176  }
177  }
178  return M;
179 }
180 
183 {
184  if (getRows() != m.getRows()) {
186  "Bad size during vpColVector (%dx1) and vpColVector "
187  "(%dx1) substraction",
188  getRows(), m.getRows()));
189  }
190  vpColVector v(rowNum);
191 
192  for (unsigned int i = 0; i < rowNum; i++)
193  v[i] = (*this)[i] - m[i];
194  return v;
195 }
196 
210 vpColVector::vpColVector(const vpColVector &v, unsigned int r, unsigned int nrows) : vpArray2D<double>(nrows, 1)
211 {
212  init(v, r, nrows);
213 }
214 
252 void vpColVector::init(const vpColVector &v, unsigned int r, unsigned int nrows)
253 {
254  unsigned int rnrows = r + nrows;
255 
256  if (rnrows > v.getRows())
257  throw(vpException(vpException::dimensionError, "Bad row dimension (%d > %d) used to initialize vpColVector", rnrows,
258  v.getRows()));
259  resize(nrows, false);
260 
261  if (this->rowPtrs == NULL) // Fix coverity scan: explicit null dereferenced
262  return; // Nothing to do
263  for (unsigned int i = r; i < rnrows; i++)
264  (*this)[i - r] = v[i];
265 }
266 
268 {
269  for (unsigned int i = 0; i < v.size(); i++)
270  (*this)[i] = v[i];
271 }
272 
274 {
275  for (unsigned int i = 0; i < p.size(); i++)
276  (*this)[i] = p[i];
277 }
278 
280 {
281  for (unsigned int i = 0; i < v.size(); i++)
282  (*this)[i] = v[i];
283 }
284 
286 vpColVector::vpColVector(const vpMatrix &M, unsigned int j) : vpArray2D<double>(M.getRows(), 1)
287 {
288  for (unsigned int i = 0; i < M.getCols(); i++)
289  (*this)[i] = M[i][j];
290 }
291 
299 {
300  if (M.getCols() != 1) {
301  throw(vpException(vpException::dimensionError, "Cannot construct a (%dx1) row vector from a (%dx%d) matrix",
302  M.getRows(), M.getRows(), M.getCols()));
303  }
304 
305  for (unsigned int i = 0; i < M.getRows(); i++)
306  (*this)[i] = M[i][0];
307 }
308 
312 vpColVector::vpColVector(const std::vector<double> &v) : vpArray2D<double>((unsigned int)v.size(), 1)
313 {
314  for (unsigned int i = 0; i < v.size(); i++)
315  (*this)[i] = v[i];
316 }
320 vpColVector::vpColVector(const std::vector<float> &v) : vpArray2D<double>((unsigned int)v.size(), 1)
321 {
322  for (unsigned int i = 0; i < v.size(); i++)
323  (*this)[i] = (double)(v[i]);
324 }
325 
326 #ifdef VISP_HAVE_CPP11_COMPATIBILITY
328 {
329  rowNum = v.rowNum;
330  colNum = v.colNum;
331  rowPtrs = v.rowPtrs;
332  dsize = v.dsize;
333  data = v.data;
334 
335  v.rowNum = 0;
336  v.colNum = 0;
337  v.rowPtrs = NULL;
338  v.dsize = 0;
339  v.data = NULL;
340 }
341 #endif
342 
354 {
355  vpColVector A;
356  A.resize(rowNum, false);
357 
358  double *vd = A.data;
359  double *d = data;
360 
361  for (unsigned int i = 0; i < rowNum; i++)
362  *(vd++) = -(*d++);
363 
364  return A;
365 }
366 
388 {
389  vpColVector v(rowNum);
390 
391  double *vd = v.data;
392  double *d = data;
393 
394  for (unsigned int i = 0; i < rowNum; i++)
395  *(vd++) = (*d++) * x;
396  return v;
397 }
398 
418 {
419  for (unsigned int i = 0; i < rowNum; i++)
420  (*this)[i] *= x;
421  return (*this);
422 }
423 
442 {
443  for (unsigned int i = 0; i < rowNum; i++)
444  (*this)[i] /= x;
445  return (*this);
446 }
447 
468 {
469  vpColVector v(rowNum);
470 
471  double *vd = v.data;
472  double *d = data;
473 
474  for (unsigned int i = 0; i < rowNum; i++)
475  *(vd++) = (*d++) / x;
476  return v;
477 }
478 
485 {
486  if (M.getCols() != 1) {
487  throw(vpException(vpException::dimensionError, "Cannot transform a (%dx%d) matrix into a column vector",
488  M.getRows(), M.getCols()));
489  }
490 
491  resize(M.getRows(), false);
492  memcpy(data, M.data, rowNum * sizeof(double));
493 
494  return (*this);
495 }
496 
500 vpColVector &vpColVector::operator=(const std::vector<double> &v)
501 {
502  resize((unsigned int)v.size(), false);
503  for (unsigned int i = 0; i < v.size(); i++)
504  (*this)[i] = v[i];
505  return *this;
506 }
510 vpColVector &vpColVector::operator=(const std::vector<float> &v)
511 {
512  resize((unsigned int)v.size(), false);
513  for (unsigned int i = 0; i < v.size(); i++)
514  (*this)[i] = (float)v[i];
515  return *this;
516 }
517 
519 {
520  unsigned int k = v.rowNum;
521  if (rowNum != k) {
522  resize(k, false);
523  }
524 
525  memcpy(data, v.data, rowNum * sizeof(double));
526  return *this;
527 }
528 
533 {
534  unsigned int k = tv.getRows();
535  if (rowNum != k) {
536  resize(k, false);
537  }
538 
539  memcpy(data, tv.data, rowNum * sizeof(double));
540  return *this;
541 }
546 {
547  unsigned int k = rv.getRows();
548  if (rowNum != k) {
549  resize(k, false);
550  }
551 
552  memcpy(data, rv.data, rowNum * sizeof(double));
553  return *this;
554 }
559 {
560  unsigned int k = p.getRows();
561  if (rowNum != k) {
562  resize(k, false);
563  }
564 
565  memcpy(data, p.data, rowNum * sizeof(double));
566  return *this;
567 }
568 
590 {
591  *this = v;
592  return *this;
593 }
594 
620 {
621  for (unsigned int i = 0; i < rowNum; i++) {
622  for (unsigned int j = 0; j < colNum; j++) {
623  rowPtrs[i][j] = *x++;
624  }
625  }
626  return *this;
627 }
628 
631 {
632  double *d = data;
633 
634  for (unsigned int i = 0; i < rowNum; i++)
635  *(d++) = x;
636  return *this;
637 }
638 
643 std::vector<double> vpColVector::toStdVector()
644 {
645  std::vector<double> v(this->size());
646 
647  for (unsigned int i = 0; i < this->size(); i++)
648  v[i] = data[i];
649  return v;
650 }
651 
652 #ifdef VISP_HAVE_CPP11_COMPATIBILITY
654 {
655  if (this != &other) {
656  free(data);
657  free(rowPtrs);
658 
659  rowNum = other.rowNum;
660  colNum = other.colNum;
661  rowPtrs = other.rowPtrs;
662  dsize = other.dsize;
663  data = other.data;
664 
665  other.rowNum = 0;
666  other.colNum = 0;
667  other.rowPtrs = NULL;
668  other.dsize = 0;
669  other.data = NULL;
670  }
671 
672  return *this;
673 }
674 #endif
675 
676 bool vpColVector::operator==(const vpColVector &v) const {
677  if (rowNum != v.rowNum ||
678  colNum != v.colNum /* should not happen */)
679  return false;
680 
681  for (unsigned int i = 0; i < rowNum; i++) {
682  if (!vpMath::equal(data[i], v.data[i], std::numeric_limits<double>::epsilon()))
683  return false;
684  }
685 
686  return true;
687 }
688 
689 bool vpColVector::operator!=(const vpColVector &v) const {
690  return !(*this == v);
691 }
692 
697 {
698  vpRowVector v(rowNum);
699  memcpy(v.data, data, rowNum * sizeof(double));
700  return v;
701 }
702 
707 vpRowVector vpColVector::transpose() const { return t(); }
708 
713 void vpColVector::transpose(vpRowVector &v) const { v = t(); }
714 
719 vpColVector operator*(const double &x, const vpColVector &v)
720 {
721  vpColVector vout;
722  vout = v * x;
723  return vout;
724 }
725 
733 double vpColVector::dotProd(const vpColVector &a, const vpColVector &b)
734 {
735  if (a.data == NULL) {
736  throw(vpException(vpException::fatalError, "Cannot compute the dot product: first vector empty"));
737  }
738  if (b.data == NULL) {
739  throw(vpException(vpException::fatalError, "Cannot compute the dot product: second vector empty"));
740  }
741  if (a.size() != b.size()) {
743  "Cannot compute the dot product between column vectors "
744  "with different dimensions (%d) and (%d)",
745  a.size(), b.size()));
746  }
747 
748  double *ad = a.data;
749  double *bd = b.data;
750 
751  double c = 0;
752  for (unsigned int i = 0; i < a.getRows(); i++)
753  c += *(ad++) * *(bd++);
754  // vpMatrix c = (a.t() * b);
755  // return c[0][0];
756  return c;
757 }
758 
767 {
768  x = x / sqrt(x.sumSquare());
769 
770  return x;
771 }
772 
781 {
782 
783  double sum_square = sumSquare();
784 
785  // if (sum != 0.0)
786  if (std::fabs(sum_square) > std::numeric_limits<double>::epsilon())
787  *this /= sqrt(sum_square);
788 
789  // If sum = 0, we have a nul vector. So we return just.
790  return *this;
791 }
792 
822 {
823  if (v.data == NULL) {
824  throw(vpException(vpException::fatalError, "Cannot sort content of column vector: vector empty"));
825  }
826  vpColVector tab;
827  tab = v;
828  unsigned int nb_permutation = 1;
829  unsigned int i = 0;
830  while (nb_permutation != 0) {
831  nb_permutation = 0;
832  for (unsigned int j = v.getRows() - 1; j >= i + 1; j--) {
833  if ((tab[j] > tab[j - 1])) {
834  double tmp = tab[j];
835  tab[j] = tab[j - 1];
836  tab[j - 1] = tmp;
837  nb_permutation++;
838  }
839  }
840  i++;
841  }
842 
843  return tab;
844 }
845 
874 {
875  if (v.data == NULL) {
876  throw(vpException(vpException::fatalError, "Cannot sort content of column vector: vector empty"));
877  }
878  vpColVector tab;
879  tab = v;
880  unsigned int nb_permutation = 1;
881  unsigned int i = 0;
882  while (nb_permutation != 0) {
883  nb_permutation = 0;
884  for (unsigned int j = v.getRows() - 1; j >= i + 1; j--) {
885  if ((tab[j] < tab[j - 1])) {
886  double tmp = tab[j];
887  tab[j] = tab[j - 1];
888  tab[j - 1] = tmp;
889  nb_permutation++;
890  }
891  }
892  i++;
893  }
894 
895  return tab;
896 }
897 
914 void vpColVector::stack(double d)
915 {
916  this->resize(rowNum + 1, false);
917  (*this)[rowNum - 1] = d;
918 }
919 
939 void vpColVector::stack(const vpColVector &v) { *this = vpColVector::stack(*this, v); }
940 
960 {
961  vpColVector C;
962  vpColVector::stack(A, B, C);
963  return C;
964 }
965 
985 {
986  unsigned int nrA = A.getRows();
987  unsigned int nrB = B.getRows();
988 
989  if (nrA == 0 && nrB == 0) {
990  C.resize(0);
991  return;
992  }
993 
994  if (nrB == 0) {
995  C = A;
996  return;
997  }
998 
999  if (nrA == 0) {
1000  C = B;
1001  return;
1002  }
1003 
1004  // General case
1005  C.resize(nrA + nrB, false);
1006 
1007  for (unsigned int i = 0; i < nrA; i++)
1008  C[i] = A[i];
1009 
1010  for (unsigned int i = 0; i < nrB; i++)
1011  C[nrA + i] = B[i];
1012 }
1013 
1018 {
1019  if (v.data == NULL || v.size() == 0) {
1020  throw(vpException(vpException::dimensionError, "Cannot compute column vector mean: vector empty"));
1021  }
1022 
1023  // Use directly sum() function
1024  double mean = v.sum();
1025 
1026  // Old code used
1027  // double *vd = v.data;
1028  // for (unsigned int i=0 ; i < v.getRows() ; i++)
1029  // mean += *(vd++);
1030 
1031  return mean / v.getRows();
1032 }
1033 
1038 {
1039  if (v.data == NULL || v.size() == 0) {
1040  throw(vpException(vpException::dimensionError, "Cannot compute column vector median: vector empty"));
1041  }
1042 
1043  std::vector<double> vectorOfDoubles(v.data, v.data + v.rowNum);
1044 
1045  return vpMath::getMedian(vectorOfDoubles);
1046 }
1047 
1051 double vpColVector::stdev(const vpColVector &v, const bool useBesselCorrection)
1052 {
1053  if (v.data == NULL || v.size() == 0) {
1054  throw(vpException(vpException::dimensionError, "Cannot compute column vector stdev: vector empty"));
1055  }
1056 
1057  double mean_value = mean(v);
1058  double sum_squared_diff = 0.0;
1059  unsigned int i = 0;
1060 
1061 #if VISP_HAVE_SSE2
1062  if (vpCPUFeatures::checkSSE2()) {
1063  __m128d v_sub, v_mul, v_sum = _mm_setzero_pd();
1064  __m128d v_mean = _mm_set1_pd(mean_value);
1065 
1066  if (v.getRows() >= 4) {
1067  for (; i <= v.getRows() - 4; i += 4) {
1068  v_sub = _mm_sub_pd(_mm_loadu_pd(v.data + i), v_mean);
1069  v_mul = _mm_mul_pd(v_sub, v_sub);
1070  v_sum = _mm_add_pd(v_mul, v_sum);
1071 
1072  v_sub = _mm_sub_pd(_mm_loadu_pd(v.data + i + 2), v_mean);
1073  v_mul = _mm_mul_pd(v_sub, v_sub);
1074  v_sum = _mm_add_pd(v_mul, v_sum);
1075  }
1076  }
1077 
1078  double res[2];
1079  _mm_storeu_pd(res, v_sum);
1080 
1081  sum_squared_diff = res[0] + res[1];
1082  }
1083 // Old code used before SSE
1084 //#else
1085 // for(unsigned int i = 0; i < v.size(); i++) {
1086 // sum_squared_diff += (v[i]-mean_value) * (v[i]-mean_value);
1087 // }
1088 #endif
1089 
1090  for (; i < v.getRows(); i++) {
1091  sum_squared_diff += (v[i] - mean_value) * (v[i] - mean_value);
1092  }
1093 
1094  double divisor = (double)v.size();
1095  if (useBesselCorrection && v.size() > 1) {
1096  divisor = divisor - 1;
1097  }
1098 
1099  return std::sqrt(sum_squared_diff / divisor);
1100 }
1101 
1117 {
1118  vpMatrix M;
1119  if (v.getRows() != 3) {
1120  throw(vpException(vpException::dimensionError, "Cannot compute skew vector of a non 3-dimention vector (%d)",
1121  v.getRows()));
1122  }
1123 
1124  M.resize(3, 3, false, false);
1125  M[0][0] = 0;
1126  M[0][1] = -v[2];
1127  M[0][2] = v[1];
1128  M[1][0] = v[2];
1129  M[1][1] = 0;
1130  M[1][2] = -v[0];
1131  M[2][0] = -v[1];
1132  M[2][1] = v[0];
1133  M[2][2] = 0;
1134 
1135  return M;
1136 }
1137 
1149 {
1150  if (a.getRows() != 3 || b.getRows() != 3) {
1152  "Cannot compute the cross product between column "
1153  "vector with dimension %d and %d",
1154  a.getRows(), b.getRows()));
1155  }
1156 
1157  return vpColVector::skew(a) * b;
1158 }
1159 
1168 vpMatrix vpColVector::reshape(const unsigned int &nrows, const unsigned int &ncols)
1169 {
1170  vpMatrix M(nrows, ncols);
1171  reshape(M, nrows, ncols);
1172  return M;
1173 }
1174 
1230 void vpColVector::reshape(vpMatrix &M, const unsigned int &nrows, const unsigned int &ncols)
1231 {
1232  if (dsize != nrows * ncols) {
1233  throw(vpException(vpException::dimensionError, "Cannot reshape (%dx1) column vector in (%dx%d) matrix", rowNum,
1234  M.getRows(), M.getCols()));
1235  }
1236  if ((M.getRows() != nrows) || (M.getCols() != ncols))
1237  M.resize(nrows, ncols, false, false);
1238 
1239  for (unsigned int j = 0; j < ncols; j++)
1240  for (unsigned int i = 0; i < nrows; i++)
1241  M[i][j] = data[j * nrows + i];
1242 }
1243 
1276 void vpColVector::insert(unsigned int i, const vpColVector &v)
1277 {
1278  if (i + v.size() > this->size())
1279  throw(vpException(vpException::dimensionError, "Unable to insert a column vector"));
1280 
1281  if (data != NULL && v.data != NULL && v.rowNum > 0) {
1282  memcpy(data + i, v.data, sizeof(double) * v.rowNum);
1283  }
1284 }
1285 
1305 int vpColVector::print(std::ostream &s, unsigned int length, char const *intro) const
1306 {
1307  typedef std::string::size_type size_type;
1308 
1309  unsigned int m = getRows();
1310  unsigned int n = 1;
1311 
1312  std::vector<std::string> values(m * n);
1313  std::ostringstream oss;
1314  std::ostringstream ossFixed;
1315  std::ios_base::fmtflags original_flags = oss.flags();
1316 
1317  // ossFixed <<std::fixed;
1318  ossFixed.setf(std::ios::fixed, std::ios::floatfield);
1319 
1320  size_type maxBefore = 0; // the length of the integral part
1321  size_type maxAfter = 0; // number of decimals plus
1322  // one place for the decimal point
1323  for (unsigned int i = 0; i < m; ++i) {
1324  oss.str("");
1325  oss << (*this)[i];
1326  if (oss.str().find("e") != std::string::npos) {
1327  ossFixed.str("");
1328  ossFixed << (*this)[i];
1329  oss.str(ossFixed.str());
1330  }
1331 
1332  values[i] = oss.str();
1333  size_type thislen = values[i].size();
1334  size_type p = values[i].find('.');
1335 
1336  if (p == std::string::npos) {
1337  maxBefore = vpMath::maximum(maxBefore, thislen);
1338  // maxAfter remains the same
1339  } else {
1340  maxBefore = vpMath::maximum(maxBefore, p);
1341  maxAfter = vpMath::maximum(maxAfter, thislen - p - 1);
1342  }
1343  }
1344 
1345  size_type totalLength = length;
1346  // increase totalLength according to maxBefore
1347  totalLength = vpMath::maximum(totalLength, maxBefore);
1348  // decrease maxAfter according to totalLength
1349  maxAfter = (std::min)(maxAfter, totalLength - maxBefore);
1350  if (maxAfter == 1)
1351  maxAfter = 0;
1352 
1353  // the following line is useful for debugging
1354  // std::cerr <<totalLength <<" " <<maxBefore <<" " <<maxAfter <<"\n";
1355 
1356  if (intro)
1357  s << intro;
1358  s << "[" << m << "," << n << "]=\n";
1359 
1360  for (unsigned int i = 0; i < m; i++) {
1361  s << " ";
1362  size_type p = values[i].find('.');
1363  s.setf(std::ios::right, std::ios::adjustfield);
1364  s.width((std::streamsize)maxBefore);
1365  s << values[i].substr(0, p).c_str();
1366 
1367  if (maxAfter > 0) {
1368  s.setf(std::ios::left, std::ios::adjustfield);
1369  if (p != std::string::npos) {
1370  s.width((std::streamsize)maxAfter);
1371  s << values[i].substr(p, maxAfter).c_str();
1372  } else {
1373  assert(maxAfter > 1);
1374  s.width((std::streamsize)maxAfter);
1375  s << ".0";
1376  }
1377  }
1378 
1379  s << ' ';
1380 
1381  s << std::endl;
1382  }
1383 
1384  s.flags(original_flags); // restore s to standard state
1385 
1386  return (int)(maxBefore + maxAfter);
1387 }
1388 
1394 double vpColVector::sum() const
1395 {
1396  double sum = 0.0;
1397  unsigned int i = 0;
1398 
1399 #if VISP_HAVE_SSE2
1400  if (vpCPUFeatures::checkSSE2()) {
1401  __m128d v_sum1 = _mm_setzero_pd(), v_sum2 = _mm_setzero_pd(), v_sum;
1402 
1403  if (rowNum >= 4) {
1404  for (; i <= rowNum - 4; i += 4) {
1405  v_sum1 = _mm_add_pd(_mm_loadu_pd(data + i), v_sum1);
1406  v_sum2 = _mm_add_pd(_mm_loadu_pd(data + i + 2), v_sum2);
1407  }
1408  }
1409 
1410  v_sum = _mm_add_pd(v_sum1, v_sum2);
1411 
1412  double res[2];
1413  _mm_storeu_pd(res, v_sum);
1414 
1415  sum = res[0] + res[1];
1416  }
1417 // Old code used before SSE
1418 //#else
1419 // for (unsigned int i=0;i<rowNum;i++) {
1420 // sum += rowPtrs[i][0];
1421 // }
1422 #endif
1423 
1424  for (; i < rowNum; i++) {
1425  sum += (*this)[i];
1426  }
1427 
1428  return sum;
1429 }
1430 
1438 {
1439  double sum_square = 0.0;
1440  unsigned int i = 0;
1441 
1442 #if VISP_HAVE_SSE2
1443  if (vpCPUFeatures::checkSSE2()) {
1444  __m128d v_mul1, v_mul2;
1445  __m128d v_sum = _mm_setzero_pd();
1446 
1447  if (rowNum >= 4) {
1448  for (; i <= rowNum - 4; i += 4) {
1449  v_mul1 = _mm_mul_pd(_mm_loadu_pd(data + i), _mm_loadu_pd(data + i));
1450  v_mul2 = _mm_mul_pd(_mm_loadu_pd(data + i + 2), _mm_loadu_pd(data + i + 2));
1451 
1452  v_sum = _mm_add_pd(v_mul1, v_sum);
1453  v_sum = _mm_add_pd(v_mul2, v_sum);
1454  }
1455  }
1456 
1457  double res[2];
1458  _mm_storeu_pd(res, v_sum);
1459 
1460  sum_square = res[0] + res[1];
1461  }
1462 // Old code used before SSE
1463 //#else
1464 // for (unsigned int i=0;i<rowNum;i++) {
1465 // double x=rowPtrs[i][0];
1466 // sum_square += x*x;
1467 // }
1468 #endif
1469 
1470  for (; i < rowNum; i++) {
1471  sum_square += (*this)[i] * (*this)[i];
1472  }
1473 
1474  return sum_square;
1475 }
1476 
1484 {
1485  // Use directly sumSquare() function
1486  double norm = sumSquare();
1487 
1488  // Old code used
1489  // for (unsigned int i=0;i<dsize;i++) {
1490  // double x = *(data +i); norm += x*x;
1491  // }
1492 
1493  return sqrt(norm);
1494 }
1495 
1503 {
1504  if (v.getRows() != rowNum || v.getCols() != colNum) {
1505  throw(vpException(vpException::dimensionError, "Hadamard product: bad dimensions!"));
1506  }
1507 
1508  vpColVector out;
1509  out.resize(rowNum, false);
1510 
1511  unsigned int i = 0;
1512 
1513 #if VISP_HAVE_SSE2
1514  if (vpCPUFeatures::checkSSE2() && dsize >= 2) {
1515  for (; i <= dsize - 2; i += 2) {
1516  __m128d vout = _mm_mul_pd(_mm_loadu_pd(data + i), _mm_loadu_pd(v.data + i));
1517  _mm_storeu_pd(out.data + i, vout);
1518  }
1519  }
1520 #endif
1521 
1522  for (; i < dsize; i++) {
1523  out.data[i] = data[i] * v.data[i];
1524  }
1525 
1526  return out;
1527 }
1528 
1541 {
1542  double norm = 0.0;
1543  for (unsigned int i = 0; i < rowNum; i++) {
1544  double x = fabs((*this)[i]);
1545  if (x > norm) {
1546  norm = x;
1547  }
1548  }
1549  return norm;
1550 }
1551 
1580 std::ostream &vpColVector::cppPrint(std::ostream &os, const std::string &matrixName, bool octet) const
1581 {
1582  os << "vpColVector " << matrixName << " (" << this->getRows() << "); " << std::endl;
1583 
1584  for (unsigned int i = 0; i < this->getRows(); ++i) {
1585 
1586  if (!octet) {
1587  os << matrixName << "[" << i << "] = " << (*this)[i] << "; " << std::endl;
1588  } else {
1589  for (unsigned int k = 0; k < sizeof(double); ++k) {
1590  os << "((unsigned char*)&(" << matrixName << "[" << i << "]) )[" << k << "] = 0x" << std::hex
1591  << (unsigned int)((unsigned char *)&((*this)[i]))[k] << "; " << std::endl;
1592  }
1593  }
1594  }
1595  std::cout << std::endl;
1596  return os;
1597 };
1598 
1625 std::ostream &vpColVector::csvPrint(std::ostream &os) const
1626 {
1627  for (unsigned int i = 0; i < this->getRows(); ++i) {
1628  os << (*this)[i];
1629 
1630  os << std::endl;
1631  }
1632  return os;
1633 };
1634 
1660 std::ostream &vpColVector::maplePrint(std::ostream &os) const
1661 {
1662  os << "([ " << std::endl;
1663  for (unsigned int i = 0; i < this->getRows(); ++i) {
1664  os << "[";
1665  os << (*this)[i] << ", ";
1666  os << "]," << std::endl;
1667  }
1668  os << "])" << std::endl;
1669  return os;
1670 };
1671 
1708 std::ostream &vpColVector::matlabPrint(std::ostream &os) const
1709 {
1710  os << "[ ";
1711  for (unsigned int i = 0; i < this->getRows(); ++i) {
1712  os << (*this)[i] << ", ";
1713  if (this->getRows() != i + 1) {
1714  os << ";" << std::endl;
1715  } else {
1716  os << "]" << std::endl;
1717  }
1718  }
1719  return os;
1720 };
1721 
1722 #if defined(VISP_BUILD_DEPRECATED_FUNCTIONS)
1723 
1737 void vpColVector::insert(const vpColVector &v, const unsigned int r, const unsigned int c)
1738 {
1739  (void)c;
1740  insert(r, v);
1741 }
1742 #endif // defined(VISP_BUILD_DEPRECATED_FUNCTIONS)
bool operator==(const vpColVector &v) const
Comparison operator.
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:104
Implementation of a generic rotation vector.
vp_deprecated void init()
Definition: vpColVector.h:311
static vpColVector invSort(const vpColVector &v)
static vpColVector sort(const vpColVector &v)
static double stdev(const vpColVector &v, const bool useBesselCorrection=false)
double euclideanNorm() const
Implementation of row vector and the associated operations.
Definition: vpRowVector.h:72
std::ostream & matlabPrint(std::ostream &os) const
static double getMedian(const std::vector< double > &v)
Definition: vpMath.cpp:222
vpColVector operator*(const double &x, const vpColVector &v)
vpColVector operator/(const double x) const
void resize(const unsigned int nrows, const unsigned int ncols, const bool flagNullify=true, const bool recopy_=true)
Definition: vpArray2D.h:171
static bool equal(double x, double y, double s=0.001)
Definition: vpMath.h:290
vpColVector operator+(const vpColVector &v) const
Operator that allows to add two column vectors.
Definition: vpColVector.cpp:67
error that can be emited by ViSP classes.
Definition: vpException.h:71
std::ostream & csvPrint(std::ostream &os) const
double * data
Address of the first element of the data array.
Definition: vpArray2D.h:84
Implementation of a generic 2D array used as vase class of matrices and vectors.
Definition: vpArray2D.h:70
vpColVector & operator*=(double x)
unsigned int size() const
Return the number of elements of the 2D array.
Definition: vpArray2D.h:158
unsigned int getCols() const
Definition: vpArray2D.h:146
vpColVector operator-() const
static double median(const vpColVector &v)
static Type maximum(const Type &a, const Type &b)
Definition: vpMath.h:137
vpColVector & operator/=(double x)
int print(std::ostream &s, unsigned int length, char const *intro=0) const
bool operator!=(const vpColVector &v) const
unsigned int rowNum
Number of rows in the array.
Definition: vpArray2D.h:74
double infinityNorm() const
vpColVector & normalize()
vpColVector hadamard(const vpColVector &v) const
std::ostream & maplePrint(std::ostream &os) const
VISP_EXPORT bool checkSSE2()
vpRowVector t() const
std::ostream & cppPrint(std::ostream &os, const std::string &matrixName="A", bool octet=false) const
vpColVector & operator<<(const vpColVector &v)
static double mean(const vpColVector &v)
unsigned int getRows() const
Definition: vpArray2D.h:156
vpColVector & operator-=(vpColVector v)
Operator that allows to substract two column vectors.
unsigned int colNum
Number of columns in the array.
Definition: vpArray2D.h:76
void insert(unsigned int i, const vpColVector &v)
double sumSquare() const
std::vector< double > toStdVector()
vpColVector & operator=(const vpColVector &v)
Copy operator. Allow operation such as A = v.
vpColVector & operator+=(vpColVector v)
Operator that allows to add two column vectors.
Implementation of column vector and the associated operations.
Definition: vpColVector.h:72
static double dotProd(const vpColVector &a, const vpColVector &b)
void stack(double d)
vpRowVector transpose() const
Implementation of a pose vector and operations on poses.
Definition: vpPoseVector.h:92
vpColVector()
Basic constructor that creates an empty 0-size column vector.
Definition: vpColVector.h:78
static vpMatrix skew(const vpColVector &v)
unsigned int dsize
Current array size (rowNum * colNum)
Definition: vpArray2D.h:80
void reshape(vpMatrix &M, const unsigned int &nrows, const unsigned int &ncols)
static vpColVector crossProd(const vpColVector &a, const vpColVector &b)
double operator*(const vpColVector &x) const
Class that consider the case of a translation vector.
double ** rowPtrs
Address of the first element of each rows.
Definition: vpArray2D.h:78
double sum() const
void resize(const unsigned int i, const bool flagNullify=true)
Definition: vpColVector.h:244