Bidirectional Path Tracing 4 - A tweak of the filter

In Veach97 10.3.1 each pixel value $I_{j}$ in the image is defined as a weighted average

$$I_{j} = \int\int_{D}h_{j}(u, v)I(u, v)dudv$$

where $D$ is the image region (in world space, not pixel space), and $h_{j}$ is the filter function (for example, a 2x2 box filter) for pixel j. When the filter function satisfies the criteria that $\int\int_{D}h_{j}(u, v)dudv=1$, we filtered the measurement equation integration result $I$ across lens and film area to one particular pixel $I_{j}$ . This is a bit counter intuitive for me at the beginning: The measurement equation measure up a value related to the whole lens and film, and pixel value (what image rendering is actually looking for) is a post process product from the measurement value (by multiply importance $W_{e}(x, \omega)$ and filter $h_{j}(u, v)$ ) .

Since measurement equation contains integration across the whole film ( the $D$ region), the samples we cast should be thrown across the whole film (For example, a 640x480 image with 4 samples per pixel, we throw out 640x480x4 samples to get the image value for pixel (123, 456). It doesn't sound make sense at first, but thinking further, all the other samples that pixel (123, 456) doesn't take get accumulated to other pixels, so we are not doing wasted work after all) The estimator for this equation would then be:

$$E(I_{j}) = 1/N \sum_{i=1}^{N}h_{j}(u_{i}, v_{i})I(u_{i}, v_{i}) / pdf_{A_{film}}(u_{i}, v_{i}) = D/N \sum_{i=1}^{N}h_{j}(u_{i}, v_{i})I(u_{i}, v_{i})$$

In the path tracing world that shooting ray from camera, we usually know which samples will contribute to a certain pixel (continue the 2x2 box filter example, the final (123, 456) pixel value would only be affect by the ray shooting across range (122, 455)->(124, 457) )  In fact, the filter equation usually used for path tracing looks like:

$$E(I_{j}) = \frac{\sum_{i=1}^{N}h_{j}(u_{i}, v_{i})I(u_{i}, v_{i})}{\sum_{i=1}^{N}h_{j}(u_{i}, v_{i})}$$

The coding logic is roughly like: for each sample we got a radiance value, we add value $h_{j}(u_{i}, v_{i})I(u_{i}, v_{i})$ across region this sample centered kernel can cover to the covered pixels, and add $h_{j}(u_{i}, v_{i})$ to a separate weight value associated to that pixel. At the end, we normalize each pixel with their associated weight. This works for camera based tracing method since, again, we know which samples would contribute to which pixel. This doesn't work for light tracing method since, for each ray shoot out from light, we don't know which pixel it will slam on film after bouncing around the scene.

As a result, Veach97 10.3.4.3 propose that using equation $E(I_{j}) = \frac{\sum_{i=1}^{N}h_{j}(u_{i}, v_{i})I(u_{i}, v_{i})}{ \sum_{i=1}^{N}h_{j}(u_{i}, v_{i})}$ for path tracing method, and use $E(I_{j})= D/N \sum_{i=1}^{N}h_{j}(u_{i}, v_{i})I(u_{i}, v_{i})$ for light tracing method. That is, for light tracing method, we throw out total N samples, whenever a sample hit the film and lens, we added the sample contribution $h_{j}(u_{i}, v_{i})I(u_{i}, v_{i})$ to nearby affected pixels. At the end, we multiply the whole image film by factor $D/N$ where $D$ is the film area in world space.

