Geometry & Recognition

In computational geometry, an alpha shape, or α-shape, is a family of piecewise linear simple curves in the Euclidean plane associated with the shape of a finite set of points. https://en.wikipedia.org/wiki/Alpha_shape

class AlphaShape {
    double alpha;
    std::set<std::pair<int,int>> Border; //pair;
public:
    AlphaShape() { }
    AlphaShape(CDC *pDC, std::vector<CfPt>& Q, const double radius = 50) {
        if (Q.size() < 2) return;
        alpha = radius;
        const double alphaSq = alpha * alpha;
        const double diameter = 2 * alpha;
        for (int i = 0; i < Q.size(); ++i) {
            const CfPt &p = Q[i];
            for (int j = 0; j < i; ++j) {
                const CfPt &q = Q[j];
                if (q == p) continue;
                const double edist = p.Dist(q); 
                // 두 점의 사잇거리가 지름보다 크면 edge가 안됨;
                if (edist > diameter) continue;                     

                // 두 점을 원주에 포함하는 반지름 alpha인 두 원을 찾음;
                // note that c1 == c2 if edist == diameter;
                double scale = sqrt(alphaSq - (edist/2)*(edist/2));
                CfPt n = (q - p) * scale / edist;
                CfPt mid = (q + p) / 2;
                CfPt c1 = mid + n.Rot90(); // CfPt(mid.x - ny, mid.y + nx);
                CfPt c2 = mid - n.Rot90(); // CfPt(mid.x + ny, mid.y - nx);

                // alpha shape의 변이 되기 위해서는 두 원 중 하나는
                // 다른 점들을 포함하지 않아야 함;
                bool c1_empty = true, c2_empty = true;
                for (int k = Q.size(); (k-->0) && (c1_empty || c2_empty);) {
                    if (Q[k] == p|| Q[k]== q) continue; //i==k || j==k?
                    if (c1.DistSq(Q[k]) < alphaSq) c1_empty = false;
                    if (c2.DistSq(Q[k]) < alphaSq) c2_empty = false;            
                }

                if (c1_empty || c2_empty) {                  
                    if (c1_empty) drawCircle(pDC, c1);
                    if (c2_empty) drawCircle(pDC, c2);
                    // draw edge;
                    drawEdge(pDC, p, q);
                    //Border.insert(std::make_pair(i,j));
                }
            }
        }
    } 
    void drawEdge(CDC *pDC, const CfPt &p, const CfPt &q) {
        CPen pen(PS_SOLID,2,RGB(255,0,255)); //magenta;
        CPen *pOld = pDC->SelectObject(&pen);
        pDC->MoveTo(p); pDC->LineTo(q);
        pDC->SelectObject(pOld);
    }
    void drawCircle(CDC *pDC, const CfPt &c) {
        CPen pen(PS_SOLID, 1, RGB(0,200,200));
        CPen *pOld = pDC->SelectObject(&pen);
        pDC->SelectObject(GetStockObject(NULL_BRUSH));
        pDC->Ellipse(c.x-alpha, c.y-alpha, c.x+alpha, c.y+alpha);
        pDC->SelectObject(pOld);		
    }
    std::set<std::pair<int,int>> getBorder() {
        return Border;
    }
};

주어진 관측값 $\{x_i, y_i\}$을 spline $y=S(x)$을 이용해서 모델링을 하자. 주어진 관측위치 $x_i$에서 spline은 관측값에 충분히 가까이 지나야 하고, 더블어 관측과 관측 사이에서는 너무 급격히 변하지 않는 곡선이어야 한다. 곡선의 변동은 곡률에 의해서 결정되므로 급격히 변하지 않기 위해서는 spline의 이차 미분값의 크기가 되도록이면 작아한다. 따라서 구하려고 하는 spline은 다음을 최소화시키는 조건에서 찾을 수 있을 것이다.

$$L= \alpha \sum_i \left( \frac{y_i - S(x_i)}{\sigma_i} \right)^2 + \beta \int_{x_0}^{x_n} dx ( S''(x) )^2 \\  S=\underset{}{\text{argmin}}~L$$

