Image Matting: Knockout method

Ref: https://grail.cs.washington.edu/projects/digital-matting/image-matting/  

세번째 이미지 = alpha 값이 적용된 전경(Gandalf)

// knockout method;
// trimap을 만들 때 전경/배경/중간 영역(경계) 구분;
#define BACKGROUND 0    // 배경 마스크 색상;
#define FOREGROUND 255  // 전경 마스크 색상;
#define UNKNOWN    128  // 경계 영역의 생상;
#define EDGEMARK   3    // 전경과 배경의 경계선을 표시하는 색상
// 주어진 픽셀에서 경계점까지 떨어진 거리(제곱)에 따라 경계점을 
// 정렬하기 위해 필요함;
struct PtDist {
    int x, y, d;
    PtDist() {}
    PtDist(int x, int y): x(x), y(y), d(0) {}
};
struct RGBTRIO {
    double rgb[3];  // b-g-r
};
static 
bool compare(const PtDist& a, const PtDist& b) {
    // 거리 제곱의 ascending order로 정렬;
    return a.d < b.d;
}
// 최소거리의 DISTANCE_THRESHOLD배 이내에 있는 경계 픽셀만을 이용해서 
// alpha값을 추정한다:
#define DISTANCE_THRESHOLD  (5.0)
static
RGBTRIO estimateColor(int x0, int y0, CRaster& raster, 
                   std::vector<PtDist>& pts) {
    for (int i = pts.size(); i-->0;) {
        int x = pts[i].x - x0; 
        int y = pts[i].y - y0;
        pts[i].d = x * x + y * y; 
    }
    std::sort(pts.begin(), pts.end(), compare);
    const double minDist = sqrt(double(pts[0].d));
    ASSERT(minDist > 0);
    const double maxDist = DISTANCE_THRESHOLD * minDist;
    const double inv_diff = 1.0/ (maxDist - minDist);
    double wsum = 0, rsum = 0, gsum = 0, bsum = 0, d;
    for (int i = 0; i < pts.size(); i++) { // sorted data;
        if ((d = sqrt(double(pts[i].d))) < maxDist) {
            double w = (maxDist - d) * inv_diff;
            DWORD q = raster.GetPixel0(pts[i].x, pts[i].y);
            bsum += w * ((q      ) & 0xFF); // blue;
            gsum += w * ((q >>  8) & 0xFF); // green;
            rsum += w * ((q >> 16) & 0xFF); // red;
            wsum += w;
        }
        else break;
    }
    if (wsum == 0) wsum = 1;
    RGBTRIO Q;
    Q.rgb[0] = (bsum / wsum);
    Q.rgb[1] = (gsum / wsum);
    Q.rgb[2] = (rsum / wsum);
    return Q;
}
#define AddVec(A,B,C)	{C[0]=A[0]+B[0];C[1]=A[1]+B[1];C[2]=A[2]+B[2];}
#define SubVec(A,B,C)	{C[0]=A[0]-B[0];C[1]=A[1]-B[1];C[2]=A[2]-B[2];}
#define NormSq(A)       (A[0]*A[0]+A[1]*A[1]+A[2]*A[2])
#define DotProd(A,B)	(A[0]*B[0]+A[1]*B[1]+A[2]*B[2])
#define ScaleVec(A,s)   {A[0]*=s;A[1]*=s;A[2]*=s;}
static
double estimateAlpha(double F[3], double B[3], double C[3]) {
    double N[3], V[3], Bp[3];
    SubVec(F, B, N);						//N=F-B;
    double nn = NormSq(N);					//nn=N.N;
    if (nn == 0) return 1;	
    SubVec(C, B, V);						//V=C-B;
    double proj = DotProd (N, V) / nn;		//proj=(V.N)/(N.N);
    ScaleVec(N, proj);						//N=N(V.N)/(N.N);
    SubVec(C, N, Bp);						//Bp=C-N(V.N)/(N.N);
    //
    SubVec(F, Bp, N);						//N=F-Bp;
    nn = NormSq(N);							//nn=N.N;
    if (nn == 0) return 1;
    SubVec(C, Bp, V);						//V=C-Bp;
    proj = DotProd(N, V) / nn;				//proj=(N.V)/(N.N);

    return proj < 0 ? 0: proj > 1 ? 1: proj;
}
void FindMatt(CRaster& raster, CRaster& trimap, CRaster& alphaMap) {
    ASSERT((trimap.GetBPP()==8) && (raster.GetBPP()==24));
    const CSize sz = raster.GetSize();
    alphaMap.SetDimensions(sz, 32);// Bitmap with alpha-channel;
    std::vector<PtDist> ptFore = binaryTracing(trimap, FOREGROUND);
    std::vector<PtDist> ptBack = binaryTracing(trimap, BACKGROUND);
    
    double C[3];
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)trimap.GetLinePtr(y);
        BYTE *q = (BYTE *)raster.GetLinePtr(y);
        BYTE *r = (BYTE *)alphaMap.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            if (UNKNOWN == *p) {// unknown pixel(0=background,255=foreground);
                RGBTRIO Qf = estimateColor(x, y, raster, ptFore);
                RGBTRIO Qb = estimateColor(x, y, raster, ptBack);
                C[0] = q[0]; C[1] = q[1]; C[2] = q[2];
                double alpha = estimateAlpha(Qf.rgb, Qb.rgb, C);
                // if ( alpha >= 0.5)
                r[0] = (int)max(0, min(255, Qf.rgb[0]));
                r[1] = (int)max(0, min(255, Qf.rgb[1]));
                r[2] = (int)max(0, min(255, Qf.rgb[2]));
                r[3] = int(alpha * 255);       //[0:1]->[0:255]
            } else if (FOREGROUND == *p) {
                r[0] = q[0];
                r[1] = q[1];
                r[2] = q[2];
                r[3] = 0xFF;
            }
            p++; q += 3; r += 4;
        }
    }
}