Lehtinen, J., Soler, C., Moon, B., Yoon, S.-E., Zwicker, M., Jarosz, W., … Rousselle, F. (2015). Recent Advances in Adaptive Sampling and Reconstruction for Monte Carlo Rendering. Computer Graphics Forum, 34(2), 667–681. https://doi.org/10.1111/cgf.12592
Abstract
그리디 알고리즘을 이용해서 지역적으로 최적화된 가우시안 필터 크기를 찾고 그 결과에 따라 샘플링 정도를 다르게 해보았다.
Introduction
Monte Carlo rendering 을 위해서 샘플링을 최적화된 방식으로 하고, 그 결과를 적절한 필터를 이용해서 정돈하는 것은 매우 중요하다. 우리는 이 문제를 MSE (mean squared error)을 이용해서 필터와 샘플링의 optimization 문제라고 생각하였다. greedy 알고리즘에서 두 개의 반복 구문을 사용하였다. 첫째는 MS에러를 최소화 하는 가우시안 필터의 크기를 각 픽셀마다 찾는 것이고, 둘쨰는 그 결과를 이용하여 샘플링 정도를 다르게 하는 것이다. 우리의 결과를 ground truth와 비교해 보면 그 차이가 무시할 만큼 작다는 것을 알 수 있다.
Previous Work
Algorithm Overview
Monte Carlo samples를 위한 정해진 샘플수에서 MSE를 최소화 시키는 필터와 샘플링 분배를 찾는 것이 목표이다. 근데 필터라고 하면 너무 종류가 많으니깐 잠재적인 필터들의 집합을 각 픽셀마다 부여한다. 즉 필터의 종류를 선택하는 것은 이번 문제의 핵심이 아니다. 그림에서 보듯, initial set of samples에서 출발하여 픽셀의 MSE를 최소화하는 필터를 찾는다. 둘쨰로 선택되어진 필터에 기반하여 MSE를 더 줄일 수 있는 샘플링 기술을 찾는다. 이 반복은 특정한 기준을 만족시킬경우 종료된다. 필터가 작은 스케일을 가질경우 작은 bias(편향)을 가지지만 큰 nosie 가 발생한다. (Blur를 생각하면 편함)
Filter Selection
Incremental MSE Minimization
여기서 MSE 는 [math]\displaystyle{ bias^2 + variance }[/math]으로 정의한다. 이때 대다수의 픽셀에서 필터가 작은(fine)것에서 큰(coarse) 것으로 이동함에 따라 픽셀 bias는 증가하고 분산은 감소하는 경향을 보인다. 여기서 분산은 픽셀이 주위픽셀과 얼마나 다른 정도인지를 이야기 하고, bias 는 실제 예측치와 얼마나 떨어진 값인지를 이야기 한다. 점진적 증가를 가정한다면, MSE의 변화량은 [math]\displaystyle{ MSE[c]-MSE[f]=Bias[c]^2-Bias[f]^2+Var[c]-Var[f] }[/math]로 표현된다. 즉 우리가 할일은 fine 에서 coase로 필터를 조정하면서 MSE의 변화량를 추적하다 그 변화량이 양수가 되면 멈추면 된다.
Quadratic Approximation
여기서는 어떻게 사전 ground truth 에 대한 지식 없이 bias 를 계산할 것인지에 대한 문제를 다룬다. 우선 c와 f는 각각 coarse 그리고 fine 필터의 축약형이다. r은 각 필터의 스케일을 말한다. 각 필터의 스케일이란 각 필터가 어느 픽셀만큼 뒤에 0으로 수렴되는 지를 결정한다. 예를 들어 가우시안 필터에서 스케일은 convolution에 사용된 Gaussian 필터의 표준편차에 해당하며, σ=0일 때 f0(x,y) = f(x,y)로 정의된다. 이때 Silverman 1986에 따르면 스케일과 편향은 다음과 같은 관계가 있다.
- [math]\displaystyle{ Bias[c]=\frac{r_c^2}{r_f^2}Bias[f] }[/math]
이를 이용해서 Bias가 측정된 값 - 실제 값임을 이용해 식을 변형시키면, 실제 값이 없어지면서 다음과 같은 수식이 도출된다.
- [math]\displaystyle{ Bias[c]^2-Bias[f]^2\approx \frac{r_c^2+r_f^2}{r_c^2-r_f^2}(c-f)^2 }[/math]
결국 최종 에러의 변화량은 다음과 같이 정의된다.
- [math]\displaystyle{ S\approx \frac{r_f^2+r_c^2}{r_c^2-r_f^2}(c-f)^2 + Var[c]-Var[f] }[/math]
(원 논문은 부호가 잘 못 되어 있다.) S를 여기선 scale selector 라고 하겠다. S가 양수이면 f로 필터링 하고 그 외의 경우에는 다음 coarse scale로 이동한다. 이때 Var는 실제값에 대한 정보가 필요하지 않음으로 (여기선 필터가 적용되는 필셀들의 평균만이 들어감으로) 근사항이 필요하지 않다.
Estimation from Noisy Data
c, f 그리고 Var는 실제 이미지를 필터에 통과했을때의 값으로 정의되어 있기때문에, Monte Carlo Rendering를 사용하는 우리로써는 바로 적용시킬 수는 없다. 그렇지만 linear 필터들은 사실 몬테 카를로 샘플들의 가중합이라는 사실에 접근하면, c와 f는 각 픽셀을 필터에 적용시킨 값이라는 사실을 알 수 있다. (너무 당연한 이야기. 우리가 적용시키려는 샘플은 당연히 필터 거친 픽셀 칼라이겠다.) 굳이 수식으로 표현하자면 s를 몬테카를로 샘플링에서 얻은 각 픽셀 샘플들이라고 하면 다음과 같이 정의된다.
- [math]\displaystyle{ f=\sum_{i=1...k}^{ }w_is_i }[/math]
또한 분산은 다음과 같이 가중치의 제곱과 각각 샘플들의 Var로 구성되어 진다. 독립된 두 확률변수의 합의 분산은 각각의 분산의 합과 같다는 성질이 있다. 또한 분산에 상수를 곱하면 제곱으로 더해짐으로 아래와 같이 정해지는 것이 자연스럽다.
- [math]\displaystyle{ Var[f]=\sum_{i=1...k}(w_i)^2Var[s_i] }[/math]
하지만 Var[s]는 몬테카를로 샘플링을 통해서 얻은 것으로 실제값이 될 수 없다. 따라서 우리는 Var를 구하기 위해서 경험적으로 얻은 샘플들의 집합 P에 대해서만 구할 수 밖에 없다. 여기서 s바는 P에 있는 샘플들의 평균이고 |P|는 P에 있는 원소들의 개수이다. 즉 통계적 추정을 사용하여 분산을 구하게 된다. 각 통계적 추적에 사용되는 원소들의 집합 P는 같은 Pixel에 있는 모든 sample들의 집합이다.
- [math]\displaystyle{ Var[s_i]\approx\frac{1}{|P|-1}\sum_{j\in P}^{ }(s_j-\bar{s})^2 }[/math]
Emprical Anlysis
1D function (선을 위아래로 변형시켜서 만들어지는 함수)를 균등한 샘플들로 샘플링해보자. f와 c는 정해진 스케일 (여기선 2)로 스케일이 조절된다. 1D function 이 0인 상수함수이고 샘플들이 노이즈로부터 얻어졌을 경우, 편향은 0 (노이즈가 0+- 로 주어지기 때문에)이고 분산은 필터가 커질수록 작아지기 때문에 이상적인 상황에서는 반드시 coarse 필터를 골라야 한다. 하지만 scale selector는 종종 잘못된 선택을하기도 한다. 왜냐하면 샘플링이 노이즈로부터 이루어 지기 때문이다. (이때 bias와 variance는 서로 독립이다. 왜냐하면 bias는 필터들의 가중합이고 variance는 다시 픽셀에서만의 분산으로 재 정의되기 때문이다.) 즉 error rate는 bias + variance 그래프에서 0과 1사이의 값을 선택할 확률이 된다. error rate 는 gamma 를 통해서 조절되며, 특수한 관계가 있다. gamma를 bias 에 곱해서 bias에 편향을 준다. 이를 통해서 error rate가 클수록 좀더 큰 필터를 사용하도록 유도한다.
error rate 과 다른 변수들의 상관관계
우선 우리는 error rate가 주어진 경우 (얼만큼의 잘못된 오차를 허용할 것인가?), 그것을 조절하기 위해서 bias와 variance의 가중합을 고려해야 한다. error rate와 다른 변수들의 상관관계를 경험적으로 살펴보면 다음과 같다.
- Sample Variance: 상관 없음
- Filter Radius: 큰 r은 살짝 작은 error rate을 만들어 낸다.
- Number of Samples per Pixel: error rate는 크게는 상수값을 유지한다. 그러나 sample 수가 작으면 error rate 는 증가한다.
결국 경험적으로 error rate가 주어졌을 경우 그것을 이루어 내기 위한 bias 에 곱해지는 가중치는 다음과 같이 정의된다. [math]\displaystyle{ z(r)=-log(1-1.9r)^(1/\sqrt{2}) }[/math] 이 근사는 0.4미만에서만 통용되지만 0.4이상의 error는 무의미 하다고 생각되어 지기 때문에 유용하다.
Post-Processing the Filter Selection
error rate 에 의해서 잘못된 선택을 하게 되면 이미지에 기분나쁜 스파이크로 표현이 된다. 따라서 이러한 outlier들을 제거하기 위한 작업이 필요하다. 이에 따르면 filter들을 선택할 확률은 연속적인 흐름으로 나타나는 경우가 많다. 따라서 인접한 픽셀이 불 연속적인 radius를 사용한다면 최대한 연속적인 radius를 사용하도록 조절해준다. 이 방식으로 하면 디테일한 부분을 놓칠 수 있자만, 보다 안정적인 결과를 얻을 수 있다. 이는 현재 픽셀을 얼마만큼 크기의 가우시안 필터를 적용시킬 것인지 정하는 맵 (벡터라든지 어레이로 정의되어 있을 것임)을 다시 추가적인 가우시안 필터를 통과시킴으로써 이룰 수 있다. 현재 픽셀의 스텝만큼 현재 픽셀은 제외하고 (현재 픽셀이 너무 기여하는 바가 커짐으로) 가우시안 필터를 적용시킴으로써 이러한 목표를 달성한다.
Sample Distribution
추가적으로 MSE를 최소화 하기위해서 새로운 샘플들을 배치하는 작업이다. 우선 예측된 MSE의 변화량들을 선택한 필터에 도달할 때 까지 모두 더한다. 이 값을 선택된 scale + epsilon으로 나누어서 상대적인 MSE값을 구한다. 여기서 epsilon을 더하는 이유는 매우 어두운 영역에서 너무 확대 해석하는 것을 막기 위해서이다. 상대적인 MSE는 필터에 기여한 샘플수에 역으로 비례하기 때문에 n개의 샘플을 주어진 픽셀에 더하는 것은 MSE를 [math]\displaystyle{ \frac{n_s}{n+n_s} }[/math]만큼 줄이게 된다. 결국 상대적인 MSE는 [math]\displaystyle{ \frac{n}{n+n_s} }[/math]만큼 줄이게 된다. 이에 따라서 그들의 잠재적인 에러 감소에 따른 prioiry queue를 가지고서, 큐에서 m개의 픽셀을 선택해서 n개의 샘플만큼 픽셀들의 필터 범위에 뿌려준다. 추가적으로 1개의 샘플을 선택한 픽셀 위치에 더해준다. 이런식으로 MSE를 최대한 감소시키는 픽셀들만을 선택해서 샘플링 할 수 있다. 예를 들어 픽셀당 n개의 샘플을 8개의 iteration에서 요구한다고 해보자. 이미지가 M개의 픽셀들로 구성되어 있고 미리 k개의 샘플로 initialization 작업을 할경우 각 반복마다 N = M(n-k)/8개의 샘플을 요구한다. 여기서 initialization 작업은 미리 몇개의 샘플들로 추가 적인 sampling을 위한 작업을 준비하는 작업을 말한다. 그리고 iteration은 샘플링을 하고 필터를 통과시키는 것을 몇번 반복했는지를 말한다. 즉 m = N/n이 된다. 즉 각 반복마다 m개의 픽셀에 N개의 샘플을 하도록 유도된다. 이러한 N개의 샘플은 priority queue 에 의해서 MSE를 최소화하는 것이 선택되도록 된다.
implementation
PBRT로 머시기 머시기
Conclusion
만약 빛이 너무 안비치는 어두운 영역의 경우에는 샘플링이 제대로 이루어 지지 않는 문제가 발생한다. 또한 엣지 근처에서 발생한 노이즈는 필터 할 수 없다는 문제가 있다.
질문
- 왜 MSE 는 [math]\displaystyle{ bias^2 + variance }[/math]으로 정의되었는가? https://en.wikipedia.org/wiki/Mean_squared_error 참고
- Silverman 1986은 왜 저렇게 Bias간의 관계를 설정했는가?
- Relative MSE 와 MSE 의 차이점은 무었인가?
코드
void Denoiser::setColor(Bitmap *source) {
input = source; result = new Bitmap(size);
}
void Denoiser::initialize(int imageWidht, int imageHeight, int sampleCount) {
this->sampleCount = sampleCount;
size = Vector2i(imageWidht, imageHeight);
samples = new std::vector<std::vector<Color3f>*>();
samples->reserve(imageWidht * imageHeight);
for(int i=0;i<imageWidht * imageHeight;i++) {
samples->push_back(new std::vector<Color3f>());
samples->at(i)->reserve(sampleCount);
}
variances = new std::vector<float>();
variances->reserve(imageWidht * imageHeight);
for(int i=0;i<imageHeight * imageWidht;i++) {
variances->push_back(0.f);
}
filterScales = new std::vector<int>();
filterScales->reserve(imageWidht * imageHeight);
for(int i=0;i<imageHeight * imageWidht;i++) {
filterScales->push_back(0);
}
float lowSampleCount = 1-1.f/sampleCount;
float gammaResult = -log(1-pow((1.9*gamma),1/sqrt(2)));
biasFactor = lowSampleCount * gammaResult;
}
void Denoiser::calculate() {
std::cout<<"Start Adaptive Denoising... "<<endl;
clock_t start, end;
start = clock();
//initailizing the filters
filters = new std::vector<float*>();
filterPdfs = new std::vector<float>();
filters->reserve(DENOISER_FILTER_NUM);
filterPdfs->reserve(DENOISER_FILTER_NUM);
getVarianceFinal();
printf(" [gamma %.2f]\n", gamma);
for(int i=0;i<DENOISER_FILTER_NUM;i++) {
float* temp = new float[DENOISER_FILTER_SIZE * DENOISER_FILTER_SIZE];
float tempPdf = 0;
float stddev = initialStddev * pow(stddevStep, i);
for(int x = 0;x<DENOISER_FILTER_SIZE;x++) {
for(int y=0;y<DENOISER_FILTER_SIZE;y++) {
float value = eval(x - DENOISER_FILTER_SIZE/2, y - DENOISER_FILTER_SIZE/2, stddev, DENOISER_FILTER_SIZE/2);
temp[x + y*DENOISER_FILTER_SIZE] = value;
tempPdf += value;
}
}
printf(" [filter %d] --standard deviation %f--\n", i + 1, stddev);
filters->push_back(temp);
filterPdfs->push_back(tempPdf);
}
for(int y=0;y<size.y();y++) {
for(int x=0;x<size.x();x++) {
calculateMSE(x,y,0);
}
}
filterSelection();
std::string filterPath = "denoiserData.ppm";
std::ofstream out(filterPath);
out<< "P3\n" << size.x() << " " << size.y() << "\n255\n";
for(int y=0;y<size.y();y++) {
for(int x=0;x<size.x();x++) {
//pushToFilter(x, y, 0);
float returnVal = filterScales->at(x + y*size.x());
int ir = int(255.99 * (1-(returnVal/DENOISER_FILTER_NUM)));
out<< ir << " " << ir << " " << ir << "\n";
pushToFilter(x, y);
}
}
out.close();
end = clock();
printf("Done! Took %.3fsec\n",(float)(end-start)/CLOCKS_PER_SEC);
}
Color3f Denoiser::pushToFilter(int x, int y, int step) {
if(x < DENOISER_FILTER_SIZE/2 || x >= size.x()-DENOISER_FILTER_SIZE/2) {
result->coeffRef(y, x) = input->coeffRef(y, x);
return input->coeffRef(y, x);
}
if(y < DENOISER_FILTER_SIZE/2 || y >= size.y()-DENOISER_FILTER_SIZE/2) {
result->coeffRef(y, x) = input->coeffRef(y, x);
return input->coeffRef(y, x);
}
result->coeffRef(y, x) = Color3f(0.f);
if(step == -1)
step = filterScales->at(x + y * size.x());
for(int i=0;i<DENOISER_FILTER_SIZE;i++) {
for(int j=0;j<DENOISER_FILTER_SIZE;j++) {
result->coeffRef(y,x) += filters->at(step)[i + DENOISER_FILTER_SIZE * j]*input->coeffRef(y+j-DENOISER_FILTER_SIZE/2,x+i-DENOISER_FILTER_SIZE/2)/filterPdfs->at(step);
}
}
return result->coeffRef(y, x);
}
int Denoiser::calculateMSE(int x, int y, int step) {
if(step >= DENOISER_FILTER_NUM - 1) {
filterScales->at(x + size.x() * y) = step;
return step + 1;
}
if(x < DENOISER_FILTER_SIZE/2 || x >= size.x()-DENOISER_FILTER_SIZE/2) {
filterScales->at(x + size.x() * y) = step;
return 0;
}
if(y < DENOISER_FILTER_SIZE/2 || y >= size.y()-DENOISER_FILTER_SIZE/2) {
filterScales->at(x + size.x() * y) = step;
return 0;
}
int returnValue = step;
float bias = 0.f;
float variance = 0.f;
float coarseScale = initialStddev * pow(stddevStep, step + 1);
float fineScale = initialStddev * pow(stddevStep , step);
float doubleCS = coarseScale * coarseScale;
float doubleFS = fineScale * fineScale;
Color3f fineResult = pushToFilter(x, y, step);
Color3f coarseResult = pushToFilter(x, y, step + 1);
float cminusf = pow(fineResult.x() - coarseResult.x(),2) +
pow(fineResult.y() - coarseResult.y(),2) +
pow(fineResult.z() - coarseResult.z(),2);
bias = ((doubleCS + doubleFS) / (doubleCS - doubleFS)) * cminusf;
bias *= biasFactor;
for(int i=0;i<DENOISER_FILTER_SIZE;i++) {
for(int j=0;j<DENOISER_FILTER_SIZE;j++) {
int index = (y + j - DENOISER_FILTER_SIZE/2) * size.x() + x + i - DENOISER_FILTER_SIZE/2;
variance += (filters->at(step + 1)[i + DENOISER_FILTER_SIZE * j] - filters->at(step)[i + DENOISER_FILTER_SIZE * j]) * variances->at(index);
}
}
if(bias + variance <= 0) {
returnValue = calculateMSE(x, y, step + 1);
}
else {
if(step == 0) {
filterScales->at(x + size.x() * y) = step;
return 0;
}
filterScales->at(x + size.x() * y) = step;
}
return returnValue + 1;
}
void Denoiser::setVariance(int x, int y, Color3f value) {
samples->at(x + y * size.x())->push_back(value);
}
void Denoiser::getVarianceFinal() {
for(int y=0;y<size.y();y++) {
for(int x=0;x<size.x();x++) {
float tempSum = 0.f;
for(int i=0;i<sampleCount;i++) {
tempSum += pow(samples->at(y * size.x() + x)->at(i).x() - input->coeffRef(y, x).x(),2) +
pow(samples->at(y * size.x() + x)->at(i).y() - input->coeffRef(y, x).y(),2) +
pow(samples->at(y * size.x() + x)->at(i).z() - input->coeffRef(y, x).z(),2);
}
variances->at(y * size.x() + x) = tempSum / (sampleCount-1);
}
}
}
void Denoiser::filterSelection() {
std::vector<int>* tempFilterScales;
tempFilterScales = new std::vector<int>();
tempFilterScales->reserve(size.x() * size.y());
for(int i=0;i<size.x() * size.y();i++) {
tempFilterScales->push_back(0);
}
for(int x = 0;x<size.x();x++) {
for(int y=0;y<size.y();y++) {
if(x < DENOISER_FILTER_SIZE/2 || x >= size.x()-DENOISER_FILTER_SIZE/2) {
tempFilterScales->at(x + y * size.x()) = filterScales->at(x + y * size.x());
continue;
}
if(y < DENOISER_FILTER_SIZE/2 || y >= size.y()-DENOISER_FILTER_SIZE/2) {
tempFilterScales->at(x + y * size.x()) = filterScales->at(x + y * size.x());
continue;
}
int step = filterScales->at(x + y * size.x()) + 1;
float tempStep = 0.f;
if(step > DENOISER_FILTER_NUM - 1)
step -= 1;
for(int i=0;i<DENOISER_FILTER_SIZE;i++) {
for(int j=0;j<DENOISER_FILTER_SIZE;j++) {
if(i == DENOISER_FILTER_SIZE/2 && j == DENOISER_FILTER_SIZE/2)
continue;
int index = (y + j - DENOISER_FILTER_SIZE/2) * size.x() + x + i - DENOISER_FILTER_SIZE/2;
tempStep += filters->at(step)[i + DENOISER_FILTER_SIZE * j]*filterScales->at(index);
}
}
tempFilterScales->at(x + y * size.x()) = (int)tempStep/(filterPdfs->at(step) - filters->at(step)[DENOISER_FILTER_SIZE + DENOISER_FILTER_SIZE*DENOISER_FILTER_SIZE/2]);
}
}
delete filterScales;
filterScales = tempFilterScales;
}
float Denoiser::eval(int x, int y, float stddev, float radius) const {
float invTwoStddev = 1.f/(stddev * stddev); float a = INV_TWOPI * 1 * invTwoStddev; float b = -(x*x + y*y) * 0.5f * invTwoStddev; return a * exp(b);
}
Bitmap* Denoiser::getResult() {
return result;
}