Image Smoothing Filters
1D→2D · RGB · 8/16bit

🧩 BoundaryMode별 count 계산 차이

Rectangular 평균에서 count(분모)는 경계 모드에 따라 다르게 계산됩니다. 이 차이가 경계 픽셀의 밝기에 직접 영향.

🔄 전체 처리 흐름 — 단계별 분석

Rectangular 필터의 처리 과정을 4단계로 분해합니다. 모든 필터 중 가장 단순하지만, 경계 처리와 8/16-bit 분기를 이해하는 데 중요합니다.

⚡ Interior 빠른 경로 — 8-bit vs 16-bit 구현 차이

16-bit 파이프라인은 Interior(경계에서 r 이상 떨어진) 픽셀에 대해 GetIndex1D() 호출을 완전히 생략하는 빠른 경로를 가집니다. 8-bit는 모든 픽셀에 GetIndex1D()를 호출합니다.

// ApplyRect8: Adaptive 분기만 별도,
// 나머지는 모두 GetIndex1D() 경로

if (mode == BoundaryMode.Adaptive)
{
    // 윈도우 축소 → 직접 루프
    int left = Math.Min(r, x);
    int right = Math.Min(r, width-1-x);
    // ...
}
else
{
    // 모든 픽셀에 GetIndex1D() 호출
    for (wy) {
        int ny = GetIndex1D(y+wy, height, mode);
        for (wx) {
            int nx = GetIndex1D(x+wx, width, mode);
            // ← 내부 픽셀도 이 검사를 통과
        }
    }
}

int fullCount = (2*r+1) * (2*r+1);
bool yInterior = y >= r && y < height-r;

if (mode == Adaptive
    && !(yInterior && x >= r && x < width-r))
{
    // 경계 Adaptive: 윈도우 축소
}
else if (yInterior && x >= r && x < width-r)
{
    // ★ Interior 빠른 경로:
    // GetIndex1D() 호출 0회!
    double sB = 0;
    for (int wy = -r; wy <= r; wy++)
        for (int wx = -r; wx <= r; wx++)
            sB += planes.B[y+wy, x+wx];
    outB[y,x] = sB / fullCount;
    // ↑ 사전 계산된 fullCount 재사용
}
else
{
    // 경계: GetIndex1D() 경로
}

⚡ 성능 특성 — Rectangular 필터

📊 denom(분모) 계산 — 경로별 차이

이항 가중 평균의 분모는 경계 모드에 따라 다르게 계산됩니다. 이 차이가 경계 픽셀의 밝기에 직접 영향.

📐 Adaptive 경계 — 축소 커널 재생성 상세

Adaptive 모드에서 경계 픽셀은 축소된 윈도우에 맞는 새로운 Binomial 계수를 생성합니다. GetBinomial1D(winW/H) 캐싱으로 동일 크기 재계산 없음.

// Adaptive 분기 핵심
int left  = Math.Min(r, x);
int right = Math.Min(r, width - 1 - x);
int winW = left + right + 1;  // 축소된 너비

var cx = GetBinomial1D(winW); // 새 Binomial 계수
var cy = GetBinomial1D(winH);

double sumWx = 0, sumWy = 0;
for (i) sumWx += cx[i];  // 축소 계수 합
for (i) sumWy += cy[i];
double denom = sumWx * sumWy; // 2D 외적 합

// 예: 모서리(x=0,y=0) r=2
// winW=1+2+1=3 → cx=[1,2,1] 합=4
// winH=1+2+1=3 → cy=[1,2,1] 합=4
// denom = 4×4 = 16 (원본 5×5: 16²=256)

🪣 Step 3 — Weighted Median 알고리즘 : Bucket 방식 (8-bit)

일반 중앙값은 값을 정렬해야 합니다(O(n log n)). 8-bit 이미지의 경우 픽셀값이 0~255 범위이므로, 크기 256짜리 버킷 배열에 가중치를 누적한 뒤 누적합이 전체의 절반을 넘는 지점을 찾는 것으로 O(n + 256) = O(n)에 해결합니다.