여기서 $\sigma_i$는 관측오차에 해당한다. $L$의 두 항은 서로 반대되는 경향으로 spline을 변화시키는데, 첫번째 항이 기여가 큰 경우는 spline이 주어진 관측값의 보간 곡선에 가까워지므로 곡선이 급격히 변할 수 있다. 두번째 항의 기여가 커지는 경우는 spline은 주어진 관측값을 통과하지는 못하지만 충분히 부드러워지게 만든다. 그리고 극단적으로 이 값만 기여하는 경우는 spline은 직선이 되게 된다. 두 항의 상대적인 가중치를 설정하기 위해서 $\beta = 1-\alpha$로 놓을 수 있지만, 첫번째 항과 두 번째 항의 물리적인 차원이 다르므로 단순 비교는 정확하지 않다. 첫번쨰 항은 차원이 없는 반면에 두 번째 항은 거리의 역수 차원을 가지므로 비교를 하기 위해서는 적당한 거리 척도를 곱해  $(1-\alpha)\times\text{distance scale}$로 하는 것이 맞다. 보통 $\alpha=1$로 잡고, $\beta$를 $0$(interpolation에 가까워짐)에서부터 점점 큰 값으로 (직선에 가까워짐)으로 변화시키면 된다.

측정값 사이 구간에서 $n-1$개의 삼차 spline

$$ S(x)=\{ S_i (x) = a_ix^3 + b_i x^2+ c_ix +d_i ,~x_i  \le x \le x_{i+1} | i=0,1,2,...,n-1\} $$

을 이용해서 smoothing spline을 구하자. 이 경우 한 구간에서 spline의 미분은 인접 spline 2차 계수의 선형보간으로 표현할 수 있다.

$$S_i''(x)=   (1-x/h_i) b_i + (x/h_i) b_{i+1}, ~~h_i = x_{i+1}-x_i$$

따라서 $L$은 다음과 같이 표현된다.

$$L = \alpha \sum_{i=0}^n \left(\frac{y_i - d_i}{\sigma_i} \right)^2 + \beta \sum_{i=0}^{n-1} \frac{4h_i}{3} (b_i^2 +b_i b_{i+1} +b_{i+1}^2 )  $$

따라서 smoothing spline을 찾는 문제는 $L$을 최소화 시키는 $b_i, d_i$을 찾는 문제로 환원되었다. 나머지 계수 $c_i$는 spline의 연속성, $a_i$는 이차 도함수의 연속성을 요구하면 구할 수 있다.

$$ a_i = \frac{b_{i+1} - b_i}{3h_i} \\ c_i = \frac{d_{i+1}-d_i}{h_i} -\frac{1}{3} (b_{i+1}-2 b_i) h_i $$

여기에 1차 도함수의 연속성을 요구하면 $b_i$와 $d_i$사이의 recursion relation을 구할 수 있으므로 최종적으로 $b_i$ 또는 $d_i$을 구하는 문제로 환원이 된다.

$$ b_{i-1} h_{i-1} + 2 b_i ( h_{i-1} + h_i  ) + b_{i+1} h_i = \frac{3}{h_i} ( d_{i+1} - d_i ) - \frac{3}{h_{i-1}}(d_i - d_{i-1}) $$

계속...

Ref: https://monoceros.physics.muni.cz/~jancely/NM/Texty/Numerika/CubicSmoothingSpline.pdf

