Geometry & Recognition

Mathematica에서 shooting method를 이용해서 비선형 미분방정식의 해를 구하자. 

\begin{align}f' &= \frac{1}{x} (1-a) f \\ a'&=- \frac{x}{2} f^2 (f^2-1) \end{align}

경계조건은 $f(0)=0$, $f(\infty)\to 0$, $a(0)=0$을 만족시켜야 한다. $x$가 증가하면 $f(x)$는 빠르게 $1$로 수렴하므로 오른쪽 경계조건은 $f(15)\to1$로 잡아도 충분하다. $x=0$에서의 apparent singularity를 Mathematica가 처리할 수 있도록 방정식을 변형시켜주어야 한다.

In relaxation methods ODEs are replaced by approximate finite difference equations (FDEs) on a grid or mesh of points that spans the domain of interest.

알고리즘 및 사용 예제 관련문서들:

NumericalRecipiesInC/c17-3.pdf

NumericalRecipiesInC/c17-4.pdf

Numerical Recipes 에서 얻을 수 있는 필요한 코드들:

1) solvde() : main routine; 

2) bksub ()  : supplementary routine

3) pinvs ()  : supplementary routine

4) red ()  : supplementary routine

5) difeq ()  : user supplied  routine in which differential equations, boundary conditions, and their jacobians are defined. called by solvde().

주의점: 경계조건 설정에서 주의해야 한다;

1. Derivative Based Algorithms:

• Thresholded absolute gradient: $$ f  =\sum_x \sum _y |I (x+1, y) - I(x, y)| $$  where   $|I(x+1, y) - I(x, y)| > \text{threshold};$
• Squared gradient: $$f=\sum_x \sum_y ( I(x+1, y)- I(x, y)) ^2 $$ where $(I(x+1, y) - I(x, y))^2 > \text{threshold};$
• Brenner gradient: $$f=\sum_x \sum_y ( I(x+2, y)- I(x, y)) ^2 $$ where $(I(x+2, y) - I(x, y))^2 > \text{threshold};$
• Tenenbaum gradient: $$f = \sum_x \sum_y (\text{SobelX}^2 (x, y)+ \text{SobelY}^2(x,y));$$
• Energy Laplace: $$f = \sum_x \sum_y ( (L * I) (x, y))^2$$  where $$L =\left[\begin{array}{ccc} -1& -4 & -1 \\  -4 & 20 & -4 \\ -1 & -4 & -1 \end{array}  \right]$$ 
• Sum of modified Laplace: $$f = \sum_x \sum_y (\text{LaplaceX}(x, y)| + |\text{LaplaceY}(x, y)|) $$
• Sum of squared Gaussian derivatives $$f = \frac{1}{\text{total pixels}} \sum_x \sum_y ((\text{Gaussian derivative X}(x,y)) ^2 + (\text{Gaussian derivative Y}(x,y))^2 )$$ where $\sigma  = d/(2\sqrt{3})$, $d$ = dimension of the smallest feature;

2. Statistical Algorithms:

•  Variance Measure ; $$f = \frac{1}{\text{total pixels}}\sum _x \sum _y ( I  (x, y) - \overline{I})^2;$$ where $\overline{I}$ is the mean of image.
• Normalized Variance Measure ; $$f = \frac{1}{\text{total pixels}\times \overline{I}} \sum_x \sum_y ( I(x, y)- \overline{I})^2;$$
• Auto-correlation Measure: $$f = \sum_x\sum_y I(x, y) I(x+1, y) - \sum _x \sum_y I(x, y) I(x+2, y);$$
• Standard Deviation-based Correlation: $$f = \sum_x \sum_y I(x, y) I(x+1, y) - \overline{I}^2(x, y)\times \text{total pixels};$$

3. Histogram-based Algorithms:

• Range Algorithm: $$f =\text{max}\{i| \text{histogram}(i)>0\} - \text{min}\{ i| \text{histogram}(i)>0\};$$
• Entropy Algorithm: $$f=-\sum_{i=0}^{255} p(i) \log_{2} p(i)$$ where $p(i)= \text{histogram}(i)/\text{total pixels};$

4. Other Algorithms:

• Threshold Contents: $$f=\sum_x \sum_y I(x, y)$$ where $I(x,y) \ge \text{threshold};$
• Thresholded Pixel Count: $$f=\sum_x\sum_y \theta(\text{threshold}- I(x, y));$$ where $\theta(x)$ is the step-function.
• Image Power: $$f= \sum_{x}\sum_{y} I^2(x, y)$$ where $I(x, y) \ge\text{threshold};$

Ref: Dynamic evaluation of autofocusing for automated microscopic analysis of blood smear and pap smear, J. Microscopy,227, 15(2007).                           

영상에서 추출한 객체(object)의 경계선은 객체의 모양에 대한 많은 정보를 담고 있지만, 어떤 경우에는 불필요하게 많은 정보가 될 수도 있다. 이 경우 원래의 경계를 충분히 닮은 다각형으로 간략화하여서 불필요한 정보를 제거할 수 있다. 이 다각형으로의  근사가  물체의 경계를 얼마나 제대로 표현할 수 있는지에 대한 기준이 있어야 한다. $N$개의 점 $\{(x_i, y_i) | i = 1,... ,N\}$으로 표현된 경계가 있다고 하자. 가장 단순한 근사는 가장 멀리 떨어진 두 점($A, B$)으로 만들어진 선분일 것이다(선분의 길이가 대략적인 물체의 크기를 주고, 선분 방향이 물체의 orientation을 알려준다). 이 선분을 기준으로 다시 가장 멀리 떨어진 점($C$)을 찾아서 일정 거리 이상 떨어져 있으면 꼭짓점으로 추가하여 두 개의 선분으로 분할한다. 두 개의 선분에 대해서 동일한 분할 작업을 재귀적으로 반복하므로 좀 더 정밀하게 물체를 표현하는 다각형을 찾을 수 있다. 선분에서 가장 멀리 떨어진 지점이 기준 거리 이내이면 분할을 멈춘다. 닫힌 다각형 근사일 때는 시작점 = 끝점으로 조건을 부여하면 된다.

** open polyline으로 우리나라 경계를 단순화시키는 예: top-right 쪽이 시작과 끝이므로 simplication 과정에서 변화지 않는다.

// distance square from P to the line segment AB ;

double pt2SegDist2(CPoint P, CPoint A, CPoint B) ;

// P에서 두 점 A, B을 연결하는 선분까지 거리
double pt2SegDist2(CPoint P, CPoint A, CPoint B) {
    double dx = B.x - A.x, dy = B.y - A.y ;
    double lenAB2 = dx * dx + dy * dy ;
    double du = P.x - A.x, dv = P.y - A.y ;
    if (lenAB2 == 0.) return du * du + dv * dv;
    double dot = dx * du + dy * dv ; 
    if (dot <= 0.) return du * du + dv * dv;
    else if (dot >= lenAB2) {
        du = P.x - B.x; dv = P.y - B.y;
        return du * du + dv * dv;
    } else {
        double slash = du * dy - dv * dx ;
        return slash * slash / lenAB2;
    }
};

//  recursive Douglas-Peucker 알고리즘;

void DouglasPeucker(double tol, std::vector<CPoint>& vtx, int start, int end, 
                   std::vector<bool>& mark) {
    if (end <= start + 1) return; //종료 조건;
    // vtx[start] to vtx[end]을 잇는 선분을 기준으로 분할점을 찾음;
    int break_id = start;			// 선분에서 가장 먼 꼭짓점;
    double maxdist2 = 0;
    const double tolsq = tol * tol;		
    for (int i = start + 1; i < end; i++) {
        double dist2 = pt2SegDist2(vtx[i], vtx[start], vtx[end]);
        if (dist2 <=  maxdist2) continue;
        break_id = i;
        maxdist2 = dist2;
    } 
    if (maxdist2 > tolsq) {		// 가장 먼 꼭짓점까지 거리가 임계값을 넘으면 => 분할;
        mark[break] = true;		// vtx[break_id]를 유지;
        // vtx[break_id] 기준으로 좌/우 polyline을 분할시도;
        DouglasPeucker(tol, vtx, start, break_id, mark);
        DouglasPeucker(tol, vtx, break_id, end, mark);
    }
};