// SmoothingConductor.cs — WeightedMedianDouble8() 전체 코드
private static byte WeightedMedianDouble8(
    byte[] values, double[] weights,
    double[] bucket, int count)
{
    if (count <= 0) return 0;

    Array.Clear(bucket, 0, 256);          // ① 버킷 초기화

    double total = 0;
    for (int i = 0; i < count; i++)
    {
        double w = weights[i];
        bucket[values[i]] += w;           // ② 픽셀값 → 버킷에 가중치 누적
        total += w;
    }

    if (total <= 0) return 0;

    double half = total / 2.0;
    double eps  = total * 1e-12;         // 부동소수점 허용 오차
    double acc  = 0;

    for (int i = 0; i < 256; i++)        // ③ 0→255 순서 탐색
    {
        double w = bucket[i];
        if (w > 0)
        {
            acc += w;

            if (acc > half + eps)              // 중앙값 발견!
                return (byte)i;

            if (acc >= half - eps)             // ④ TieBreak 보간
            {
                for (int k = i + 1; k < 256; k++)
                {
                    if (bucket[k] > 0)
                        return (byte)((i + k + 1) >> 1);
                }
                return (byte)i;
            }
        }
    }
    return 255;
}

🎯 실제 사용 사례 — Binomial Weighted Median Filter

Binomial Weighted Median은 엣지와 원본 구조를 보존하면서 임펄스·랜덤 노이즈만 선택적으로 제거하는 것이 핵심 강점이며, 이 특성은 다양한 산업 및 연구 분야에서 실질적인 이점을 제공합니다.

🔄 전체 처리 흐름 — 단계별 분석

Gaussian Weighted Median은 Binomial Median과 완전히 동일한 파이프라인을 가집니다. 유일한 차이는 Step 1의 가중치 생성 함수입니다.

📐 sigmaFactor 영향 분석 — 가우시안 폭이 중앙값에 미치는 효과

sigmaFactor(σF)는 가우시안 커널의 유효 범위를 결정합니다. σF가 클수록 중앙에 집중, σF가 작을수록 넓게 분산됩니다. 이것이 중앙값 결과에 미치는 영향:

⚡ 성능 최적화 포인트 — Gaussian Median

① 원리의 차이 — 가중치 분포 형태

두 필터의 가중치를 같은 r에서 시각적으로 비교합니다. 슬라이더로 반경을 조절해 어떻게 달라지는지 확인하세요.

② 이미지 결과물 차이 — 시나리오별 분석

실제 픽셀값 예시로 두 필터가 어떻게 다른 결과를 내는지 계산합니다. (픽셀값 0~255 기준)

③ 장단점 종합 비교

⑤ 커널 크기(r) 증가 시 — Gaussian이 왜곡되고 Binomial이 원본에 가까운 이유

Weighted Median의 핵심 조건은 단 하나입니다: 중앙 픽셀의 누적 가중치가 50% 미만이어야 주변 픽셀들이 노이즈를 억제합니다. 50% 이상이면 가중 중앙값 탐색이 노이즈 픽셀에서 멈춰 노이즈가 그대로 출력됩니다.

double half = total / 2.0;    // 이항 계수=정수값 → 정확
double eps = total * 1e-12;  // 범용 코드 공유 (Binomial엔 불필요)
if (acc > half + eps) return i; // 정수값끼리 비교 → 결과 정확

double half = total / 2.0;   // 누적 오차 존재
double eps = total * 1e-12;  // 보정값 도입 필요
if (acc > half + eps) return i; // 근사 비교

⑥ 선택 가이드 — 언제 어떤 필터를?

⑦ 수식 나란히 비교 — 같은 r에서 각 픽셀의 기여도