struct Spline {
    double x, y;
    double a, b, c, d;
    double eval(double t) {
        return d + (c + (b + a*t)*t)*t;
    }
    void toBezier(const double h, CPoint Q[4]) {
        double xp[4], yp[4];
        toBezier(h, xp, yp);
        for (int i = 0; i < 4; i++) 
            Q[i] = CPoint(int(xp[i]+.5), int(yp[i]+.5));
    }
    void toBezier(const double h, double xp[4], double yp[4]) {
        xp[0] = x; 
        xp[1] = x + h/3;
        xp[2] = x + 2*h/3;
        xp[3] = x + h;
        yp[0] = d;
        yp[1] = d + c * h / 3;
        yp[2] = d + (2 * c + b * h) * h / 3;
        yp[3] = d + (c + (b + a * h) * h) * h; 	
    }
};
void solveFiveDiag(std::vector<double> &u,
                   std::vector<double> &v,
                   std::vector<double> &w,
                   std::vector<double> &q) 
{
    const int n = u.size();
    u[0] = 0; //u[1]은 외부에서;
    v[1] = v[1]/u[1];
    w[1] = w[1]/u[1];
    for (int j = 2; j < n; j++) {
        u[j] = u[j] - u[j-2] * w[j-2] * w[j-2] - u[j-1] * v[j-1] * v[j-1];
        v[j] = (v[j] - u[j-1]*v[j-1]*w[j-1])/u[j];
        w[j] = w[j] / u[j];
    }
    //forward substitution
    q[1] = q[1] - v[0] * q[0];
    for (int j = 2; j < n; j++)
        q[j] = q[j] - v[j-1] * q[j-1] - w[j-2] * q[j-2];

    for (int j = 1; j < n; j++) 
        q[j] = q[j] / u[j];

    // back substitution
    q[n] = 0;
    // q[n-1]은 외부에서 정해짐;
    for (int j = n-2; j > 0; j--) 
        q[j] = q[j] - v[j] * q[j+1] - w[j]*q[j+2];
}
std::vector<Spline> SmoothingSpline1(
    std::vector<double> &x, //n+1
    std::vector<double> &y, //n+1
    std::vector<double> &sigma, //n+1
    double alpha, double beta)
{
    const int n = x.size()-1;
    std::vector<Spline> S(x.size());
    for (int i = x.size(); i-->0;)
        S[i].x = x[i], S[i].y = y[i];

    std::vector<double> h(n), r(x.size()), q(x.size()), u(n), v(n), w(n);
    h[0] = x[1] - x[0];
    r[0] = 3 / h[0];
    r[n] = 0;
    for (int i = 1; i < n; i++) {
        h[i] = x[i+1] - x[i];
        r[i] = 3/h[i];
        q[i] = (y[i+1] - y[i])*r[i] - (y[i] - y[i-1])*r[i-1];
    }
    //
    double mu = 2 * beta / (3 * alpha);
    for (int i = 1; i < n; i++) {
         double ui = SQR(r[i-1]) * sigma[i-1]
               + SQR(r[i-1] + r[i]) * sigma[i] + SQR(r[i]) * sigma[i+1];
         u[i] = mu * ui + 2 * (x[i+1] - x[i-1]);
         double vi = -(r[i-1] + r[i]) * r[i] * sigma[i] 
                - r[i] * (r[i] + r[i+1]) * sigma[i+1];
         v[i] = mu * vi + h[i];
         w[i] = mu * r[i] * r[i+1] * sigma[i+1];
    }
    solveFiveDiag(u, v, w, q);
    S[0].d = S[0].y - mu * r[0] * q[1] * sigma[0];
    double S1d = S[1].y - mu * ((-r[0] - r[1]) * q[1] + r[1] * q[2]) * sigma[1];
    S[0].a = q[1] / (3 * h[0]);
    S[0].b = 0;
    S[0].c = (S1d - S[0].d) / h[0] - q[1] * h[0] / 3;
    
    for (int j = 1; j < n; j++) {
        S[j].a = (q[j+1] - q[j]) / (3 * h[j]);
        S[j].b = q[j];
        S[j].c = (q[j] + q[j-1]) * h[j-1] + S[j-1].c;
        S[j].d = r[j-1] * q[j-1] - (r[j-1] + r[j]) * q[j] + r[j] * q[j+1];
        S[j].d = y[j] - mu * S[j].d * sigma[j];
    }
    return S;
}

주어진 점집합의 Convex hull을 찾는 알고리즘은 보통 $n \log(n)$ 정도의 time complexity를 가지고 있다. 그러나 점집합의 사이즈가 매우 커지는 경우에는 degeneracy로 인해 optimal한 알고리즘을 사용할 때 매우 주의를 기울려야 한다. 영상처리에서 convex hull을 사용할 때 근사적인 알고리즘만 있어도 충분한 경우가 많이 있는데 이런 경우에 이용할 수 있도록 알고리즘을 수정하자. Convex hull을 만들기 위해서는 점집한을 일정한 순서로 정렬을 하는 과정이 있는데 이 과정이 알고리즘의 가장 복잡한 부분이기도 하다. 이 정렬 과정을 근사적으로 수행하기 위해서 점집합을 수평(또는 수직)으로 일정한 간격으로 나누어진 구간으로 분류를 하자. 이 과정은 주어진 점집합을 한 번 순환만 하면 완료된다. 그런 다음 각 구간에서 수직방향으로(처음 수직방향으로 구간을 나누었으면 수평방향) 가장 바깥에 있는 두 점을 그 구간의 대표점으로 잡아 convex hull을 구하면 원 집합의 근사적인 convex hull이 된다. 이 구간 대표점들은 이미 정렬이 되어 있으므로 대표점 수(~구간 수)에 비례하는 시간 복잡도 이내에 알고리즘이 수행된다. 그리고 구간의 수를 점점 늘릴수록 정확한 convex hull에 근접하게 된다.  아래 그림은 500,000개 점의 근사적인 convex hull을 보여준다(bins=128)