// driver routine;

int polySimplify(double tol, std::vector <CPoint>& V, std::vector <CPoint>& decimated) {
    if (V.size() < 2) return 0;    
    std::vector<bool> mark(V.size(), false); // 꼭짓점 마킹용 버퍼;
    // 알고리즘 호출 전에 너무 붙어있는 꼭짓점을 제거하면 보다 빠르게 동작;
    // 폐곡선이 아닌 경우 처음,끝은 marking; 
    // 폐곡선은 가장 멀리 떨어진 두점중 하나만 해도 충분;
    mark.front() = mark.back() = true; 
    DouglasPeucker( tol, V, 0, V.size()-1, mark );
    // save marked vertices;
    decimated.clear();    
    for (int i = 0; i < V.size(); i++)
        if (mark[i]) decimated.push_back(V[i]);
    return decimated.size();
}

nonrecursive version:

// nonrecursive version;
int DouglasPeucker_nonrec(double tol, std::vector<CPoint>& points,
                            std::vector<CPoint>& decimated) {
    if (points.size() <= 2) return;
    double tol2 = tol * tol; //squaring tolerance;
    // 미리 가까이 충분히 가까이 있는 꼭지점은 제거하면 더 빨라짐;
    std::vector<CPoint> vt;
    vt.push_back(points.front()) ;  // start at the beginning;
    for (int i = 1, pv = 0; i < points.size(); i++) {
        int dx = points[i].x - points[pv].x, dy = points[i].y - points[pv].y;
        double dist2 = dx * dx + dy * dy ;
        if (dist2 < tol2) continue;
        vt.push_back(points[i]); 
        pv = i; // 다시 거리를 재는 기준점;
    }
    if (vt.back()!=points.back()) vt.push_back(points.back());  // finish at the end
    points.swap(vt);
    //
    std::vector<bool> mark(points.size(), false); // vertices marking
    std::vector<int> stack(points.size());  // 분할된 라인 처음-끝 저장;
    // 폐곡선이 아닌 경우 처음,끝은 marking; 
    // 폐곡선은 가장 멀리 떨어진 두점 중 하나만 해도 충분
    mark.front() = mark.back() = true;
    int top = -1;
    // 분할 해야할 폴리라인의 양끝점을 스택에 넣음;
    stack[++top] = 0;
    stack[++top] = points.size() - 1;
    while (top >= 0) {
        // 분할을 시도할 구간의 두 끝점을 꺼냄;
        int end = stack[top--];
        int start = stack[top--];
        // 최대로 멀리 떨어진 점;
        double max_distance2 = 0;
        int break_id = -1;
        for (int i = start + 1; i < end; i++) {
            double dsq = pt2SegDist2(points[i], points[start], points[end]);
            if (dsq > max_distance2) {
                max_distance2 = dsq;
                break_id = i;
            }
        }
        // 주어진 임계값 이상 떨어졌으면 꼭지점을 유지하도록 표시;
        if (max_distance2 > tol2) {
            mark[break_id] = true;       
            // 주어진 꼭지점을 기준으로 양쪽을 다시 검사하기 위해 스택에 쌓음;
            stack[++top] = start;
            stack[++top] = break_id;
            stack[++top] = break_id;
            stack[++top] = end;
        }
    }
    // marking 된 꼭지점을 출력;
    decimated.clear();
    for (int i = 0; i < points.size(); i++)
        if (mark[i]) decimated.push_back(points[i]);
    return decimated.size();
}