// r=2: [1,4,6,4,1] / 16
// 중앙 기여: 6/16 = 37.5%
// 최외곽: 1/16 = 6.25%
// 비율: 6× 차이

// r=2, σF=6: σ=5/6≈0.83
// 중앙 기여(1D): ≈47.9% (σF=6 기준)
// 최외곽(1D): ≈2.7% (σF=6 기준)
// σF 변화로 비율 연속 조절 가능

📐 σ(sigma) 계산 — sigmaFactor의 역할

sigmaFactor(σF)는 가우시안 커널의 폭을 제어합니다. σF가 클수록 σ가 작아져 중앙에 집중, σF가 작을수록 σ가 커져 넓게 분산.

// 코드에서의 구현 — ApplyGaussian8() 시작 부분
int windowSize = checked(2 * r + 1);
double sigma = windowSize / sigmaFactor;  // σ = (2r+1) / σF
var g1D = ComputeGaussian1D(windowSize, sigma);

// Adaptive 경계: 축소된 윈도우에 맞는 σ를 재계산
double sigmaLocalX = winW / sigmaFactor; // winW < windowSize
double sigmaLocalY = winH / sigmaFactor;
var gx = ComputeGaussian1D(winW, sigmaLocalX);
var gy = ComputeGaussian1D(winH, sigmaLocalY);

🔄 전체 처리 흐름 — 단계별 분석

Gaussian Average 필터의 처리 과정을 Binomial Average와 비교하며 4단계로 분해합니다.

📊 BoundaryMode별 denom(분모) 계산 경로

Gaussian Average에서 분모(denom)는 경계 모드에 따라 다르게 계산됩니다. Binomial Average와 동일한 패턴입니다.

📐 Adaptive 경계 — 축소 가우시안 커널 재생성

Adaptive 모드에서 경계 픽셀은 축소된 윈도우에 맞는 새로운 가우시안 커널을 생성합니다. σ도 축소된 윈도우에 비례하여 재계산됩니다.

// ApplyGaussian8 — Adaptive 분기
int left  = Math.Min(r, x);
int right = Math.Min(r, width - 1 - x);
int winW  = left + right + 1;

// ★ σ도 축소된 윈도우에 비례하여 재계산
double sigmaLocalX = winW / sigmaFactor;
double sigmaLocalY = winH / sigmaFactor;
var gx = ComputeGaussian1D(winW, sigmaLocalX);
var gy = ComputeGaussian1D(winH, sigmaLocalY);

// 예: 모서리(x=0, y=0) r=2
// winW=3, σX=3/6=0.5 → 더 뾰족한 가우시안
// winH=3, σY=3/6=0.5 → 마찬가지

// denom 재계산 (합=1 정규화이므로 ~1.0)
double denom = 0;
for (yy) for (xx) denom += gy[yy] * gx[xx];
if (denom <= 0) denom = 1;

🔬 전체 처리 흐름 — derivOrder에 따른 분기

SG 필터는 derivOrder=0(스무딩)과 derivOrder>0(그래디언트)에서 처리 경로가 완전히 다릅니다. 스무딩은 2-pass, 그래디언트는 12-pass(채널×4)입니다.

// ApplySavitzkyGolaySeparable8 — totalPasses 계산
int totalPasses;
if (isSg)
    totalPasses = (derivOrder == 0) ? 2 : 12;  // 스무딩 2 / 그래디언트 12
else
    totalPasses = 1;

var proxy = new ProgressProxy(progress, checked(totalPasses * height));
// ↑ 총 진행 단위 = pass 수 × 이미지 높이

🔧 ComputeSgCoefficientsViaQR() 핵심 코드

일반적인 (AᵀA)⁻¹Aᵀ 정규 방정식이 아닌 Householder QR 경로를 사용합니다. AᵀA의 조건수 = A의 조건수²이므로, 높은 polyOrder에서 수치적 안정성을 위해 QR이 필수적입니다.

// 작업 복사본 W = A
var W = new double[windowSize, cols];
var vecs = new double[cols][];

