Visual Servoing Platform  version 3.2.1 under development (2019-05-26)
vpImageTools.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  * Image tools.
33  *
34  * Authors:
35  * Fabien Spindler
36  *
37  *****************************************************************************/
38 
39 #include <visp3/core/vpCPUFeatures.h>
40 #include <visp3/core/vpImageConvert.h>
41 #include <visp3/core/vpImageTools.h>
42 
43 #if defined __SSE2__ || defined _M_X64 || (defined _M_IX86_FP && _M_IX86_FP >= 2)
44 #include <emmintrin.h>
45 #define VISP_HAVE_SSE2 1
46 
47 #if defined __SSE3__ || (defined _MSC_VER && _MSC_VER >= 1500)
48 #include <pmmintrin.h>
49 #define VISP_HAVE_SSE3 1
50 #endif
51 #if defined __SSSE3__ || (defined _MSC_VER && _MSC_VER >= 1500)
52 #include <tmmintrin.h>
53 #define VISP_HAVE_SSSE3 1
54 #endif
55 #endif
56 
113 void vpImageTools::changeLUT(vpImage<unsigned char> &I, unsigned char A, unsigned char A_star, unsigned char B,
114  unsigned char B_star)
115 {
116  // Test if input values are valid
117  if (B <= A) {
118  vpERROR_TRACE("Bad gray levels");
120  }
121  unsigned char v;
122 
123  double factor = (double)(B_star - A_star) / (double)(B - A);
124 
125  for (unsigned int i = 0; i < I.getHeight(); i++)
126  for (unsigned int j = 0; j < I.getWidth(); j++) {
127  v = I[i][j];
128 
129  if (v <= A)
130  I[i][j] = A_star;
131  else if (v >= B)
132  I[i][j] = B_star;
133  else
134  I[i][j] = (unsigned char)(A_star + factor * (v - A));
135  }
136 }
137 
151  vpImage<unsigned char> &Idiff)
152 {
153  if ((I1.getHeight() != I2.getHeight()) || (I1.getWidth() != I2.getWidth())) {
154  throw(vpException(vpException::dimensionError, "The two images have not the same size"));
155  }
156 
157  if ((I1.getHeight() != Idiff.getHeight()) || (I1.getWidth() != Idiff.getWidth()))
158  Idiff.resize(I1.getHeight(), I1.getWidth());
159 
160  bool checkSSSE3 = vpCPUFeatures::checkSSSE3();
161 #if !VISP_HAVE_SSSE3
162  checkSSSE3 = false;
163 #endif
164 
165  unsigned int i = 0;
166  if (checkSSSE3) {
167 #if VISP_HAVE_SSSE3
168  if (I1.getSize() >= 16) {
169  const __m128i mask1 = _mm_set_epi8(-1, 14, -1, 12, -1, 10, -1, 8, -1, 6, -1, 4, -1, 2, -1, 0);
170  const __m128i mask2 = _mm_set_epi8(-1, 15, -1, 13, -1, 11, -1, 9, -1, 7, -1, 5, -1, 3, -1, 1);
171 
172  const __m128i mask_out2 = _mm_set_epi8(14, -1, 12, -1, 10, -1, 8, -1, 6, -1, 4, -1, 2, -1, 0, -1);
173 
174  for (; i <= I1.getSize()-16; i+= 16) {
175  const __m128i vdata1 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(I1.bitmap + i));
176  const __m128i vdata2 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(I2.bitmap + i));
177 
178  __m128i vdata1_reorg = _mm_shuffle_epi8(vdata1, mask1);
179  __m128i vdata2_reorg = _mm_shuffle_epi8(vdata2, mask1);
180 
181  const __m128i vshift = _mm_set1_epi16(128);
182  __m128i vdata_diff = _mm_add_epi16(_mm_sub_epi16(vdata1_reorg, vdata2_reorg), vshift);
183 
184  const __m128i v255 = _mm_set1_epi16(255);
185  const __m128i vzero = _mm_setzero_si128();
186  const __m128i vdata_diff_min_max1 = _mm_max_epi16(_mm_min_epi16(vdata_diff, v255), vzero);
187 
188  vdata1_reorg = _mm_shuffle_epi8(vdata1, mask2);
189  vdata2_reorg = _mm_shuffle_epi8(vdata2, mask2);
190 
191  vdata_diff = _mm_add_epi16(_mm_sub_epi16(vdata1_reorg, vdata2_reorg), vshift);
192  const __m128i vdata_diff_min_max2 = _mm_max_epi16(_mm_min_epi16(vdata_diff, v255), vzero);
193 
194  _mm_storeu_si128(reinterpret_cast<__m128i *>(Idiff.bitmap + i), _mm_or_si128(_mm_shuffle_epi8(vdata_diff_min_max1, mask1),
195  _mm_shuffle_epi8(vdata_diff_min_max2, mask_out2)));
196  }
197  }
198 #endif
199  }
200 
201  for (; i < I1.getSize(); i++) {
202  int diff = I1.bitmap[i] - I2.bitmap[i] + 128;
203  Idiff.bitmap[i] = static_cast<unsigned char>(vpMath::maximum(vpMath::minimum(diff, 255), 0));
204  }
205 }
206 
221 {
222  if ((I1.getHeight() != I2.getHeight()) || (I1.getWidth() != I2.getWidth())) {
223  throw(vpException(vpException::dimensionError, "Cannot compute image difference. The two images "
224  "(%ux%u) and (%ux%u) have not the same size",
225  I1.getWidth(), I1.getHeight(), I2.getWidth(), I2.getHeight()));
226  }
227 
228  if ((I1.getHeight() != Idiff.getHeight()) || (I1.getWidth() != Idiff.getWidth()))
229  Idiff.resize(I1.getHeight(), I1.getWidth());
230 
231  bool checkSSSE3 = vpCPUFeatures::checkSSSE3();
232 #if !VISP_HAVE_SSSE3
233  checkSSSE3 = false;
234 #endif
235 
236  unsigned int i = 0;
237  if (checkSSSE3) {
238 #if VISP_HAVE_SSSE3
239  if (I1.getSize() >= 4) {
240  const __m128i mask1 = _mm_set_epi8(-1, 14, -1, 12, -1, 10, -1, 8, -1, 6, -1, 4, -1, 2, -1, 0);
241  const __m128i mask2 = _mm_set_epi8(-1, 15, -1, 13, -1, 11, -1, 9, -1, 7, -1, 5, -1, 3, -1, 1);
242 
243  const __m128i mask_out2 = _mm_set_epi8(14, -1, 12, -1, 10, -1, 8, -1, 6, -1, 4, -1, 2, -1, 0, -1);
244 
245  for (; i <= I1.getSize()-4; i+= 4) {
246  const __m128i vdata1 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(I1.bitmap + i));
247  const __m128i vdata2 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(I2.bitmap + i));
248 
249  __m128i vdata1_reorg = _mm_shuffle_epi8(vdata1, mask1);
250  __m128i vdata2_reorg = _mm_shuffle_epi8(vdata2, mask1);
251 
252  const __m128i vshift = _mm_set1_epi16(128);
253  __m128i vdata_diff = _mm_add_epi16(_mm_sub_epi16(vdata1_reorg, vdata2_reorg), vshift);
254 
255  const __m128i v255 = _mm_set1_epi16(255);
256  const __m128i vzero = _mm_setzero_si128();
257  const __m128i vdata_diff_min_max1 = _mm_max_epi16(_mm_min_epi16(vdata_diff, v255), vzero);
258 
259  vdata1_reorg = _mm_shuffle_epi8(vdata1, mask2);
260  vdata2_reorg = _mm_shuffle_epi8(vdata2, mask2);
261 
262  vdata_diff = _mm_add_epi16(_mm_sub_epi16(vdata1_reorg, vdata2_reorg), vshift);
263  const __m128i vdata_diff_min_max2 = _mm_max_epi16(_mm_min_epi16(vdata_diff, v255), vzero);
264 
265  _mm_storeu_si128(reinterpret_cast<__m128i *>(Idiff.bitmap + i), _mm_or_si128(_mm_shuffle_epi8(vdata_diff_min_max1, mask1),
266  _mm_shuffle_epi8(vdata_diff_min_max2, mask_out2)));
267  }
268  }
269 #endif
270  }
271 
272  for (; i < I1.getSize(); i++) {
273  int diffR = I1.bitmap[i].R - I2.bitmap[i].R + 128;
274  int diffG = I1.bitmap[i].G - I2.bitmap[i].G + 128;
275  int diffB = I1.bitmap[i].B - I2.bitmap[i].B + 128;
276  int diffA = I1.bitmap[i].A - I2.bitmap[i].A + 128;
277  Idiff.bitmap[i].R = static_cast<unsigned char>(vpMath::maximum(vpMath::minimum(diffR, 255), 0));
278  Idiff.bitmap[i].G = static_cast<unsigned char>(vpMath::maximum(vpMath::minimum(diffG, 255), 0));
279  Idiff.bitmap[i].B = static_cast<unsigned char>(vpMath::maximum(vpMath::minimum(diffB, 255), 0));
280  Idiff.bitmap[i].A = static_cast<unsigned char>(vpMath::maximum(vpMath::minimum(diffA, 255), 0));
281  }
282 }
283 
295  vpImage<unsigned char> &Idiff)
296 {
297  if ((I1.getHeight() != I2.getHeight()) || (I1.getWidth() != I2.getWidth())) {
298  throw(vpException(vpException::dimensionError, "The two images do not have the same size"));
299  }
300 
301  if ((I1.getHeight() != Idiff.getHeight()) || (I1.getWidth() != Idiff.getWidth()))
302  Idiff.resize(I1.getHeight(), I1.getWidth());
303 
304  unsigned int n = I1.getHeight() * I1.getWidth();
305  for (unsigned int b = 0; b < n; b++) {
306  int diff = I1.bitmap[b] - I2.bitmap[b];
307  Idiff.bitmap[b] = static_cast<unsigned char>(vpMath::abs(diff));
308  }
309 }
310 
319 {
320  if ((I1.getHeight() != I2.getHeight()) || (I1.getWidth() != I2.getWidth())) {
321  throw(vpException(vpException::dimensionError, "The two images do not have the same size"));
322  }
323 
324  if ((I1.getHeight() != Idiff.getHeight()) || (I1.getWidth() != Idiff.getWidth()))
325  Idiff.resize(I1.getHeight(), I1.getWidth());
326 
327  unsigned int n = I1.getHeight() * I1.getWidth();
328  for (unsigned int b = 0; b < n; b++) {
329  Idiff.bitmap[b] = vpMath::abs(I1.bitmap[b] - I2.bitmap[b]);
330  }
331 }
332 
347 {
348  if ((I1.getHeight() != I2.getHeight()) || (I1.getWidth() != I2.getWidth())) {
349  throw(vpException(vpException::dimensionError, "The two images do not have the same size"));
350  }
351 
352  if ((I1.getHeight() != Idiff.getHeight()) || (I1.getWidth() != Idiff.getWidth()))
353  Idiff.resize(I1.getHeight(), I1.getWidth());
354 
355  unsigned int n = I1.getHeight() * I1.getWidth();
356  for (unsigned int b = 0; b < n; b++) {
357  int diffR = I1.bitmap[b].R - I2.bitmap[b].R;
358  int diffG = I1.bitmap[b].G - I2.bitmap[b].G;
359  int diffB = I1.bitmap[b].B - I2.bitmap[b].B;
360  // int diffA = I1.bitmap[b].A - I2.bitmap[b].A;
361  Idiff.bitmap[b].R = static_cast<unsigned char>(vpMath::abs(diffR));
362  Idiff.bitmap[b].G = static_cast<unsigned char>(vpMath::abs(diffG));
363  Idiff.bitmap[b].B = static_cast<unsigned char>(vpMath::abs(diffB));
364  // Idiff.bitmap[b].A = diffA;
365  Idiff.bitmap[b].A = 0;
366  }
367 }
368 
379  vpImage<unsigned char> &Ires, const bool saturate)
380 {
381  if ((I1.getHeight() != I2.getHeight()) || (I1.getWidth() != I2.getWidth())) {
382  throw(vpException(vpException::dimensionError, "The two images do not have the same size"));
383  }
384 
385  if ((I1.getHeight() != Ires.getHeight()) || (I1.getWidth() != Ires.getWidth())) {
386  Ires.resize(I1.getHeight(), I1.getWidth());
387  }
388 
389  unsigned char *ptr_I1 = I1.bitmap;
390  unsigned char *ptr_I2 = I2.bitmap;
391  unsigned char *ptr_Ires = Ires.bitmap;
392  unsigned int cpt = 0;
393 
394 #if VISP_HAVE_SSE2
395  if (vpCPUFeatures::checkSSE2() && Ires.getSize() >= 16) {
396  for (; cpt <= Ires.getSize() - 16; cpt += 16, ptr_I1 += 16, ptr_I2 += 16, ptr_Ires += 16) {
397  const __m128i v1 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(ptr_I1));
398  const __m128i v2 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(ptr_I2));
399  const __m128i vres = saturate ? _mm_adds_epu8(v1, v2) : _mm_add_epi8(v1, v2);
400 
401  _mm_storeu_si128(reinterpret_cast<__m128i *>(ptr_Ires), vres);
402  }
403  }
404 #endif
405 
406  for (; cpt < Ires.getSize(); cpt++, ++ptr_I1, ++ptr_I2, ++ptr_Ires) {
407  *ptr_Ires = saturate ? vpMath::saturate<unsigned char>((short int)*ptr_I1 + (short int)*ptr_I2) : *ptr_I1 + *ptr_I2;
408  }
409 }
410 
421  vpImage<unsigned char> &Ires, const bool saturate)
422 {
423  if ((I1.getHeight() != I2.getHeight()) || (I1.getWidth() != I2.getWidth())) {
424  throw(vpException(vpException::dimensionError, "The two images do not have the same size"));
425  }
426 
427  if ((I1.getHeight() != Ires.getHeight()) || (I1.getWidth() != Ires.getWidth())) {
428  Ires.resize(I1.getHeight(), I1.getWidth());
429  }
430 
431  unsigned char *ptr_I1 = I1.bitmap;
432  unsigned char *ptr_I2 = I2.bitmap;
433  unsigned char *ptr_Ires = Ires.bitmap;
434  unsigned int cpt = 0;
435 
436 #if VISP_HAVE_SSE2
437  if (vpCPUFeatures::checkSSE2() && Ires.getSize() >= 16) {
438  for (; cpt <= Ires.getSize() - 16; cpt += 16, ptr_I1 += 16, ptr_I2 += 16, ptr_Ires += 16) {
439  const __m128i v1 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(ptr_I1));
440  const __m128i v2 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(ptr_I2));
441  const __m128i vres = saturate ? _mm_subs_epu8(v1, v2) : _mm_sub_epi8(v1, v2);
442 
443  _mm_storeu_si128(reinterpret_cast<__m128i *>(ptr_Ires), vres);
444  }
445  }
446 #endif
447 
448  for (; cpt < Ires.getSize(); cpt++, ++ptr_I1, ++ptr_I2, ++ptr_Ires) {
449  *ptr_Ires = saturate ?
450  vpMath::saturate<unsigned char>(static_cast<short int>(*ptr_I1) - static_cast<short int>(*ptr_I2)) :
451  *ptr_I1 - *ptr_I2;
452  }
453 }
454 
466 void vpImageTools::initUndistortMap(const vpCameraParameters &cam, unsigned int width, unsigned int height,
467  vpArray2D<int> &mapU, vpArray2D<int> &mapV,
468  vpArray2D<float> &mapDu, vpArray2D<float> &mapDv)
469 {
470  mapU.resize(height, width, false, false);
471  mapV.resize(height, width, false, false);
472  mapDu.resize(height, width, false, false);
473  mapDv.resize(height, width, false, false);
474 
475  float u0 = static_cast<float>(cam.get_u0());
476  float v0 = static_cast<float>(cam.get_v0());
477  float px = static_cast<float>(cam.get_px());
478  float py = static_cast<float>(cam.get_py());
479  float kud = static_cast<float>(cam.get_kud());
480 
481  if (std::fabs(static_cast<double>(kud)) <= std::numeric_limits<double>::epsilon()) {
482  // There is no need to undistort the image
483  for (unsigned int i = 0; i < height; i++) {
484  for (unsigned int j = 0; j < width; j++) {
485  mapU[i][j] = static_cast<int>(j);
486  mapV[i][j] = static_cast<int>(i);
487  mapDu[i][j] = 0;
488  mapDv[i][j] = 0;
489  }
490  }
491 
492  return;
493  }
494 
495  float invpx = 1.0f / px;
496  float invpy = 1.0f / py;
497 
498  float kud_px2 = kud * invpx * invpx;
499  float kud_py2 = kud * invpy * invpy;
500 
501  for (unsigned int v = 0; v < height; v++) {
502  float deltav = v - v0;
503  float fr1 = 1.0f + kud_py2 * deltav * deltav;
504 
505  for (unsigned int u = 0; u < width; u++) {
506  // computation of u,v : corresponding pixel coordinates in I.
507  float deltau = u - u0;
508  float fr2 = fr1 + kud_px2 * deltau * deltau;
509 
510  float u_float = deltau * fr2 + u0;
511  float v_float = deltav * fr2 + v0;
512 
513  int u_round = static_cast<int>(u_float);
514  int v_round = static_cast<int>(v_float);
515 
516  mapU[v][u] = u_round;
517  mapV[v][u] = v_round;
518 
519  mapDu[v][u] = u_float - u_round;
520  mapDv[v][u] = v_float - v_round;
521  }
522  }
523 }
524 
537 {
538  if (I.getSize() == 0) {
539  std::cerr << "Error, input image is empty." << std::endl;
540  return;
541  }
542 
543  II.resize(I.getHeight() + 1, I.getWidth() + 1, 0.0);
544  IIsq.resize(I.getHeight() + 1, I.getWidth() + 1, 0.0);
545 
546  for (unsigned int i = 1; i < II.getHeight(); i++) {
547  for (unsigned int j = 1; j < II.getWidth(); j++) {
548  II[i][j] = I[i - 1][j - 1] + II[i - 1][j] + II[i][j - 1] - II[i - 1][j - 1];
549  IIsq[i][j] = vpMath::sqr(I[i - 1][j - 1]) + IIsq[i - 1][j] + IIsq[i][j - 1] - IIsq[i - 1][j - 1];
550  }
551  }
552 }
553 
562 #if VISP_HAVE_SSE2
563  const bool useOptimized)
564 #else
565  const bool)
566 #endif
567 {
568  if ((I1.getHeight() != I2.getHeight()) || (I1.getWidth() != I2.getWidth())) {
569  throw vpException(vpException::dimensionError, "Error: in vpImageTools::normalizedCorrelation(): "
570  "image dimension mismatch between I1=%ux%u and I2=%ux%u",
571  I1.getHeight(), I1.getWidth(), I2.getHeight(), I2.getWidth());
572  }
573 
574  const double a = I1.getMeanValue();
575  const double b = I2.getMeanValue();
576 
577  double ab = 0.0;
578  double a2 = 0.0;
579  double b2 = 0.0;
580 
581  unsigned int cpt = 0;
582 
583 #if VISP_HAVE_SSE2
584  if (vpCPUFeatures::checkSSE2() && I1.getSize() >= 2 && useOptimized) {
585  const double *ptr_I1 = I1.bitmap;
586  const double *ptr_I2 = I2.bitmap;
587 
588  const __m128d v_mean_a = _mm_set1_pd(a);
589  const __m128d v_mean_b = _mm_set1_pd(b);
590  __m128d v_ab = _mm_setzero_pd();
591  __m128d v_a2 = _mm_setzero_pd();
592  __m128d v_b2 = _mm_setzero_pd();
593 
594  for (; cpt <= I1.getSize() - 2; cpt += 2, ptr_I1 += 2, ptr_I2 += 2) {
595  const __m128d v1 = _mm_loadu_pd(ptr_I1);
596  const __m128d v2 = _mm_loadu_pd(ptr_I2);
597  const __m128d norm_a = _mm_sub_pd(v1, v_mean_a);
598  const __m128d norm_b = _mm_sub_pd(v2, v_mean_b);
599  v_ab = _mm_add_pd(v_ab, _mm_mul_pd(norm_a, norm_b));
600  v_a2 = _mm_add_pd(v_a2, _mm_mul_pd(norm_a, norm_a));
601  v_b2 = _mm_add_pd(v_b2, _mm_mul_pd(norm_b, norm_b));
602  }
603 
604  double v_res_ab[2], v_res_a2[2], v_res_b2[2];
605  _mm_storeu_pd(v_res_ab, v_ab);
606  _mm_storeu_pd(v_res_a2, v_a2);
607  _mm_storeu_pd(v_res_b2, v_b2);
608 
609  ab = v_res_ab[0] + v_res_ab[1];
610  a2 = v_res_a2[0] + v_res_a2[1];
611  b2 = v_res_b2[0] + v_res_b2[1];
612  }
613 #endif
614 
615  for (; cpt < I1.getSize(); cpt++) {
616  ab += (I1.bitmap[cpt] - a) * (I2.bitmap[cpt] - b);
617  a2 += vpMath::sqr(I1.bitmap[cpt] - a);
618  b2 += vpMath::sqr(I2.bitmap[cpt] - b);
619  }
620 
621  return ab / sqrt(a2 * b2);
622 }
623 
632 {
633  unsigned int height = I.getHeight(), width = I.getWidth();
634  V.resize(width); // resize and nullify
635 
636  for (unsigned int i = 0; i < height; ++i)
637  for (unsigned int j = 0; j < width; ++j)
638  V[j] += I[i][j];
639  for (unsigned int j = 0; j < width; ++j)
640  V[j] /= height;
641 }
642 
648 {
649  double s = I.getSum();
650  for (unsigned int i = 0; i < I.getHeight(); ++i)
651  for (unsigned int j = 0; j < I.getWidth(); ++j)
652  I(i, j, I(i, j) / s);
653 }
654 
663  const vpImageInterpolationType &method)
664 {
665  switch (method) {
667  return I(vpMath::round(point.get_i()), vpMath::round(point.get_j()));
668  case INTERPOLATION_LINEAR: {
669  int x1 = (int)floor(point.get_i());
670  int x2 = (int)ceil(point.get_i());
671  int y1 = (int)floor(point.get_j());
672  int y2 = (int)ceil(point.get_j());
673  double v1, v2;
674  if (x1 == x2) {
675  v1 = I(x1, y1);
676  v2 = I(x1, y2);
677  } else {
678  v1 = (x2 - point.get_i()) * I(x1, y1) + (point.get_i() - x1) * I(x2, y1);
679  v2 = (x2 - point.get_i()) * I(x1, y2) + (point.get_i() - x1) * I(x2, y2);
680  }
681  if (y1 == y2)
682  return v1;
683  return (y2 - point.get_j()) * v1 + (point.get_j() - y1) * v2;
684  }
685  case INTERPOLATION_CUBIC: {
687  "vpImageTools::interpolate(): bi-cubic interpolation is not implemented.");
688  }
689  default: {
690  throw vpException(vpException::notImplementedError, "vpImageTools::interpolate(): invalid interpolation type");
691  }
692  }
693 }
694 
702 {
703  unsigned int x_d = vpMath::round(r.getHeight());
704  unsigned int y_d = vpMath::round(r.getWidth());
705  double x1 = r.getTopLeft().get_i();
706  double y1 = r.getTopLeft().get_j();
707  double t = r.getOrientation();
708  Dst.resize(x_d, y_d);
709  for (unsigned int x = 0; x < x_d; ++x) {
710  for (unsigned int y = 0; y < y_d; ++y) {
711  Dst(x, y,
712  (unsigned char)interpolate(Src, vpImagePoint(x1 + x * cos(t) + y * sin(t), y1 - x * sin(t) + y * cos(t)),
714  }
715  }
716 }
717 
725 {
726  unsigned int x_d = vpMath::round(r.getHeight());
727  unsigned int y_d = vpMath::round(r.getWidth());
728  double x1 = r.getTopLeft().get_i();
729  double y1 = r.getTopLeft().get_j();
730  double t = r.getOrientation();
731  Dst.resize(x_d, y_d);
732  for (unsigned int x = 0; x < x_d; ++x) {
733  for (unsigned int y = 0; y < y_d; ++y) {
734  Dst(x, y, interpolate(Src, vpImagePoint(x1 + x * cos(t) + y * sin(t), y1 - x * sin(t) + y * cos(t)),
736  }
737  }
738 }
739 
755  vpImage<double> &I_score, const unsigned int step_u, const unsigned int step_v,
756  const bool useOptimized)
757 {
758  if (I.getSize() == 0) {
759  std::cerr << "Error, input image is empty." << std::endl;
760  return;
761  }
762 
763  if (I_tpl.getSize() == 0) {
764  std::cerr << "Error, template image is empty." << std::endl;
765  return;
766  }
767 
768  if (I_tpl.getHeight() > I.getHeight() || I_tpl.getWidth() > I.getWidth()) {
769  std::cerr << "Error, template image is bigger than input image." << std::endl;
770  return;
771  }
772 
773  vpImage<double> I_double, I_tpl_double;
774  vpImageConvert::convert(I, I_double);
775  vpImageConvert::convert(I_tpl, I_tpl_double);
776 
777  const unsigned int height_tpl = I_tpl.getHeight(), width_tpl = I_tpl.getWidth();
778  I_score.resize(I.getHeight() - height_tpl, I.getWidth() - width_tpl, 0.0);
779 
780  if (useOptimized) {
781  vpImage<double> II, IIsq;
782  integralImage(I, II, IIsq);
783 
784  vpImage<double> II_tpl, IIsq_tpl;
785  integralImage(I_tpl, II_tpl, IIsq_tpl);
786 
787  // zero-mean template image
788  const double sum2 = (II_tpl[height_tpl][width_tpl] + II_tpl[0][0] - II_tpl[0][width_tpl] - II_tpl[height_tpl][0]);
789  const double mean2 = sum2 / I_tpl.getSize();
790  for (unsigned int cpt = 0; cpt < I_tpl_double.getSize(); cpt++) {
791  I_tpl_double.bitmap[cpt] -= mean2;
792  }
793 
794 #if defined _OPENMP && _OPENMP >= 200711 // OpenMP 3.1
795 #pragma omp parallel for schedule(dynamic)
796  for (unsigned int i = 0; i < I.getHeight() - height_tpl; i += step_v) {
797  for (unsigned int j = 0; j < I.getWidth() - width_tpl; j += step_u) {
798  I_score[i][j] = normalizedCorrelation(I_double, I_tpl_double, II, IIsq, II_tpl, IIsq_tpl, i, j);
799  }
800  }
801 #else
802  // error C3016: 'i': index variable in OpenMP 'for' statement must have signed integral type
803  int end = (int)((I.getHeight() - height_tpl) / step_v) + 1;
804  std::vector<unsigned int> vec_step_v((size_t)end);
805  for (unsigned int cpt = 0, idx = 0; cpt < I.getHeight() - height_tpl; cpt += step_v, idx++) {
806  vec_step_v[(size_t)idx] = cpt;
807  }
808 #if defined _OPENMP // only to disable warning: ignoring #pragma omp parallel [-Wunknown-pragmas]
809 #pragma omp parallel for schedule(dynamic)
810 #endif
811  for (int cpt = 0; cpt < end; cpt++) {
812  for (unsigned int j = 0; j < I.getWidth() - width_tpl; j += step_u) {
813  I_score[vec_step_v[cpt]][j] =
814  normalizedCorrelation(I_double, I_tpl_double, II, IIsq, II_tpl, IIsq_tpl, vec_step_v[cpt], j);
815  }
816  }
817 #endif
818  } else {
819  vpImage<double> I_cur;
820 
821  for (unsigned int i = 0; i < I.getHeight() - height_tpl; i += step_v) {
822  for (unsigned int j = 0; j < I.getWidth() - width_tpl; j += step_u) {
823  vpRect roi(vpImagePoint(i, j), vpImagePoint(i + height_tpl - 1, j + width_tpl - 1));
824  vpImageTools::crop(I_double, roi, I_cur);
825 
826  I_score[i][j] = vpImageTools::normalizedCorrelation(I_cur, I_tpl_double, useOptimized);
827  }
828  }
829  }
830 }
831 
832 // Reference:
833 // http://blog.demofox.org/2015/08/15/resizing-images-with-bicubic-interpolation/
834 // t is a value that goes from 0 to 1 to interpolate in a C1 continuous way
835 // across uniformly sampled data points. when t is 0, this will return B.
836 // When t is 1, this will return C. In between values will return an
837 // interpolation between B and C. A and B are used to calculate the slopes at
838 // the edges.
839 float vpImageTools::cubicHermite(const float A, const float B, const float C, const float D, const float t)
840 {
841  float a = (-A + 3.0f * B - 3.0f * C + D) / 2.0f;
842  float b = A + 2.0f * C - (5.0f * B + D) / 2.0f;
843  float c = (-A + C) / 2.0f;
844  float d = B;
845 
846  return a * t * t * t + b * t * t + c * t + d;
847 }
848 
849 float vpImageTools::lerp(const float A, const float B, const float t) {
850  return A * (1.0f - t) + B * t;
851 }
852 
853 int64_t vpImageTools::lerp2(int64_t A, int64_t B, int64_t t, int64_t t_1) {
854  return A * t_1 + B * t;
855 }
856 
858  const vpImage<double> &II, const vpImage<double> &IIsq,
859  const vpImage<double> &II_tpl, const vpImage<double> &IIsq_tpl,
860  const unsigned int i0, const unsigned int j0)
861 {
862  double ab = 0.0;
863 #if VISP_HAVE_SSE2
864  bool use_sse_version = true;
865  if (vpCPUFeatures::checkSSE2() && I2.getWidth() >= 2) {
866  const double *ptr_I1 = I1.bitmap;
867  const double *ptr_I2 = I2.bitmap;
868 
869  __m128d v_ab = _mm_setzero_pd();
870 
871  for (unsigned int i = 0; i < I2.getHeight(); i++) {
872  unsigned int j = 0;
873  ptr_I1 = &I1.bitmap[(i0 + i) * I1.getWidth() + j0];
874 
875  for (; j <= I2.getWidth() - 2; j += 2, ptr_I1 += 2, ptr_I2 += 2) {
876  const __m128d v1 = _mm_loadu_pd(ptr_I1);
877  const __m128d v2 = _mm_loadu_pd(ptr_I2);
878  v_ab = _mm_add_pd(v_ab, _mm_mul_pd(v1, v2));
879  }
880 
881  for (; j < I2.getWidth(); j++) {
882  ab += (I1[i0 + i][j0 + j]) * I2[i][j];
883  }
884  }
885 
886  double v_res_ab[2];
887  _mm_storeu_pd(v_res_ab, v_ab);
888 
889  ab += v_res_ab[0] + v_res_ab[1];
890  } else {
891  use_sse_version = false;
892  }
893 #else
894  bool use_sse_version = false;
895 #endif
896 
897  if (!use_sse_version) {
898  for (unsigned int i = 0; i < I2.getHeight(); i++) {
899  for (unsigned int j = 0; j < I2.getWidth(); j++) {
900  ab += (I1[i0 + i][j0 + j]) * I2[i][j];
901  }
902  }
903  }
904 
905  const unsigned int height_tpl = I2.getHeight(), width_tpl = I2.getWidth();
906  const double sum1 =
907  (II[i0 + height_tpl][j0 + width_tpl] + II[i0][j0] - II[i0][j0 + width_tpl] - II[i0 + height_tpl][j0]);
908  const double sum2 = (II_tpl[height_tpl][width_tpl] + II_tpl[0][0] - II_tpl[0][width_tpl] - II_tpl[height_tpl][0]);
909 
910  double a2 = ((IIsq[i0 + I2.getHeight()][j0 + I2.getWidth()] + IIsq[i0][j0] - IIsq[i0][j0 + I2.getWidth()] -
911  IIsq[i0 + I2.getHeight()][j0]) -
912  (1.0 / I2.getSize()) * vpMath::sqr(sum1));
913 
914  double b2 = ((IIsq_tpl[I2.getHeight()][I2.getWidth()] + IIsq_tpl[0][0] - IIsq_tpl[0][I2.getWidth()] -
915  IIsq_tpl[I2.getHeight()][0]) -
916  (1.0 / I2.getSize()) * vpMath::sqr(sum2));
917  return ab / sqrt(a2 * b2);
918 }
919 
931  const vpArray2D<float> &mapDu, const vpArray2D<float> &mapDv, vpImage<unsigned char> &Iundist)
932 {
933  Iundist.resize(I.getHeight(), I.getWidth());
934 
935 #if defined _OPENMP // only to disable warning: ignoring #pragma omp parallel [-Wunknown-pragmas]
936 #pragma omp parallel for schedule(dynamic)
937 #endif
938  for (int i_ = 0; i_ < static_cast<int>(I.getHeight()); i_++) {
939  const unsigned int i = static_cast<unsigned int>(i_);
940  for (unsigned int j = 0; j < I.getWidth(); j++) {
941 
942  int u_round = mapU[i][j];
943  int v_round = mapV[i][j];
944 
945  float du = mapDu[i][j];
946  float dv = mapDv[i][j];
947 
948  if (0 <= u_round && 0 <= v_round && u_round < static_cast<int>(I.getWidth()) - 1
949  && v_round < static_cast<int>(I.getHeight()) - 1) {
950  // process interpolation
951  float col0 = lerp(I[v_round][u_round], I[v_round][u_round + 1], du);
952  float col1 = lerp(I[v_round + 1][u_round], I[v_round + 1][u_round + 1], du);
953  float value = lerp(col0, col1, dv);
954 
955  Iundist[i][j] = static_cast<unsigned char>(value);
956  } else {
957  Iundist[i][j] = 0;
958  }
959  }
960  }
961 }
962 
973 void vpImageTools::remap(const vpImage<vpRGBa> &I, const vpArray2D<int> &mapU, const vpArray2D<int> &mapV,
974  const vpArray2D<float> &mapDu, const vpArray2D<float> &mapDv, vpImage<vpRGBa> &Iundist)
975 {
976  Iundist.resize(I.getHeight(), I.getWidth());
977 
978  bool checkSSE2 = vpCPUFeatures::checkSSE2();
979 #if !VISP_HAVE_SSE2
980  checkSSE2 = false;
981 #endif
982 
983  if (checkSSE2) {
984 #if defined VISP_HAVE_SSE2
985 #if defined _OPENMP // only to disable warning: ignoring #pragma omp parallel [-Wunknown-pragmas]
986 #pragma omp parallel for schedule(dynamic)
987 #endif
988  for (int i_ = 0; i_ < static_cast<int>(I.getHeight()); i_++) {
989  const unsigned int i = static_cast<unsigned int>(i_);
990  for (unsigned int j = 0; j < I.getWidth(); j++) {
991 
992  int u_round = mapU[i][j];
993  int v_round = mapV[i][j];
994 
995  const __m128 vdu = _mm_set1_ps(mapDu[i][j]);
996  const __m128 vdv = _mm_set1_ps(mapDv[i][j]);
997 
998  if (0 <= u_round && 0 <= v_round && u_round < static_cast<int>(I.getWidth()) - 1
999  && v_round < static_cast<int>(I.getHeight()) - 1) {
1000  #define VLERP(va, vb, vt) _mm_add_ps(va, _mm_mul_ps(_mm_sub_ps(vb, va), vt));
1001 
1002  // process interpolation
1003  const __m128 vdata1 =
1004  _mm_set_ps(static_cast<float>(I[v_round][u_round].A), static_cast<float>(I[v_round][u_round].B),
1005  static_cast<float>(I[v_round][u_round].G), static_cast<float>(I[v_round][u_round].R));
1006 
1007  const __m128 vdata2 =
1008  _mm_set_ps(static_cast<float>(I[v_round][u_round + 1].A), static_cast<float>(I[v_round][u_round + 1].B),
1009  static_cast<float>(I[v_round][u_round + 1].G), static_cast<float>(I[v_round][u_round + 1].R));
1010 
1011  const __m128 vdata3 =
1012  _mm_set_ps(static_cast<float>(I[v_round + 1][u_round].A), static_cast<float>(I[v_round + 1][u_round].B),
1013  static_cast<float>(I[v_round + 1][u_round].G), static_cast<float>(I[v_round + 1][u_round].R));
1014 
1015  const __m128 vdata4 = _mm_set_ps(
1016  static_cast<float>(I[v_round + 1][u_round + 1].A), static_cast<float>(I[v_round + 1][u_round + 1].B),
1017  static_cast<float>(I[v_round + 1][u_round + 1].G), static_cast<float>(I[v_round + 1][u_round + 1].R));
1018 
1019  const __m128 vcol0 = VLERP(vdata1, vdata2, vdu);
1020  const __m128 vcol1 = VLERP(vdata3, vdata4, vdu);
1021  const __m128 vvalue = VLERP(vcol0, vcol1, vdv);
1022 
1023  #undef VLERP
1024 
1025  float values[4];
1026  _mm_storeu_ps(values, vvalue);
1027  Iundist[i][j].R = static_cast<unsigned char>(values[0]);
1028  Iundist[i][j].G = static_cast<unsigned char>(values[1]);
1029  Iundist[i][j].B = static_cast<unsigned char>(values[2]);
1030  Iundist[i][j].A = static_cast<unsigned char>(values[3]);
1031  } else {
1032  Iundist[i][j] = 0;
1033  }
1034  }
1035  }
1036 #endif
1037  } else {
1038 #if defined _OPENMP // only to disable warning: ignoring #pragma omp parallel [-Wunknown-pragmas]
1039 #pragma omp parallel for schedule(dynamic)
1040 #endif
1041  for (int i_ = 0; i_ < static_cast<int>(I.getHeight()); i_++) {
1042  const unsigned int i = static_cast<unsigned int>(i_);
1043  for (unsigned int j = 0; j < I.getWidth(); j++) {
1044 
1045  int u_round = mapU[i][j];
1046  int v_round = mapV[i][j];
1047 
1048  float du = mapDu[i][j];
1049  float dv = mapDv[i][j];
1050 
1051  if (0 <= u_round && 0 <= v_round && u_round < static_cast<int>(I.getWidth()) - 1
1052  && v_round < static_cast<int>(I.getHeight()) - 1) {
1053  // process interpolation
1054  float col0 = lerp(I[v_round][u_round].R, I[v_round][u_round + 1].R, du);
1055  float col1 = lerp(I[v_round + 1][u_round].G, I[v_round + 1][u_round + 1].G, du);
1056  float value = lerp(col0, col1, dv);
1057 
1058  Iundist[i][j].R = static_cast<unsigned char>(value);
1059 
1060  col0 = lerp(I[v_round][u_round].G, I[v_round][u_round + 1].G, du);
1061  col1 = lerp(I[v_round + 1][u_round].G, I[v_round + 1][u_round + 1].G, du);
1062  value = lerp(col0, col1, dv);
1063 
1064  Iundist[i][j].G = static_cast<unsigned char>(value);
1065 
1066  col0 = lerp(I[v_round][u_round].B, I[v_round][u_round + 1].B, du);
1067  col1 = lerp(I[v_round + 1][u_round].B, I[v_round + 1][u_round + 1].B, du);
1068  value = lerp(col0, col1, dv);
1069 
1070  Iundist[i][j].B = static_cast<unsigned char>(value);
1071 
1072  col0 = lerp(I[v_round][u_round].A, I[v_round][u_round + 1].A, du);
1073  col1 = lerp(I[v_round + 1][u_round].A, I[v_round + 1][u_round + 1].A, du);
1074  value = lerp(col0, col1, dv);
1075 
1076  Iundist[i][j].A = static_cast<unsigned char>(value);
1077  } else {
1078  Iundist[i][j] = 0;
1079  }
1080  }
1081  }
1082  }
1083 }
static void extract(const vpImage< unsigned char > &Src, vpImage< unsigned char > &Dst, const vpRectOriented &r)
static void imageDifferenceAbsolute(const vpImage< unsigned char > &I1, const vpImage< unsigned char > &I2, vpImage< unsigned char > &Idiff)
double get_u0() const
double get_i() const
Definition: vpImagePoint.h:204
unsigned int getWidth() const
Definition: vpImage.h:237
static void convert(const vpImage< unsigned char > &src, vpImage< vpRGBa > &dest)
unsigned char B
Blue component.
Definition: vpRGBa.h:150
Implementation of row vector and the associated operations.
Definition: vpRowVector.h:115
Type * bitmap
points toward the bitmap
Definition: vpImage.h:134
#define vpERROR_TRACE
Definition: vpDebug.h:393
void resize(const unsigned int nrows, const unsigned int ncols, const bool flagNullify=true, const bool recopy_=true)
Definition: vpArray2D.h:305
error that can be emited by ViSP classes.
Definition: vpException.h:71
Implementation of a generic 2D array used as base class for matrices and vectors. ...
Definition: vpArray2D.h:131
double get_py() const
double getHeight() const
Get the rectangle height.
Error that can be emited by the vpImage class and its derivates.
static void normalize(vpImage< double > &I)
unsigned char G
Green component.
Definition: vpRGBa.h:149
static void imageDifference(const vpImage< unsigned char > &I1, const vpImage< unsigned char > &I2, vpImage< unsigned char > &Idiff)
double getWidth() const
Get the rectangle width.
static int round(const double x)
Definition: vpMath.h:241
double get_j() const
Definition: vpImagePoint.h:215
static double interpolate(const vpImage< unsigned char > &I, const vpImagePoint &point, const vpImageInterpolationType &method=INTERPOLATION_NEAREST)
static Type abs(const Type &x)
Definition: vpMath.h:158
double getSum() const
Definition: vpImage.h:1609
static Type maximum(const Type &a, const Type &b)
Definition: vpMath.h:143
double getOrientation() const
Get the rectangle orientation (rad).
static void remap(const vpImage< unsigned char > &I, const vpArray2D< int > &mapU, const vpArray2D< int > &mapV, const vpArray2D< float > &mapDu, const vpArray2D< float > &mapDv, vpImage< unsigned char > &Iundist)
double get_v0() const
unsigned int getSize() const
Definition: vpImage.h:218
static double sqr(double x)
Definition: vpMath.h:114
VISP_EXPORT bool checkSSE2()
unsigned char A
Additionnal component.
Definition: vpRGBa.h:151
Generic class defining intrinsic camera parameters.
void resize(const unsigned int h, const unsigned int w)
resize the image : Image initialization
Definition: vpImage.h:863
static Type minimum(const Type &a, const Type &b)
Definition: vpMath.h:151
VISP_EXPORT bool checkSSSE3()
double get_px() const
void resize(const unsigned int i, const bool flagNullify=true)
Definition: vpRowVector.h:268
double get_kud() const
static void initUndistortMap(const vpCameraParameters &cam, unsigned int width, unsigned int height, vpArray2D< int > &mapU, vpArray2D< int > &mapV, vpArray2D< float > &mapDu, vpArray2D< float > &mapDv)
static void columnMean(const vpImage< double > &I, vpRowVector &result)
static void integralImage(const vpImage< unsigned char > &I, vpImage< double > &II, vpImage< double > &IIsq)
static void imageAdd(const vpImage< unsigned char > &I1, const vpImage< unsigned char > &I2, vpImage< unsigned char > &Ires, const bool saturate=false)
static void crop(const vpImage< Type > &I, double roi_top, double roi_left, unsigned int roi_height, unsigned int roi_width, vpImage< Type > &crop, unsigned int v_scale=1, unsigned int h_scale=1)
Definition: vpImageTools.h:282
vpImagePoint getTopLeft() const
Get the top-left corner.
static void imageSubtract(const vpImage< unsigned char > &I1, const vpImage< unsigned char > &I2, vpImage< unsigned char > &Ires, const bool saturate=false)
unsigned char R
Red component.
Definition: vpRGBa.h:148
Type getMeanValue() const
Return the mean value of the bitmap.
Definition: vpImage.h:966
unsigned int getHeight() const
Definition: vpImage.h:179
Defines a rectangle in the plane.
Definition: vpRect.h:78
Class that defines a 2D point in an image. This class is useful for image processing and stores only ...
Definition: vpImagePoint.h:88
static void changeLUT(vpImage< unsigned char > &I, unsigned char A, unsigned char newA, unsigned char B, unsigned char newB)
Defines an oriented rectangle in the plane.
static void templateMatching(const vpImage< unsigned char > &I, const vpImage< unsigned char > &I_tpl, vpImage< double > &I_score, const unsigned int step_u, const unsigned int step_v, const bool useOptimized=true)
static double normalizedCorrelation(const vpImage< double > &I1, const vpImage< double > &I2, const bool useOptimized=true)