40 #include <visp3/core/vpException.h> 41 #include <visp3/tt_mi/vpTemplateTrackerMI.h> 42 #include <visp3/tt_mi/vpTemplateTrackerMIBSpline.h> 106 Pt =
new double[
Ncb];
107 Pr =
new double[
Ncb];
142 Pt =
new double[
Ncb];
143 Pr =
new double[
Ncb];
154 unsigned int Ncb_ = (
unsigned int)
Ncb;
155 unsigned int Nc_ = (
unsigned int)
Nc;
156 unsigned int influBspline_ = (
unsigned int)
influBspline;
158 memset(
Prt, 0, Ncb_ * Ncb_ *
sizeof(
double));
159 memset(
PrtD, 0, Nc_ * Nc_ * influBspline_ *
sizeof(
double));
162 Warp->computeCoeff(tp);
163 for (
unsigned int point = 0; point <
templateSize; point++) {
169 Warp->computeDenom(
X1, tp);
175 if ((i2 >= 0) && (j2 >= 0) && (i2 < I.
getHeight() - 1) && (j2 < I.
getWidth() - 1)) {
184 int cr = (int)((IW * (
Nc - 1)) / 255.);
185 int ct = (int)((Tij * (
Nc - 1)) / 255.);
186 double er = (IW * (
Nc - 1)) / 255. - cr;
187 double et = ((double)Tij * (
Nc - 1)) / 255. - ct;
191 vpTemplateTrackerMIBSpline::PutPVBsplineD(
PrtD, cr, er, ct, et,
Nc, 1.,
bspline);
198 for (
int r = 0; r <
Nc; r++)
199 for (
int t = 0; t <
Nc; t++) {
204 Prt[r2 *
Ncb + t2] += *pt;
212 for (
unsigned int r = 0; r < Ncb_; r++)
213 for (
unsigned int t = 0; t < Ncb_; t++)
215 Prt[r * Ncb_ + t] =
Prt[r * Ncb_ + t] / Nbpoint;
217 memset(
Pr, 0, Ncb_ *
sizeof(
double));
218 for (
unsigned int r = 0; r < Ncb_; r++) {
219 for (
unsigned int t = 0; t < Ncb_; t++)
220 Pr[r] +=
Prt[r * Ncb_ + t];
224 memset(
Pt, 0, Ncb_ *
sizeof(
double));
225 for (
unsigned int t = 0; t < Ncb_; t++) {
226 for (
unsigned int r = 0; r < Ncb_; r++)
227 Pt[t] +=
Prt[r * Ncb_ + t];
229 for (
unsigned int r = 0; r < Ncb_; r++)
231 if (std::fabs(
Pr[r]) > std::numeric_limits<double>::epsilon())
232 MI -=
Pr[r] * log(
Pr[r]);
234 for (
unsigned int t = 0; t < Ncb_; t++)
236 if (std::fabs(
Pt[t]) > std::numeric_limits<double>::epsilon())
237 MI -=
Pt[t] * log(
Pt[t]);
239 for (
unsigned int r = 0; r < Ncb_; r++)
240 for (
unsigned int t = 0; t < Ncb_; t++)
242 if (std::fabs(
Prt[r * Ncb_ + t]) > std::numeric_limits<double>::epsilon())
243 MI +=
Prt[r * Ncb_ + t] * log(
Prt[r * Ncb_ + t]);
260 double Prt_[256][256];
262 memset(Pr_, 0, 256 *
sizeof(
double));
263 memset(Pt_, 0, 256 *
sizeof(
double));
264 memset(Prt_, 0, 256 * 256 *
sizeof(
double));
266 for (
unsigned int point = 0; point <
templateSize; point++) {
272 Warp->computeDenom(
X1, tp);
277 if ((i2 >= 0) && (j2 >= 0) && (i2 < I.
getHeight() - 1) && (j2 < I.
getWidth() - 1)) {
281 IW = I[(int)i2][(
int)j2];
287 Prt_[(int)Tij][(
int)IW]++;
291 for (
int i = 0; i < 256; i++) {
294 for (
int j = 0; j < 256; j++)
295 Prt_[i][j] /= Nbpoint;
298 for (
int r = 0; r < 256; r++)
300 if (std::fabs(Pr_[r]) > std::numeric_limits<double>::epsilon())
301 MI -= Pr_[r] * log(Pr_[r]);
303 for (
int t = 0; t < 256; t++)
305 if (std::fabs(Pt_[t]) > std::numeric_limits<double>::epsilon())
306 MI -= Pt_[t] * log(Pt_[t]);
309 for (
int r = 0; r < 256; r++)
310 for (
int t = 0; t < 256; t++)
312 if (std::fabs(Prt_[r][t]) > std::numeric_limits<double>::epsilon())
313 denom -= (Prt_[r][t] * log(Prt_[r][t]));
316 if (std::fabs(denom) > std::numeric_limits<double>::epsilon())
351 unsigned int Nc_ = (
unsigned int)
Nc;
352 unsigned int Ncb_ = (
unsigned int)
Ncb;
353 unsigned int bspline_ = (
unsigned int)
bspline;
355 for (
unsigned int r = 0; r < Nc_; r++) {
356 for (
unsigned int t = 0; t < Nc_; t++) {
357 for (
unsigned int r2 = 0; r2 < bspline_; r2++) {
358 for (
unsigned int t2 = 0; t2 < bspline_; t2++) {
359 Prt[((r2 + r) * Ncb_ + (t2 + t))] += *pt++;
360 for (
unsigned int ip = 0; ip <
nbParam; ip++) {
361 dPrt[((r2 + r) * Ncb_ + (t2 + t)) * nbParam + ip] += *pt++;
362 for (
unsigned int it = 0; it <
nbParam; it++) {
363 d2Prt[((r2 + r) * Ncb_ + (t2 + t)) * nbParam * nbParam + ip * nbParam + it] += *pt++;
420 unsigned int indd, indd2;
422 for (
volatile int i = 0; i <
Ncb *
Ncb; i++) {
423 Prt[i] =
Prt[i] / nbpoint;
424 for (
unsigned int j = 0; j <
nbParam; j++) {
427 for (
unsigned int k = 0; k <
nbParam; k++) {
437 unsigned int Ncb_ = (
unsigned int)
Ncb;
440 memset(
Pr, 0, Ncb_ *
sizeof(
double));
441 for (
unsigned int r = 0; r < Ncb_; r++) {
442 for (
unsigned int t = 0; t < Ncb_; t++)
443 Pr[r] +=
Prt[r * Ncb_ + t];
447 memset(
Pt, 0, Ncb_ *
sizeof(
double));
448 for (
unsigned int t = 0; t < Ncb_; t++) {
449 for (
unsigned int r = 0; r < Ncb_; r++)
450 Pt[t] +=
Prt[r * Ncb_ + t];
455 for (
unsigned int r = 0; r < Ncb_; r++) {
457 if (std::fabs(
Pr[r]) > std::numeric_limits<double>::epsilon()) {
459 MI -=
Pr[r] * log(
Pr[r]);
464 for (
unsigned int t = 0; t < Ncb_; t++) {
466 if (std::fabs(
Pt[t]) > std::numeric_limits<double>::epsilon()) {
468 MI -=
Pt[t] * log(
Pt[t]);
472 for (
unsigned int r = 0; r < Ncb_; r++)
473 for (
unsigned int t = 0; t < Ncb_; t++)
475 if (std::fabs(
Prt[r * Ncb_ + t]) > std::numeric_limits<double>::epsilon())
476 MI +=
Prt[r * Ncb_ + t] * log(
Prt[r * Ncb_ + t]);
481 double seuilevitinf = 1e-200;
484 unsigned int Ncb_ = (
unsigned int)
Ncb;
485 for (
unsigned int t = 0; t < Ncb_; t++) {
487 if (
Pt[t] > seuilevitinf) {
488 for (
unsigned int r = 0; r < Ncb_; r++) {
490 if (
Prt[r * Ncb_ + t] > seuilevitinf) {
491 for (
unsigned int it = 0; it <
nbParam; it++)
494 dtemp = 1. + log(
Prt[r * Ncb_ + t] /
Pt[t]);
495 for (
unsigned int it = 0; it <
nbParam; it++)
496 for (
unsigned int jt = 0; jt <
nbParam; jt++)
499 d2Prt[(r * Ncb_ + t) * nbParam * nbParam + it * nbParam + jt] * dtemp;
501 Hessian[it][jt] +=
d2Prt[(r * Ncb_ + t) * nbParam * nbParam + it * nbParam + jt] * dtemp;
517 double seuilevitinf = 1e-200;
519 double u = 0, v = 0, B = 0;
520 double *du =
new double[
nbParam];
521 double *dv =
new double[
nbParam];
522 double *A =
new double[
nbParam];
523 double *dB =
new double[
nbParam];
524 double **d2u =
new double *[
nbParam];
525 double **d2v =
new double *[
nbParam];
526 double **dA =
new double *[
nbParam];
527 for (
unsigned int i = 0; i <
nbParam; i++) {
533 memset(du, 0, nbParam *
sizeof(
double));
534 memset(dv, 0, nbParam *
sizeof(
double));
535 memset(A, 0, nbParam *
sizeof(
double));
536 memset(dB, 0, nbParam *
sizeof(
double));
537 memset(
dprtemp, 0, nbParam *
sizeof(
double));
539 for (
unsigned int it = 0; it <
nbParam; it++) {
540 for (
unsigned int jt = 0; jt <
nbParam; jt++) {
541 memset(d2u[it], 0, nbParam *
sizeof(
double));
542 memset(d2v[it], 0, nbParam *
sizeof(
double));
543 memset(dA[it], 0, nbParam *
sizeof(
double));
547 unsigned int Ncb_ = (
unsigned int)
Ncb;
549 for (
unsigned int t = 0; t < Ncb_; t++) {
551 if (
Pt[t] > seuilevitinf) {
552 for (
unsigned int r = 0; r < Ncb_; r++) {
554 if (
Prt[r * Ncb_ + t] > seuilevitinf) {
555 for (
unsigned int it = 0; it <
nbParam; it++) {
561 u +=
Prt[r * Ncb_ + t] * log(
Pt[t] *
Pr[r]);
563 v +=
Prt[r * Ncb_ + t] * log(
Prt[r * Ncb_ + t]);
565 for (
unsigned int it = 0; it <
nbParam; it++) {
569 dv[it] +=
dprtemp[it] * (1 + log(
Prt[r * Ncb_ + t]));
571 for (
unsigned int it = 0; it <
nbParam; it++) {
572 for (
unsigned int jt = 0; jt <
nbParam; jt++) {
573 d2u[it][jt] +=
d2Prt[(r * Ncb_ + t) * nbParam * nbParam + it * nbParam + jt] * log(
Pt[t] *
Pr[r]) +
576 d2Prt[(r * Ncb_ + t) * nbParam * nbParam + it * nbParam + jt] * (1 + log(
Prt[r * Ncb_ + t])) +
586 for (
unsigned int it = 0; it <
nbParam; it++) {
588 A[it] = du[it] * v - u * dv[it];
590 dB[it] = 2 * v * dv[it];
593 for (
unsigned int jt = 0; jt <
nbParam; jt++) {
595 dA[it][jt] = d2u[it][jt] * v - d2v[it][jt] * u;
597 Hessian[it][jt] = (dA[it][jt] * B - A[it] * dB[it]) / (B * B);
605 for (
unsigned int i = 0; i <
nbParam; i++) {
622 double seuilevitinf = 1e-200;
624 unsigned int Ncb_ = (
unsigned int)
Ncb;
626 for (
unsigned int t = 0; t < Ncb_; t++) {
628 if (
Pt[t] > seuilevitinf) {
629 for (
unsigned int r = 0; r < Ncb_; r++) {
630 if (
Prt[r * Ncb_ + t] > seuilevitinf)
633 for (
unsigned int it = 0; it <
nbParam; it++)
636 dtemp = 1. + log(
Prt[r * Ncb_ + t] /
Pt[t]);
638 for (
unsigned int it = 0; it <
nbParam; it++)
647 unsigned int Ncb_ = (
unsigned int)
Ncb;
648 unsigned int Nc_ = (
unsigned int)
Nc;
649 unsigned int influBspline_ = (
unsigned int)
influBspline;
651 memset(
Prt, 0, Ncb_ * Ncb_ *
sizeof(
double));
652 memset(
dPrt, 0, Ncb_ * Ncb_ *
nbParam *
sizeof(
double));
665 unsigned int tNcb = (
unsigned int)(nc + bspline_);
666 unsigned int tinfluBspline = (
unsigned int)(bspline_ * bspline_);
667 unsigned int nc_ = (
unsigned int)nc;
668 double *tPrtD =
new double[nc_ * nc_ * tinfluBspline];
669 double *tPrt =
new double[tNcb * tNcb];
670 double *tPr =
new double[tNcb];
671 double *tPt =
new double[tNcb];
674 volatile int Nbpoint = 0;
680 memset(tPrt, 0, tNcb * tNcb *
sizeof(
double));
681 memset(tPrtD, 0, nc_ * nc_ * tinfluBspline *
sizeof(
double));
684 Warp->computeCoeff(tp);
685 for (
unsigned int point = 0; point <
templateSize; point++) {
691 Warp->computeDenom(
X1, tp);
697 if ((i2 >= 0) && (j2 >= 0) && (i2 < I.
getHeight() - 1) && (j2 < I.
getWidth()) - 1) {
706 int cr = (int)((IW * (nc - 1)) / 255.);
707 int ct = (int)((Tij * (nc - 1)) / 255.);
708 double er = (IW * (nc - 1)) / 255. - cr;
709 double et = ((double)Tij * (nc - 1)) / 255. - ct;
713 vpTemplateTrackerMIBSpline::PutPVBsplineD(tPrtD, cr, er, ct, et, nc, 1., bspline_);
717 int tNcb_ = (int)tNcb;
718 int tinfluBspline_ = (int)tinfluBspline;
719 for (
int r = 0; r < nc; r++)
720 for (
int t = 0; t < nc; t++) {
721 for (
volatile int i = 0; i < tinfluBspline_; i++) {
723 r2 = r + i / bspline_;
724 t2 = t + i % bspline_;
725 tPrt[r2 * tNcb_ + t2] += *pt;
739 for (
unsigned int r = 0; r < tNcb; r++)
740 for (
unsigned int t = 0; t < tNcb; t++)
742 tPrt[r * tNcb + t] = tPrt[r * tNcb + t] / Nbpoint;
744 memset(tPr, 0, tNcb *
sizeof(
double));
745 for (
unsigned int r = 0; r < tNcb; r++) {
746 for (
unsigned int t = 0; t < tNcb; t++)
747 tPr[r] += tPrt[r * tNcb + t];
751 memset(tPt, 0, (
size_t)(tNcb *
sizeof(
double)));
752 for (
unsigned int t = 0; t < tNcb; t++) {
753 for (
unsigned int r = 0; r < tNcb; r++)
754 tPt[t] += tPrt[r * tNcb + t];
756 for (
unsigned int r = 0; r < tNcb; r++)
758 if (std::fabs(tPr[r]) > std::numeric_limits<double>::epsilon())
759 MI -= tPr[r] * log(tPr[r]);
761 for (
unsigned int t = 0; t < tNcb; t++)
763 if (std::fabs(tPt[t]) > std::numeric_limits<double>::epsilon())
764 MI -= tPt[t] * log(tPt[t]);
766 for (
unsigned int r = 0; r < tNcb; r++)
767 for (
unsigned int t = 0; t < tNcb; t++)
769 if (std::fabs(tPrt[r * tNcb + t]) > std::numeric_limits<double>::epsilon())
770 MI += tPrt[r * tNcb + t] * log(tPrt[r * tNcb + t]);
789 volatile int Nbpoint = 0;
790 unsigned int Tij, IW;
797 Warp->computeCoeff(tp);
798 for (
unsigned int point = 0; point <
templateSize; point++) {
804 Warp->computeDenom(
X1, tp);
810 if ((i2 >= 0) && (j2 >= 0) && (i2 < I.
getHeight() - 1) && (j2 < I.
getWidth()) - 1) {
815 IW = (
unsigned int)I.
getValue(i2, j2);
817 IW = (
unsigned int)GaussI.
getValue(i2, j2);
828 Prt256 = Prt256 / Nbpoint;
829 Pr256 = Pr256 / Nbpoint;
830 Pt256 = Pt256 / Nbpoint;
834 for (
unsigned int t = 0; t < 256; t++) {
835 for (
unsigned int r = 0; r < 256; r++) {
837 if (std::fabs(Prt256[r][t]) > std::numeric_limits<double>::epsilon())
838 MI += Prt256[r][t] * log(Prt256[r][t]);
841 if (std::fabs(Pt256[t]) > std::numeric_limits<double>::epsilon())
842 MI += -Pt256[t] * log(Pt256[t]);
844 if (std::fabs(Pr256[t]) > std::numeric_limits<double>::epsilon())
845 MI += -Pr256[t] * log(Pr256[t]);
vpTemplateTrackerMI()
Default constructor.
Implementation of a matrix and operations on matrices.
double NMI_postEstimation
void computeHessien(vpMatrix &H)
unsigned int getWidth() const
vpMatrix covarianceMatrix
vpTemplateTrackerPoint * ptTemplate
virtual void warpX(const int &i, const int &j, double &i2, double &j2, const vpColVector &ParamM)=0
void resize(const unsigned int nrows, const unsigned int ncols, const bool flagNullify=true, const bool recopy_=true)
error that can be emited by ViSP classes.
vpHessienApproximationType ApproxHessian
double getMI256(const vpImage< unsigned char > &I, const vpColVector &tp)
unsigned int templateSize
Type getValue(unsigned int i, unsigned int j) const
Error that can be emited by the vpTracker class and its derivates.
double getCost(const vpImage< unsigned char > &I, const vpColVector &tp)
void computeMI(double &MI)
double getNormalizedCost(const vpImage< unsigned char > &I, const vpColVector &tp)
Implementation of column vector and the associated operations.
unsigned int getHeight() const
static void filter(const vpImage< double > &I, vpImage< double > &Iu, vpImage< double > &Iv, const vpMatrix &M, const bool convolve=false)
void computeHessienNormalized(vpMatrix &H)
void setBspline(const vpBsplineType &newbs)
virtual ~vpTemplateTrackerMI()
vpTemplateTrackerWarp * Warp
void computeProba(int &nbpoint)
vpHessienType hessianComputation
void resize(const unsigned int i, const bool flagNullify=true)