for (int k = 0; k < cols; k++)
{
    int len = windowSize - k;
    var v = new double[len];
    for (int i = 0; i < len; i++)
        v[i] = W[k + i, k];

    double norm = Math.Sqrt(sigma);
    // 부호 선택: catastrophic cancellation 방지
    double alpha = v[0] >= 0
                 ? -norm : norm;
    v[0] -= alpha;

    // v 정규화
    double vnorm2 = 0;
    for (i) vnorm2 += v[i] * v[i];
    double inv = 1.0 / Math.Sqrt(vnorm2);
    for (i) v[i] *= inv;
    vecs[k] = v;

    // 나머지 열 반사:
    // W[k:,j] -= 2·v·(vᵀ·W[k:,j])
    for (int j = k; j < cols; j++)
    {
        double dot = 0;
        for (i) dot += v[i] * W[k+i,j];
        dot *= 2.0;
        for (i) W[k+i,j] -= dot * v[i];
    }
}
// W[0..cols-1, 0..cols-1] → R (상삼각)

var z = new double[cols];
for (int i = 0; i < cols; i++)
{
    double rhs = (i == derivOrder)
                 ? 1.0 : 0.0;
    for (int j = 0; j < i; j++)
        rhs -= W[j, i] * z[j];
    // Rᵀ[i,j] = R[j,i] = W[j,i]
    z[i] = rhs / W[i, i];
}

var h = new double[windowSize];
for (j < cols) h[j] = z[j];
// h[cols..] 은 0 유지

for (int k = cols-1; k >= 0; k--)
{
    var v = vecs[k];
    int len = windowSize - k;
    double dot = 0;
    for (i) dot += v[i] * h[k+i];
    dot *= 2.0;
    for (i) h[k+i] -= dot * v[i];
}

// STEP 5: 정규화 / 스케일링
if (derivOrder == 0) {
    double s = Σh[i];
    for (i) h[i] /= s; // 합=1 보존
} else {
    double scale =
        FactorialAsDouble(derivOrder)
        / Math.Pow(delta, derivOrder);
    for (i) h[i] *= scale;
}

🗃️ 비대칭 SG 계수 캐싱 — GetAsymmetricSg / GetAsymmetricSgDeriv

경계 픽셀에서 좌/우(또는 상/하) 가용 폭이 다르므로 비대칭 Vandermonde 행렬로 계수를 재계산합니다. 동일한 (left, right, polyOrder) 조합은 Dictionary에 캐싱됩니다.

// 스무딩 계수 캐시
private static readonly
  Dictionary<Tuple<int,int,int>, double[]>
    _sgAsymSmoothCache;

private static double[] GetAsymmetricSg(
    int left, int right, int polyOrder)
{
    var key = Tuple.Create(
        left, right, polyOrder);
    lock (_sgAsymSmoothCacheLock) {
        if (_sgAsymSmoothCache
            .TryGetValue(key, out var c))
            return c; // 캐시 적중
        c = ComputeSGAsymmetric(
            left, right, polyOrder,
            derivOrder:0, delta:1.0);
        _sgAsymSmoothCache[key] = c;
        return c;
    }
}

// 미분 계수 캐시 (4-tuple 키)
private static readonly
  Dictionary<Tuple<int,int,int,int>, double[]>
    _sgAsymDerivCache;

private static double[] GetAsymmetricSgDeriv(
    int left, int right,
    int polyOrder, int derivOrder)
{
    var key = Tuple.Create(
        left, right, polyOrder, derivOrder);
    lock (_sgAsymDerivCacheLock) {
        if (cache.TryGetValue(key, out var c))
            return c;
        c = ComputeSGAsymmetric(
            left, right, polyOrder,
            derivOrder, delta:1.0);
        cache[key] = c;
        return c;
    }
}

📐 derivOrder > 0 — 그래디언트 크기(Gradient Magnitude) 계산

