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

JS 시뮬레이터 vs 실제 C# 구현 - 영역별 일치도

🧩 BoundaryMode 별 count 계산 차이 Rectangular 기준

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
// left = Math.Min(r = 2, x = 0) = 0, right = Math.Min(r = 2, w - 1) = 2
// winW = 0 + 2 + 1 = 3 → cx = [1, 2, 1] 합 = 4
// winH = 0 + 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);               // 1. 버킷 초기화

    double total = 0;
    for (int i = 0; i < count; i++)
    {
        double w = weights[i];
        bucket[values[i]] += w;                // 2. 픽셀 값 → 버킷에 가중치 누적
        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++)              // 3. 0 → 255 순서 탐색
    {
        double w = bucket[i];
        if (w > 0)
        {
            acc += w;

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

            if (acc >= half - eps)             // 4. 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 Weighted Median 과 완전히 동일한 파이프라인을 가집니다. 유일한 차이는 Step 1 의 가중치 생성 함수입니다.

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

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

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

1. 원리의 차이 - 가중치 분포 형태

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

2. 이미지 결과물 차이 - 시나리오별 분석

실제 픽셀 값 예시로 두 필터가 어떻게 다른 결과를 도출하는지 확인하세요. (픽셀 값 0 ~ 255 기준)

3. 장단점 종합 비교

5. 커널 크기 (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; // 근사 비교

6. 선택 가이드 - 언제 어떤 필터를?

7. 수식 나란히 비교 - 같은 r 에서 각 픽셀의 기여도

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

// r = 2, σF = 6.0 : σ = 5 / 6.0 ≈ 0.83
// 중앙 기여(1D) : ≈ 47.9% (σF = 6.0 기준)
// 최외곽(1D) : ≈ 2.7% (σF = 6.0 기준)
// σ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 (Smoothing) 과 derivOrder > 0 (Gradient) 에서 처리 경로가 완전히 다릅니다. Smoothing 은 2-pass, 그래디언트는 12-pass (채널 × 4) 입니다.

// ApplySavitzkyGolaySeparable8 - totalPasses 계산 (RGB 8-bit 파이프라인)
int totalPasses;
if (isSg)
totalPasses = (derivOrder == 0) ? 2 : 12;  // Smoothing 2 / Gradient 12 (채널 3 × 4-pass)
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 에 캐싱됩니다.

// Smoothing 계수 캐시
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 방향 Smoothing
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 방향 Smoothing
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 는 범위 밖 픽셀에 대해 처리 방식을 결정합니다.

🔬 심층 분석 : 대칭적 축소(Symmetric Shrinking) vs AdaptiveMask

"AdaptiveMask 대신 커널 창을 항상 대칭으로 줄이면 완벽한 Zero-phase 를 달성할 수 있지 않을까?" 라는 질문에 대한 심층 분석입니다. 결론부터 말하면 1D 에서는 합리적 절충이 될 수 있지만, 2D 에서는 구조적으로 실용성이 붕괴됩니다.

// x = 1, r = 3 → wx = -3, -2, -1, 0, 1, 2, 3
//              skip skip skip ✓  ✓  ✓  ✓ → 4개 유효
for (int wx = -r; wx <= r; wx++) {
    int nx = x + wx;
    if ((uint) nx >= (uint) width) continue; // OOB 건너뜀
    double w = coeff1D[wx + r]; // 전체 커널 가중치 유지
    // ...
}

// x = 1, r = 3 → symR = min(1, width - 2) = 1
//              wx = -1, 0, 1 → 3 개 유효 (항상 대칭)
int symR = Math.Min(Math.Min(r, x),
                    width - 1 - x);
for (int wx = -symR; wx <= symR; wx++) {
    double w = coeff1D[wx + r]; // 중심부 가중치만
    int nx = x + wx;
    // ...
}