// cross(P1-P0, P2-P0);
// = := colinear;
// < := <p0,p1,p2> CW order;
// > := <p0,p1,p2> CCW order;
static
int ccw(const CPoint &P0, const CPoint &P1, const CPoint &P2 ) {
    return (P1.x - P0.x)*(P2.y - P0.y) - (P2.x - P0.x)*(P1.y - P0.y);
}
std::vector<CPoint> approximateHull(const std::vector<CPoint>& Q, int bins) {
    if (Q.size() < 3) return std::vector<CPoint> ();
    if (bins < 0) bins = max(1, Q.size()/10);
    int Lmin = 0, Lmax = 0, Rmin = 0, Rmax = 0;
    int left = Q[0].x, right = left;
    for (int i = 1; i < Q.size(); i++) {
        if (Q[i].x <= left) {
            if (Q[i].x < left) {
                Lmin = Lmax = i;
                left = Q[i].x;
            } 
            else 
                if (Q[i].y < Q[Lmin].y) Lmin = i;
                else if (Q[i].y > Q[Lmax].y) Lmax = i;
        }
        if (Q[i].x >= right) {
            if (Q[i].x > right) {
                Rmin = Rmax = i;
                right = Q[i].x;
            }
            else 
                if (Q[i].y < Q[Rmin].y) Rmin = i;
                else if (Q[i].y > Q[Rmax].y) Rmax = i;
        }
    }
    if (left == right) return std::vector<CPoint> (); //vertical line;
    //
    std::vector<int> slotYhigh(bins + 2, -1);
    std::vector<int> slotYlow(bins + 2, -1);
    slotYlow.front() = Lmin; 
    slotYhigh.front() = Lmax;
    slotYlow.back()  = Rmin;
    slotYhigh.back()  = Rmax;
    //
    const CPoint &A = Q[Lmin];
    const CPoint &B = Q[Rmin];
    const CPoint &C = Q[Lmax];
    const CPoint &D = Q[Rmax];
    for (int i = 0; i < Q.size(); i++) {
        if (Q[i].x == right || Q[i].x == left) continue;
        if (ccw(A, B, Q[i]) < 0) { // below line(A,B);
            int b = 1 + (bins * (Q[i].x - left)) / (right - left);
            if (slotYlow[b]==-1) slotYlow[b] = i;
            else if (Q[i].y < Q[slotYlow[b]].y) slotYlow[b] = i;
            continue ;
        }
        if (ccw(C, D, Q[i]) > 0) {// above line(C,D);
            int b = 1 + (bins * (Q[i].x - left)) / (right - left);
            if (slotYhigh[b]==-1) slotYhigh[b] = i;
            else if (Q[i].y > Q[slotYhigh[b]].y) slotYhigh[b] = i;
            continue;
        }
    }
    std::vector<CPoint> hull(Q.size());
    int s = -1;
    for (int i = 0; i <= bins + 1; i++) {
        if (slotYlow[i]==-1) continue;
        while (s > 0) 
            if (ccw(hull[s-1], hull[s], Q[slotYlow[i]]) > 0) break;
            else s--;
        
        hull[++s] = Q[slotYlow[i]];
    };
    if (Rmax != Rmin) hull[++s] = Q[Rmax];
    //
    int s0 = s ;
    for (int i = bins; i >= 0; i--) {
        if (slotYhigh[i]==-1) continue;
        while (s > s0) 
            if (ccw(hull[s-1], hull[s], Q[slotYhigh[i]]) > 0) break;
            else s--;
        
        hull[++s] = Q[slotYhigh[i]];
    }
    if (hull[s] == hull[0]) s--;
    hull.resize(s+1);
    return hull;
}

