이미지에 Poisson Noise 넣기

평균이 $\lambda$인 Poisson 분포를 가지는 random number를 만드는 레시피는 \([0,1]\)에서 추출된 균일한 random number 수열이 $u_1, u_2, u_3,...$로 주어질 때, $u_1  u_2 u_3....u_{k+1}\le  e^{-\lambda}$을 처음 만족하는 $k$를 찾으면 됨이 Knuth에 의해서 알려졌다.

int PoissonValue(double lambda) { // lambda = mean value
    double L = exp(-lambda);
    int k = 0;
    double p = 1;
    do {
        k++;
        // generate uniform random number u in [0,1] and let p ← p × u.
        p *= double(rand()) / RAND_MAX;
    } while (p > L);
    return (k - 1);
}

디지털 영상에서 각 픽셀은 CCD 또는 CMOS 이미지 센서가 일정한 노출 시간 동안 검출한 광자 수에 비례하는 값을 나타낸다. 광자의 도착은 Poisson 과정을 따르므로, 각 픽셀에서 검출되는 광자 수는 평균이 $\lambda$인 Poisson 분포를 따른다고 모델링할 수 있다. 실제 영상에서는 $\lambda$를 알 수 없으므로 관측된 픽셀값을 $\lambda$의 추정값으로 간주한다. 따라서 Poisson 노이즈를 추가하려면 각 픽셀값 \(I(x,y)\)를 평균으로 하는 Poisson 분포에서 난수를 생성하여 원래 픽셀값을 그 난수로 대체하면 된다.

double AddPoissonNoise(CRaster& raster, CRaster& noised) {
    if (raster.IsEmpty()) return -1;
    CSize sz = raster.GetSize();
    int channel = raster.GetBPP() / 8;
    noised = raster;
    srand(unsigned(time(0)));
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        BYTE *q = (BYTE *)noised.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++)
            for (int i = 0; i < channel; i++) {
            	int a = PoissonValue(*p++);
            	*q++ = a > 0xFF ? 0xFF: a;
            }
    }
    return PSNR(raster, noised);
}
double PSNR(CRaster& raster, CRaster& noised) {
    ASSERT(raster.GetBPP()  == noised.GetBPP());
    ASSERT(raster.GetSize() == noised.GetSize());
    CSize sz = raster.GetSize();
    int channel = raster.GetBPP() / 8;
    double mse = 0; // mean-squared error;
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        BYTE *q = (BYTE *)noised.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++)
            for (int i = 0; i < channel; i++) {
                int a = *p++ - *q++;
                mse +=  a * a;
            }
    }
    mse /= (channel * sz.cx * sz.cy);
    double decibel = 10. * log10(255 * 255 / mse);
    return decibel;
};

 

아래 예제는 \( \cos(t) \)로 표현되는 신호에 Poisson noise를 추가하는 과정을 보여준다. Poisson 분포는 평균값이 0 이상인 경우에만 정의되므로, 음수 값을 갖는 신호에는 직접 적용할 수 없다. 따라서 먼저 원신호에 1을 더하여 \(\cos(t)+1\)을 사용하면 신호의 범위가 \([0,2]\)가 되어 모든 값이 0 이상이 된다.

그러나 \(\cos(t)+1\)은 실제 event의 평균 개수가 아니라 상대적인 신호의 세기를 나타낸다. 따라서 Poisson 노이즈를 현실적으로 모사하려면 이를 실제 평균 event 수에 대응시키는 스케일링이 필요하다. 예를 들어 신호의 최댓값이 평균 50개의 event에 해당한다고 가정하면, 먼저 신호를 50배 하여 \(\lambda(t)=50\bigl(\cos(t)+1\bigr)\)을 각 시점에서의 Poisson 분포의 평균값으로 사용한다. 그런 다음 평균이 \(\lambda(t)\)인 Poisson 난수를 생성하고, 생성된 값을 다시 50으로 나누어 원래의 신호 스케일로 복원한다. 즉,

\[y(t)=\frac{\operatorname{Poisson}\left(50(\cos(t)+1)\right)}{50}\]과 같이 계산하면 된다.

이 과정으로 생성된 신호는 평균적으로 원래의 \(\cos(t)+1\) 신호를 유지하면서도, 실제 광자 계수(photon counting) 과정에서 발생하는 shot noise와 동일한 통계적 특성을 갖는 Poisson 노이즈가 추가된다. 또한 스케일링 계수는 노이즈의 크기를 결정하는 중요한 역할을 한다. 스케일링 계수가 클수록 평균 event 수가 증가하므로 상대적인 노이즈의 크기(표준편차/평균)는 \(\frac{1}{\sqrt{\lambda}}\)에 비례하여 감소한다. 반대로 스케일링 계수가 작을수록 평균 event 수가 감소하여 Poisson 노이즈가 더욱 크게 나타난다. 따라서 스케일링 계수는 실제 센서에서 측정되는 평균 광자 수 또는 계수(count)의 크기를 반영하는 파라미터로 볼 수 있다.

void PoissonNoiseEx1D() {
	const double scale = 50.;
    FILE *fp = fopen("poisson.txt", "w");
    for (double t = 0; t < 10; t += 0.01) {
        double y = scale * (cos(t) + 1.);
        fprintf(fp, "%f %f %f\n", t, y, (double)PoissonValue(y) / scale);
    }
    fclose(fp);
}

Poisson 잡음은 random noise처럼 단순히 더해지는 잡음과 달리 신호값이 클 때 잡음도 더 커짐을 알 수 있다. 그러나 SNR은 신호값이 클수록 더 좋아진다.