ViSP  2.10.0
vpMomentObject.cpp
1 /****************************************************************************
2  *
3  * $Id: vpMomentObject.cpp 5301 2015-02-10 16:36:52Z mbakthav $
4  *
5  * This file is part of the ViSP software.
6  * Copyright (C) 2005 - 2014 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  * Object input structure used by moments.
36  *
37  * Authors:
38  * Filip Novotny
39  *
40  *****************************************************************************/
41 
42 #include <visp/vpMomentBasic.h>
43 #include <visp/vpMomentObject.h>
44 #include <visp/vpCameraParameters.h>
45 #include <visp/vpPixelMeterConversion.h>
46 #include <visp/vpConfig.h>
47 #include <stdexcept>
48 
49 #include <cmath>
50 #include <limits>
51 
52 #ifdef VISP_HAVE_OPENMP
53 #include <omp.h>
54 #endif
55 #include <cassert>
56 
66 double vpMomentObject::calc_mom_polygon(unsigned int p, unsigned int q, const std::vector<vpPoint>& points){
67  unsigned int i,k,l;
68  double den,mom,s;
69  double x_k;
70  double x_p_k;
71  double y_l;
72  double y_q_l;
73 
74  den = static_cast<double>( (p+q+2)*(p+q+1)*vpMath::comb(p+q,p) );
75 
76  mom = 0.0;
77  for (i=1;i<=points.size()-1;i++)
78  {
79  s = 0.0;
80  x_k=1;
81  for (k=0;k<=p;k++)
82  {
83  y_l=1;
84  x_p_k = pow(points[i-1].get_x(), (int)(p-k));
85  for (l=0;l<=q;l++)
86  {
87  y_q_l=pow(points[i-1].get_y(), (int)(q-l));
88 
89  s += static_cast<double>(
90  vpMath::comb(k+l,l)
91  *vpMath::comb(p+q-k-l,q-l)
92  *x_k
93  *x_p_k
94  *y_l
95  *y_q_l );
96 
97  y_l*=points[i].get_y();
98 
99  }
100  x_k*=points[i].get_x();
101 
102  }
103 
104  s *= ((points[i-1].get_x())*(points[i].get_y())-(points[i].get_x())*(points[i-1].get_y()));
105  mom += s;
106  }
107  mom /= den;
108  return(mom);
109 }
110 
124 void vpMomentObject::cacheValues(std::vector<double>& cache,double x, double y){
125  cache[0]=1;
126 
127  for(register unsigned int i=1;i<order;i++)
128  cache[i]=cache[i-1]*x;
129 
130  for(register unsigned int j=order;j<order*order;j+=order)
131  cache[j]=cache[j-order]*y;
132 
133  for(register unsigned int j=1;j<order;j++){
134  for(register unsigned int i=1;i<order-j;i++){
135  cache[j*order+i] = cache[j*order]*cache[i];
136  }
137  }
138 }
139 
144 void vpMomentObject::cacheValues(std::vector<double>& cache,double x, double y, double IntensityNormalized) {
145 
146  cache[0]=IntensityNormalized;
147 
148  double invIntensityNormalized = 0.;
149  if (std::fabs(IntensityNormalized)>=std::numeric_limits<double>::epsilon())
150  invIntensityNormalized = 1.0/IntensityNormalized;
151 
152  for(register unsigned int i=1;i<order;i++)
153  cache[i]=cache[i-1]*x;
154 
155  for(register unsigned int j=order;j<order*order;j+=order)
156  cache[j]=cache[j-order]*y;
157 
158  for(register unsigned int j=1;j<order;j++){
159  for(register unsigned int i=1;i<order-j;i++){
160  cache[j*order+i] = cache[j*order]*cache[i]*invIntensityNormalized;
161  }
162  }
163 }
164 
165 
238 void vpMomentObject::fromVector(std::vector<vpPoint>& points){
240  if(
241  std::abs(points.rbegin()->get_x()-points.begin()->get_x())>std::numeric_limits<double>::epsilon() ||
242  std::abs(points.rbegin()->get_y()-points.begin()->get_y())>std::numeric_limits<double>::epsilon()
243  ){
244  points.resize(points.size()+1);
245  points[points.size()-1] = points[0];
246  }
247  for(register unsigned int j=0;j<order*order;j++){
248  values[j]=calc_mom_polygon(j%order,j/order,points);
249  }
250  }else{
251  std::vector<double> cache(order*order,0.);
252  values.assign(order*order,0);
253  for(register unsigned int i=0;i<points.size();i++){
254  cacheValues(cache,points[i].get_x(),points[i].get_y());
255  for(register unsigned int j=0;j<order;j++){
256  for(register unsigned int k=0;k<order-j;k++){
257  values[j*order+k]+=cache[j*order+k];
258  }
259  }
260  }
261  }
262 }
263 
293 void vpMomentObject::fromImage(const vpImage<unsigned char>& image, unsigned char threshold, const vpCameraParameters& cam){
294 #ifdef VISP_HAVE_OPENMP
295  #pragma omp parallel shared(threshold)
296  {
297  std::vector<double> curvals(order*order);
298  curvals.assign(order*order,0.);
299  unsigned int i_, j_;
300 
301  #pragma omp for nowait//automatically organize loop counter between threads
302  for(int i=0;i<(int)image.getCols();i++){
303  for(int j=0;j<(int)image.getRows();j++){
304  i_ = static_cast<unsigned int>(i);
305  j_ = static_cast<unsigned int>(j);
306  if(image[j_][i_]>threshold){
307  double x=0;
308  double y=0;
310 
311  double xval=1.;
312  double yval=1.;
313  for(register unsigned int k=0;k<order;k++){
314  xval=1.;
315  for(register unsigned int l=0;l<order-k;l++){
316  curvals[(k*order+l)]+=(xval*yval);
317  xval*=x;
318  }
319  yval*=y;
320  }
321  }
322  }
323  }
324 
325  #pragma omp master //only set this variable in master thread
326  {
327  values.assign(order*order, 0.);
328  }
329 
330  #pragma omp barrier
331 
332  for(register unsigned int k=0;k<order;k++){
333  for(register unsigned int l=0;l<order-k;l++){
334  #pragma omp atomic
335  values[k*order+l]+= curvals[k*order+l];
336  }
337  }
338 
339  }
340 #else
341  std::vector<double> cache(order*order,0.);
342  values.assign(order*order,0);
343  for(register unsigned int i=0;i<image.getCols();i++){
344  for(register unsigned int j=0;j<image.getRows();j++){
345  if(image[j][i]>threshold){
346  double x=0;
347  double y=0;
349  cacheValues(cache,x,y);
350  for(register unsigned int k=0;k<order;k++){
351  for(register unsigned int l=0;l<order-k;l++){
352  values[k*order+l]+=cache[k*order+l];
353  }
354  }
355  }
356  }
357  }
358 #endif
359 
360  //Normalisation equivalent to sampling interval/pixel size delX x delY
361  double norm_factor = 1./(cam.get_px()*cam.get_py());
362  for (std::vector<double>::iterator it = values.begin(); it!=values.end(); it++) {
363  *it = (*it) * norm_factor;
364  }
365 }
366 
378  vpCameraImgBckGrndType bg_type, bool normalize_with_pix_size)
379 {
380  std::vector<double> cache(order*order,0.);
381  values.assign(order*order,0);
382 
383  // (x,y) - Pixel co-ordinates in metres
384  double x=0;
385  double y=0;
386  //for indexing into cache[] and values[]
387  unsigned int idx = 0;
388  unsigned int kidx = 0;
389 
390  double intensity = 0;
391  double intensity_white = 0;
392 
393  //double Imax = static_cast<double>(image.getMaxValue());
394  double Imax = 255.; // To check the effect of gray level change. ISR Coimbra
395 
396  double iscale = 1.0;
397  if (flg_normalize_intensity) // This makes the image a probability density function
398  iscale = 1.0/Imax;
399 
400  if (bg_type == vpMomentObject::WHITE) {
402  for(register unsigned int j=0;j<image.getRows();j++){
403  for(register unsigned int i=0;i<image.getCols();i++){
404  x = 0;
405  y = 0;
406  intensity = (double)(image[j][i])*iscale;
407  intensity_white = 1. - intensity;
408 
410  cacheValues(cache,x,y, intensity_white); // Modify 'cache' which has x^p*y^q to x^p*y^q*(1 - I(x,y))
411 
412  // Copy to "values"
413  for(register unsigned int k=0;k<order;k++){
414  kidx = k*order;
415  for(register unsigned int l=0;l<order-k;l++){
416  idx = kidx+l;
417  values[idx]+= cache[idx];
418  }
419  }
420  }
421  }
422  }
423  else {
425  for(register unsigned int j=0;j<image.getRows();j++){
426  for(register unsigned int i=0;i<image.getCols();i++){
427  x = 0;
428  y = 0;
429  intensity = (double)(image[j][i])*iscale;
431 
432  // Cache values for fast moment calculation
433  cacheValues(cache,x,y, intensity); // Modify 'cache' which has x^p*y^q to x^p*y^q*I(x,y)
434 
435  // Copy to moments array 'values'
436  for(register unsigned int k=0;k<order;k++){
437  kidx = k*order;
438  for(register unsigned int l=0;l<order-k;l++){
439  idx = kidx+l;
440  values[idx]+= cache[idx];
441  }
442  }
443 
444  }
445  }
446  }
447 
448  if (normalize_with_pix_size){
449  // Normalisation equivalent to sampling interval/pixel size delX x delY
450  double norm_factor = 1./(cam.get_px()*cam.get_py());
451  for (std::vector<double>::iterator it = values.begin(); it!=values.end(); it++) {
452  *it = (*it) * norm_factor;
453  }
454  }
455 }
456 
461 void
462 vpMomentObject::init(unsigned int orderinp) {
463  order = orderinp + 1;
465  flg_normalize_intensity = true; // By default, the intensity values are normalized
466  values.resize((order+1)*(order+1));
467  values.assign((order+1)*(order+1),0);
468 }
469 
473 void
475  order = objin.getOrder()+1;
476  type = objin.getType();
478  values.resize(objin.values.size());
479  values = objin.values;
480 }
481 
495 vpMomentObject::vpMomentObject(unsigned int max_order)
496  : flg_normalize_intensity(true), order(max_order+1), type(vpMomentObject::DENSE_FULL_OBJECT),
497  values()
498 {
499  init(max_order);
500 }
501 
506  : flg_normalize_intensity(true), order(1), type(vpMomentObject::DENSE_FULL_OBJECT),
507  values()
508 {
509  init(srcobj);
510 }
511 
531 const std::vector<double>& vpMomentObject::get() const {
532  return values;
533 
534 }
535 
542 double vpMomentObject::get(unsigned int i, unsigned int j) const {
543  assert(i+j<=getOrder());
544  if(i+j>=order) throw vpException(vpException::badValue,"The requested value has not been computed, you should specify a higher order.");
545 
546  return values[j*order+i];
547 }
548 
555 void vpMomentObject::set(unsigned int i, unsigned int j, const double& value_ij){
556  assert(i+j<=getOrder());
557  if(i+j>=order) throw vpException(vpException::badValue,"The requested value cannot be set, you should specify a higher order for the moment object.");
558  values[j*order+i] = value_ij;
559 }
560 
576 VISP_EXPORT std::ostream & operator<<(std::ostream & os, const vpMomentObject& m){
577  for(unsigned int i = 0;i<m.values.size();i++){
578 
579  if(i%(m.order)==0)
580  os << std::endl;
581 
582  if((i%(m.order)+i/(m.order))<m.order)
583  os << m.values[i] ;
584  else
585  os << "x";
586 
587  os << "\t";
588 
589  }
590 
591  return os;
592 }
593 
598 void
599 vpMomentObject::printWithIndices(const vpMomentObject& momobj, std::ostream& os) {
600  std::vector<double> moment = momobj.get();
601  os << std::endl <<"Order of vpMomentObject: "<<momobj.getOrder()<<std::endl;
602  // Print out values. This is same as printing using operator <<
603  for(unsigned int k=0; k<=momobj.getOrder(); k++) {
604  for(unsigned int l=0; l<(momobj.getOrder()+1)-k; l++){
605  os << "m[" << l << "," << k << "] = " << moment[k*(momobj.getOrder()+1)+ l] << "\t";
606  }
607  os << std::endl;
608  }
609  os <<std::endl;
610 }
611 
638 vpMatrix
640  std::vector<double> moment = momobj.get();
641  unsigned int order = momobj.getOrder();
642  vpMatrix M(order+1, order+1);
643  for(unsigned int k=0; k<=order; k++) {
644  for(unsigned int l=0; l<(order+1)-k; l++){
645  M[l][k] = moment[k*(order+1)+ l];
646  }
647  }
648  return M;
649 }
650 
659 // deliberate empty
660 }
unsigned int getCols() const
Definition: vpImage.h:180
Definition of the vpMatrix class.
Definition: vpMatrix.h:98
vpMomentObject(unsigned int order)
static vpMatrix convertTovpMatrix(const vpMomentObject &momobj)
static void printWithIndices(const vpMomentObject &momobj, std::ostream &os)
error that can be emited by ViSP classes.
Definition: vpException.h:76
Class for generic objects.
static void convertPoint(const vpCameraParameters &cam, const double &u, const double &v, double &x, double &y)
Point coordinates conversion from pixel coordinates to normalized coordinates in meter...
const std::vector< double > & get() const
double get_py() const
virtual ~vpMomentObject()
void fromImage(const vpImage< unsigned char > &image, unsigned char threshold, const vpCameraParameters &cam)
void init(unsigned int orderinp)
unsigned int getRows() const
Definition: vpImage.h:171
void cacheValues(std::vector< double > &cache, double x, double y)
void set(unsigned int i, unsigned int j, const double &value_ij)
Generic class defining intrinsic camera parameters.
bool flg_normalize_intensity
unsigned int order
void fromVector(std::vector< vpPoint > &points)
double get_px() const
std::vector< double > values
vpObjectType getType() const
vpObjectType type
static long double comb(unsigned int n, unsigned int p)
Definition: vpMath.h:213
unsigned int getOrder() const