derivOrder>0이면 X/Y 방향 미분을 합성하여 엣지 강도 맵을 생성합니다. 채널당 4-pass × 3채널 = 12-pass.

// ComputeGradientMagnitudeForChannel() — 채널당 4-pass
var tmp = new double[h, w];
var gx  = new double[h, w];

// Pass 1: X 방향 미분
Parallel.For(0, h, y => {
    for (x) tmp[y,x] = Convolve1D_X(src, y, x, w, r, coeffDeriv, mode, polyOrder, derivOrder);
    proxy?.StepRows(1);
});
// Pass 2: Y 방향 스무딩
Parallel.For(0, h, y => {
    for (x) gx[y,x] = Convolve1D_Y(tmp, y, x, h, r, coeffSmooth, mode, polyOrder, 0);
    proxy?.StepRows(1);
});
// Pass 3: X 방향 스무딩
Parallel.For(0, h, y => {
    for (x) tmp[y,x] = Convolve1D_X(src, y, x, w, r, coeffSmooth, mode, polyOrder, 0);
    proxy?.StepRows(1);
});
// Pass 4: Y 방향 미분 + 합성
Parallel.For(0, h, y => {
    for (x) {
        double gyVal = Convolve1D_Y(tmp, y, x, h, r, coeffDeriv, mode, polyOrder, derivOrder);
        double gxVal = gx[y,x];
        dst[y,x] = Math.Sqrt(gxVal*gxVal + gyVal*gyVal);
    }
    proxy?.StepRows(1);
});

📊 SG 필터 vs 다른 필터 — 특성 비교

🔍 왜 경계 처리가 필요한가? — 문제 원리

반경 r인 커널을 이미지 가장자리 픽셀에 적용하면, 커널의 일부가 이미지 밖(-1, -2, ... 또는 w, w+1, ...) 좌표를 참조합니다. 이 범위 밖 픽셀에 어떤 값을 할당하느냐가 BoundaryMode입니다.

모드 선택 가이드

📊 필터별 경계 처리 방식 차이

동일한 BoundaryMode라도 각 필터가 경계를 처리하는 방식은 다릅니다. 특히 Adaptive 모드에서 가중치 재계산 방식이 필터마다 크게 다릅니다.

🔬 SG 전용 — 비대칭 경계 계수의 특수성

SG 필터는 경계에서 다른 필터와 근본적으로 다른 접근을 합니다. 단순 축소가 아닌 비대칭 Vandermonde 행렬로 QR 분해를 수행합니다.

// Convolve1D_X — Adaptive/ValidOnly 경계 경로
int left  = Math.Min(r, x);
int right = Math.Min(r, width - 1 - x);
int localWindow = left + right + 1;

// polyOrder/derivOrder 자동 축소 (과적합 방지)
int localPoly  = Math.Min(polyOrder, localWindow - 1);
int localDeriv = Math.Min(derivOrder, localPoly);

// 비대칭 SG 계수 (Dictionary 캐싱)
double[] h = (localDeriv == 0)
    ? GetAsymmetricSg(left, right, localPoly)           // 스무딩 캐시
    : GetAsymmetricSgDeriv(left, right, localPoly, localDeriv); // 미분 캐시

// 비대칭 계수로 합성곱
double sum = 0;
int start = x - left;
for (int i = 0; i < h.Length; i++)
    sum += h[i] * src[y, start + i];
return sum;

📝 GetBoundaryMethodText() — 사용자 표시용 텍스트

public static string GetBoundaryMethodText(BoundaryMode mode) {
    switch (mode) {
        case Symmetric:    return "Symmetric (Mirror)";
        case Replicate:    return "Replicate (Edge Clamp)";
        case ZeroPad:      return "Zero Padding";
        case Adaptive:     return "Adaptive (Window Shrink)";
        case AdaptiveMask: return "Adaptive Mask (Skip Invalid)";
        case ValidOnly:    return "Valid Only";
    }
}