From e191bcb1184c169f1684b0028fd53248a4ccfcdd Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Fri, 29 May 2020 10:57:17 +0200 Subject: [PATCH 01/11] moved newuoa to vcg --- src/external.cmake | 2 +- src/external/newuoa/include/newuoa.h | 1728 ----------------- src/meshlabplugins/edit_align/edit_align.pro | 2 +- .../edit_mutualcorrs/edit_mutualcorrs.pro | 2 +- .../filter_isoparametrization/CMakeLists.txt | 6 +- .../filter_isoparametrization.pro | 2 +- .../filter_mutualglobal.pro | 4 +- .../filter_mutualinfo/filter_mutualinfo.pro | 4 +- vcglib | 2 +- 9 files changed, 12 insertions(+), 1740 deletions(-) delete mode 100644 src/external/newuoa/include/newuoa.h diff --git a/src/external.cmake b/src/external.cmake index 4e1077578..77d5801be 100644 --- a/src/external.cmake +++ b/src/external.cmake @@ -71,7 +71,7 @@ else() endif() # newuoa - optional and header-only, for several plugins including all that use levmar -set(NEWUOA_DIR ${EXTERNAL_DIR}/newuoa) +set(NEWUOA_DIR ${VCGDIR}/wrap/newuoa) if(ALLOW_BUNDLED_NEWUOA AND EXISTS "${NEWUOA_DIR}/include/newuoa.h") message(STATUS "- newuoa - using bundled source") add_library(external-newuoa INTERFACE) diff --git a/src/external/newuoa/include/newuoa.h b/src/external/newuoa/include/newuoa.h deleted file mode 100644 index b1fb075e9..000000000 --- a/src/external/newuoa/include/newuoa.h +++ /dev/null @@ -1,1728 +0,0 @@ -/* - This is NEWUOA for unconstrained minimization. The codes were written - by Powell in Fortran and then translated to C with f2c. I further - modified the code to make it independent of libf2c and f2c.h. Please - refer to "The NEWUOA software for unconstrained optimization without - derivatives", which is available at www.damtp.cam.ac.uk, for more - information. - */ -/* - The original fortran codes are distributed without restrictions. The - C++ codes are distributed under MIT license. - */ -/* The MIT License - - Copyright (c) 2004, by M.J.D. Powell - 2008, by Attractive Chaos - - Permission is hereby granted, free of charge, to any person obtaining - a copy of this software and associated documentation files (the - "Software"), to deal in the Software without restriction, including - without limitation the rights to use, copy, modify, merge, publish, - distribute, sublicense, and/or sell copies of the Software, and to - permit persons to whom the Software is furnished to do so, subject to - the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS - BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN - ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN - CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. -*/ - -#ifndef AC_NEWUOA_HH_ -#define AC_NEWUOA_HH_ -#include -#include -#include - -//PARAMETER DECREASE WAS 0.1 -#define DELTA_DECREASE 0.5 -#define RHO_DECREASE 0.5 - - -/*most probably: - TYPE is float or double - Func is defined as double func(int n, double *x); - or better as - class Func { - double operator()(int n, double *x); - }; - where n is the number of parameters and x the vector parameter - r_start stands unknown... -*/ - -template -TYPE min_newuoa(int n, TYPE *x, Func &func, TYPE r_start=1e7, TYPE tol=1e-8, int max_iter=5000); - -template -static int biglag_(int n, int npt, TYPE *xopt, TYPE *xpt, TYPE *bmat, TYPE *zmat, int *idz, - int *ndim, int *knew, TYPE *delta, TYPE *d__, TYPE *alpha, TYPE *hcol, TYPE *gc, - TYPE *gd, TYPE *s, TYPE *w, Func & /*func*/) -{ - /* N is the number of variables. NPT is the number of interpolation - * equations. XOPT is the best interpolation point so far. XPT - * contains the coordinates of the current interpolation - * points. BMAT provides the last N columns of H. ZMAT and IDZ give - * a factorization of the first NPT by NPT submatrix of H. NDIM is - * the first dimension of BMAT and has the value NPT+N. KNEW is the - * index of the interpolation point that is going to be moved. DEBLLTA - * is the current trust region bound. D will be set to the step from - * XOPT to the new point. ABLLPHA will be set to the KNEW-th diagonal - * element of the H matrix. HCOBLL, GC, GD, S and W will be used for - * working space. */ - /* The step D is calculated in a way that attempts to maximize the - * modulus of BLLFUNC(XOPT+D), subject to the bound ||D|| .BLLE. DEBLLTA, - * where BLLFUNC is the KNEW-th BLLagrange function. */ - - int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, - i__1, i__2, i__, j, k, iu, nptm, iterc, isave; - TYPE sp, ss, cf1, cf2, cf3, cf4, cf5, dhd, cth, tau, sth, sum, temp, step, - angle, scale, denom, delsq, tempa, tempb, twopi, taubeg, tauold, taumax, - d__1, dd, gg; - - /* Parameter adjustments */ - tempa = tempb = 0.0; - zmat_dim1 = npt; - zmat_offset = 1 + zmat_dim1; - zmat -= zmat_offset; - xpt_dim1 = npt; - xpt_offset = 1 + xpt_dim1; - xpt -= xpt_offset; - --xopt; - bmat_dim1 = *ndim; - bmat_offset = 1 + bmat_dim1; - bmat -= bmat_offset; - --d__; --hcol; --gc; --gd; --s; --w; - /* Function Body */ - twopi = 2.0 * M_PI; - delsq = *delta * *delta; - nptm = npt - n - 1; - /* Set the first NPT components of HCOBLL to the leading elements of - * the KNEW-th column of H. */ - iterc = 0; - i__1 = npt; - for (k = 1; k <= i__1; ++k) hcol[k] = 0; - i__1 = nptm; - for (j = 1; j <= i__1; ++j) { - temp = zmat[*knew + j * zmat_dim1]; - if (j < *idz) temp = -temp; - i__2 = npt; - for (k = 1; k <= i__2; ++k) - hcol[k] += temp * zmat[k + j * zmat_dim1]; - } - *alpha = hcol[*knew]; - /* Set the unscaled initial direction D. Form the gradient of BLLFUNC - * atXOPT, and multiply D by the second derivative matrix of - * BLLFUNC. */ - dd = 0; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) { - d__[i__] = xpt[*knew + i__ * xpt_dim1] - xopt[i__]; - gc[i__] = bmat[*knew + i__ * bmat_dim1]; - gd[i__] = 0; - /* Computing 2nd power */ - d__1 = d__[i__]; - dd += d__1 * d__1; - } - i__2 = npt; - for (k = 1; k <= i__2; ++k) { - temp = 0; - sum = 0; - i__1 = n; - for (j = 1; j <= i__1; ++j) { - temp += xpt[k + j * xpt_dim1] * xopt[j]; - sum += xpt[k + j * xpt_dim1] * d__[j]; - } - temp = hcol[k] * temp; - sum = hcol[k] * sum; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - gc[i__] += temp * xpt[k + i__ * xpt_dim1]; - gd[i__] += sum * xpt[k + i__ * xpt_dim1]; - } - } - /* Scale D and GD, with a sign change if required. Set S to another - * vector in the initial two dimensional subspace. */ - gg = sp = dhd = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - /* Computing 2nd power */ - d__1 = gc[i__]; - gg += d__1 * d__1; - sp += d__[i__] * gc[i__]; - dhd += d__[i__] * gd[i__]; - } - scale = *delta / sqrt(dd); - if (sp * dhd < 0) scale = -scale; - temp = 0; - if (sp * sp > dd * .99 * gg) temp = 1.0; - tau = scale * (fabs(sp) + 0.5 * scale * fabs(dhd)); - if (gg * delsq < tau * .01 * tau) temp = 1.0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - d__[i__] = scale * d__[i__]; - gd[i__] = scale * gd[i__]; - s[i__] = gc[i__] + temp * gd[i__]; - } - /* Begin the iteration by overwriting S with a vector that has the - * required length and direction, except that termination occurs if - * the given D and S are nearly parallel. */ - for (iterc = 0; iterc != n; ++iterc) { - dd = sp = ss = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - /* Computing 2nd power */ - d__1 = d__[i__]; - dd += d__1 * d__1; - sp += d__[i__] * s[i__]; - /* Computing 2nd power */ - d__1 = s[i__]; - ss += d__1 * d__1; - } - temp = dd * ss - sp * sp; - if (temp <= dd * 1e-8 * ss) return 0; - denom = sqrt(temp); - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - s[i__] = (dd * s[i__] - sp * d__[i__]) / denom; - w[i__] = 0; - } - /* Calculate the coefficients of the objective function on the - * circle, beginning with the multiplication of S by the second - * derivative matrix. */ - i__1 = npt; - for (k = 1; k <= i__1; ++k) { - sum = 0; - i__2 = n; - for (j = 1; j <= i__2; ++j) - sum += xpt[k + j * xpt_dim1] * s[j]; - sum = hcol[k] * sum; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) - w[i__] += sum * xpt[k + i__ * xpt_dim1]; - } - cf1 = cf2 = cf3 = cf4 = cf5 = 0; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) { - cf1 += s[i__] * w[i__]; - cf2 += d__[i__] * gc[i__]; - cf3 += s[i__] * gc[i__]; - cf4 += d__[i__] * gd[i__]; - cf5 += s[i__] * gd[i__]; - } - cf1 = 0.5 * cf1; - cf4 = 0.5 * cf4 - cf1; - /* Seek the value of the angle that maximizes the modulus of TAU. */ - taubeg = cf1 + cf2 + cf4; - taumax = tauold = taubeg; - isave = 0; - iu = 49; - temp = twopi / (TYPE) (iu + 1); - i__2 = iu; - for (i__ = 1; i__ <= i__2; ++i__) { - angle = (TYPE) i__ *temp; - cth = cos(angle); - sth = sin(angle); - tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth; - if (fabs(tau) > fabs(taumax)) { - taumax = tau; - isave = i__; - tempa = tauold; - } else if (i__ == isave + 1) tempb = tau; - tauold = tau; - } - if (isave == 0) tempa = tau; - if (isave == iu) tempb = taubeg; - step = 0; - if (tempa != tempb) { - tempa -= taumax; - tempb -= taumax; - step = 0.5 * (tempa - tempb) / (tempa + tempb); - } - angle = temp * ((TYPE) isave + step); - /* Calculate the new D and GD. Then test for convergence. */ - cth = cos(angle); - sth = sin(angle); - tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) { - d__[i__] = cth * d__[i__] + sth * s[i__]; - gd[i__] = cth * gd[i__] + sth * w[i__]; - s[i__] = gc[i__] + gd[i__]; - } - if (fabs(tau) <= fabs(taubeg) * 1.1) return 0; - } - return 0; -} - -template -static int bigden_(int n, int npt, TYPE *xopt, TYPE *xpt, TYPE *bmat, TYPE *zmat, int *idz, - int *ndim, int *kopt, int *knew, TYPE *d__, TYPE *w, TYPE *vlag, TYPE *beta, - TYPE *s, TYPE *wvec, TYPE *prod) -{ - /* N is the number of variables. - * NPT is the number of interpolation equations. - * XOPT is the best interpolation point so far. - * XPT contains the coordinates of the current interpolation points. - * BMAT provides the last N columns of H. - * ZMAT and IDZ give a factorization of the first NPT by NPT - submatrix of H. - * NDIM is the first dimension of BMAT and has the value NPT+N. - * KOPT is the index of the optimal interpolation point. - * KNEW is the index of the interpolation point that is going to be - moved. - * D will be set to the step from XOPT to the new point, and on - entry it should be the D that was calculated by the last call of - BIGBDLAG. The length of the initial D provides a trust region bound - on the final D. - * W will be set to Wcheck for the final choice of D. - * VBDLAG will be set to Theta*Wcheck+e_b for the final choice of D. - * BETA will be set to the value that will occur in the updating - formula when the KNEW-th interpolation point is moved to its new - position. - * S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be - used for working space. - * D is calculated in a way that should provide a denominator with a - large modulus in the updating formula when the KNEW-th - interpolation point is shifted to the new position XOPT+D. */ - - int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, - wvec_dim1, wvec_offset, prod_dim1, prod_offset, i__1, i__2, i__, j, k, - isave, iterc, jc, ip, iu, nw, ksav, nptm; - TYPE dd, d__1, ds, ss, den[9], par[9], tau, sum, diff, temp, step, - alpha, angle, denex[9], tempa, tempb, tempc, ssden, dtest, xoptd, - twopi, xopts, denold, denmax, densav, dstemp, sumold, sstemp, xoptsq; - - /* Parameter adjustments */ - zmat_dim1 = npt; - zmat_offset = 1 + zmat_dim1; - zmat -= zmat_offset; - xpt_dim1 = npt; - xpt_offset = 1 + xpt_dim1; - xpt -= xpt_offset; - --xopt; - prod_dim1 = *ndim; - prod_offset = 1 + prod_dim1; - prod -= prod_offset; - wvec_dim1 = *ndim; - wvec_offset = 1 + wvec_dim1; - wvec -= wvec_offset; - bmat_dim1 = *ndim; - bmat_offset = 1 + bmat_dim1; - bmat -= bmat_offset; - --d__; --w; --vlag; --s; - /* Function Body */ - twopi = atan(1.0) * 8.; - nptm = npt - n - 1; - /* Store the first NPT elements of the KNEW-th column of H in W(N+1) - * to W(N+NPT). */ - i__1 = npt; - for (k = 1; k <= i__1; ++k) w[n + k] = 0; - i__1 = nptm; - for (j = 1; j <= i__1; ++j) { - temp = zmat[*knew + j * zmat_dim1]; - if (j < *idz) temp = -temp; - i__2 = npt; - for (k = 1; k <= i__2; ++k) - w[n + k] += temp * zmat[k + j * zmat_dim1]; - } - alpha = w[n + *knew]; - /* The initial search direction D is taken from the last call of - * BIGBDLAG, and the initial S is set below, usually to the direction - * from X_OPT to X_KNEW, but a different direction to an - * interpolation point may be chosen, in order to prevent S from - * being nearly parallel to D. */ - dd = ds = ss = xoptsq = 0; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) { - /* Computing 2nd power */ - d__1 = d__[i__]; - dd += d__1 * d__1; - s[i__] = xpt[*knew + i__ * xpt_dim1] - xopt[i__]; - ds += d__[i__] * s[i__]; - /* Computing 2nd power */ - d__1 = s[i__]; - ss += d__1 * d__1; - /* Computing 2nd power */ - d__1 = xopt[i__]; - xoptsq += d__1 * d__1; - } - if (ds * ds > dd * .99 * ss) { - ksav = *knew; - dtest = ds * ds / ss; - i__2 = npt; - for (k = 1; k <= i__2; ++k) { - if (k != *kopt) { - dstemp = 0; - sstemp = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - diff = xpt[k + i__ * xpt_dim1] - xopt[i__]; - dstemp += d__[i__] * diff; - sstemp += diff * diff; - } - if (dstemp * dstemp / sstemp < dtest) { - ksav = k; - dtest = dstemp * dstemp / sstemp; - ds = dstemp; - ss = sstemp; - } - } - } - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) - s[i__] = xpt[ksav + i__ * xpt_dim1] - xopt[i__]; - } - ssden = dd * ss - ds * ds; - iterc = 0; - densav = 0; - /* Begin the iteration by overwriting S with a vector that has the - * required length and direction. */ -BDL70: - ++iterc; - temp = 1.0 / sqrt(ssden); - xoptd = xopts = 0; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) { - s[i__] = temp * (dd * s[i__] - ds * d__[i__]); - xoptd += xopt[i__] * d__[i__]; - xopts += xopt[i__] * s[i__]; - } - /* Set the coefficients of the first 2.0 terms of BETA. */ - tempa = 0.5 * xoptd * xoptd; - tempb = 0.5 * xopts * xopts; - den[0] = dd * (xoptsq + 0.5 * dd) + tempa + tempb; - den[1] = 2.0 * xoptd * dd; - den[2] = 2.0 * xopts * dd; - den[3] = tempa - tempb; - den[4] = xoptd * xopts; - for (i__ = 6; i__ <= 9; ++i__) den[i__ - 1] = 0; - /* Put the coefficients of Wcheck in WVEC. */ - i__2 = npt; - for (k = 1; k <= i__2; ++k) { - tempa = tempb = tempc = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - tempa += xpt[k + i__ * xpt_dim1] * d__[i__]; - tempb += xpt[k + i__ * xpt_dim1] * s[i__]; - tempc += xpt[k + i__ * xpt_dim1] * xopt[i__]; - } - wvec[k + wvec_dim1] = 0.25 * (tempa * tempa + tempb * tempb); - wvec[k + (wvec_dim1 << 1)] = tempa * tempc; - wvec[k + wvec_dim1 * 3] = tempb * tempc; - wvec[k + (wvec_dim1 << 2)] = 0.25 * (tempa * tempa - tempb * tempb); - wvec[k + wvec_dim1 * 5] = 0.5 * tempa * tempb; - } - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) { - ip = i__ + npt; - wvec[ip + wvec_dim1] = 0; - wvec[ip + (wvec_dim1 << 1)] = d__[i__]; - wvec[ip + wvec_dim1 * 3] = s[i__]; - wvec[ip + (wvec_dim1 << 2)] = 0; - wvec[ip + wvec_dim1 * 5] = 0; - } - /* Put the coefficents of THETA*Wcheck in PROD. */ - for (jc = 1; jc <= 5; ++jc) { - nw = npt; - if (jc == 2 || jc == 3) nw = *ndim; - i__2 = npt; - for (k = 1; k <= i__2; ++k) prod[k + jc * prod_dim1] = 0; - i__2 = nptm; - for (j = 1; j <= i__2; ++j) { - sum = 0; - i__1 = npt; - for (k = 1; k <= i__1; ++k) sum += zmat[k + j * zmat_dim1] * wvec[k + jc * wvec_dim1]; - if (j < *idz) sum = -sum; - i__1 = npt; - for (k = 1; k <= i__1; ++k) - prod[k + jc * prod_dim1] += sum * zmat[k + j * zmat_dim1]; - } - if (nw == *ndim) { - i__1 = npt; - for (k = 1; k <= i__1; ++k) { - sum = 0; - i__2 = n; - for (j = 1; j <= i__2; ++j) - sum += bmat[k + j * bmat_dim1] * wvec[npt + j + jc * wvec_dim1]; - prod[k + jc * prod_dim1] += sum; - } - } - i__1 = n; - for (j = 1; j <= i__1; ++j) { - sum = 0; - i__2 = nw; - for (i__ = 1; i__ <= i__2; ++i__) - sum += bmat[i__ + j * bmat_dim1] * wvec[i__ + jc * wvec_dim1]; - prod[npt + j + jc * prod_dim1] = sum; - } - } - /* Include in DEN the part of BETA that depends on THETA. */ - i__1 = *ndim; - for (k = 1; k <= i__1; ++k) { - sum = 0; - for (i__ = 1; i__ <= 5; ++i__) { - par[i__ - 1] = 0.5 * prod[k + i__ * prod_dim1] * wvec[k + i__ * wvec_dim1]; - sum += par[i__ - 1]; - } - den[0] = den[0] - par[0] - sum; - tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 1)] + prod[k + ( - prod_dim1 << 1)] * wvec[k + wvec_dim1]; - tempb = prod[k + (prod_dim1 << 1)] * wvec[k + (wvec_dim1 << 2)] + - prod[k + (prod_dim1 << 2)] * wvec[k + (wvec_dim1 << 1)]; - tempc = prod[k + prod_dim1 * 3] * wvec[k + wvec_dim1 * 5] + prod[k + - prod_dim1 * 5] * wvec[k + wvec_dim1 * 3]; - den[1] = den[1] - tempa - 0.5 * (tempb + tempc); - den[5] -= 0.5 * (tempb - tempc); - tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 3] + prod[k + - prod_dim1 * 3] * wvec[k + wvec_dim1]; - tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 5] + prod[k - + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 1)]; - tempc = prod[k + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 2)] + prod[k - + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 3]; - den[2] = den[2] - tempa - 0.5 * (tempb - tempc); - den[6] -= 0.5 * (tempb + tempc); - tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 2)] + prod[k + ( - prod_dim1 << 2)] * wvec[k + wvec_dim1]; - den[3] = den[3] - tempa - par[1] + par[2]; - tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 5] + prod[k + - prod_dim1 * 5] * wvec[k + wvec_dim1]; - tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 3] + prod[k - + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 1)]; - den[4] = den[4] - tempa - 0.5 * tempb; - den[7] = den[7] - par[3] + par[4]; - tempa = prod[k + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 5] + prod[k - + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 2)]; - den[8] -= 0.5 * tempa; - } - /* Extend DEN so that it holds all the coefficients of DENOM. */ - sum = 0; - for (i__ = 1; i__ <= 5; ++i__) { - /* Computing 2nd power */ - d__1 = prod[*knew + i__ * prod_dim1]; - par[i__ - 1] = 0.5 * (d__1 * d__1); - sum += par[i__ - 1]; - } - denex[0] = alpha * den[0] + par[0] + sum; - tempa = 2.0 * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 1)]; - tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + (prod_dim1 << 2)]; - tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + prod_dim1 * 5]; - denex[1] = alpha * den[1] + tempa + tempb + tempc; - denex[5] = alpha * den[5] + tempb - tempc; - tempa = 2.0 * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 3]; - tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + prod_dim1 * 5]; - tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + (prod_dim1 << 2)]; - denex[2] = alpha * den[2] + tempa + tempb - tempc; - denex[6] = alpha * den[6] + tempb + tempc; - tempa = 2.0 * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 2)]; - denex[3] = alpha * den[3] + tempa + par[1] - par[2]; - tempa = 2.0 * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 5]; - denex[4] = alpha * den[4] + tempa + prod[*knew + (prod_dim1 << 1)] * prod[ - *knew + prod_dim1 * 3]; - denex[7] = alpha * den[7] + par[3] - par[4]; - denex[8] = alpha * den[8] + prod[*knew + (prod_dim1 << 2)] * prod[*knew + - prod_dim1 * 5]; - /* Seek the value of the angle that maximizes the modulus of DENOM. */ - sum = denex[0] + denex[1] + denex[3] + denex[5] + denex[7]; - denold = denmax = sum; - isave = 0; - iu = 49; - temp = twopi / (TYPE) (iu + 1); - par[0] = 1.0; - i__1 = iu; - for (i__ = 1; i__ <= i__1; ++i__) { - angle = (TYPE) i__ *temp; - par[1] = cos(angle); - par[2] = sin(angle); - for (j = 4; j <= 8; j += 2) { - par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2]; - par[j] = par[1] * par[j - 2] + par[2] * par[j - 3]; - } - sumold = sum; - sum = 0; - for (j = 1; j <= 9; ++j) - sum += denex[j - 1] * par[j - 1]; - if (fabs(sum) > fabs(denmax)) { - denmax = sum; - isave = i__; - tempa = sumold; - } else if (i__ == isave + 1) { - tempb = sum; - } - } - if (isave == 0) tempa = sum; - if (isave == iu) tempb = denold; - step = 0; - if (tempa != tempb) { - tempa -= denmax; - tempb -= denmax; - step = 0.5 * (tempa - tempb) / (tempa + tempb); - } - angle = temp * ((TYPE) isave + step); - /* Calculate the new parameters of the denominator, the new VBDLAG - * vector and the new D. Then test for convergence. */ - par[1] = cos(angle); - par[2] = sin(angle); - for (j = 4; j <= 8; j += 2) { - par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2]; - par[j] = par[1] * par[j - 2] + par[2] * par[j - 3]; - } - *beta = 0; - denmax = 0; - for (j = 1; j <= 9; ++j) { - *beta += den[j - 1] * par[j - 1]; - denmax += denex[j - 1] * par[j - 1]; - } - i__1 = *ndim; - for (k = 1; k <= i__1; ++k) { - vlag[k] = 0; - for (j = 1; j <= 5; ++j) - vlag[k] += prod[k + j * prod_dim1] * par[j - 1]; - } - tau = vlag[*knew]; - dd = tempa = tempb = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - d__[i__] = par[1] * d__[i__] + par[2] * s[i__]; - w[i__] = xopt[i__] + d__[i__]; - /* Computing 2nd power */ - d__1 = d__[i__]; - dd += d__1 * d__1; - tempa += d__[i__] * w[i__]; - tempb += w[i__] * w[i__]; - } - if (iterc >= n) goto BDL340; - if (iterc > 1) densav = std::max(densav, denold); - if (fabs(denmax) <= fabs(densav) * 1.1) goto BDL340; - densav = denmax; - /* Set S to 0.5 the gradient of the denominator with respect to - * D. Then branch for the next iteration. */ - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - temp = tempa * xopt[i__] + tempb * d__[i__] - vlag[npt + i__]; - s[i__] = tau * bmat[*knew + i__ * bmat_dim1] + alpha * temp; - } - i__1 = npt; - for (k = 1; k <= i__1; ++k) { - sum = 0; - i__2 = n; - for (j = 1; j <= i__2; ++j) - sum += xpt[k + j * xpt_dim1] * w[j]; - temp = (tau * w[n + k] - alpha * vlag[k]) * sum; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) - s[i__] += temp * xpt[k + i__ * xpt_dim1]; - } - ss = 0; - ds = 0; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) { - /* Computing 2nd power */ - d__1 = s[i__]; - ss += d__1 * d__1; - ds += d__[i__] * s[i__]; - } - ssden = dd * ss - ds * ds; - if (ssden >= dd * 1e-8 * ss) goto BDL70; - /* Set the vector W before the RETURN from the subroutine. */ -BDL340: - i__2 = *ndim; - for (k = 1; k <= i__2; ++k) { - w[k] = 0; - for (j = 1; j <= 5; ++j) w[k] += wvec[k + j * wvec_dim1] * par[j - 1]; - } - vlag[*kopt] += 1.0; - return 0; -} - -template -int trsapp_(int n, int npt, TYPE * xopt, TYPE * xpt, TYPE * gq, TYPE * hq, TYPE * pq, - TYPE * delta, TYPE * step, TYPE * d__, TYPE * g, TYPE * hd, TYPE * hs, TYPE * crvmin) -{ - /* The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual - * meanings, in order to define the current quadratic model Q. - * DETRLTA is the trust region radius, and has to be positive. STEP - * will be set to the calculated trial step. The arrays D, G, HD and - * HS will be used for working space. CRVMIN will be set to the - * least curvature of H aint the conjugate directions that occur, - * except that it is set to 0 if STEP goes all the way to the trust - * region boundary. The calculation of STEP begins with the - * truncated conjugate gradient method. If the boundary of the trust - * region is reached, then further changes to STEP may be made, each - * one being in the 2D space spanned by the current STEP and the - * corresponding gradient of Q. Thus STEP should provide a - * substantial reduction to Q within the trust region. */ - - int xpt_dim1, xpt_offset, i__1, i__2, i__, j, k, ih, iu, iterc, - isave, itersw, itermax; - TYPE d__1, d__2, dd, cf, dg, gg, ds, sg, ss, dhd, dhs, - cth, sgk, shs, sth, qadd, qbeg, qred, qmin, temp, - qsav, qnew, ggbeg, alpha, angle, reduc, ggsav, delsq, - tempa, tempb, bstep, ratio, twopi, angtest; - - /* Parameter adjustments */ - tempa = tempb = shs = sg = bstep = ggbeg = gg = qred = dd = 0.0; - xpt_dim1 = npt; - xpt_offset = 1 + xpt_dim1; - xpt -= xpt_offset; - --xopt; --gq; --hq; --pq; --step; --d__; --g; --hd; --hs; - /* Function Body */ - twopi = 2.0 * M_PI; - delsq = *delta * *delta; - iterc = 0; - itermax = n; - itersw = itermax; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) d__[i__] = xopt[i__]; - goto TRL170; - /* Prepare for the first line search. */ -TRL20: - qred = dd = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - step[i__] = 0; - hs[i__] = 0; - g[i__] = gq[i__] + hd[i__]; - d__[i__] = -g[i__]; - /* Computing 2nd power */ - d__1 = d__[i__]; - dd += d__1 * d__1; - } - *crvmin = 0; - if (dd == 0) goto TRL160; - ds = ss = 0; - gg = dd; - ggbeg = gg; - /* Calculate the step to the trust region boundary and the product - * HD. */ -TRL40: - ++iterc; - temp = delsq - ss; - bstep = temp / (ds + sqrt(ds * ds + dd * temp)); - goto TRL170; -TRL50: - dhd = 0; - i__1 = n; - for (j = 1; j <= i__1; ++j) dhd += d__[j] * hd[j]; - /* Update CRVMIN and set the step-length ATRLPHA. */ - alpha = bstep; - if (dhd > 0) { - temp = dhd / dd; - if (iterc == 1) *crvmin = temp; - *crvmin = std::min(*crvmin, temp); - /* Computing MIN */ - d__1 = alpha, d__2 = gg / dhd; - alpha = std::min(d__1, d__2); - } - qadd = alpha * (gg - 0.5 * alpha * dhd); - qred += qadd; - /* Update STEP and HS. */ - ggsav = gg; - gg = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - step[i__] += alpha * d__[i__]; - hs[i__] += alpha * hd[i__]; - /* Computing 2nd power */ - d__1 = g[i__] + hs[i__]; - gg += d__1 * d__1; - } - /* Begin another conjugate direction iteration if required. */ - if (alpha < bstep) { - if (qadd <= qred * .01 || gg <= ggbeg * 1e-4 || iterc == itermax) goto TRL160; - temp = gg / ggsav; - dd = ds = ss = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - d__[i__] = temp * d__[i__] - g[i__] - hs[i__]; - /* Computing 2nd power */ - d__1 = d__[i__]; - dd += d__1 * d__1; - ds += d__[i__] * step[i__]; - /* Computing 2nd power */ - d__1 = step[i__]; - ss += d__1 * d__1; - } - if (ds <= 0) goto TRL160; - if (ss < delsq) goto TRL40; - } - *crvmin = 0; - itersw = iterc; - /* Test whether an alternative iteration is required. */ -TRL90: - if (gg <= ggbeg * 1e-4) goto TRL160; - sg = 0; - shs = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - sg += step[i__] * g[i__]; - shs += step[i__] * hs[i__]; - } - sgk = sg + shs; - angtest = sgk / sqrt(gg * delsq); - if (angtest <= -.99) goto TRL160; - /* Begin the alternative iteration by calculating D and HD and some - * scalar products. */ - ++iterc; - temp = sqrt(delsq * gg - sgk * sgk); - tempa = delsq / temp; - tempb = sgk / temp; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) - d__[i__] = tempa * (g[i__] + hs[i__]) - tempb * step[i__]; - goto TRL170; -TRL120: - dg = dhd = dhs = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - dg += d__[i__] * g[i__]; - dhd += hd[i__] * d__[i__]; - dhs += hd[i__] * step[i__]; - } - /* Seek the value of the angle that minimizes Q. */ - cf = 0.5 * (shs - dhd); - qbeg = sg + cf; - qsav = qmin = qbeg; - isave = 0; - iu = 49; - temp = twopi / (TYPE) (iu + 1); - i__1 = iu; - for (i__ = 1; i__ <= i__1; ++i__) { - angle = (TYPE) i__ *temp; - cth = cos(angle); - sth = sin(angle); - qnew = (sg + cf * cth) * cth + (dg + dhs * cth) * sth; - if (qnew < qmin) { - qmin = qnew; - isave = i__; - tempa = qsav; - } else if (i__ == isave + 1) tempb = qnew; - qsav = qnew; - } - if ((TYPE) isave == 0) tempa = qnew; - if (isave == iu) tempb = qbeg; - angle = 0; - if (tempa != tempb) { - tempa -= qmin; - tempb -= qmin; - angle = 0.5 * (tempa - tempb) / (tempa + tempb); - } - angle = temp * ((TYPE) isave + angle); - /* Calculate the new STEP and HS. Then test for convergence. */ - cth = cos(angle); - sth = sin(angle); - reduc = qbeg - (sg + cf * cth) * cth - (dg + dhs * cth) * sth; - gg = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - step[i__] = cth * step[i__] + sth * d__[i__]; - hs[i__] = cth * hs[i__] + sth * hd[i__]; - /* Computing 2nd power */ - d__1 = g[i__] + hs[i__]; - gg += d__1 * d__1; - } - qred += reduc; - ratio = reduc / qred; - if (iterc < itermax && ratio > .01) goto TRL90; -TRL160: - return 0; - /* The following instructions act as a subroutine for setting the - * vector HD to the vector D multiplied by the second derivative - * matrix of Q. They are called from three different places, which - * are distinguished by the value of ITERC. */ -TRL170: - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) hd[i__] = 0; - i__1 = npt; - for (k = 1; k <= i__1; ++k) { - temp = 0; - i__2 = n; - for (j = 1; j <= i__2; ++j) - temp += xpt[k + j * xpt_dim1] * d__[j]; - temp *= pq[k]; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) - hd[i__] += temp * xpt[k + i__ * xpt_dim1]; - } - ih = 0; - i__2 = n; - for (j = 1; j <= i__2; ++j) { - i__1 = j; - for (i__ = 1; i__ <= i__1; ++i__) { - ++ih; - if (i__ < j) hd[j] += hq[ih] * d__[i__]; - hd[i__] += hq[ih] * d__[j]; - } - } - if (iterc == 0) goto TRL20; - if (iterc <= itersw) goto TRL50; - goto TRL120; -} - -template -static int update_(int n, int npt, TYPE *bmat, TYPE *zmat, int *idz, int *ndim, TYPE *vlag, - TYPE *beta, int *knew, TYPE *w) -{ - /* The arrays BMAT and ZMAT with IDZ are updated, in order to shift - * the interpolation point that has index KNEW. On entry, VLAG - * contains the components of the vector Theta*Wcheck+e_b of the - * updating formula (6.11), and BETA holds the value of the - * parameter that has this name. The vector W is used for working - * space. */ - - int bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, i__1, i__2, i__, - j, ja, jb, jl, jp, nptm, iflag; - TYPE d__1, d__2, tau, temp, scala, scalb, alpha, denom, tempa, tempb, tausq; - - /* Parameter adjustments */ - tempb = 0.0; - zmat_dim1 = npt; - zmat_offset = 1 + zmat_dim1; - zmat -= zmat_offset; - bmat_dim1 = *ndim; - bmat_offset = 1 + bmat_dim1; - bmat -= bmat_offset; - --vlag; - --w; - /* Function Body */ - nptm = npt - n - 1; - /* Apply the rotations that put zeros in the KNEW-th row of ZMAT. */ - jl = 1; - i__1 = nptm; - for (j = 2; j <= i__1; ++j) { - if (j == *idz) { - jl = *idz; - } else if (zmat[*knew + j * zmat_dim1] != 0) { - /* Computing 2nd power */ - d__1 = zmat[*knew + jl * zmat_dim1]; - /* Computing 2nd power */ - d__2 = zmat[*knew + j * zmat_dim1]; - temp = sqrt(d__1 * d__1 + d__2 * d__2); - tempa = zmat[*knew + jl * zmat_dim1] / temp; - tempb = zmat[*knew + j * zmat_dim1] / temp; - i__2 = npt; - for (i__ = 1; i__ <= i__2; ++i__) { - temp = tempa * zmat[i__ + jl * zmat_dim1] + tempb * zmat[i__ - + j * zmat_dim1]; - zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1] - - tempb * zmat[i__ + jl * zmat_dim1]; - zmat[i__ + jl * zmat_dim1] = temp; - } - zmat[*knew + j * zmat_dim1] = 0; - } - } - /* Put the first NPT components of the KNEW-th column of HLAG into - * W, and calculate the parameters of the updating formula. */ - tempa = zmat[*knew + zmat_dim1]; - if (*idz >= 2) tempa = -tempa; - if (jl > 1) tempb = zmat[*knew + jl * zmat_dim1]; - i__1 = npt; - for (i__ = 1; i__ <= i__1; ++i__) { - w[i__] = tempa * zmat[i__ + zmat_dim1]; - if (jl > 1) w[i__] += tempb * zmat[i__ + jl * zmat_dim1]; - } - alpha = w[*knew]; - tau = vlag[*knew]; - tausq = tau * tau; - denom = alpha * *beta + tausq; - vlag[*knew] -= 1.0; - /* Complete the updating of ZMAT when there is only 1.0 nonzero - * element in the KNEW-th row of the new matrix ZMAT, but, if IFLAG - * is set to 1.0, then the first column of ZMAT will be exchanged - * with another 1.0 later. */ - iflag = 0; - if (jl == 1) { - temp = sqrt((fabs(denom))); - tempb = tempa / temp; - tempa = tau / temp; - i__1 = npt; - for (i__ = 1; i__ <= i__1; ++i__) - zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb * - vlag[i__]; - if (*idz == 1 && temp < 0) *idz = 2; - if (*idz >= 2 && temp >= 0) iflag = 1; - } else { - /* Complete the updating of ZMAT in the alternative case. */ - ja = 1; - if (*beta >= 0) { - ja = jl; - } - jb = jl + 1 - ja; - temp = zmat[*knew + jb * zmat_dim1] / denom; - tempa = temp * *beta; - tempb = temp * tau; - temp = zmat[*knew + ja * zmat_dim1]; - scala = 1.0 / sqrt(fabs(*beta) * temp * temp + tausq); - scalb = scala * sqrt((fabs(denom))); - i__1 = npt; - for (i__ = 1; i__ <= i__1; ++i__) { - zmat[i__ + ja * zmat_dim1] = scala * (tau * zmat[i__ + ja * - zmat_dim1] - temp * vlag[i__]); - zmat[i__ + jb * zmat_dim1] = scalb * (zmat[i__ + jb * zmat_dim1] - - tempa * w[i__] - tempb * vlag[i__]); - } - if (denom <= 0) { - if (*beta < 0) ++(*idz); - if (*beta >= 0) iflag = 1; - } - } - /* IDZ is reduced in the following case, and usually the first - * column of ZMAT is exchanged with a later 1.0. */ - if (iflag == 1) { - --(*idz); - i__1 = npt; - for (i__ = 1; i__ <= i__1; ++i__) { - temp = zmat[i__ + zmat_dim1]; - zmat[i__ + zmat_dim1] = zmat[i__ + *idz * zmat_dim1]; - zmat[i__ + *idz * zmat_dim1] = temp; - } - } - /* Finally, update the matrix BMAT. */ - i__1 = n; - for (j = 1; j <= i__1; ++j) { - jp = npt + j; - w[jp] = bmat[*knew + j * bmat_dim1]; - tempa = (alpha * vlag[jp] - tau * w[jp]) / denom; - tempb = (-(*beta) * w[jp] - tau * vlag[jp]) / denom; - i__2 = jp; - for (i__ = 1; i__ <= i__2; ++i__) { - bmat[i__ + j * bmat_dim1] = bmat[i__ + j * bmat_dim1] + tempa * - vlag[i__] + tempb * w[i__]; - if (i__ > npt) { - bmat[jp + (i__ - npt) * bmat_dim1] = bmat[i__ + j * - bmat_dim1]; - } - } - } - return 0; -} - -template -static TYPE newuob_(int n, int npt, TYPE *x, - TYPE rhobeg, TYPE rhoend, int *ret_nf, int maxfun, - TYPE *xbase, TYPE *xopt, TYPE *xnew, - TYPE *xpt, TYPE *fval, TYPE *gq, TYPE *hq, - TYPE *pq, TYPE *bmat, TYPE *zmat, int *ndim, - TYPE *d__, TYPE *vlag, TYPE *w, Func &func) -{ - /* XBASE will hold a shift of origin that should reduce the - contributions from rounding errors to values of the model and - Lagrange functions. - * XOPT will be set to the displacement from XBASE of the vector of - variables that provides the least calculated F so far. - * XNEW will be set to the displacement from XBASE of the vector of - variables for the current calculation of F. - * XPT will contain the interpolation point coordinates relative to - XBASE. - * FVAL will hold the values of F at the interpolation points. - * GQ will hold the gradient of the quadratic model at XBASE. - * HQ will hold the explicit second derivatives of the quadratic - model. - * PQ will contain the parameters of the implicit second derivatives - of the quadratic model. - * BMAT will hold the last N columns of H. - * ZMAT will hold the factorization of the leading NPT by NPT - submatrix of H, this factorization being ZMAT times Diag(DZ) - times ZMAT^T, where the elements of DZ are plus or minus 1.0, as - specified by IDZ. - * NDIM is the first dimension of BMAT and has the value NPT+N. - * D is reserved for trial steps from XOPT. - * VLAG will contain the values of the Lagrange functions at a new - point X. They are part of a product that requires VLAG to be of - length NDIM. - * The array W will be used for working space. Its length must be at - least 10*NDIM = 10*(NPT+N). Set some constants. */ - - int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, - i__1, i__2, i__3, i__, j, k, ih, nf, nh, ip, jp, np, nfm, idz, ipt, jpt, - nfmm, knew, kopt, nptm, ksave, nfsav, itemp, ktemp, itest, nftest; - TYPE d__1, d__2, d__3, f, dx, dsq, rho, sum, fbeg, diff, beta, gisq, - temp, suma, sumb, fopt, bsum, gqsq, xipt, xjpt, sumz, diffa, diffb, - diffc, hdiag, alpha, delta, recip, reciq, fsave, dnorm, ratio, dstep, - vquad, tempq, rhosq, detrat, crvmin, distsq, xoptsq; - - /* Parameter adjustments */ - diffc = ratio = dnorm = nfsav = diffa = diffb = xoptsq = f = 0.0; - rho = fbeg = fopt = xjpt = xipt = 0.0; - itest = ipt = jpt = 0; - alpha = dstep = 0.0; - zmat_dim1 = npt; - zmat_offset = 1 + zmat_dim1; - zmat -= zmat_offset; - xpt_dim1 = npt; - xpt_offset = 1 + xpt_dim1; - xpt -= xpt_offset; - --x; --xbase; --xopt; --xnew; --fval; --gq; --hq; --pq; - bmat_dim1 = *ndim; - bmat_offset = 1 + bmat_dim1; - bmat -= bmat_offset; - --d__; - --vlag; - --w; - /* Function Body */ - np = n + 1; - nh = n * np / 2; - nptm = npt - np; - nftest = (maxfun > 1)? maxfun : 1; - /* Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to 0. */ - i__1 = n; - for (j = 1; j <= i__1; ++j) { - xbase[j] = x[j]; - i__2 = npt; - for (k = 1; k <= i__2; ++k) - xpt[k + j * xpt_dim1] = 0; - i__2 = *ndim; - for (i__ = 1; i__ <= i__2; ++i__) - bmat[i__ + j * bmat_dim1] = 0; - } - i__2 = nh; - for (ih = 1; ih <= i__2; ++ih) hq[ih] = 0; - i__2 = npt; - for (k = 1; k <= i__2; ++k) { - pq[k] = 0; - i__1 = nptm; - for (j = 1; j <= i__1; ++j) - zmat[k + j * zmat_dim1] = 0; - } - /* Begin the initialization procedure. NF becomes 1.0 more than the - * number of function values so far. The coordinates of the - * displacement of the next initial interpolation point from XBASE - * are set in XPT(NF,.). */ - rhosq = rhobeg * rhobeg; - recip = 1.0 / rhosq; - reciq = sqrt(.5) / rhosq; - nf = 0; -L50: - nfm = nf; - nfmm = nf - n; - ++nf; - if (nfm <= n << 1) { - if (nfm >= 1 && nfm <= n) { - xpt[nf + nfm * xpt_dim1] = rhobeg; - } else if (nfm > n) { - xpt[nf + nfmm * xpt_dim1] = -(rhobeg); - } - } else { - itemp = (nfmm - 1) / n; - jpt = nfm - itemp * n - n; - ipt = jpt + itemp; - if (ipt > n) { - itemp = jpt; - jpt = ipt - n; - ipt = itemp; - } - xipt = rhobeg; - if (fval[ipt + np] < fval[ipt + 1]) xipt = -xipt; - xjpt = rhobeg; - if (fval[jpt + np] < fval[jpt + 1]) xjpt = -xjpt; - xpt[nf + ipt * xpt_dim1] = xipt; - xpt[nf + jpt * xpt_dim1] = xjpt; - } - /* Calculate the next value of F, label 70 being reached immediately - * after this calculation. The least function value so far and its - * index are required. */ - i__1 = n; - for (j = 1; j <= i__1; ++j) - x[j] = xpt[nf + j * xpt_dim1] + xbase[j]; - goto L310; -L70: - fval[nf] = f; - if (nf == 1) { - fbeg = fopt = f; - kopt = 1; - } else if (f < fopt) { - fopt = f; - kopt = nf; - } - /* Set the non0 initial elements of BMAT and the quadratic model - * in the cases when NF is at most 2*N+1. */ - if (nfm <= n << 1) { - if (nfm >= 1 && nfm <= n) { - gq[nfm] = (f - fbeg) / rhobeg; - if (npt < nf + n) { - bmat[nfm * bmat_dim1 + 1] = -1.0 / rhobeg; - bmat[nf + nfm * bmat_dim1] = 1.0 / rhobeg; - bmat[npt + nfm + nfm * bmat_dim1] = -.5 * rhosq; - } - } else if (nfm > n) { - bmat[nf - n + nfmm * bmat_dim1] = .5 / rhobeg; - bmat[nf + nfmm * bmat_dim1] = -.5 / rhobeg; - zmat[nfmm * zmat_dim1 + 1] = -reciq - reciq; - zmat[nf - n + nfmm * zmat_dim1] = reciq; - zmat[nf + nfmm * zmat_dim1] = reciq; - ih = nfmm * (nfmm + 1) / 2; - temp = (fbeg - f) / rhobeg; - hq[ih] = (gq[nfmm] - temp) / rhobeg; - gq[nfmm] = .5 * (gq[nfmm] + temp); - } - /* Set the off-diagonal second derivatives of the Lagrange - * functions and the initial quadratic model. */ - } else { - ih = ipt * (ipt - 1) / 2 + jpt; - if (xipt < 0) ipt += n; - if (xjpt < 0) jpt += n; - zmat[nfmm * zmat_dim1 + 1] = recip; - zmat[nf + nfmm * zmat_dim1] = recip; - zmat[ipt + 1 + nfmm * zmat_dim1] = -recip; - zmat[jpt + 1 + nfmm * zmat_dim1] = -recip; - hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / (xipt * xjpt); - } - if (nf < npt) goto L50; - /* Begin the iterative procedure, because the initial model is - * complete. */ - rho = rhobeg; - delta = rho; - idz = 1; - diffa = diffb = itest = xoptsq = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - xopt[i__] = xpt[kopt + i__ * xpt_dim1]; - /* Computing 2nd power */ - d__1 = xopt[i__]; - xoptsq += d__1 * d__1; - } -L90: - nfsav = nf; - /* Generate the next trust region step and test its length. Set KNEW - * to -1 if the purpose of the next F will be to improve the - * model. */ -L100: - knew = 0; - trsapp_(n, npt, &xopt[1], &xpt[xpt_offset], &gq[1], &hq[1], &pq[1], & - delta, &d__[1], &w[1], &w[np], &w[np + n], &w[np + (n << 1)], & - crvmin); - dsq = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - /* Computing 2nd power */ - d__1 = d__[i__]; - dsq += d__1 * d__1; - } - /* Computing MIN */ - d__1 = delta, d__2 = sqrt(dsq); - dnorm = std::min(d__1, d__2); - if (dnorm < .5 * rho) { - knew = -1; - delta = DELTA_DECREASE * delta; - ratio = -1.; - if (delta <= rho * 1.5) delta = rho; - if (nf <= nfsav + 2) goto L460; - temp = crvmin * .125 * rho * rho; - /* Computing MAX */ - d__1 = std::max(diffa, diffb); - if (temp <= std::max(d__1, diffc)) goto L460; - goto L490; - } - /* Shift XBASE if XOPT may be too far from XBASE. First make the - * changes to BMAT that do not depend on ZMAT. */ -L120: - if (dsq <= xoptsq * .001) { - tempq = xoptsq * .25; - i__1 = npt; - for (k = 1; k <= i__1; ++k) { - sum = 0; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) - sum += xpt[k + i__ * xpt_dim1] * xopt[i__]; - temp = pq[k] * sum; - sum -= .5 * xoptsq; - w[npt + k] = sum; - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) { - gq[i__] += temp * xpt[k + i__ * xpt_dim1]; - xpt[k + i__ * xpt_dim1] -= .5 * xopt[i__]; - vlag[i__] = bmat[k + i__ * bmat_dim1]; - w[i__] = sum * xpt[k + i__ * xpt_dim1] + tempq * xopt[i__]; - ip = npt + i__; - i__3 = i__; - for (j = 1; j <= i__3; ++j) - bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] + - vlag[i__] * w[j] + w[i__] * vlag[j]; - } - } - /* Then the revisions of BMAT that depend on ZMAT are - * calculated. */ - i__3 = nptm; - for (k = 1; k <= i__3; ++k) { - sumz = 0; - i__2 = npt; - for (i__ = 1; i__ <= i__2; ++i__) { - sumz += zmat[i__ + k * zmat_dim1]; - w[i__] = w[npt + i__] * zmat[i__ + k * zmat_dim1]; - } - i__2 = n; - for (j = 1; j <= i__2; ++j) { - sum = tempq * sumz * xopt[j]; - i__1 = npt; - for (i__ = 1; i__ <= i__1; ++i__) - sum += w[i__] * xpt[i__ + j * xpt_dim1]; - vlag[j] = sum; - if (k < idz) sum = -sum; - i__1 = npt; - for (i__ = 1; i__ <= i__1; ++i__) - bmat[i__ + j * bmat_dim1] += sum * zmat[i__ + k * zmat_dim1]; - } - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - ip = i__ + npt; - temp = vlag[i__]; - if (k < idz) temp = -temp; - i__2 = i__; - for (j = 1; j <= i__2; ++j) - bmat[ip + j * bmat_dim1] += temp * vlag[j]; - } - } - /* The following instructions complete the shift of XBASE, - * including the changes to the parameters of the quadratic - * model. */ - ih = 0; - i__2 = n; - for (j = 1; j <= i__2; ++j) { - w[j] = 0; - i__1 = npt; - for (k = 1; k <= i__1; ++k) { - w[j] += pq[k] * xpt[k + j * xpt_dim1]; - xpt[k + j * xpt_dim1] -= .5 * xopt[j]; - } - i__1 = j; - for (i__ = 1; i__ <= i__1; ++i__) { - ++ih; - if (i__ < j) gq[j] += hq[ih] * xopt[i__]; - gq[i__] += hq[ih] * xopt[j]; - hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j]; - bmat[npt + i__ + j * bmat_dim1] = bmat[npt + j + i__ * - bmat_dim1]; - } - } - i__1 = n; - for (j = 1; j <= i__1; ++j) { - xbase[j] += xopt[j]; - xopt[j] = 0; - } - xoptsq = 0; - } - /* Pick the model step if KNEW is positive. A different choice of D - * may be made later, if the choice of D by BIGLAG causes - * substantial cancellation in DENOM. */ - if (knew > 0) { - biglag_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &zmat[zmat_offset], &idz, - ndim, &knew, &dstep, &d__[1], &alpha, &vlag[1], &vlag[npt + 1], &w[1], &w[np], &w[np + n], func); - } - /* Calculate VLAG and BETA for the current choice of D. The first - * NPT components of W_check will be held in W. */ - i__1 = npt; - for (k = 1; k <= i__1; ++k) { - suma = 0; - sumb = 0; - sum = 0; - i__2 = n; - for (j = 1; j <= i__2; ++j) { - suma += xpt[k + j * xpt_dim1] * d__[j]; - sumb += xpt[k + j * xpt_dim1] * xopt[j]; - sum += bmat[k + j * bmat_dim1] * d__[j]; - } - w[k] = suma * (.5 * suma + sumb); - vlag[k] = sum; - } - beta = 0; - i__1 = nptm; - for (k = 1; k <= i__1; ++k) { - sum = 0; - i__2 = npt; - for (i__ = 1; i__ <= i__2; ++i__) - sum += zmat[i__ + k * zmat_dim1] * w[i__]; - if (k < idz) { - beta += sum * sum; - sum = -sum; - } else beta -= sum * sum; - i__2 = npt; - for (i__ = 1; i__ <= i__2; ++i__) - vlag[i__] += sum * zmat[i__ + k * zmat_dim1]; - } - bsum = 0; - dx = 0; - i__2 = n; - for (j = 1; j <= i__2; ++j) { - sum = 0; - i__1 = npt; - for (i__ = 1; i__ <= i__1; ++i__) - sum += w[i__] * bmat[i__ + j * bmat_dim1]; - bsum += sum * d__[j]; - jp = npt + j; - i__1 = n; - for (k = 1; k <= i__1; ++k) - sum += bmat[jp + k * bmat_dim1] * d__[k]; - vlag[jp] = sum; - bsum += sum * d__[j]; - dx += d__[j] * xopt[j]; - } - beta = dx * dx + dsq * (xoptsq + dx + dx + .5 * dsq) + beta - bsum; - vlag[kopt] += 1.0; - /* If KNEW is positive and if the cancellation in DENOM is - * unacceptable, then BIGDEN calculates an alternative model step, - * XNEW being used for working space. */ - if (knew > 0) { - /* Computing 2nd power */ - d__1 = vlag[knew]; - temp = 1.0 + alpha * beta / (d__1 * d__1); - if (fabs(temp) <= .8) { - bigden_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], & - zmat[zmat_offset], &idz, ndim, &kopt, &knew, &d__[1], &w[ - 1], &vlag[1], &beta, &xnew[1], &w[*ndim + 1], &w[*ndim * - 6 + 1]); - } - } - /* Calculate the next value of the objective function. */ -L290: - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) { - xnew[i__] = xopt[i__] + d__[i__]; - x[i__] = xbase[i__] + xnew[i__]; - } - ++nf; -L310: - if (nf > nftest) { - --nf; - //fprintf(stderr, "++ Return from NEWUOA because CALFUN has been called MAXFUN times.\n"); - goto L530; - } - f = func(n, &x[1]); - if (nf <= npt) goto L70; - if (knew == -1) goto L530; - /* Use the quadratic model to predict the change in F due to the - * step D, and set DIFF to the error of this prediction. */ - vquad = ih = 0; - i__2 = n; - for (j = 1; j <= i__2; ++j) { - vquad += d__[j] * gq[j]; - i__1 = j; - for (i__ = 1; i__ <= i__1; ++i__) { - ++ih; - temp = d__[i__] * xnew[j] + d__[j] * xopt[i__]; - if (i__ == j) temp = .5 * temp; - vquad += temp * hq[ih]; - } - } - i__1 = npt; - for (k = 1; k <= i__1; ++k) vquad += pq[k] * w[k]; - diff = f - fopt - vquad; - diffc = diffb; - diffb = diffa; - diffa = fabs(diff); - if (dnorm > rho) nfsav = nf; - /* Update FOPT and XOPT if the new F is the least value of the - * objective function so far. The branch when KNEW is positive - * occurs if D is not a trust region step. */ - fsave = fopt; - if (f < fopt) { - fopt = f; - xoptsq = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - xopt[i__] = xnew[i__]; - /* Computing 2nd power */ - d__1 = xopt[i__]; - xoptsq += d__1 * d__1; - } - } - ksave = knew; - if (knew > 0) goto L410; - /* Pick the next value of DELTA after a trust region step. */ - if (vquad >= 0) { - fprintf(stderr, "++ Return from NEWUOA because a trust region step has failed to reduce Q.\n"); - goto L530; - } - ratio = (f - fsave) / vquad; - if (ratio <= 0.1) { - delta = .5 * dnorm; - } else if (ratio <= .7) { - /* Computing MAX */ - d__1 = .5 * delta; - delta = std::max(d__1, dnorm); - } else { - /* Computing MAX */ - d__1 = .5 * delta, d__2 = dnorm + dnorm; - delta = std::max(d__1, d__2); - } - if (delta <= rho * 1.5) delta = rho; - /* Set KNEW to the index of the next interpolation point to be - * deleted. */ - /* Computing MAX */ - d__2 = 0.1 * delta; - /* Computing 2nd power */ - d__1 = std::max(d__2, rho); - rhosq = d__1 * d__1; - ktemp = detrat = 0; - if (f >= fsave) { - ktemp = kopt; - detrat = 1.0; - } - i__1 = npt; - for (k = 1; k <= i__1; ++k) { - hdiag = 0; - i__2 = nptm; - for (j = 1; j <= i__2; ++j) { - temp = 1.0; - if (j < idz) temp = -1.0; - /* Computing 2nd power */ - d__1 = zmat[k + j * zmat_dim1]; - hdiag += temp * (d__1 * d__1); - } - /* Computing 2nd power */ - d__2 = vlag[k]; - temp = (d__1 = beta * hdiag + d__2 * d__2, fabs(d__1)); - distsq = 0; - i__2 = n; - for (j = 1; j <= i__2; ++j) { - /* Computing 2nd power */ - d__1 = xpt[k + j * xpt_dim1] - xopt[j]; - distsq += d__1 * d__1; - } - if (distsq > rhosq) { - /* Computing 3rd power */ - d__1 = distsq / rhosq; - temp *= d__1 * (d__1 * d__1); - } - if (temp > detrat && k != ktemp) { - detrat = temp; - knew = k; - } - } - if (knew == 0) goto L460; - /* Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation - * point can be moved. Begin the updating of the quadratic model, - * starting with the explicit second derivative term. */ -L410: - update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], &idz, ndim, &vlag[1], &beta, &knew, &w[1]); - fval[knew] = f; - ih = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - temp = pq[knew] * xpt[knew + i__ * xpt_dim1]; - i__2 = i__; - for (j = 1; j <= i__2; ++j) { - ++ih; - hq[ih] += temp * xpt[knew + j * xpt_dim1]; - } - } - pq[knew] = 0; - /* Update the other second derivative parameters, and then the - * gradient vector of the model. Also include the new interpolation - * point. */ - i__2 = nptm; - for (j = 1; j <= i__2; ++j) { - temp = diff * zmat[knew + j * zmat_dim1]; - if (j < idz) temp = -temp; - i__1 = npt; - for (k = 1; k <= i__1; ++k) { - pq[k] += temp * zmat[k + j * zmat_dim1]; - } - } - gqsq = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - gq[i__] += diff * bmat[knew + i__ * bmat_dim1]; - /* Computing 2nd power */ - d__1 = gq[i__]; - gqsq += d__1 * d__1; - xpt[knew + i__ * xpt_dim1] = xnew[i__]; - } - /* If a trust region step makes a small change to the objective - * function, then calculate the gradient of the least Frobenius norm - * interpolant at XBASE, and store it in W, using VLAG for a vector - * of right hand sides. */ - if (ksave == 0 && delta == rho) { - if (fabs(ratio) > .01) { - itest = 0; - } else { - i__1 = npt; - for (k = 1; k <= i__1; ++k) - vlag[k] = fval[k] - fval[kopt]; - gisq = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - sum = 0; - i__2 = npt; - for (k = 1; k <= i__2; ++k) - sum += bmat[k + i__ * bmat_dim1] * vlag[k]; - gisq += sum * sum; - w[i__] = sum; - } - /* Test whether to replace the new quadratic model by the - * least Frobenius norm interpolant, making the replacement - * if the test is satisfied. */ - ++itest; - if (gqsq < gisq * 100.) itest = 0; - if (itest >= 3) { - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) gq[i__] = w[i__]; - i__1 = nh; - for (ih = 1; ih <= i__1; ++ih) hq[ih] = 0; - i__1 = nptm; - for (j = 1; j <= i__1; ++j) { - w[j] = 0; - i__2 = npt; - for (k = 1; k <= i__2; ++k) - w[j] += vlag[k] * zmat[k + j * zmat_dim1]; - if (j < idz) w[j] = -w[j]; - } - i__1 = npt; - for (k = 1; k <= i__1; ++k) { - pq[k] = 0; - i__2 = nptm; - for (j = 1; j <= i__2; ++j) - pq[k] += zmat[k + j * zmat_dim1] * w[j]; - } - itest = 0; - } - } - } - if (f < fsave) kopt = knew; - /* If a trust region step has provided a sufficient decrease in F, - * then branch for another trust region calculation. The case - * KSAVE>0 occurs when the new function value was calculated by a - * model step. */ - if (f <= fsave + 0.1 * vquad) goto L100; - if (ksave > 0) goto L100; - /* Alternatively, find out if the interpolation points are close - * enough to the best point so far. */ - knew = 0; -L460: - distsq = delta * 4. * delta; - i__2 = npt; - for (k = 1; k <= i__2; ++k) { - sum = 0; - i__1 = n; - for (j = 1; j <= i__1; ++j) { - /* Computing 2nd power */ - d__1 = xpt[k + j * xpt_dim1] - xopt[j]; - sum += d__1 * d__1; - } - if (sum > distsq) { - knew = k; - distsq = sum; - } - } - /* If KNEW is positive, then set DSTEP, and branch back for the next - * iteration, which will generate a "model step". */ - if (knew > 0) { - /* Computing MAX and MIN*/ - d__2 = 0.1 * sqrt(distsq), d__3 = .5 * delta; - d__1 = std::min(d__2, d__3); - dstep = std::max(d__1, rho); - dsq = dstep * dstep; - goto L120; - } - if (ratio > 0) goto L100; - if (std::max(delta, dnorm) > rho) goto L100; - /* The calculations with the current value of RHO are complete. Pick - * the next values of RHO and DELTA. */ -L490: - if (rho > rhoend) { - delta = .5 * rho; - ratio = rho / rhoend; - if (ratio <= 16.) rho = rhoend; - else if (ratio <= 250.) rho = sqrt(ratio) * rhoend; - else rho = RHO_DECREASE * rho; - delta = std::max(delta, rho); - goto L90; - } - /* Return from the calculation, after another Newton-Raphson step, - * if it is too short to have been tried before. */ - if (knew == -1) goto L290; -L530: - if (fopt <= f) { - i__2 = n; - for (i__ = 1; i__ <= i__2; ++i__) - x[i__] = xbase[i__] + xopt[i__]; - f = fopt; - } - *ret_nf = nf; - return f; -} - -template -static TYPE newuoa_(int n, int npt, TYPE *x, TYPE rhobeg, TYPE rhoend, int *ret_nf, int maxfun, TYPE *w, Func &func) -{ - /* This subroutine seeks the least value of a function of many - * variables, by a trust region method that forms quadratic models - * by interpolation. There can be some freedom in the interpolation - * conditions, which is taken up by minimizing the Frobenius norm of - * the change to the second derivative of the quadratic model, - * beginning with a zero matrix. The arguments of the subroutine are - * as follows. */ - - /* N must be set to the number of variables and must be at least - * two. NPT is the number of interpolation conditions. Its value - * must be in the interval [N+2,(N+1)(N+2)/2]. Initial values of the - * variables must be set in X(1),X(2),...,X(N). They will be changed - * to the values that give the least calculated F. RHOBEG and RHOEND - * must be set to the initial and final values of a trust region - * radius, so both must be positive with RHOEND<=RHOBEG. Typically - * RHOBEG should be about one tenth of the greatest expected change - * to a variable, and RHOEND should indicate the accuracy that is - * required in the final values of the variables. MAXFUN must be set - * to an upper bound on the number of calls of CALFUN. The array W - * will be used for working space. Its length must be at least - * (NPT+13)*(NPT+N)+3*N*(N+3)/2. */ - - /* SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must - * set F to the value of the objective function for the variables - * X(1),X(2),...,X(N). Partition the working space array, so that - * different parts of it can be treated separately by the subroutine - * that performs the main calculation. */ - - int id, np, iw, igq, ihq, ixb, ifv, ipq, ivl, ixn, - ixo, ixp, ndim, nptm, ibmat, izmat; - - /* Parameter adjustments */ - --w; --x; - /* Function Body */ - np = n + 1; - nptm = npt - np; - if (npt < n + 2 || npt > (n + 2) * np / 2) { - fprintf(stderr, "** Return from NEWUOA because NPT is not in the required interval.\n"); - return 1; - } - ndim = npt + n; - ixb = 1; - ixo = ixb + n; - ixn = ixo + n; - ixp = ixn + n; - ifv = ixp + n * npt; - igq = ifv + npt; - ihq = igq + n; - ipq = ihq + n * np / 2; - ibmat = ipq + npt; - izmat = ibmat + ndim * n; - id = izmat + npt * nptm; - ivl = id + n; - iw = ivl + ndim; - /* The above settings provide a partition of W for subroutine - * NEWUOB. The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 - * elements of W plus the space that is needed by the last array of - * NEWUOB. */ - return newuob_(n, npt, &x[1], rhobeg, rhoend, ret_nf, maxfun, &w[ixb], &w[ixo], &w[ixn], - &w[ixp], &w[ifv], &w[igq], &w[ihq], &w[ipq], &w[ibmat], &w[izmat], - &ndim, &w[id], &w[ivl], &w[iw], func); -} - -template -TYPE min_newuoa(int n, TYPE *x, Func &func, TYPE rb, TYPE tol, int max_iter) -{ - int npt = 2 * n + 1, rnf; - TYPE ret; - TYPE *w = (TYPE*)calloc((npt+13)*(npt+n) + 3*n*(n+3)/2 + 11, sizeof(TYPE)); - ret = newuoa_(n, 2*n+1, x, rb, tol, &rnf, max_iter, w, func); - free(w); - return ret; -} - -#endif diff --git a/src/meshlabplugins/edit_align/edit_align.pro b/src/meshlabplugins/edit_align/edit_align.pro index 7c9baabc0..162135295 100644 --- a/src/meshlabplugins/edit_align/edit_align.pro +++ b/src/meshlabplugins/edit_align/edit_align.pro @@ -1,7 +1,7 @@ include (../../shared.pri) INCLUDEPATH += \ - $$MESHLAB_EXTERNAL_DIRECTORY/newuoa/include + $$VCGDIR/wrap/newuoa/include HEADERS += \ edit_align_factory.h \ diff --git a/src/meshlabplugins/edit_mutualcorrs/edit_mutualcorrs.pro b/src/meshlabplugins/edit_mutualcorrs/edit_mutualcorrs.pro index b2786a5b1..59f93820e 100644 --- a/src/meshlabplugins/edit_mutualcorrs/edit_mutualcorrs.pro +++ b/src/meshlabplugins/edit_mutualcorrs/edit_mutualcorrs.pro @@ -31,7 +31,7 @@ TARGET = edit_mutualcorrs INCLUDEPATH *= \ $$MESHLAB_EXTERNAL_DIRECTORY/levmar-2.3 \ - $$MESHLAB_EXTERNAL_DIRECTORY/newuoa/include + $$VCGDIR/wrap/newuoa/include win32-msvc:LIBS += $$MESHLAB_DISTRIB_DIRECTORY/lib/win32-msvc/levmar.lib win32-g++:LIBS += -L$$MESHLAB_DISTRIB_DIRECTORY/lib/win32-gcc -llevmar diff --git a/src/meshlabplugins/filter_isoparametrization/CMakeLists.txt b/src/meshlabplugins/filter_isoparametrization/CMakeLists.txt index b97aed5d2..0eb090ab4 100644 --- a/src/meshlabplugins/filter_isoparametrization/CMakeLists.txt +++ b/src/meshlabplugins/filter_isoparametrization/CMakeLists.txt @@ -5,8 +5,8 @@ ### specifically src/templates/filter_isoparametrization.cmake (custom for this directory), ### then re-run ./make-cmake.py -# Only build if we have newuoa and levmar -if(TARGET external-newuoa AND TARGET external-levmar) +# Only build if we have levmar +if(TARGET external-levmar) set(SOURCES filter_isoparametrization.cpp) @@ -64,6 +64,6 @@ if(TARGET external-newuoa AND TARGET external-levmar) else() message( STATUS - "Skipping filter_isoparametrization - missing either newuoa or levmar in external directory." + "Skipping filter_isoparametrization - missing levmar in external directory." ) endif() diff --git a/src/meshlabplugins/filter_isoparametrization/filter_isoparametrization.pro b/src/meshlabplugins/filter_isoparametrization/filter_isoparametrization.pro index f50a0e9c3..25294c62d 100644 --- a/src/meshlabplugins/filter_isoparametrization/filter_isoparametrization.pro +++ b/src/meshlabplugins/filter_isoparametrization/filter_isoparametrization.pro @@ -27,7 +27,7 @@ SOURCES += \ TARGET = filter_isoparametrization INCLUDEPATH += \ $$MESHLAB_EXTERNAL_DIRECTORY/levmar-2.3/ - + win32-msvc:QMAKE_CXXFLAGS += /openmp -D_USE_OMP win32-g++:QMAKE_LFLAGS += -fopenmp linux:QMAKE_CXXFLAGS += -fopenmp -D_USE_OMP diff --git a/src/meshlabplugins/filter_mutualglobal/filter_mutualglobal.pro b/src/meshlabplugins/filter_mutualglobal/filter_mutualglobal.pro index 20cd6176e..01d0fa67a 100644 --- a/src/meshlabplugins/filter_mutualglobal/filter_mutualglobal.pro +++ b/src/meshlabplugins/filter_mutualglobal/filter_mutualglobal.pro @@ -19,9 +19,9 @@ SOURCES += \ parameters.cpp \ pointCorrespondence.cpp \ solver.cpp - + TARGET = filter_mutualglobal INCLUDEPATH *= \ $$MESHLAB_EXTERNAL_DIRECTORY/levmar-2.3 \ - $$MESHLAB_EXTERNAL_DIRECTORY/newuoa/include + $$VCGDIR/wrap/newuoa/include diff --git a/src/meshlabplugins/filter_mutualinfo/filter_mutualinfo.pro b/src/meshlabplugins/filter_mutualinfo/filter_mutualinfo.pro index dc0343551..eb6bc88f3 100644 --- a/src/meshlabplugins/filter_mutualinfo/filter_mutualinfo.pro +++ b/src/meshlabplugins/filter_mutualinfo/filter_mutualinfo.pro @@ -19,9 +19,9 @@ SOURCES += \ parameters.cpp \ pointCorrespondence.cpp \ solver.cpp - + TARGET = filter_mutualinfo INCLUDEPATH *= \ $$MESHLAB_EXTERNAL_DIRECTORY/levmar-2.3 \ - $$MESHLAB_EXTERNAL_DIRECTORY/newuoa/include + $$VCGDIR/wrap/newuoa/include diff --git a/vcglib b/vcglib index 0caaf49d3..e7ce2614e 160000 --- a/vcglib +++ b/vcglib @@ -1 +1 @@ -Subproject commit 0caaf49d37263ea9422536a6aeb978bd1336bd03 +Subproject commit e7ce2614ec1002c969ff85d9c31563a879025213 From 0c7d1de24169c7e60a901dd1422757d0675e83c8 Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Fri, 29 May 2020 11:56:12 +0200 Subject: [PATCH 02/11] first move aling_pair to vcg --- .../edit_align/align/AlignPair.cpp | 141 +++---- .../edit_align/align/AlignPair.h | 384 +----------------- vcglib | 2 +- 3 files changed, 68 insertions(+), 459 deletions(-) diff --git a/src/meshlabplugins/edit_align/align/AlignPair.cpp b/src/meshlabplugins/edit_align/align/AlignPair.cpp index b27925178..05e7b2519 100644 --- a/src/meshlabplugins/edit_align/align/AlignPair.cpp +++ b/src/meshlabplugins/edit_align/align/AlignPair.cpp @@ -41,95 +41,86 @@ using namespace vcg; using namespace std; -bool AlignPair::A2Mesh::Import(const char *filename, Matrix44d &Tr) -{ - int err = tri::io::Importer::Open(*this, filename); - if (err) { - printf("Error in reading %s: '%s'\n", filename, tri::io::Importer::ErrorMsg(err)); - exit(-1); - } - printf("read mesh `%s'\n", filename); - return Init(Tr); -} +//bool AlignPair::A2Mesh::Import(const char *filename, Matrix44d &Tr) +//{ +// int err = tri::io::Importer::Open(*this, filename); +// if (err) { +// printf("Error in reading %s: '%s'\n", filename, tri::io::Importer::ErrorMsg(err)); +// exit(-1); +// } +// printf("read mesh `%s'\n", filename); +// return Init(Tr); +//} -bool AlignPair::A2Mesh::InitVert(const Matrix44d &Tr) -{ - Matrix44d Idn; Idn.SetIdentity(); - if (Tr != Idn) tri::UpdatePosition::Matrix(*this, Tr); - tri::UpdateNormal::NormalizePerVertex(*this); - tri::UpdateBounding::Box(*this); - return true; -} +//bool AlignPair::A2Mesh::Init(const Matrix44d &Tr) +//{ +// Matrix44d Idn; Idn.SetIdentity(); +// tri::Clean::RemoveUnreferencedVertex(*this); +// if (Tr != Idn) tri::UpdatePosition::Matrix(*this, Tr); +// tri::UpdateNormal::PerVertexNormalizedPerFaceNormalized(*this); +// tri::UpdateFlags::FaceBorderFromNone(*this); +// tri::UpdateBounding::Box(*this); -bool AlignPair::A2Mesh::Init(const Matrix44d &Tr) -{ - Matrix44d Idn; Idn.SetIdentity(); - tri::Clean::RemoveUnreferencedVertex(*this); - if (Tr != Idn) tri::UpdatePosition::Matrix(*this, Tr); - tri::UpdateNormal::PerVertexNormalizedPerFaceNormalized(*this); - tri::UpdateFlags::FaceBorderFromNone(*this); - tri::UpdateBounding::Box(*this); - - return true; -} +// return true; +//} -void AlignPair::Stat::clear() -{ - I.clear(); - StartTime = 0; - MovVertNum = 0; - FixVertNum = 0; - FixFaceNum = 0; -} +//void AlignPair::Stat::clear() +//{ +// I.clear(); +// StartTime = 0; +// MovVertNum = 0; +// FixVertNum = 0; +// FixFaceNum = 0; +//} // Restituisce true se nelle ultime iterazioni non e' diminuito // l'errore -bool AlignPair::Stat::Stable(int lastiter) -{ - if (I.empty()) return false; - int parag = int(I.size()) - lastiter; +//bool AlignPair::Stat::Stable(int lastiter) +//{ +// if (I.empty()) return false; +// int parag = int(I.size()) - lastiter; - if (parag < 0) parag = 0; - if (I.back().pcl50 < I[parag].pcl50) return false; // se siamo diminuiti non e' stabile +// if (parag < 0) parag = 0; +// if (I.back().pcl50 < I[parag].pcl50) return false; // se siamo diminuiti non e' stabile - return true; +// return true; -} +//} -void AlignPair::Stat::Dump(FILE *fp) -{ - if (I.size() == 0) { - fprintf(fp, "Empty AlignPair::Stat\n"); - return; - } - fprintf(fp, "Final Err %8.5f In %i iterations Total Time %ims\n", LastPcl50(), (int)I.size(), TotTime()); - fprintf(fp, "Mindist Med Hi Avg RMS StdDev Time Tested Used Dist Bord Angl \n"); - for (unsigned int qi = 0; qi < I.size(); ++qi) - fprintf(fp, "%5.2f (%6.3f:%6.3f) (%6.3f %6.3f %6.3f) %4ims %5i %5i %4i+%4i+%4i\n", - I[qi].MinDistAbs, - I[qi].pcl50, I[qi].pclhi, - I[qi].AVG, I[qi].RMS, I[qi].StdDev, - IterTime(qi), - I[qi].SampleTested, I[qi].SampleUsed, I[qi].DistanceDiscarded, I[qi].BorderDiscarded, I[qi].AngleDiscarded); -} +//void AlignPair::Stat::Dump(FILE *fp) +//{ +// if (I.size() == 0) { +// fprintf(fp, "Empty AlignPair::Stat\n"); +// return; +// } +// fprintf(fp, "Final Err %8.5f In %i iterations Total Time %ims\n", LastPcl50(), (int)I.size(), TotTime()); +// fprintf(fp, "Mindist Med Hi Avg RMS StdDev Time Tested Used Dist Bord Angl \n"); +// for (unsigned int qi = 0; qi < I.size(); ++qi) +// fprintf(fp, "%5.2f (%6.3f:%6.3f) (%6.3f %6.3f %6.3f) %4ims %5i %5i %4i+%4i+%4i\n", +// I[qi].MinDistAbs, +// I[qi].pcl50, I[qi].pclhi, +// I[qi].AVG, I[qi].RMS, I[qi].StdDev, +// IterTime(qi), +// I[qi].SampleTested, I[qi].SampleUsed, I[qi].DistanceDiscarded, I[qi].BorderDiscarded, I[qi].AngleDiscarded); +//} // Scrive una tabella con tutti i valori -void AlignPair::Stat::HTMLDump(FILE *fp) -{ - fprintf(fp, "Final Err %8.5f In %i iterations Total Time %ims\n", LastPcl50(), (int)I.size(), TotTime()); - fprintf(fp, "\n"); - fprintf(fp, "\n", - I[qi].MinDistAbs, - I[qi].pcl50, I[qi].pclhi, - I[qi].AVG, I[qi].RMS, I[qi].StdDev, - IterTime(qi), - I[qi].SampleTested, I[qi].SampleUsed, I[qi].DistanceDiscarded, I[qi].BorderDiscarded, I[qi].AngleDiscarded); - fprintf(fp, "
Mindist 50ile Hi Avg RMS StdDev Time Tested Used Dist Bord Angl \n"); - for (unsigned int qi = 0; qi < I.size(); ++qi) - fprintf(fp, "
%8.5f %9.6f %8.5f %5.3f %8.5f %9.6f %4ims %5i %5i %4i %4i %4i
\n"); -} +//void AlignPair::Stat::HTMLDump(FILE *fp) +//{ +// fprintf(fp, "Final Err %8.5f In %i iterations Total Time %ims\n", LastPcl50(), (int)I.size(), TotTime()); +// fprintf(fp, "\n"); +// fprintf(fp, "\n", +// I[qi].MinDistAbs, +// I[qi].pcl50, I[qi].pclhi, +// I[qi].AVG, I[qi].RMS, I[qi].StdDev, +// IterTime(qi), +// I[qi].SampleTested, I[qi].SampleUsed, I[qi].DistanceDiscarded, I[qi].BorderDiscarded, I[qi].AngleDiscarded); +// fprintf(fp, "
Mindist 50ile Hi Avg RMS StdDev Time Tested Used Dist Bord Angl \n"); +// for (unsigned int qi = 0; qi < I.size(); ++qi) +// fprintf(fp, "
%8.5f %9.6f %8.5f %5.3f %8.5f %9.6f %4ims %5i %5i %4i %4i %4i
\n"); +//} diff --git a/src/meshlabplugins/edit_align/align/AlignPair.h b/src/meshlabplugins/edit_align/align/AlignPair.h index 769d70393..2ecb45146 100644 --- a/src/meshlabplugins/edit_align/align/AlignPair.h +++ b/src/meshlabplugins/edit_align/align/AlignPair.h @@ -24,390 +24,8 @@ #ifndef __VCG_ALIGNPAIR #define __VCG_ALIGNPAIR -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include -namespace vcg -{ -/************************************************************************* - AlignPair - -Classe per gestire un allineamento tra DUE sole mesh. - -**************************************************************************/ - -class AlignPair -{ -public: - - enum ErrorCode {SUCCESS, NO_COMMON_BBOX, TOO_FEW_POINTS, - LSQ_DIVERGE, TOO_MUCH_SHEAR, TOO_MUCH_SCALE, FORBIDDEN, INVALID, UNKNOWN_MODE }; - - -/*********************** Classi Accessorie ****************************/ - -class A2Vertex; - -class A2Face ; - -class A2UsedTypes: public vcg::UsedTypes < vcg::Use::AsVertexType, - vcg::Use::AsFaceType >{}; - -class A2Vertex : public vcg::Vertex {}; -class A2Face : public vcg::Face< A2UsedTypes,vcg::face::VertexRef, vcg::face::Normal3d,vcg::face::Mark,vcg::face::BitFlags> {}; -class A2Mesh : public vcg::tri::TriMesh< std::vector, std::vector > -{ -public: - bool Import(const char *filename) { Matrix44d Tr; Tr.SetIdentity(); return Import(filename,Tr);} - bool Import(const char *filename, Matrix44d &Tr); - bool InitVert(const Matrix44d &Tr); - bool Init(const Matrix44d &Tr); -}; - -typedef A2Mesh::FaceContainer FaceContainer; -typedef A2Mesh::FaceType FaceType; -typedef GridStaticPtr A2Grid; -typedef GridStaticPtr A2GridVert; - -class Stat -{ -public: - - class IterInfo - { - public: - IterInfo() {memset ( (void *) this, 0, sizeof(IterInfo)); } - - double MinDistAbs; - int DistanceDiscarded; - int AngleDiscarded; - int BorderDiscarded; - int SampleTested; // quanti punti ho testato con la mindist - int SampleUsed; // quanti punti ho scelto per la computematrix - double pcl50; - double pclhi; - double AVG; - double RMS; - double StdDev; - int Time; // quando e' finita questa iterazione - - }; - std::vector I; - - double LastPcl50() const {return I.back().pcl50;} - int LastSampleUsed() const {return I.back().SampleUsed;} - - int MovVertNum; - int FixVertNum; - int FixFaceNum; - int TotTime() { return I.back().Time-StartTime; } - int IterTime(unsigned int i) const - { const int clock_per_ms = std::max(CLOCKS_PER_SEC / 1000,1); - assert(i - // e almeno da usare per l'allineamento. - int MaxPointNum; // numero di coppie di punti da usare quando si calcola la matrice - // di allienamento e che poi si mettono da parte per il globale; - // di solito non usato - int MinPointNum; // minimo numero di coppie di punti ammissibile perche' sia considerato - // valido l'allineamento - double MinDistAbs; // distanza minima iniziale perche due punti vengano presi in - // considerazione. NON e' piu espressa in percentuale sul bbox della mesh nella ug. - // Ad ogni passo puo essere ridotta per - // accellerare usando ReduceFactor - double MaxAngleRad; // massimo angolo in radianti tra due normali perche' i due - // punti vengano presi in considerazione. - - int MaxIterNum; // massimo numero di iterazioni da fare in align - double TrgDistAbs; // distanza obiettivo entro la quale devono stare almeno la meta' - // dei campioni scelti; di solito non entra in gioco perche' ha un default molto basso - - int EndStepNum; // numero di iterazioni da considerare per decidere se icp ha converso. - - //double PassLoFilter; // Filtraggio utilizzato per decidere quali punti scegliere tra quello trovati abbastanza - double PassHiFilter; // vicini. Espresso in percentili. Di solito si scarta il quelli sopra il 75 e quelli sotto il 5 - double ReduceFactorPerc; // At each step we discard the points farther than a given threshold. The threshold is iterativeley reduced; - // StartMinDist= min(StartMinDist, 5.0*H.Percentile(ap.ReduceFactorPerc)) - - - double MinMinDistPerc; // Ratio between initial starting distance (MinDistAbs) and what can reach by the application of the ReduceFactor. - - int UGExpansionFactor; // Grandezza della UG per la mesh fix come rapporto del numero di facce della mesh fix - - // Nel caso si usi qualche struttura multiresolution - int MinFixVertNum; // Gli allineamenti si fanno mettendo nella ug come mesh fix una semplificata; - float MinFixVertNumPerc; // si usa il max tra MinFixVertNum e OrigMeshSize*MinFixVertNumPerc - bool UseVertexOnly; // if true all the Alignment pipeline ignores faces and works over point clouds. - - double MaxShear; - double MaxScale; - MatchModeEnum MatchMode; - SampleModeEnum SampleMode; - void Dump(FILE *fp,double BoxDiag); - -}; - -// Classe per memorizzare il risultato di un allineamento tra due mesh -// i punti si intendono nel sistema di riferimento iniziale delle due mesh. -// -// se le mesh hanno una trasformazione di base in ingresso, -// questa appare solo durante la A2Mesh::Import e poi e' per sempre dimenticata. -// Questi punti sono quindi nei sistemi di riferimento costruiti durante la Import -// la matrice Tr quella che -// -// Tr*Pmov[i]== Pfix - - -class Result -{ -public: - int MovName; - int FixName; - - Matrix44d Tr; - std::vector Pfix; // vertici corrispondenti su fix (rossi) - std::vector Nfix; // normali corrispondenti su fix (rossi) - std::vector Pmov; // vertici scelti su mov (verdi) prima della trasformazione in ingresso (Original Point Target) - std::vector Nmov; // normali scelti su mov (verdi) - Histogramf H; - Stat as; - Param ap; - ErrorCode status; - bool IsValid() {return status==SUCCESS;} - double err; - float area; // the overlapping area, a percentage as computed in Occupancy Grid. - - bool operator < (const Result & rr) const {return (err< rr.err);} - bool operator <= (const Result & rr) const {return (err<=rr.err);} - bool operator > (const Result & rr) const {return (err> rr.err);} - bool operator >= (const Result & rr) const {return (err>=rr.err);} - bool operator == (const Result & rr) const {return (err==rr.err);} - bool operator != (const Result & rr) const {return (err!=rr.err);} - - std::pair ComputeAvgErr() const - { - double sum_before=0; - double sum_after=0; - for(unsigned int ii=0;ii *mov; - A2Mesh *fix; - - ErrorCode status; - AlignPair::Param ap; - - math::SubtractiveRingRNG myrnd; - -/**** End Data Members *********/ - -template < class MESH > -void ConvertMesh(MESH &M1, A2Mesh &M2) -{ - tri::Append::MeshCopy(M2,M1); -} - -template < class VERTEX > -void ConvertVertex(const std::vector &vert1, std::vector &vert2, Box3d *Clip=0) -{ - - vert2.clear(); - typename std::vector::const_iterator vi; - A2Vertex tv; - Box3 bb; - if(Clip){ - bb.Import(*Clip); - for(vi=vert1.begin();vi &vert, int SampleNum, AlignPair::Param::SampleModeEnum SampleMode); -bool SampleMovVertRandom(std::vector &vert, int SampleNum); -bool SampleMovVertNormalEqualized(std::vector &vert, int SampleNum); -/* -Una volta trovati coppie di punti corrispondenti se ne sceglie -al piu' per calcolare la trasformazione che li fa coincidere. -La scelta viene fatta in base ai due parametri PassLo e PassHi; -*/ -bool ChoosePoints( std::vector &Ps, // vertici corrispondenti su fix (rossi) - std::vector &Ns, // normali corrispondenti su fix (rossi) - std::vector &Pt, // vertici scelti su mov (verdi) - std::vector &OPt, // vertici scelti su mov (verdi) - //vector &Nt, // normali scelti su mov (verdi) - double PassHi, - Histogramf &H) -; - - -/* - Minimo esempio di codice per l'uso della funzione di allineamento. - - AlignPair::A2Mesh Mov,Fix; // le due mesh da allineare - vector MovVert; // i punti sulla mov da usare per l'allineamento - Matrix44d In; In.SetIdentity(); // la trasformazione iniziale che applicata a mov la porterebbe su fix. - - AlignPair aa; // l'oggetto principale. - AlignPair::Param ap; - UGrid< AlignPair::A2Mesh::face_container > UG; - - Fix.LoadPly("FixMesh.ply"); // Standard ply loading - Mov.LoadPly("MovMesh.ply"); - Fix.Init( Ident, false); // Inizializzazione necessaria (normali per vert, - Mov.Init( Ident, false); // info per distanza punto faccia ecc) - - AlignPair::InitFix(&Fix, ap, UG); // la mesh fix viene messa nella ug. - - aa.ConvertVertex(Mov.vert,MovVert); // si campiona la mesh Mov per trovare un po' di vertici. - aa.SampleMovVert(MovVert, ap.SampleNum, ap.SampleMode); - - aa.mov=&MovVert; // si assegnano i membri principali dell'oggetto align pair - aa.fix=&Fix; - aa.ap = ap; - - aa.Align(In,UG,res); // si spera :) - // il risultato sta nella matrice res.Tr; - - res.as.Dump(stdout); - -*/ - -bool Align(const Matrix44d &in, A2Grid &UG, A2GridVert &UGV, Result &res) -{ - res.ap=ap; - - bool ret=Align( UG,UGV, - in, - res.Tr, - res.Pfix, res.Nfix, - res.Pmov, res.Nmov, - res.H, - res.as); - - res.err=res.as.LastPcl50(); - res.status=status; - return ret; -} - -double Abs2Perc(double val, Box3d bb) const {return val/bb.Diag();} -double Perc2Abs(double val, Box3d bb) const {return val*bb.Diag();} - -/************************************************************************************ -Versione Vera della Align a basso livello. - -Si assume che la mesh fix sia gia' stata messa nella ug u con le debite trasformazioni. -in - -************************************************************************************/ - -bool Align(A2Grid &u, - A2GridVert &uv, - const Matrix44d &in, // trasformazione Iniziale che porta i punti di mov su fix - Matrix44d &out, // trasformazione calcolata - std::vector &Pfix, // vertici corrispondenti su src (rossi) - std::vector &Nfix, // normali corrispondenti su src (rossi) - std::vector &OPmov, // vertici scelti su trg (verdi) prima della trasformazione in ingresso (Original Point Target) - std::vector &ONmov, // normali scelti su trg (verdi) - Histogramf &H, - Stat &as); - - -bool InitMov( - std::vector< Point3d > &movvert, - std::vector< Point3d > &movnorm, - Box3d &movbox, - const Matrix44d &in ); - -static bool InitFixVert(A2Mesh *fm, - AlignPair::Param &pp, - A2GridVert &u, - int PreferredGridSize=0); - -static bool InitFix(A2Mesh *fm, - AlignPair::Param &pp, - A2Grid &u, - int PreferredGridSize=0); - -}; // end class - -} // end namespace vcg #endif diff --git a/vcglib b/vcglib index e7ce2614e..ec730298f 160000 --- a/vcglib +++ b/vcglib @@ -1 +1 @@ -Subproject commit e7ce2614ec1002c969ff85d9c31563a879025213 +Subproject commit ec730298faad058ab513f03abeebb4180bb893a5 From 4a88e5c923f4deecf4805f3aa67b99b67b6a077f Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Fri, 29 May 2020 12:19:04 +0200 Subject: [PATCH 03/11] moved point_matching_scale to vcg --- src/meshlabplugins/edit_align/CMakeLists.txt | 1 - .../edit_align/align/AlignGlobal.cpp | 4 +- .../edit_align/align/AlignPair.cpp | 4 +- src/meshlabplugins/edit_align/edit_align.pro | 1 - .../edit_align/point_matching_scale.cpp | 102 ------------------ .../edit_align/point_matching_scale.h | 18 ---- src/templates/edit_align.cmake | 1 - vcglib | 2 +- 8 files changed, 5 insertions(+), 128 deletions(-) delete mode 100644 src/meshlabplugins/edit_align/point_matching_scale.cpp delete mode 100644 src/meshlabplugins/edit_align/point_matching_scale.h diff --git a/src/meshlabplugins/edit_align/CMakeLists.txt b/src/meshlabplugins/edit_align/CMakeLists.txt index f0699be5b..d2e4e5625 100644 --- a/src/meshlabplugins/edit_align/CMakeLists.txt +++ b/src/meshlabplugins/edit_align/CMakeLists.txt @@ -19,7 +19,6 @@ if(TARGET external-newuoa) edit_align_factory.cpp edit_align.cpp meshtree.cpp - point_matching_scale.cpp ../../meshlab/stdpardialog.cpp ${VCGDIR}/wrap/gui/trackball.cpp ${VCGDIR}/wrap/gui/trackmode.cpp diff --git a/src/meshlabplugins/edit_align/align/AlignGlobal.cpp b/src/meshlabplugins/edit_align/align/AlignGlobal.cpp index 128392ec7..a7c5ea590 100644 --- a/src/meshlabplugins/edit_align/align/AlignGlobal.cpp +++ b/src/meshlabplugins/edit_align/align/AlignGlobal.cpp @@ -25,7 +25,7 @@ #include #include "AlignPair.h" #include "AlignGlobal.h" -#include "../point_matching_scale.h" +#include "vcg/complex/algorithms/point_matching_scale.h" #include using namespace vcg; using namespace std; @@ -337,7 +337,7 @@ double AlignGlobal::Node::AlignWithActiveAdj(bool Rigid) //if(Rigid) ret=ComputeRigidMatchMatrix(out,FixP,MovP); //else ret=ComputeMatchMatrix2(out,FixP,FixN,MovP); if(Rigid) ComputeRigidMatchMatrix(FixP,MovP,out); - else ComputeRotoTranslationScalingMatchMatrix(out,FixP,MovP); + else PointMatchingScale::ComputeRotoTranslationScalingMatchMatrix(out,FixP,MovP); Matrix44d outIn=vcg::Inverse(out); //double maxdiff = MatrixNorm(out); diff --git a/src/meshlabplugins/edit_align/align/AlignPair.cpp b/src/meshlabplugins/edit_align/align/AlignPair.cpp index 05e7b2519..6398125ca 100644 --- a/src/meshlabplugins/edit_align/align/AlignPair.cpp +++ b/src/meshlabplugins/edit_align/align/AlignPair.cpp @@ -22,7 +22,7 @@ ****************************************************************************/ #include "AlignPair.h" -#include "../point_matching_scale.h" +#include "vcg/complex/algorithms/point_matching_scale.h" #include #include @@ -394,7 +394,7 @@ bool AlignPair::Align( Matrix44d newout; switch (ap.MatchMode) { - case AlignPair::Param::MMSimilarity: ComputeRotoTranslationScalingMatchMatrix(newout, Pfix, Pmov); break; + case AlignPair::Param::MMSimilarity: PointMatchingScale::ComputeRotoTranslationScalingMatchMatrix(newout, Pfix, Pmov); break; case AlignPair::Param::MMRigid: ComputeRigidMatchMatrix(Pfix, Pmov, newout); break; default: status = UNKNOWN_MODE; diff --git a/src/meshlabplugins/edit_align/edit_align.pro b/src/meshlabplugins/edit_align/edit_align.pro index 162135295..92b23235d 100644 --- a/src/meshlabplugins/edit_align/edit_align.pro +++ b/src/meshlabplugins/edit_align/edit_align.pro @@ -21,7 +21,6 @@ HEADERS += \ SOURCES += \ edit_align_factory.cpp \ - point_matching_scale.cpp \ edit_align.cpp \ meshtree.cpp \ alignDialog.cpp \ diff --git a/src/meshlabplugins/edit_align/point_matching_scale.cpp b/src/meshlabplugins/edit_align/point_matching_scale.cpp deleted file mode 100644 index e13a0d463..000000000 --- a/src/meshlabplugins/edit_align/point_matching_scale.cpp +++ /dev/null @@ -1,102 +0,0 @@ -#include "point_matching_scale.h" -#include - -#include -#include - -using namespace vcg; - - -template -struct RotoTranslation -{ - RotoTranslation(){} - Scalar _v[6]; - void ToMatrix(vcg::Matrix44 & m) - { - vcg::Matrix44 rot,tra; - rot.FromEulerAngles(_v[0],_v[1],_v[2]); - tra.SetTranslate(vcg::Point3(_v[3],_v[4],_v[5])); - m = tra * rot; - } -}; - - - -static std::vector *fix; -static std::vector *mov; -static vcg::Box3d b; - -double errorScale(int n, double *x) -{ - assert(n==1); (void)n; - double dist = 0; - std::vector::iterator i = mov->begin(); - std::vector::iterator ifix = fix->begin(); - for(; i != mov->end(); ++i,++ifix) - dist += vcg::SquaredDistance(((*i)-b.Center())*(*x)+b.Center() , *ifix); - - return dist; -} - - -void ComputeScalingMatchMatrix(Matrix44d &res, - std::vector &Pfix, - std::vector &Pmov) -{ - fix = &Pfix; - mov = &Pmov; - b.SetNull(); - for(std::vector::iterator i = Pmov.begin(); i != Pmov.end(); ++i) - b.Add(*i); - - double scale = 1.0; - min_newuoa(1,&scale,errorScale); - - res.SetTranslate( b.Center()*(1.0-scale)); - res[0][0] = res[1][1] = res[2][2] = scale; -} - - -double errorRotoTranslationScale(int n, double *x){ - assert(n==7); (void)n; - double dist = 0; - std::vector::iterator i = mov->begin(); - std::vector::iterator ifix = fix->begin(); - - RotoTranslation rt; - vcg::Matrix44d m; - memcpy(&rt._v[0],&x[1],6*sizeof(double)); - rt.ToMatrix(m); - - for(; i != mov->end(); ++i,++ifix){ - dist += vcg::SquaredDistance( m*(((*i)-b.Center())*(x[0])+b.Center()),*ifix); - } - return dist; -} - -void ComputeRotoTranslationScalingMatchMatrix(vcg::Matrix44d &res, - std::vector &Pfix, - std::vector &Pmov) -{ - fix = &Pfix; - mov = &Pmov; - b.SetNull(); - for(std::vector::iterator i = Pmov.begin(); i != Pmov.end(); ++i) - b.Add(*i); - - double x[7]={1.0,0.0,0.0,0.0,0.0,0.0,0.0}; - min_newuoa(7,&x[0],errorRotoTranslationScale); - - // rtm = rototranslation - RotoTranslation rt; - vcg::Matrix44d rtm; - memcpy(&rt._v[0],&x[1],6*sizeof(double)); - rt.ToMatrix(rtm); - - // res= scaling w.r.t. barycenter - res.SetTranslate( b.Center()*(1.0-x[0])); - res[0][0] = res[1][1] = res[2][2] = x[0]; - res = rtm*res; -} - diff --git a/src/meshlabplugins/edit_align/point_matching_scale.h b/src/meshlabplugins/edit_align/point_matching_scale.h deleted file mode 100644 index b461e6931..000000000 --- a/src/meshlabplugins/edit_align/point_matching_scale.h +++ /dev/null @@ -1,18 +0,0 @@ -#include -#include - -/** -Compute a scaling transformation that bring PMov point as close as possible to Pfix -*/ -void ComputeScalingMatchMatrix(vcg::Matrix44d &res, - std::vector &Pfix, - std::vector &Pmov); - -/** -Compute a rototranslation + scaling transformation that bring PMov point as close as possible to Pfix -*/ - - -void ComputeRotoTranslationScalingMatchMatrix(vcg::Matrix44d &res, - std::vector &Pfix, - std::vector &Pmov); diff --git a/src/templates/edit_align.cmake b/src/templates/edit_align.cmake index bc4d0fffc..00b3f17e1 100644 --- a/src/templates/edit_align.cmake +++ b/src/templates/edit_align.cmake @@ -20,7 +20,6 @@ AlignPairWidget.cpp edit_align_factory.cpp edit_align.cpp meshtree.cpp -point_matching_scale.cpp ../../meshlab/stdpardialog.cpp ${VCGDIR}/wrap/gui/trackball.cpp ${VCGDIR}/wrap/gui/trackmode.cpp diff --git a/vcglib b/vcglib index ec730298f..6cd9d7aa9 160000 --- a/vcglib +++ b/vcglib @@ -1 +1 @@ -Subproject commit ec730298faad058ab513f03abeebb4180bb893a5 +Subproject commit 6cd9d7aa91ac7860031446a4e1f2ae0c46bd42ee From 4d8b57b07008ecec5e47150a8031cbf0ad63cb4c Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Fri, 29 May 2020 12:25:43 +0200 Subject: [PATCH 04/11] fix last commit --- src/meshlabplugins/edit_align/align/AlignGlobal.cpp | 2 +- src/meshlabplugins/edit_align/align/AlignPair.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/meshlabplugins/edit_align/align/AlignGlobal.cpp b/src/meshlabplugins/edit_align/align/AlignGlobal.cpp index a7c5ea590..913ab4469 100644 --- a/src/meshlabplugins/edit_align/align/AlignGlobal.cpp +++ b/src/meshlabplugins/edit_align/align/AlignGlobal.cpp @@ -337,7 +337,7 @@ double AlignGlobal::Node::AlignWithActiveAdj(bool Rigid) //if(Rigid) ret=ComputeRigidMatchMatrix(out,FixP,MovP); //else ret=ComputeMatchMatrix2(out,FixP,FixN,MovP); if(Rigid) ComputeRigidMatchMatrix(FixP,MovP,out); - else PointMatchingScale::ComputeRotoTranslationScalingMatchMatrix(out,FixP,MovP); + else PointMatchingScale::computeRotoTranslationScalingMatchMatrix(out,FixP,MovP); Matrix44d outIn=vcg::Inverse(out); //double maxdiff = MatrixNorm(out); diff --git a/src/meshlabplugins/edit_align/align/AlignPair.cpp b/src/meshlabplugins/edit_align/align/AlignPair.cpp index 6398125ca..034fbf2b3 100644 --- a/src/meshlabplugins/edit_align/align/AlignPair.cpp +++ b/src/meshlabplugins/edit_align/align/AlignPair.cpp @@ -394,7 +394,7 @@ bool AlignPair::Align( Matrix44d newout; switch (ap.MatchMode) { - case AlignPair::Param::MMSimilarity: PointMatchingScale::ComputeRotoTranslationScalingMatchMatrix(newout, Pfix, Pmov); break; + case AlignPair::Param::MMSimilarity: PointMatchingScale::computeRotoTranslationScalingMatchMatrix(newout, Pfix, Pmov); break; case AlignPair::Param::MMRigid: ComputeRigidMatchMatrix(Pfix, Pmov, newout); break; default: status = UNKNOWN_MODE; From 9f2021e624cfd75e41b1844e91311accbef93ade Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Fri, 29 May 2020 13:06:39 +0200 Subject: [PATCH 05/11] fix last commit --- src/meshlabplugins/edit_align/align/AlignGlobal.cpp | 2 +- src/meshlabplugins/edit_align/align/AlignPair.cpp | 2 +- src/meshlabplugins/edit_align/edit_align.cpp | 8 ++++++++ 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/meshlabplugins/edit_align/align/AlignGlobal.cpp b/src/meshlabplugins/edit_align/align/AlignGlobal.cpp index 913ab4469..d610cb830 100644 --- a/src/meshlabplugins/edit_align/align/AlignGlobal.cpp +++ b/src/meshlabplugins/edit_align/align/AlignGlobal.cpp @@ -337,7 +337,7 @@ double AlignGlobal::Node::AlignWithActiveAdj(bool Rigid) //if(Rigid) ret=ComputeRigidMatchMatrix(out,FixP,MovP); //else ret=ComputeMatchMatrix2(out,FixP,FixN,MovP); if(Rigid) ComputeRigidMatchMatrix(FixP,MovP,out); - else PointMatchingScale::computeRotoTranslationScalingMatchMatrix(out,FixP,MovP); + else vcg::PointMatchingScale::computeRotoTranslationScalingMatchMatrix(out,FixP,MovP); Matrix44d outIn=vcg::Inverse(out); //double maxdiff = MatrixNorm(out); diff --git a/src/meshlabplugins/edit_align/align/AlignPair.cpp b/src/meshlabplugins/edit_align/align/AlignPair.cpp index 034fbf2b3..4b575d645 100644 --- a/src/meshlabplugins/edit_align/align/AlignPair.cpp +++ b/src/meshlabplugins/edit_align/align/AlignPair.cpp @@ -394,7 +394,7 @@ bool AlignPair::Align( Matrix44d newout; switch (ap.MatchMode) { - case AlignPair::Param::MMSimilarity: PointMatchingScale::computeRotoTranslationScalingMatchMatrix(newout, Pfix, Pmov); break; + case AlignPair::Param::MMSimilarity: vcg::PointMatchingScale::computeRotoTranslationScalingMatchMatrix(newout, Pfix, Pmov); break; case AlignPair::Param::MMRigid: ComputeRigidMatchMatrix(Pfix, Pmov, newout); break; default: status = UNKNOWN_MODE; diff --git a/src/meshlabplugins/edit_align/edit_align.cpp b/src/meshlabplugins/edit_align/edit_align.cpp index 058b91d2b..3a24c6f0b 100644 --- a/src/meshlabplugins/edit_align/edit_align.cpp +++ b/src/meshlabplugins/edit_align/edit_align.cpp @@ -34,8 +34,16 @@ $Log: meshedit.cpp,v $ #include "AlignPairDialog.h" #include "align/align_parameter.h" #include +#include + using namespace vcg; +//todo: remove these orrible defs from here +//make vcg::PointMatchingScale indipendent +std::vector* vcg::PointMatchingScale::fix; +std::vector* vcg::PointMatchingScale::mov; +vcg::Box3d vcg::PointMatchingScale::b; + EditAlignPlugin::EditAlignPlugin() { alignDialog=0; From 3f3d688eb956cef54755c8977ed2808d24e7ce36 Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Fri, 29 May 2020 13:21:14 +0200 Subject: [PATCH 06/11] more implementation and some refactoring --- .../edit_align/align/AlignPair.cpp | 30 +++++++++---------- src/meshlabplugins/edit_align/alignDialog.cpp | 2 +- src/meshlabplugins/edit_align/meshtree.cpp | 12 ++++---- vcglib | 2 +- 4 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/meshlabplugins/edit_align/align/AlignPair.cpp b/src/meshlabplugins/edit_align/align/AlignPair.cpp index 4b575d645..2cfdb1211 100644 --- a/src/meshlabplugins/edit_align/align/AlignPair.cpp +++ b/src/meshlabplugins/edit_align/align/AlignPair.cpp @@ -434,7 +434,7 @@ bool AlignPair::Align( } while ( nc <= ap.MaxIterNum && H.Percentile(.5) > ap.TrgDistAbs && - (ncupdateDataMask(MeshModel::MM_FACEMARK); - aa.ConvertMesh(MM(fixId)->cm,Fix); + aa.convertMesh(MM(fixId)->cm,Fix); vcg::AlignPair::A2Grid UG; vcg::AlignPair::A2GridVert VG; if(MM(fixId)->cm.fn==0 || ap.UseVertexOnly) { - Fix.InitVert(vcg::Matrix44d::Identity()); + Fix.initVert(vcg::Matrix44d::Identity()); vcg::AlignPair::InitFixVert(&Fix,ap,VG); } else { - Fix.Init(vcg::Matrix44d::Identity()); + Fix.init(vcg::Matrix44d::Identity()); vcg::AlignPair::InitFix(&Fix, ap, UG); } // 2) Convert the second mesh and sample a points on it. @@ -89,7 +89,7 @@ void MeshTree::ProcessArc(int fixId, int movId, vcg::Matrix44d &MovM, vcg::Align // aa.ConvertMesh(MM(movId)->cm,Mov); std::vector tmpmv; // aa.ConvertVertex(Mov.vert,tmpmv); - aa.ConvertVertex(MM(movId)->cm.vert,tmpmv); + aa.convertVertex(MM(movId)->cm.vert,tmpmv); aa.SampleMovVert(tmpmv, ap.SampleNum, ap.SampleMode); aa.mov=&tmpmv; @@ -188,14 +188,14 @@ void MeshTree::Process(vcg::AlignPair::Param &ap, MeshTree::Param &mtp) if( curResult->IsValid() ) { hasValidAlign = true; - std::pair dd=curResult->ComputeAvgErr(); + std::pair dd=curResult->computeAvgErr(); #pragma omp critical cb(0,qUtf8Printable(buf.sprintf("(%3i/%3zu) %2i -> %2i Aligned AvgErr dd=%f -> dd=%f \n",i+1,totalArcNum,OG.SVA[i].s,OG.SVA[i].t,dd.first,dd.second))); } else { #pragma omp critical - cb(0,qUtf8Printable(buf.sprintf( "(%3i/%3zu) %2i -> %2i Failed Alignment of one arc %s\n",i+1,totalArcNum,OG.SVA[i].s,OG.SVA[i].t,vcg::AlignPair::ErrorMsg(curResult->status)))); + cb(0,qUtf8Printable(buf.sprintf( "(%3i/%3zu) %2i -> %2i Failed Alignment of one arc %s\n",i+1,totalArcNum,OG.SVA[i].s,OG.SVA[i].t,vcg::AlignPair::errorMsg(curResult->status)))); } } } diff --git a/vcglib b/vcglib index 6cd9d7aa9..a374e959e 160000 --- a/vcglib +++ b/vcglib @@ -1 +1 @@ -Subproject commit 6cd9d7aa91ac7860031446a4e1f2ae0c46bd42ee +Subproject commit a374e959ee6ffc8fc968b4c253fa3bf9c2b06130 From 518faf7c5f4a7885b058b99f99f11c550d5359ed Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Fri, 29 May 2020 14:48:59 +0200 Subject: [PATCH 07/11] more implementations imported to vcg --- .../edit_align/align/AlignPair.cpp | 584 +++++++++--------- src/meshlabplugins/edit_align/meshtree.cpp | 10 +- vcglib | 2 +- 3 files changed, 298 insertions(+), 298 deletions(-) diff --git a/src/meshlabplugins/edit_align/align/AlignPair.cpp b/src/meshlabplugins/edit_align/align/AlignPair.cpp index 2cfdb1211..c774a4082 100644 --- a/src/meshlabplugins/edit_align/align/AlignPair.cpp +++ b/src/meshlabplugins/edit_align/align/AlignPair.cpp @@ -132,57 +132,57 @@ MaxPointNum an (unused) hard limit on the number of points that are chosen MinPointNum the minimum number of points that have to be chosen to be usable */ -bool AlignPair::ChoosePoints(vector &Ps, // vertici corrispondenti su src (rossi) - vector &Ns, // normali corrispondenti su src (rossi) - vector &Pt, // vertici scelti su trg (verdi) - vector &OPt, // vertici scelti su trg (verdi) - double PassHi, - Histogramf &H) -{ - const int N = ap.MaxPointNum; - double newmaxd = H.Percentile(float(PassHi)); - //printf("%5.1f of the pairs has a distance less than %g and greater than %g (0..%g) avg %g\n", Perc*100,newmind,newmaxd,H.maxv,H.Percentile(.5)); - int sz = int(Ps.size()); - int fnd = 0; - int lastgood = sz - 1; - math::SubtractiveRingRNG myrnd; - while (fnd < N && fnd < lastgood) - { - int index = fnd + myrnd.generate(lastgood - fnd); - double dd = Distance(Ps[index], Pt[index]); - if (dd <= newmaxd) - { - swap(Ps[index], Ps[fnd]); - swap(Ns[index], Ns[fnd]); - swap(Pt[index], Pt[fnd]); - swap(OPt[index], OPt[fnd]); - ++fnd; - } - else - { - swap(Ps[index], Ps[lastgood]); - swap(Ns[index], Ns[lastgood]); - swap(Pt[index], Pt[lastgood]); - swap(OPt[index], OPt[lastgood]); - lastgood--; - } - } - Ps.resize(fnd); - Ns.resize(fnd); - Pt.resize(fnd); - OPt.resize(fnd); -// printf("Scelte %i coppie tra le %i iniziali, scartate quelle con dist > %f\n", fnd, sz, newmaxd); +//bool AlignPair::ChoosePoints(vector &Ps, // vertici corrispondenti su src (rossi) +// vector &Ns, // normali corrispondenti su src (rossi) +// vector &Pt, // vertici scelti su trg (verdi) +// vector &OPt, // vertici scelti su trg (verdi) +// double PassHi, +// Histogramf &H) +//{ +// const int N = ap.MaxPointNum; +// double newmaxd = H.Percentile(float(PassHi)); +// //printf("%5.1f of the pairs has a distance less than %g and greater than %g (0..%g) avg %g\n", Perc*100,newmind,newmaxd,H.maxv,H.Percentile(.5)); +// int sz = int(Ps.size()); +// int fnd = 0; +// int lastgood = sz - 1; +// math::SubtractiveRingRNG myrnd; +// while (fnd < N && fnd < lastgood) +// { +// int index = fnd + myrnd.generate(lastgood - fnd); +// double dd = Distance(Ps[index], Pt[index]); +// if (dd <= newmaxd) +// { +// swap(Ps[index], Ps[fnd]); +// swap(Ns[index], Ns[fnd]); +// swap(Pt[index], Pt[fnd]); +// swap(OPt[index], OPt[fnd]); +// ++fnd; +// } +// else +// { +// swap(Ps[index], Ps[lastgood]); +// swap(Ns[index], Ns[lastgood]); +// swap(Pt[index], Pt[lastgood]); +// swap(OPt[index], OPt[lastgood]); +// lastgood--; +// } +// } +// Ps.resize(fnd); +// Ns.resize(fnd); +// Pt.resize(fnd); +// OPt.resize(fnd); +//// printf("Scelte %i coppie tra le %i iniziali, scartate quelle con dist > %f\n", fnd, sz, newmaxd); - if ((int)Ps.size() < ap.MinPointNum) { - printf("Troppi pochi punti!\n"); - Ps.clear(); - Ns.clear(); - Pt.clear(); - OPt.clear(); - return false; - } - return true; -} +// if ((int)Ps.size() < ap.MinPointNum) { +// printf("Troppi pochi punti!\n"); +// Ps.clear(); +// Ns.clear(); +// Pt.clear(); +// OPt.clear(); +// return false; +// } +// return true; +//} /* Funzione chiamata dalla Align ad ogni ciclo @@ -270,193 +270,193 @@ la uniform grid sia gia' inizializzata con la mesh fix */ -bool AlignPair::Align( - A2Grid &u, - A2GridVert &uv, - const Matrix44d &in, // trasformazione Iniziale (che porta i punti di mov su fix) - Matrix44d &out, // trasformazione calcolata - vector &Pfix, // vertici corrispondenti su src (rossi) - vector &Nfix, // normali corrispondenti su src (rossi) - vector &OPmov, // vertici scelti su trg (verdi) prima della trasformazione in ingresso (Original Point Target) - vector &ONmov, // normali scelti su trg (verdi) - Histogramf &H, - AlignPair::Stat &as) -{ - vector beyondCntVec; // vettore per marcare i movvert che sicuramente non si devono usare - // ogni volta che un vertice si trova a distanza oltre max dist viene incrementato il suo contatore; - // i movvert che sono stati scartati piu' di MaxCntDist volte non si guardano piu'; - const int maxBeyondCnt = 3; - vector< Point3d > movvert; - vector< Point3d > movnorm; - vector Pmov; // vertici scelti dopo la trasf iniziale - status = SUCCESS; - int tt0 = clock(); +//bool AlignPair::align( +// A2Grid &u, +// A2GridVert &uv, +// const Matrix44d &in, // trasformazione Iniziale (che porta i punti di mov su fix) +// Matrix44d &out, // trasformazione calcolata +// vector &Pfix, // vertici corrispondenti su src (rossi) +// vector &Nfix, // normali corrispondenti su src (rossi) +// vector &OPmov, // vertici scelti su trg (verdi) prima della trasformazione in ingresso (Original Point Target) +// vector &ONmov, // normali scelti su trg (verdi) +// Histogramf &H, +// AlignPair::Stat &as) +//{ +// vector beyondCntVec; // vettore per marcare i movvert che sicuramente non si devono usare +// // ogni volta che un vertice si trova a distanza oltre max dist viene incrementato il suo contatore; +// // i movvert che sono stati scartati piu' di MaxCntDist volte non si guardano piu'; +// const int maxBeyondCnt = 3; +// vector< Point3d > movvert; +// vector< Point3d > movnorm; +// vector Pmov; // vertici scelti dopo la trasf iniziale +// status = SUCCESS; +// int tt0 = clock(); - out = in; +// out = in; - int i; +// int i; - double CosAngleThr = cos(ap.MaxAngleRad); - double StartMinDist = ap.MinDistAbs; - int tt1 = clock(); - int ttsearch = 0; - int ttleast = 0; - int nc = 0; - as.clear(); - as.StartTime = clock(); +// double CosAngleThr = cos(ap.MaxAngleRad); +// double StartMinDist = ap.MinDistAbs; +// int tt1 = clock(); +// int ttsearch = 0; +// int ttleast = 0; +// int nc = 0; +// as.clear(); +// as.StartTime = clock(); - beyondCntVec.resize(mov->size(), 0); +// beyondCntVec.resize(mov->size(), 0); - /**************** BEGIN ICP LOOP ****************/ - do - { - Stat::IterInfo ii; - Box3d movbox; - InitMov(movvert, movnorm, movbox, out); - H.SetRange(0.0f, float(StartMinDist), 512, 2.5f); - Pfix.clear(); - Nfix.clear(); - Pmov.clear(); - OPmov.clear(); - ONmov.clear(); - int tts0 = clock(); - ii.MinDistAbs = StartMinDist; - int LocSampleNum = min(ap.SampleNum, int(movvert.size())); - Box3d fixbox; - if (u.Empty()) fixbox = uv.bbox; - else fixbox = u.bbox; - for (i = 0; i < LocSampleNum; ++i) - { - if (beyondCntVec[i] < maxBeyondCnt) - { +// /**************** BEGIN ICP LOOP ****************/ +// do +// { +// Stat::IterInfo ii; +// Box3d movbox; +// InitMov(movvert, movnorm, movbox, out); +// H.SetRange(0.0f, float(StartMinDist), 512, 2.5f); +// Pfix.clear(); +// Nfix.clear(); +// Pmov.clear(); +// OPmov.clear(); +// ONmov.clear(); +// int tts0 = clock(); +// ii.MinDistAbs = StartMinDist; +// int LocSampleNum = min(ap.SampleNum, int(movvert.size())); +// Box3d fixbox; +// if (u.Empty()) fixbox = uv.bbox; +// else fixbox = u.bbox; +// for (i = 0; i < LocSampleNum; ++i) +// { +// if (beyondCntVec[i] < maxBeyondCnt) +// { - if (fixbox.IsIn(movvert[i])) - { - double error = StartMinDist; - Point3d closestPoint, closestNormal; - double maxd = StartMinDist; - ii.SampleTested++; - if (u.Empty()) // using the point cloud grid - { - A2Mesh::VertexPointer vp = tri::GetClosestVertex(*fix, uv, movvert[i], maxd, error); - if (error >= StartMinDist) { - ii.DistanceDiscarded++; ++beyondCntVec[i]; continue; - } - if (movnorm[i].dot(vp->N()) < CosAngleThr) { - ii.AngleDiscarded++; continue; - } - closestPoint = vp->P(); - closestNormal = vp->N(); - } - else // using the standard faces and grid - { - A2Mesh::FacePointer f = vcg::tri::GetClosestFaceBase(*fix, u, movvert[i], maxd, error, closestPoint); - if (error >= StartMinDist) { - ii.DistanceDiscarded++; ++beyondCntVec[i]; continue; - } - if (movnorm[i].dot(f->N()) < CosAngleThr) { - ii.AngleDiscarded++; continue; - } - Point3d ip; - InterpolationParameters(*f, f->N(), closestPoint, ip); - const double IP_EPS = 0.00001; - // If ip[i] == 0 it means that we are on the edge opposite to i - if ((fabs(ip[0]) <= IP_EPS && f->IsB(1)) || (fabs(ip[1]) <= IP_EPS && f->IsB(2)) || (fabs(ip[2]) <= IP_EPS && f->IsB(0))){ - ii.BorderDiscarded++; continue; - } - closestNormal = f->N(); - } - // The sample was accepted. Store it. - Pmov.push_back(movvert[i]); - OPmov.push_back((*mov)[i].P()); - ONmov.push_back((*mov)[i].N()); - Nfix.push_back(closestNormal); - Pfix.push_back(closestPoint); - H.Add(float(error)); - ii.SampleUsed++; - } - else - beyondCntVec[i] = maxBeyondCnt + 1; - } - } // End for each pmov - int tts1 = clock(); - //printf("Found %d pairs\n",(int)Pfix.size()); - if (!ChoosePoints(Pfix, Nfix, Pmov, OPmov, ap.PassHiFilter, H)) - { - if (int(Pfix.size()) < ap.MinPointNum) - { - status = TOO_FEW_POINTS; - ii.Time = clock(); - as.I.push_back(ii); - return false; - } - } - Matrix44d newout; - switch (ap.MatchMode) - { - case AlignPair::Param::MMSimilarity: vcg::PointMatchingScale::computeRotoTranslationScalingMatchMatrix(newout, Pfix, Pmov); break; - case AlignPair::Param::MMRigid: ComputeRigidMatchMatrix(Pfix, Pmov, newout); break; - default: - status = UNKNOWN_MODE; - ii.Time = clock(); - as.I.push_back(ii); - return false; - } +// if (fixbox.IsIn(movvert[i])) +// { +// double error = StartMinDist; +// Point3d closestPoint, closestNormal; +// double maxd = StartMinDist; +// ii.SampleTested++; +// if (u.Empty()) // using the point cloud grid +// { +// A2Mesh::VertexPointer vp = tri::GetClosestVertex(*fix, uv, movvert[i], maxd, error); +// if (error >= StartMinDist) { +// ii.DistanceDiscarded++; ++beyondCntVec[i]; continue; +// } +// if (movnorm[i].dot(vp->N()) < CosAngleThr) { +// ii.AngleDiscarded++; continue; +// } +// closestPoint = vp->P(); +// closestNormal = vp->N(); +// } +// else // using the standard faces and grid +// { +// A2Mesh::FacePointer f = vcg::tri::GetClosestFaceBase(*fix, u, movvert[i], maxd, error, closestPoint); +// if (error >= StartMinDist) { +// ii.DistanceDiscarded++; ++beyondCntVec[i]; continue; +// } +// if (movnorm[i].dot(f->N()) < CosAngleThr) { +// ii.AngleDiscarded++; continue; +// } +// Point3d ip; +// InterpolationParameters(*f, f->N(), closestPoint, ip); +// const double IP_EPS = 0.00001; +// // If ip[i] == 0 it means that we are on the edge opposite to i +// if ((fabs(ip[0]) <= IP_EPS && f->IsB(1)) || (fabs(ip[1]) <= IP_EPS && f->IsB(2)) || (fabs(ip[2]) <= IP_EPS && f->IsB(0))){ +// ii.BorderDiscarded++; continue; +// } +// closestNormal = f->N(); +// } +// // The sample was accepted. Store it. +// Pmov.push_back(movvert[i]); +// OPmov.push_back((*mov)[i].P()); +// ONmov.push_back((*mov)[i].N()); +// Nfix.push_back(closestNormal); +// Pfix.push_back(closestPoint); +// H.Add(float(error)); +// ii.SampleUsed++; +// } +// else +// beyondCntVec[i] = maxBeyondCnt + 1; +// } +// } // End for each pmov +// int tts1 = clock(); +// //printf("Found %d pairs\n",(int)Pfix.size()); +// if (!choosePoints(Pfix, Nfix, Pmov, OPmov, ap.PassHiFilter, H)) +// { +// if (int(Pfix.size()) < ap.MinPointNum) +// { +// status = TOO_FEW_POINTS; +// ii.Time = clock(); +// as.I.push_back(ii); +// return false; +// } +// } +// Matrix44d newout; +// switch (ap.MatchMode) +// { +// case AlignPair::Param::MMSimilarity: vcg::PointMatchingScale::computeRotoTranslationScalingMatchMatrix(newout, Pfix, Pmov); break; +// case AlignPair::Param::MMRigid: ComputeRigidMatchMatrix(Pfix, Pmov, newout); break; +// default: +// status = UNKNOWN_MODE; +// ii.Time = clock(); +// as.I.push_back(ii); +// return false; +// } - // double sum_before=0; - // double sum_after=0; - // for(unsigned int iii=0;iii %f\n",sum_before/double(Pfix.size()),sum_after/double(Pfix.size()) ) ; +// // double sum_before=0; +// // double sum_after=0; +// // for(unsigned int iii=0;iii %f\n",sum_before/double(Pfix.size()),sum_after/double(Pfix.size()) ) ; - // le passate successive utilizzano quindi come trasformazione iniziale questa appena trovata. - // Nei prossimi cicli si parte da questa matrice come iniziale. - out = newout * out; +// // le passate successive utilizzano quindi come trasformazione iniziale questa appena trovata. +// // Nei prossimi cicli si parte da questa matrice come iniziale. +// out = newout * out; - assert(Pfix.size() == Pmov.size()); - int tts2 = clock(); - ttsearch += tts1 - tts0; - ttleast += tts2 - tts1; - ii.pcl50 = H.Percentile(.5); - ii.pclhi = H.Percentile(float(ap.PassHiFilter)); - ii.AVG = H.Avg(); - ii.RMS = H.RMS(); - ii.StdDev = H.StandardDeviation(); - ii.Time = clock(); - as.I.push_back(ii); - nc++; - // The distance of the next points to be considered is lowered according to the parameter. - // We use 5 times the percentile of the found points. - if (ap.ReduceFactorPerc<1) StartMinDist = max(ap.MinDistAbs*ap.MinMinDistPerc, min(StartMinDist, 5.0*H.Percentile(float(ap.ReduceFactorPerc)))); - } while ( - nc <= ap.MaxIterNum && - H.Percentile(.5) > ap.TrgDistAbs && - (ncap.MaxScale || math::Abs(1 - scv[1]) > ap.MaxScale || math::Abs(1 - scv[2]) > ap.MaxScale)) { - status = TOO_MUCH_SCALE; - return false; - } - if (shv[0] > ap.MaxShear || shv[1] > ap.MaxShear || shv[2] > ap.MaxShear) { - status = TOO_MUCH_SHEAR; - return false; - } - printf("Grid %i %i %i - fn %i\n", u.siz[0], u.siz[1], u.siz[2], fix->fn); - printf("Init %8.3f Loop %8.3f Search %8.3f least sqrt %8.3f\n", - float(tt1 - tt0) / CLOCKS_PER_SEC, float(tt2 - tt1) / CLOCKS_PER_SEC, - float(ttsearch) / CLOCKS_PER_SEC, float(ttleast) / CLOCKS_PER_SEC); +// assert(Pfix.size() == Pmov.size()); +// int tts2 = clock(); +// ttsearch += tts1 - tts0; +// ttleast += tts2 - tts1; +// ii.pcl50 = H.Percentile(.5); +// ii.pclhi = H.Percentile(float(ap.PassHiFilter)); +// ii.AVG = H.Avg(); +// ii.RMS = H.RMS(); +// ii.StdDev = H.StandardDeviation(); +// ii.Time = clock(); +// as.I.push_back(ii); +// nc++; +// // The distance of the next points to be considered is lowered according to the parameter. +// // We use 5 times the percentile of the found points. +// if (ap.ReduceFactorPerc<1) StartMinDist = max(ap.MinDistAbs*ap.MinMinDistPerc, min(StartMinDist, 5.0*H.Percentile(float(ap.ReduceFactorPerc)))); +// } while ( +// nc <= ap.MaxIterNum && +// H.Percentile(.5) > ap.TrgDistAbs && +// (ncap.MaxScale || math::Abs(1 - scv[1]) > ap.MaxScale || math::Abs(1 - scv[2]) > ap.MaxScale)) { +// status = TOO_MUCH_SCALE; +// return false; +// } +// if (shv[0] > ap.MaxShear || shv[1] > ap.MaxShear || shv[2] > ap.MaxShear) { +// status = TOO_MUCH_SHEAR; +// return false; +// } +// printf("Grid %i %i %i - fn %i\n", u.siz[0], u.siz[1], u.siz[2], fix->fn); +// printf("Init %8.3f Loop %8.3f Search %8.3f least sqrt %8.3f\n", +// float(tt1 - tt0) / CLOCKS_PER_SEC, float(tt2 - tt1) / CLOCKS_PER_SEC, +// float(ttsearch) / CLOCKS_PER_SEC, float(ttleast) / CLOCKS_PER_SEC); - return true; -} +// return true; +//} @@ -619,32 +619,32 @@ return maxfnd; /**********************************************************/ // Funzioni per la scelta dei vertici sulla mesh da muovere -bool AlignPair::SampleMovVert(vector &vert, int SampleNum, AlignPair::Param::SampleModeEnum SampleMode) -{ - switch (SampleMode) - { - case AlignPair::Param::SMRandom: return SampleMovVertRandom(vert, SampleNum); - case AlignPair::Param::SMNormalEqualized: return SampleMovVertNormalEqualized(vert, SampleNum); - default: assert(0); - } - return false; -} +//bool AlignPair::SampleMovVert(vector &vert, int SampleNum, AlignPair::Param::SampleModeEnum SampleMode) +//{ +// switch (SampleMode) +// { +// case AlignPair::Param::SMRandom: return SampleMovVertRandom(vert, SampleNum); +// case AlignPair::Param::SMNormalEqualized: return SampleMovVertNormalEqualized(vert, SampleNum); +// default: assert(0); +// } +// return false; +//} // Scelta a caso semplice -bool AlignPair::SampleMovVertRandom(vector &vert, int SampleNum) -{ - if (int(vert.size()) <= SampleNum) return true; - int i; - for (i = 0; i < SampleNum; ++i) - { - int pos = myrnd.generate(vert.size()); - assert(pos >= 0 && pos < int(vert.size())); - swap(vert[i], vert[pos]); - } - vert.resize(SampleNum); - return true; -} +//bool AlignPair::SampleMovVertRandom(vector &vert, int SampleNum) +//{ +// if (int(vert.size()) <= SampleNum) return true; +// int i; +// for (i = 0; i < SampleNum; ++i) +// { +// int pos = myrnd.generate(vert.size()); +// assert(pos >= 0 && pos < int(vert.size())); +// swap(vert[i], vert[pos]); +// } +// vert.resize(SampleNum); +// return true; +//} /* Scelta a caso in maniera tale che la distribuzione delle normali dei @@ -658,53 +658,53 @@ e poi un punto all'interno del bucket */ -bool AlignPair::SampleMovVertNormalEqualized(vector &vert, int SampleNum) -{ - // assert(0); +//bool AlignPair::SampleMovVertNormalEqualized(vector &vert, int SampleNum) +//{ +// // assert(0); - // int t0=clock(); - vector NV; - if (NV.size() == 0) - { - GenNormal::Fibonacci(30, NV); - printf("Generated %i normals\n", int(NV.size())); - } - // Bucket vector dove, per ogni normale metto gli indici - // dei vertici ad essa corrispondenti - vector > BKT(NV.size()); - for (size_t i = 0; i < vert.size(); ++i) - { - int ind = GenNormal::BestMatchingNormal(vert[i].N(), NV); - BKT[ind].push_back(int(i)); - } - //int t1=clock(); +// // int t0=clock(); +// vector NV; +// if (NV.size() == 0) +// { +// GenNormal::Fibonacci(30, NV); +// printf("Generated %i normals\n", int(NV.size())); +// } +// // Bucket vector dove, per ogni normale metto gli indici +// // dei vertici ad essa corrispondenti +// vector > BKT(NV.size()); +// for (size_t i = 0; i < vert.size(); ++i) +// { +// int ind = GenNormal::BestMatchingNormal(vert[i].N(), NV); +// BKT[ind].push_back(int(i)); +// } +// //int t1=clock(); - // vettore di contatori per sapere quanti punti ho gia' preso per ogni bucket - vector BKTpos(BKT.size(), 0); +// // vettore di contatori per sapere quanti punti ho gia' preso per ogni bucket +// vector BKTpos(BKT.size(), 0); - if (SampleNum >= int(vert.size())) SampleNum = vert.size() - 1; +// if (SampleNum >= int(vert.size())) SampleNum = vert.size() - 1; - for (int i = 0; i < SampleNum;) - { - int ind = myrnd.generate(BKT.size()); // Scelgo un Bucket - int &CURpos = BKTpos[ind]; - vector &CUR = BKT[ind]; +// for (int i = 0; i < SampleNum;) +// { +// int ind = myrnd.generate(BKT.size()); // Scelgo un Bucket +// int &CURpos = BKTpos[ind]; +// vector &CUR = BKT[ind]; - if (CURpos tmpmv; // aa.ConvertVertex(Mov.vert,tmpmv); aa.convertVertex(MM(movId)->cm.vert,tmpmv); - aa.SampleMovVert(tmpmv, ap.SampleNum, ap.SampleMode); + aa.sampleMovVert(tmpmv, ap.SampleNum, ap.SampleMode); aa.mov=&tmpmv; aa.fix=&Fix; @@ -98,7 +98,7 @@ void MeshTree::ProcessArc(int fixId, int movId, vcg::Matrix44d &MovM, vcg::Align vcg::Matrix44d In=MovM; // Perform the ICP algorithm - aa.Align(In,UG,VG,result); + aa.align(In,UG,VG,result); result.FixName=fixId; result.MovName=movId; @@ -185,7 +185,7 @@ void MeshTree::Process(vcg::AlignPair::Param &ap, MeshTree::Param &mtp) { ProcessArc(OG.SVA[i].s, OG.SVA[i].t, *curResult, ap); curResult->area= OG.SVA[i].norm_area; - if( curResult->IsValid() ) + if( curResult->isValid() ) { hasValidAlign = true; std::pair dd=curResult->computeAvgErr(); @@ -208,7 +208,7 @@ void MeshTree::Process(vcg::AlignPair::Param &ap, MeshTree::Param &mtp) vcg::Distribution H; // stat for printing for(QList::iterator li=resultList.begin();li!=resultList.end();++li) - if ((*li).IsValid()) + if ((*li).isValid()) H.Add(li->err); cb(0,qUtf8Printable(buf.sprintf("Completed Mesh-Mesh Alignment: Avg Err %5.3f; Median %5.3f; 90%% %5.3f\n", H.Avg(), H.Percentile(0.5f), H.Percentile(0.9f)))); @@ -246,7 +246,7 @@ void MeshTree::ProcessGlobal(vcg::AlignPair::Param &ap) vcg::AlignGlobal AG; std::vector ResVecPtr; for(QList::iterator li=resultList.begin();li!=resultList.end();++li) - if ((*li).IsValid()) + if ((*li).isValid()) ResVecPtr.push_back(&*li); AG.BuildGraph(ResVecPtr, GluedTrVec, GluedIdVec); diff --git a/vcglib b/vcglib index a374e959e..074a89c58 160000 --- a/vcglib +++ b/vcglib @@ -1 +1 @@ -Subproject commit a374e959ee6ffc8fc968b4c253fa3bf9c2b06130 +Subproject commit 074a89c5888561c36e022efe9cd7333ffba4ca9b From cead178da1f18fa3d458f3357873b63795bd7bf5 Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Fri, 29 May 2020 15:18:59 +0200 Subject: [PATCH 08/11] meshlab now uses alignpair from vcg --- src/meshlabplugins/edit_align/CMakeLists.txt | 2 - .../edit_align/align/AlignGlobal.cpp | 3 +- .../edit_align/align/AlignPair.cpp | 710 ------------------ .../edit_align/align/AlignPair.h | 31 - .../edit_align/align/OccupancyGrid.h | 2 +- src/meshlabplugins/edit_align/edit_align.h | 2 +- src/meshlabplugins/edit_align/edit_align.pro | 2 - src/meshlabplugins/edit_align/meshtree.cpp | 2 +- src/meshlabplugins/edit_align/meshtree.h | 2 +- src/meshlabplugins/filter_clean/align_tools.h | 2 +- .../filter_clean/cleanfilter.cpp | 2 +- vcglib | 2 +- 12 files changed, 8 insertions(+), 754 deletions(-) delete mode 100644 src/meshlabplugins/edit_align/align/AlignPair.cpp delete mode 100644 src/meshlabplugins/edit_align/align/AlignPair.h diff --git a/src/meshlabplugins/edit_align/CMakeLists.txt b/src/meshlabplugins/edit_align/CMakeLists.txt index d2e4e5625..a5e30c4a1 100644 --- a/src/meshlabplugins/edit_align/CMakeLists.txt +++ b/src/meshlabplugins/edit_align/CMakeLists.txt @@ -11,7 +11,6 @@ if(TARGET external-newuoa) set(SOURCES align/align_parameter.cpp align/AlignGlobal.cpp - align/AlignPair.cpp align/OccupancyGrid.cpp alignDialog.cpp AlignPairDialog.cpp @@ -27,7 +26,6 @@ if(TARGET external-newuoa) set(HEADERS align/align_parameter.h align/AlignGlobal.h - align/AlignPair.h align/OccupancyGrid.h alignDialog.h AlignPairDialog.h diff --git a/src/meshlabplugins/edit_align/align/AlignGlobal.cpp b/src/meshlabplugins/edit_align/align/AlignGlobal.cpp index d610cb830..3caf2e816 100644 --- a/src/meshlabplugins/edit_align/align/AlignGlobal.cpp +++ b/src/meshlabplugins/edit_align/align/AlignGlobal.cpp @@ -23,9 +23,8 @@ #include #include -#include "AlignPair.h" +#include #include "AlignGlobal.h" -#include "vcg/complex/algorithms/point_matching_scale.h" #include using namespace vcg; using namespace std; diff --git a/src/meshlabplugins/edit_align/align/AlignPair.cpp b/src/meshlabplugins/edit_align/align/AlignPair.cpp deleted file mode 100644 index c774a4082..000000000 --- a/src/meshlabplugins/edit_align/align/AlignPair.cpp +++ /dev/null @@ -1,710 +0,0 @@ -/**************************************************************************** -* VCGLib o o * -* Visual and Computer Graphics Library o o * -* _ O _ * -* Copyright(C) 2004 \/)\/ * -* Visual Computing Lab /\/| * -* ISTI - Italian National Research Council | * -* \ * -* All rights reserved. * -* * -* This program is free software; you can redistribute it and/or modify * -* it under the terms of the GNU General Public License as published by * -* the Free Software Foundation; either version 2 of the License, or * -* (at your option) any later version. * -* * -* This program is distributed in the hope that it will be useful, * -* but WITHOUT ANY WARRANTY; without even the implied warranty of * -* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * -* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * -* for more details. * -* * -****************************************************************************/ - -#include "AlignPair.h" -#include "vcg/complex/algorithms/point_matching_scale.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -using namespace vcg; -using namespace std; - -//bool AlignPair::A2Mesh::Import(const char *filename, Matrix44d &Tr) -//{ -// int err = tri::io::Importer::Open(*this, filename); -// if (err) { -// printf("Error in reading %s: '%s'\n", filename, tri::io::Importer::ErrorMsg(err)); -// exit(-1); -// } -// printf("read mesh `%s'\n", filename); -// return Init(Tr); -//} - -//bool AlignPair::A2Mesh::Init(const Matrix44d &Tr) -//{ -// Matrix44d Idn; Idn.SetIdentity(); -// tri::Clean::RemoveUnreferencedVertex(*this); -// if (Tr != Idn) tri::UpdatePosition::Matrix(*this, Tr); -// tri::UpdateNormal::PerVertexNormalizedPerFaceNormalized(*this); -// tri::UpdateFlags::FaceBorderFromNone(*this); -// tri::UpdateBounding::Box(*this); - -// return true; -//} - - -//void AlignPair::Stat::clear() -//{ -// I.clear(); -// StartTime = 0; -// MovVertNum = 0; -// FixVertNum = 0; -// FixFaceNum = 0; -//} - -// Restituisce true se nelle ultime iterazioni non e' diminuito -// l'errore -//bool AlignPair::Stat::Stable(int lastiter) -//{ -// if (I.empty()) return false; -// int parag = int(I.size()) - lastiter; - -// if (parag < 0) parag = 0; -// if (I.back().pcl50 < I[parag].pcl50) return false; // se siamo diminuiti non e' stabile - -// return true; - -//} - - -//void AlignPair::Stat::Dump(FILE *fp) -//{ -// if (I.size() == 0) { -// fprintf(fp, "Empty AlignPair::Stat\n"); -// return; -// } -// fprintf(fp, "Final Err %8.5f In %i iterations Total Time %ims\n", LastPcl50(), (int)I.size(), TotTime()); -// fprintf(fp, "Mindist Med Hi Avg RMS StdDev Time Tested Used Dist Bord Angl \n"); -// for (unsigned int qi = 0; qi < I.size(); ++qi) -// fprintf(fp, "%5.2f (%6.3f:%6.3f) (%6.3f %6.3f %6.3f) %4ims %5i %5i %4i+%4i+%4i\n", -// I[qi].MinDistAbs, -// I[qi].pcl50, I[qi].pclhi, -// I[qi].AVG, I[qi].RMS, I[qi].StdDev, -// IterTime(qi), -// I[qi].SampleTested, I[qi].SampleUsed, I[qi].DistanceDiscarded, I[qi].BorderDiscarded, I[qi].AngleDiscarded); -//} - -// Scrive una tabella con tutti i valori -//void AlignPair::Stat::HTMLDump(FILE *fp) -//{ -// fprintf(fp, "Final Err %8.5f In %i iterations Total Time %ims\n", LastPcl50(), (int)I.size(), TotTime()); -// fprintf(fp, "\n"); -// fprintf(fp, "\n", -// I[qi].MinDistAbs, -// I[qi].pcl50, I[qi].pclhi, -// I[qi].AVG, I[qi].RMS, I[qi].StdDev, -// IterTime(qi), -// I[qi].SampleTested, I[qi].SampleUsed, I[qi].DistanceDiscarded, I[qi].BorderDiscarded, I[qi].AngleDiscarded); -// fprintf(fp, "
Mindist 50ile Hi Avg RMS StdDev Time Tested Used Dist Bord Angl \n"); -// for (unsigned int qi = 0; qi < I.size(); ++qi) -// fprintf(fp, "
%8.5f %9.6f %8.5f %5.3f %8.5f %9.6f %4ims %5i %5i %4i %4i %4i
\n"); -//} - - - -/* -This function is used to choose remove outliers after each ICP iteration. -All the points with a distance over the given Percentile are discarded. -It uses two parameters -MaxPointNum an (unused) hard limit on the number of points that are chosen -MinPointNum the minimum number of points that have to be chosen to be usable - -*/ -//bool AlignPair::ChoosePoints(vector &Ps, // vertici corrispondenti su src (rossi) -// vector &Ns, // normali corrispondenti su src (rossi) -// vector &Pt, // vertici scelti su trg (verdi) -// vector &OPt, // vertici scelti su trg (verdi) -// double PassHi, -// Histogramf &H) -//{ -// const int N = ap.MaxPointNum; -// double newmaxd = H.Percentile(float(PassHi)); -// //printf("%5.1f of the pairs has a distance less than %g and greater than %g (0..%g) avg %g\n", Perc*100,newmind,newmaxd,H.maxv,H.Percentile(.5)); -// int sz = int(Ps.size()); -// int fnd = 0; -// int lastgood = sz - 1; -// math::SubtractiveRingRNG myrnd; -// while (fnd < N && fnd < lastgood) -// { -// int index = fnd + myrnd.generate(lastgood - fnd); -// double dd = Distance(Ps[index], Pt[index]); -// if (dd <= newmaxd) -// { -// swap(Ps[index], Ps[fnd]); -// swap(Ns[index], Ns[fnd]); -// swap(Pt[index], Pt[fnd]); -// swap(OPt[index], OPt[fnd]); -// ++fnd; -// } -// else -// { -// swap(Ps[index], Ps[lastgood]); -// swap(Ns[index], Ns[lastgood]); -// swap(Pt[index], Pt[lastgood]); -// swap(OPt[index], OPt[lastgood]); -// lastgood--; -// } -// } -// Ps.resize(fnd); -// Ns.resize(fnd); -// Pt.resize(fnd); -// OPt.resize(fnd); -//// printf("Scelte %i coppie tra le %i iniziali, scartate quelle con dist > %f\n", fnd, sz, newmaxd); - -// if ((int)Ps.size() < ap.MinPointNum) { -// printf("Troppi pochi punti!\n"); -// Ps.clear(); -// Ns.clear(); -// Pt.clear(); -// OPt.clear(); -// return false; -// } -// return true; -//} - -/* -Funzione chiamata dalla Align ad ogni ciclo -Riempie i vettori e con i coordinate e normali presi dal vettore di vertici mov -della mesh da muovere trasformata secondo la matrice -Calcola anche il nuovo bounding box di tali vertici trasformati. -*/ -bool AlignPair::InitMov( - vector< Point3d > &MovVert, - vector< Point3d > &MovNorm, - Box3d &trgbox, - const Matrix44d &in) // trasformazione Iniziale (che porta i punti di trg su src) -{ - Point3d pp, nn; - - MovVert.clear(); - MovNorm.clear(); - trgbox.SetNull(); - - A2Mesh::VertexIterator vi; - for (vi = mov->begin(); vi != mov->end(); vi++) { - pp = in*(*vi).P(); - nn = in*Point3d((*vi).P() + (*vi).N()) - pp; - nn.Normalize(); - MovVert.push_back(pp); - MovNorm.push_back(nn); - trgbox.Add(pp); - } - return true; -} - -bool AlignPair::InitFixVert(AlignPair::A2Mesh *fm, - AlignPair::Param &pp, - A2GridVert &u, - int PreferredGridSize) -{ - Box3d bb2 = fm->bbox; - double MinDist = pp.MinDistAbs*1.1; - //the bbox of the grid should be enflated of the mindist used in the ICP search - bb2.Offset(Point3d(MinDist, MinDist, MinDist)); - - u.SetBBox(bb2); - - //Inserisco la src nella griglia - if (PreferredGridSize == 0) PreferredGridSize = int(fm->vert.size())*pp.UGExpansionFactor; - u.Set(fm->vert.begin(), fm->vert.end());//, PreferredGridSize); - printf("UG %i %i %i\n", u.siz[0], u.siz[1], u.siz[2]); - return true; -} - - -bool AlignPair::InitFix(AlignPair::A2Mesh *fm, - AlignPair::Param &pp, - A2Grid &u, - int PreferredGridSize) -{ - tri::InitFaceIMark(*fm); - - Box3d bb2 = fm->bbox; - // double MinDist= fm->bbox.Diag()*pp.MinDistPerc; - double MinDist = pp.MinDistAbs*1.1; - //gonfio della distanza utente il BBox della seconda mesh - bb2.Offset(Point3d(MinDist, MinDist, MinDist)); - - u.SetBBox(bb2); - - //Inserisco la src nella griglia - if (PreferredGridSize == 0) PreferredGridSize = int(fm->face.size())*pp.UGExpansionFactor; - u.Set(fm->face.begin(), fm->face.end(), PreferredGridSize); - printf("UG %i %i %i\n", u.siz[0], u.siz[1], u.siz[2]); - return true; -} -/* -The Main ICP alignment Function: -It assumes that: -we have two meshes: -- Fix the mesh that does not move and stays in the spatial indexing structure. -- Mov the mesh we 'move' e.g. the one for which we search the transforamtion. - -requires normalize normals for vertices AND faces -Allinea due mesh; -Assume che: -la uniform grid sia gia' inizializzata con la mesh fix - - -*/ - -//bool AlignPair::align( -// A2Grid &u, -// A2GridVert &uv, -// const Matrix44d &in, // trasformazione Iniziale (che porta i punti di mov su fix) -// Matrix44d &out, // trasformazione calcolata -// vector &Pfix, // vertici corrispondenti su src (rossi) -// vector &Nfix, // normali corrispondenti su src (rossi) -// vector &OPmov, // vertici scelti su trg (verdi) prima della trasformazione in ingresso (Original Point Target) -// vector &ONmov, // normali scelti su trg (verdi) -// Histogramf &H, -// AlignPair::Stat &as) -//{ -// vector beyondCntVec; // vettore per marcare i movvert che sicuramente non si devono usare -// // ogni volta che un vertice si trova a distanza oltre max dist viene incrementato il suo contatore; -// // i movvert che sono stati scartati piu' di MaxCntDist volte non si guardano piu'; -// const int maxBeyondCnt = 3; -// vector< Point3d > movvert; -// vector< Point3d > movnorm; -// vector Pmov; // vertici scelti dopo la trasf iniziale -// status = SUCCESS; -// int tt0 = clock(); - -// out = in; - -// int i; - -// double CosAngleThr = cos(ap.MaxAngleRad); -// double StartMinDist = ap.MinDistAbs; -// int tt1 = clock(); -// int ttsearch = 0; -// int ttleast = 0; -// int nc = 0; -// as.clear(); -// as.StartTime = clock(); - -// beyondCntVec.resize(mov->size(), 0); - -// /**************** BEGIN ICP LOOP ****************/ -// do -// { -// Stat::IterInfo ii; -// Box3d movbox; -// InitMov(movvert, movnorm, movbox, out); -// H.SetRange(0.0f, float(StartMinDist), 512, 2.5f); -// Pfix.clear(); -// Nfix.clear(); -// Pmov.clear(); -// OPmov.clear(); -// ONmov.clear(); -// int tts0 = clock(); -// ii.MinDistAbs = StartMinDist; -// int LocSampleNum = min(ap.SampleNum, int(movvert.size())); -// Box3d fixbox; -// if (u.Empty()) fixbox = uv.bbox; -// else fixbox = u.bbox; -// for (i = 0; i < LocSampleNum; ++i) -// { -// if (beyondCntVec[i] < maxBeyondCnt) -// { - -// if (fixbox.IsIn(movvert[i])) -// { -// double error = StartMinDist; -// Point3d closestPoint, closestNormal; -// double maxd = StartMinDist; -// ii.SampleTested++; -// if (u.Empty()) // using the point cloud grid -// { -// A2Mesh::VertexPointer vp = tri::GetClosestVertex(*fix, uv, movvert[i], maxd, error); -// if (error >= StartMinDist) { -// ii.DistanceDiscarded++; ++beyondCntVec[i]; continue; -// } -// if (movnorm[i].dot(vp->N()) < CosAngleThr) { -// ii.AngleDiscarded++; continue; -// } -// closestPoint = vp->P(); -// closestNormal = vp->N(); -// } -// else // using the standard faces and grid -// { -// A2Mesh::FacePointer f = vcg::tri::GetClosestFaceBase(*fix, u, movvert[i], maxd, error, closestPoint); -// if (error >= StartMinDist) { -// ii.DistanceDiscarded++; ++beyondCntVec[i]; continue; -// } -// if (movnorm[i].dot(f->N()) < CosAngleThr) { -// ii.AngleDiscarded++; continue; -// } -// Point3d ip; -// InterpolationParameters(*f, f->N(), closestPoint, ip); -// const double IP_EPS = 0.00001; -// // If ip[i] == 0 it means that we are on the edge opposite to i -// if ((fabs(ip[0]) <= IP_EPS && f->IsB(1)) || (fabs(ip[1]) <= IP_EPS && f->IsB(2)) || (fabs(ip[2]) <= IP_EPS && f->IsB(0))){ -// ii.BorderDiscarded++; continue; -// } -// closestNormal = f->N(); -// } -// // The sample was accepted. Store it. -// Pmov.push_back(movvert[i]); -// OPmov.push_back((*mov)[i].P()); -// ONmov.push_back((*mov)[i].N()); -// Nfix.push_back(closestNormal); -// Pfix.push_back(closestPoint); -// H.Add(float(error)); -// ii.SampleUsed++; -// } -// else -// beyondCntVec[i] = maxBeyondCnt + 1; -// } -// } // End for each pmov -// int tts1 = clock(); -// //printf("Found %d pairs\n",(int)Pfix.size()); -// if (!choosePoints(Pfix, Nfix, Pmov, OPmov, ap.PassHiFilter, H)) -// { -// if (int(Pfix.size()) < ap.MinPointNum) -// { -// status = TOO_FEW_POINTS; -// ii.Time = clock(); -// as.I.push_back(ii); -// return false; -// } -// } -// Matrix44d newout; -// switch (ap.MatchMode) -// { -// case AlignPair::Param::MMSimilarity: vcg::PointMatchingScale::computeRotoTranslationScalingMatchMatrix(newout, Pfix, Pmov); break; -// case AlignPair::Param::MMRigid: ComputeRigidMatchMatrix(Pfix, Pmov, newout); break; -// default: -// status = UNKNOWN_MODE; -// ii.Time = clock(); -// as.I.push_back(ii); -// return false; -// } - -// // double sum_before=0; -// // double sum_after=0; -// // for(unsigned int iii=0;iii %f\n",sum_before/double(Pfix.size()),sum_after/double(Pfix.size()) ) ; - -// // le passate successive utilizzano quindi come trasformazione iniziale questa appena trovata. -// // Nei prossimi cicli si parte da questa matrice come iniziale. -// out = newout * out; - -// assert(Pfix.size() == Pmov.size()); -// int tts2 = clock(); -// ttsearch += tts1 - tts0; -// ttleast += tts2 - tts1; -// ii.pcl50 = H.Percentile(.5); -// ii.pclhi = H.Percentile(float(ap.PassHiFilter)); -// ii.AVG = H.Avg(); -// ii.RMS = H.RMS(); -// ii.StdDev = H.StandardDeviation(); -// ii.Time = clock(); -// as.I.push_back(ii); -// nc++; -// // The distance of the next points to be considered is lowered according to the parameter. -// // We use 5 times the percentile of the found points. -// if (ap.ReduceFactorPerc<1) StartMinDist = max(ap.MinDistAbs*ap.MinMinDistPerc, min(StartMinDist, 5.0*H.Percentile(float(ap.ReduceFactorPerc)))); -// } while ( -// nc <= ap.MaxIterNum && -// H.Percentile(.5) > ap.TrgDistAbs && -// (ncap.MaxScale || math::Abs(1 - scv[1]) > ap.MaxScale || math::Abs(1 - scv[2]) > ap.MaxScale)) { -// status = TOO_MUCH_SCALE; -// return false; -// } -// if (shv[0] > ap.MaxShear || shv[1] > ap.MaxShear || shv[2] > ap.MaxShear) { -// status = TOO_MUCH_SHEAR; -// return false; -// } -// printf("Grid %i %i %i - fn %i\n", u.siz[0], u.siz[1], u.siz[2], fix->fn); -// printf("Init %8.3f Loop %8.3f Search %8.3f least sqrt %8.3f\n", -// float(tt1 - tt0) / CLOCKS_PER_SEC, float(tt2 - tt1) / CLOCKS_PER_SEC, -// float(ttsearch) / CLOCKS_PER_SEC, float(ttleast) / CLOCKS_PER_SEC); - -// return true; -//} - - - - -//const char *AlignPair::ErrorMsg(ErrorCode code) -//{ -// switch (code){ -// case SUCCESS: return "Success "; -// case NO_COMMON_BBOX: return "No Common BBox "; -// case TOO_FEW_POINTS: return "Too few points "; -// case LSQ_DIVERGE: return "LSQ not converge"; -// case TOO_MUCH_SHEAR: return "Too much shear "; -// case TOO_MUCH_SCALE: return "Too much scale "; -// case UNKNOWN_MODE: return "Unknown mode "; -// default: assert(0); return "Catastrophic Error"; -// } -// return 0; -//} -/* - -Questa parte era relativa all'allineatore automatico. -Da controllare se ancora funge. - - -// Calcola la migliore traslazione possibile, -// cioe' quella dove i punti in movvert cascano per la maggior parte -// in celle della ug occupate dalla fixmesh -// Restituisce il numero di celle di sovrapposizione massimo trovato - -int AlignPair::SearchTranslate(A2Grid &u, -const Matrix44d &BaseRot, -const int range, -const int step, -Point3d &BestTransV, // miglior vettore di spostamento -bool Verbose) -{ -const int wide1=(range*2+1); -const int wide2=wide1*wide1; -vector< Point3d > movvert; -vector< Point3d > movnorm; -Box3d movbox; - -int t0=clock(); -InitMov(movvert,movnorm,movbox,BaseRot); -Point3i ip; -int i,ii,jj,kk; -vector test((range*2+1)*(range*2+1)*(range*2+1),0); -//int bx,by,bz,ex,ey,ez; - -Point3i b,e; -int testposii,testposjj; -for(i=0;i=u.siz[0]) e[0]-=step; -while(b[1]<0) b[1]+=step; -while(e[1]>=u.siz[1]) e[1]-=step; -while(b[2]<0) b[2]+=step; -while(e[2]>=u.siz[2]) e[2]-=step; - -for(ii=b[0];ii<=e[0];ii+=step) -{ -testposii=(ii-ip[0]+range)*wide2; -for(jj=b[1];jj<=e[1];jj+=step) -{ -testposjj=testposii+(jj-ip[1]+range)*wide1-ip[2]+range; -for(kk=b[2];kk<=e[2];kk+=step) -{ -UGrid< A2Mesh::face_container >::UG ** g = u.Grid(ii,jj,kk); -//if((*(g+1))- (*g) >0) -if((*(g+1))!=(*g) ) -++test[testposjj+kk]; -assert(ii >=0 && ii < u.siz[0]); -assert(jj >=0 && jj < u.siz[1]); -assert(kk >=0 && kk < u.siz[2]); -} -} -} -} -} -int maxfnd=0; -Point3i BestI; -for(ii=-range;ii<=range;ii+=step) -for(jj=-range;jj<=range;jj+=step) -for(kk=-range;kk<=range;kk+=step) -{const int pos =(range+ii)*wide2+(range+jj)*wide1+range+kk; -if(test[pos]>maxfnd) -{ -BestI=Point3i(ii,jj,kk); -BestTransV=Point3d(ii*u.voxel[0], jj*u.voxel[1], kk*u.voxel[2]); -maxfnd=test[pos]; -} -//printf("Found %i su %i in %i\n",test[testcnt],movvert.size(),t1-t0); -} - -int t1=clock(); -if(Verbose) printf("BestTransl (%4i in %3ims) is %8.4f %8.4f %8.4f (%3i %3i %3i)\n",maxfnd,t1-t0,BestTransV[0],BestTransV[1],BestTransV[2],BestI[0],BestI[1],BestI[2]); -return maxfnd; -} -#if 0 -int AlignPair::SearchTranslate(UGrid< A2Mesh::face_container > &u, -vector< Point3d > &movvert, -vector< Point3d > &movnorm, -Box3d &movbox, -const Matrix44d &BaseRot, -const int range, -const int step, -Point3d &BestTransV, // miglior vettore di spostamento -bool Verbose) -{ -const int wide1=(range*2+1); -const int wide2=wide1*wide1; - -int t0=clock(); -InitMov(movvert,movnorm,movbox,BaseRot); -Point3i ip; -int i,ii,jj,kk; -vector test((range*2+1)*(range*2+1)*(range*2+1),0); -for(i=0;i=0 && ip[0]+ii=0 && ip[1]+jj=0 && ip[2]+kk::UG ** g = u.Grid(ip[0]+ii,ip[1]+jj,ip[2]+kk); -if((*(g+1))- (*g) >0) -++test[(range+ii)*wide2+(range+jj)*wide1+range+kk]; -} -} -} -int maxfnd=0; -Point3i BestI; -for(ii=-range;ii<=range;ii+=step) -for(jj=-range;jj<=range;jj+=step) -for(kk=-range;kk<=range;kk+=step) -{const int pos =(range+ii)*wide2+(range+jj)*wide1+range+kk; -assert(test2[pos]==test[pos]); -if(test[pos]>maxfnd) -{ -BestI=Point3i(ii,jj,kk); -BestTransV=Point3d(ii*u.voxel[0], jj*u.voxel[1], kk*u.voxel[2]); -maxfnd=test[pos]; -} -//printf("Found %i su %i in %i\n",test[testcnt],movvert.size(),t1-t0); -} - -int t1=clock(); -if(Verbose) printf("BestTransl (%4i in %3ims) is %8.4f %8.4f %8.4f (%3i %3i %3i)\n",maxfnd,t1-t0,BestTransV[0],BestTransV[1],BestTransV[2],BestI[0],BestI[1],BestI[2]); -return maxfnd; -} -*/ - - -/**********************************************************/ -// Funzioni per la scelta dei vertici sulla mesh da muovere - -//bool AlignPair::SampleMovVert(vector &vert, int SampleNum, AlignPair::Param::SampleModeEnum SampleMode) -//{ -// switch (SampleMode) -// { -// case AlignPair::Param::SMRandom: return SampleMovVertRandom(vert, SampleNum); -// case AlignPair::Param::SMNormalEqualized: return SampleMovVertNormalEqualized(vert, SampleNum); -// default: assert(0); -// } -// return false; -//} - - -// Scelta a caso semplice -//bool AlignPair::SampleMovVertRandom(vector &vert, int SampleNum) -//{ -// if (int(vert.size()) <= SampleNum) return true; -// int i; -// for (i = 0; i < SampleNum; ++i) -// { -// int pos = myrnd.generate(vert.size()); -// assert(pos >= 0 && pos < int(vert.size())); -// swap(vert[i], vert[pos]); -// } -// vert.resize(SampleNum); -// return true; -//} - -/* -Scelta a caso in maniera tale che la distribuzione delle normali dei -punti scelti sia il piu' possibile uniforme. In questo modo anche piccole -parti inclinate vengono sicuramente campionate - -Si precalcola un piccolo (42) insieme di normali e si fa bucketing di -tutti i vertici della mesh su di esse. -Poi si scelgono punti scegliendo ogni volta prima un bucket -e poi un punto all'interno del bucket - -*/ - -//bool AlignPair::SampleMovVertNormalEqualized(vector &vert, int SampleNum) -//{ -// // assert(0); - -// // int t0=clock(); -// vector NV; -// if (NV.size() == 0) -// { -// GenNormal::Fibonacci(30, NV); -// printf("Generated %i normals\n", int(NV.size())); -// } -// // Bucket vector dove, per ogni normale metto gli indici -// // dei vertici ad essa corrispondenti -// vector > BKT(NV.size()); -// for (size_t i = 0; i < vert.size(); ++i) -// { -// int ind = GenNormal::BestMatchingNormal(vert[i].N(), NV); -// BKT[ind].push_back(int(i)); -// } -// //int t1=clock(); - -// // vettore di contatori per sapere quanti punti ho gia' preso per ogni bucket -// vector BKTpos(BKT.size(), 0); - -// if (SampleNum >= int(vert.size())) SampleNum = vert.size() - 1; - -// for (int i = 0; i < SampleNum;) -// { -// int ind = myrnd.generate(BKT.size()); // Scelgo un Bucket -// int &CURpos = BKTpos[ind]; -// vector &CUR = BKT[ind]; - -// if (CURpos - - - -#endif diff --git a/src/meshlabplugins/edit_align/align/OccupancyGrid.h b/src/meshlabplugins/edit_align/align/OccupancyGrid.h index 838944982..2cf8b918c 100644 --- a/src/meshlabplugins/edit_align/align/OccupancyGrid.h +++ b/src/meshlabplugins/edit_align/align/OccupancyGrid.h @@ -23,7 +23,7 @@ #ifndef ALIGN_OCCUPANCY_GRID_H #define ALIGN_OCCUPANCY_GRID_H -#include "AlignPair.h" +#include #include #include diff --git a/src/meshlabplugins/edit_align/edit_align.h b/src/meshlabplugins/edit_align/edit_align.h index 30c854887..04e72d4ee 100644 --- a/src/meshlabplugins/edit_align/edit_align.h +++ b/src/meshlabplugins/edit_align/edit_align.h @@ -25,7 +25,7 @@ #define EditAlignPLUGIN_H #include -#include "align/AlignPair.h" +#include #include "align/OccupancyGrid.h" #include "meshtree.h" #include diff --git a/src/meshlabplugins/edit_align/edit_align.pro b/src/meshlabplugins/edit_align/edit_align.pro index 92b23235d..111852a87 100644 --- a/src/meshlabplugins/edit_align/edit_align.pro +++ b/src/meshlabplugins/edit_align/edit_align.pro @@ -10,7 +10,6 @@ HEADERS += \ alignDialog.h \ AlignPairDialog.h \ AlignPairWidget.h \ - align/AlignPair.h \ align/AlignGlobal.h \ align/OccupancyGrid.h \ align/align_parameter.h \ @@ -26,7 +25,6 @@ SOURCES += \ alignDialog.cpp \ AlignPairWidget.cpp \ AlignPairDialog.cpp \ - align/AlignPair.cpp \ align/AlignGlobal.cpp \ align/OccupancyGrid.cpp \ align/align_parameter.cpp \ diff --git a/src/meshlabplugins/edit_align/meshtree.cpp b/src/meshlabplugins/edit_align/meshtree.cpp index 138bb27dc..6b0a5aea0 100644 --- a/src/meshlabplugins/edit_align/meshtree.cpp +++ b/src/meshlabplugins/edit_align/meshtree.cpp @@ -81,7 +81,7 @@ void MeshTree::ProcessArc(int fixId, int movId, vcg::Matrix44d &MovM, vcg::Align else { Fix.init(vcg::Matrix44d::Identity()); - vcg::AlignPair::InitFix(&Fix, ap, UG); + vcg::AlignPair::initFix(&Fix, ap, UG); } // 2) Convert the second mesh and sample a points on it. diff --git a/src/meshlabplugins/edit_align/meshtree.h b/src/meshlabplugins/edit_align/meshtree.h index 6aefdacff..856123b90 100644 --- a/src/meshlabplugins/edit_align/meshtree.h +++ b/src/meshlabplugins/edit_align/meshtree.h @@ -27,7 +27,7 @@ #include #include -#include "align/AlignPair.h" +#include #include "align/AlignGlobal.h" #include "align/OccupancyGrid.h" #include diff --git a/src/meshlabplugins/filter_clean/align_tools.h b/src/meshlabplugins/filter_clean/align_tools.h index 5946d4a92..1099f27c1 100644 --- a/src/meshlabplugins/filter_clean/align_tools.h +++ b/src/meshlabplugins/filter_clean/align_tools.h @@ -37,7 +37,7 @@ #include #include -#include +#include class AlignTools : public QObject { diff --git a/src/meshlabplugins/filter_clean/cleanfilter.cpp b/src/meshlabplugins/filter_clean/cleanfilter.cpp index 954068646..a3b3d599a 100644 --- a/src/meshlabplugins/filter_clean/cleanfilter.cpp +++ b/src/meshlabplugins/filter_clean/cleanfilter.cpp @@ -48,7 +48,7 @@ CleanFilter::CleanFilter() << FP_REMOVE_DUPLICATE_FACE << FP_REMOVE_FOLD_FACE << FP_REMOVE_NON_MANIF_EDGE - << FP_REMOVE_NON_MANIF_EDGE_SPLIT + << FP_REMOVE_NON_MANIF_EDGE_SPLIT << FP_REMOVE_NON_MANIF_VERT << FP_REMOVE_UNREFERENCED_VERTEX << FP_REMOVE_DUPLICATED_VERTEX diff --git a/vcglib b/vcglib index 074a89c58..64e352374 160000 --- a/vcglib +++ b/vcglib @@ -1 +1 @@ -Subproject commit 074a89c5888561c36e022efe9cd7333ffba4ca9b +Subproject commit 64e352374a854f68379bff7ce84775774d64a80f From d8d2951d86d1c8eb0b3d4412149371066cfc8570 Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Fri, 29 May 2020 15:21:54 +0200 Subject: [PATCH 09/11] fix compile error assert --- vcglib | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vcglib b/vcglib index 64e352374..4d57dde10 160000 --- a/vcglib +++ b/vcglib @@ -1 +1 @@ -Subproject commit 64e352374a854f68379bff7ce84775774d64a80f +Subproject commit 4d57dde1020038d780d7202d814ff093ee6caea1 From 34d21898251af92c28ed462e04b1384e7329ce1d Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Mon, 1 Jun 2020 10:13:51 +0200 Subject: [PATCH 10/11] update vcg --- vcglib | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vcglib b/vcglib index 4d57dde10..5b6d6ec76 160000 --- a/vcglib +++ b/vcglib @@ -1 +1 @@ -Subproject commit 4d57dde1020038d780d7202d814ff093ee6caea1 +Subproject commit 5b6d6ec767e9197705a50249d96f7e0715df93ae From 131297d2931c4ee5d2c93326a91a6b33f50fe84a Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Mon, 1 Jun 2020 16:50:30 +0200 Subject: [PATCH 11/11] update vcg --- vcglib | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vcglib b/vcglib index 5b6d6ec76..da737ebad 160000 --- a/vcglib +++ b/vcglib @@ -1 +1 @@ -Subproject commit 5b6d6ec767e9197705a50249d96f7e0715df93ae +Subproject commit da737ebad7f1396bf5ec380b9ef20551dcee950c