평면상에  주어진 점들을 보간하는 함수로써 Catmull-Rom spline을 사용할 때는 매개변수를 써서 표현하는 것보다는 $x$-좌표를 사용해서 표현하는 것이 더 편리할 수 있다. 이 경우 control point의 $x$-좌표 간격이 균일하지 않으므로 이를 고려해야 한다. 

• 주어진 구간에서 spline은 $x$의 삼차함수로  써진다:
  $$S_i (x) = a + b(x-x_i) + c(x-x_i)^2 + d(x-x_i)^3,~ x_i \le  x\le x_{i+1}$$
• 곡선은 모든 control point를 통과해야 한다:
  $$S_i (x_i) = y_i,  ~S_i(x_{i+1}) = y_{i+1}$$
• 1차 미분값은 양 이웃점의 기울기로 결정한다: 
  $$S'_i(x_i) = \frac{y_{i+1}-y_{i-1}}{x_{i+1} - x_{i-1}},~ S'_{i}(x_{i+1}) = \frac{y_{i+2}-y_i}{x_{i+2}-x_{i}}$$ 

// Catmull-Rom spline interpolating {q[i]};
void CatmullRom(const std::vector<CfPt> &q, CDC *pDC) {
    if (q.size() < 2) return;
    for (int k = 0; k < q.size()-1; k++) {
        int km1 = max(k-1, 0);
        int kp2 = min(k+2, q.size()-1);

        /* compute spline coefficients */
        double dx  = q[k+1].x - q[k].x; 
        double dy  = (q[k+1].y - q[k].y) / dx;
        double dy1 = (q[k+1].y - [km1].y) / (q[k+1].x - q[km1].x);
        double dy2 = (q[kp2].y - q[k].y) / (q[kp2].x - q[k].x);
        double a = q[k].y;
        double b = dy1;
        double c = (3 * dy - 2 * dy1 - dy2) / dx;
        double d = (-2 * dy +   dy1 + dy2) / dx / dx;
#define STEPS 20
        if (k == 0) pDC->MoveTo(q[0]);
        const double deltax = q[k+1].x - q[k].x;
        for (double j = 1; j <= STEPS; j++)  { 
            double x = j * deltax / STEPS;
            /* use Horner's rule to calculate cubic polynomial */
            double y = ((d * x + c)*x + b)*x + a;
            pDC->LineTo(CfPt(q[k].x + x, y));
        }
    }
}

주어진 점들을 감싸는 최소  면적의 직사각형은 먼저 convex hull을 구한 후 rotating caliper을 써서 구할 수 있다. 다음 코드는 반시계방향으로 정렬된 convex hull을 이용해서 minimum volume box을 구한다. 처음 좌표축 방향으로 정렬된 bounding box을 이용해서 가장 작은 각을 형성하는 convex hull의 한 변을 찾고, 그 변으로 만들어지는 bounding box를 구한다. 반시계방향으로 bounding box가 convex hull의 변을 순차적으로 접하도록 회전시키면서 최소 면적인 경우를 찾으면 된다. 이는 convex hull을 바닥에서 굴리면서 bounding box을 찾는 경우로 생각해 보면 더 쉽다. bounding box가 convex hull의 모든 변을 다 순회하면 알고리즘이 종료되므로  time complexity는 $O(N)$이다.