My previous filtering implementation use the path tracing method above, and there is one missing puzzle that need to be filled to get the light tracing method work: $\int\int_{D}h_{j}(u, v)dudv=1$ . The filter was not normalized to film world area space, for example, again the 2x2 box filter, $h_{j}(u_{i}, v_{i}) = 1$ if $(u_{i}, v_{i})$ is in the 2x2 kernel range else 0, the integration result for this filter would be 2x2 = 4 pixel area. To normalize this box filter, $h_{j}(u_{i}, v_{i})$ need to be divide by 4 pixel area so the final kernel integration become 1. I added a virtual method getNormalizedTerm for all filters to query kernel integration result in pixel space, and during the adding sample to final image space, I add in one more divide by pixel area step to normalize pixel space integration to world space. The reason I didn't put thing together is that I don't want the filter functionality depends on film size (which should able to change its size without worrying about what filter it's paired to) So here we go, the getNormalizedTerm :

	float BoxFilter::getNormalizeTerm() const {
	return 4.0f * mXWidth * mYWidth;
	}
	float TriangleFilter::getNormalizeTerm() const {
	// 4 * integrate (integrate (w - x) * (h - y) over 0-w) over 0-h =
	// w * w * h * h
	return mXWidth * mXWidth * mYWidth * mYWidth;
	}
	float GaussianFilter::getNormalizeTerm() const {
	// for now just simply use numercal approximation
	// this integration actually involves an erf function
	// that need to be approximated anyway...
	size_t step = 20;
	float deltaX = mXWidth / static_cast<float>(step);
	float deltaY = mYWidth / static_cast<float>(step);
	float result = 0.0f;
	for (size_t i = 0; i < step; ++i) {
	for (size_t j = 0; j < step; ++j) {
	result += 4.0f * deltaX *deltaY *
	gaussian(i * deltaX, mExpX) *
	gaussian(j * deltaY, mExpY);
	}
	}
	return result;
	}
	float MitchellFilter::getNormalizeTerm() const {
	return 4.0f * ((12 - 9 * mB - 6 * mC) / 4 +
	(-18 + 12 * mB +6 * mC) / 3 + (6 - 2 * mB) +
	15 * (-mB - 6 * mB) / 4 + 7 *(6 * mB + 30 * mC) / 3 +
	3 * (-12 * mB - 48 * mC) / 2 + (8 * mB + 24 * mC)) / 6.0f;
	}

view raw getNormalizedTerm.cpp hosted with ❤ by GitHub

And the related film constructor/addSample changes:

	Film::Film(...) {
	// skip some unrelated initializatio code above...
	mInvXRes = 1.0f / (float)mXRes;
	mInvYRes = 1.0f / (float)mYRes;
	// precompute filter equation as a lookup table
	// we only computer the right up corner since for all the supported
	// filters: f(x, y) = f(|x|, |y|)
	float deltaX = mFilter->getXWidth() / FILTER_TABLE_WIDTH;
	float deltaY = mFilter->getYWidth() / FILTER_TABLE_WIDTH;
	float normalizeTerm = mFilter->getNormalizeTerm();
	size_t index = 0;
	for(int y = 0; y < FILTER_TABLE_WIDTH; ++y) {
	float fy = y * deltaY;
	for(int x = 0; x < FILTER_TABLE_WIDTH; ++x) {
	float fx = x * deltaX;
	mFilterTable[index++] = mFilter->evaluate(fx, fy) /
	normalizeTerm;
	}
	}
	// skip some unrelated initializatio code below...
	}
	void Film::addSample(float imageX, float imageY, const Color& L) {
	// skip some unrelated initializatio code above...
	// transform continuous space sample to discrete space
	float dImageX = imageX - 0.5f;
	float dImageY = imageY - 0.5f;
	// calculate the pixel range covered by filter center at sample
	float xWidth = mFilter->getXWidth();
	float yWidth = mFilter->getYWidth();
	int x0 = ceilInt(dImageX - xWidth);
	int x1 = floorInt(dImageX + xWidth);
	int y0 = ceilInt(dImageY - yWidth);
	int y1 = floorInt(dImageY + yWidth);
	x0 = max(x0, mXStart);
	x1 = min(x1, mXStart + mXCount - 1);
	y0 = max(y0, mYStart);
	y1 = min(y1, mYStart + mYCount - 1);
	for(int y = y0; y <= y1; ++y) {
	float fy = fabs(FILTER_TABLE_WIDTH * (y - dImageY) / yWidth);
	int iy = min(floorInt(fy), FILTER_TABLE_WIDTH - 1);
	for(int x = x0; x <= x1; ++x) {
	float fx = fabs(FILTER_TABLE_WIDTH * (x - dImageX) / xWidth);
	int ix = min(floorInt(fx), FILTER_TABLE_WIDTH - 1);
	float weight = mFilterTable[iy * FILTER_TABLE_WIDTH + ix] *
	mInvPixelArea;
	int index = y * mXRes + x;
	mPixels[index].color += weight * L;
	mPixels[index].weight += weight;
	}
	}
	}

view raw film_update.cpp hosted with ❤ by GitHub

One thing that I feel a bit awkward is that it's actually not a easy thing to integrate Gaussian filter in finite range (surprise! consider how beautiful it converge when integrate to infinity) that turns out I use a clumsy numerical approximation to get an approximated integration result. It's not really beautiful, but it serves the purpose ;)

In the next post I would like to write down the thoughts while looking in Veach97 10.2 : the ways to form a path, and its contribution to the final measurement (which is the \alpha weight I briefly mentioned in previous post ) It will be the ground work to start writing light tracing, and push it further to bidirectional path tracing.