void MinVolumeBox(const std::vector<CfPt>& chull, CfPt mvbox[4]) {
    if (chull.size() < 3) return;
    std::vector<CfPt> edges(chull.size());
    std::vector<int>  visited(chull.size(), 0);  // initialized = 0;
    enum {BB_NONE, BB_LEFT, BB_RIGHT, BB_BOTTOM, BB_TOP};
    for (int i = 0; i < chull.size(); i++) {
        const CfPt e = chull[(i+1) % chull.size()] - chull[i];
        edges[i] = e / e.Norm();// make a unit vector;
    }
    // 먼저 좌표의 최소-최대 값을 구한다;
    // index을 순환할 때, 항상 최대-최소가 (index+chull.size)%chull.size가 되게
    double xmin = chull[0].x, xmax = xmin;
    double ymin = chull[0].y, ymax = curr;
    int left = 0, bot = 0, right = 0, top = 0 ;
    for (int i = 1; i < chull.size(); i++) {
        if (chull[i].x <= xmin)      { xmin = chull[i].x; left = i;}
        else if (chull[i].x >= xmax) { xmax = chull[i].x; right = i;}
        if (chull[i].y <= ymin)      { ymin = chull[i].y; bot = i;}
        else if (chull[i].y >= ymax) { ymax = chull[i].y; top = i;}
    }
    if (chull[0].x <= xmin)       { xmin = chull[0].x; left = 0;}
    else if ( chull[0].x >= xmax) { xmax = chull[0].x; right = 0;}
    if (chull[0].y <= ymin)       { ymin = chull[0].y; bot = 0;}
    else if (chull[0].y >= ymax)  { ymax = chull[0].y; top = 0;}

    /*____Initial mvbox Setting____*/
    CfPt center = CfPt(xmin + xmax, ymin + ymax)/2;
    CfPt ax0 = CfPt(1,0), ax1 = CfPt(0,1);	// mvbox의 unit vectors;
    CfPt eh = ax0, ev = ax1;				// working unit vectors;
    double hscale  = (xmax - xmin) / 2; 
    double vscale  = (ymax - ymin) / 2;
    double minarea = hscale * vscale;		// area / 4;
    int done = 0;
    while (!done) {
        // edges[bot]는 chull[bot+1]-chull[bot];
        // edges[bot],edges[right],edges[top],edges[left] 중에서
        // 현재의 mvbox의 변(eh, ev,-eh,-ev)과 가장 작은 각(max_cosine)을 이루는 
        // edge를 찾음;
        int iflag = BB_NONE;
        double	max_cosine = 0;
        double	cosine = eh * edges[bot];
        if (cosine > max_cosine) { max_cosine = cosine; iflag = BB_BOTTOM;}

        cosine = ev * edges[right]; ;
        if (cosine > max_cosine) { max_cosine = cosine; iflag = BB_RIGHT;}
        
        cosine = -eh * edges[top]; 
        if (cosine > max_cosine) { max_cosine = cosine; iflag = BB_TOP;}
        
        cosine = -ev * edges[left]; 
        if (cosine > max_cosine) { max_cosine = cosine; iflag = BB_LEFT;}
        
        switch (iflag) {
        case BB_BOTTOM:
            if (!visited[bot]) {// edges[bot]가 mvbox의 한변에 포함됨;
                eh = edges[bot]; ev = eh.Rot90(); 
                // 다음 변으로 회전;
                visited[bot] = 1;
                if (++bot == chull.size()) bot = 0;
            } else done = 1;
            break;
        case BB_RIGHT:
            if (!visited[right]) {// edges[right]가 mvbox의 한변에 포함됨;
                ev = edges[right]; eh = -ev.Rot90(); 
                // 다음 변으로 회전;
                visited[right] = 1;
                if (++right == chull.size()) right = 0;
            } else done = 1;
            break;
        case BB_TOP:
            if (!visited[top]) { // edges[top]이 mvbox의 한변에 포함됨;
                eh = -edges[top]; ev = eh.Rot90(); 
                // 다음 변으로 회전;
                visited[top] = 1;
                if (++top == chull.size()) top = 0;
            } else done = 1;
            break;
        case BB_LEFT:
            if (!visited[left]) {// edges[left]가 mvbox의 한변에 포함됨;
                ev = -edges[left]; eh = -ev.Rot90(); 
                // 다음 변으로 회전;
                visited[left] = 1;
                if (++left == chull.size()) left = 0;
            } else done = 1;
            break;
        case BB_NONE: // polygon이 직사각형임;
            done = 1;
            break;
        }
        // check area of mvbox;
        double len0 = eh * (chull[right] - chull[left]) / 2;
        double len1 = ev * (chull[top] - chull[bot]) / 2;
        double area = len0 * len1; //면적의 1/4 임;
        if (area < minarea) {
            minarea	= area;
            ax0	= eh; ax1 = ev;
            hscale = len0; vscale = len1;
            // box center: tt = left 기준 top-bot을 연결하는 vector/2;
            CfPt tt = (chull[top] + chull[bot]) / 2 - chull[left];
            double zz = ev * tt; 
            center = eh * len0 + ev * zz + chull[left];
        }
    }
    eh = ax0 * hscale; // recover original scale;
    ev = ax1 * vscale; // recover original scale;
    mvbox[0] = center + eh + ev;
    mvbox[1] = center - eh + ev;
    mvbox[2] = center - eh - ev;
    mvbox[3] = center + eh - ev;
}