Reconstructing a height map from a normal map

2014-08-20 Permalink

At times it happens that we have the normal map but not the height map that was used to create it. For some rendering techniques, like parallax mapping, it is the height map which is needed. In addition, having the height map is beneficial for editing the normal map.

A height map and its normal map, used in the samples below.

Recovering the height map boils down to integrating the normal map. This is what is suggested in multiple places around the web (e.g. here, here and here), but it is not explained how to do this in practice.

The problem is that the normal map at hand is typically encoded with low precision, 8-bits per channel usually, and might be compressed with a lossy algorithm afterwards. In the end, the line integrals along closed curves aren’t zero, and so the end result significantly depends on the curves along which we integrate. In fact, no matter how we integrate, we end up with significant noise in the output.

Different integration paths give different results. None are good.

Through deconvolution

Denote by f(x, y) the height map and n(x, y) = (nx, ny, 1) the normal map that was calculated by some finite difference scheme, that is by convolving the original height map with the kernels dx and dy. For the sake of boundary condition, we assume that both f and n are periodic of periods w and h. We seek a function f which minimizes:

E = 2−1 ( ǁdx ∗ f − nxǁ2 + ǁdy ∗ f − nyǁ2 )

Taking the Fourier transform F, Dx, Dy, Nx, Ny of f, dx, dy, nx, ny respectively, and applying Parseval's theorem:

E = 2−1 ( ǁDx F − Nxǁ2 + ǁDy F − Nyǁ2 )

The last expression can be minimized by each F(x,y) independently. The minimum is achieved by:

F = (|Dx|2 + |Dy|2)−1 (conj(Dx) Nx + conj(Dy) Ny)

For when |Dx|2 + |Dy|2 = 0, E is independent of F(x,y), and so we can choose any value of the later. In such cases we seek the minimum norm solution, which is achieved by setting F = 0 at these points.

Applying the inverse Fourier transform to F leads to the minimum error, minimum norm solution f.

The reconstructed height map. Running time: < 50 ms


Sample use:

static const int w = 4, h = 4;
float normal[w*h*3] = { ... };
float result[w*h];
reconstruct_height_map1(normal, w, h, result);
// ... or ...
reconstruct_height_map2(normal, w, h, result);

The reconstruct_height_map1 function assumes one sided differences:

nx(x, y) = f(x + 1, y) − f(x, y)
ny(x, y) = f(x, y + 1) − f(x, y)

The reconstruct_height_map2 function assumes central differences:

nx(x, y) = f(x + 1, y) − f(x − 1, y)
ny(x, y) = f(x, y + 1) − f(x, y − 1)

I avoid dividing here by two so that to make use of the whole [−1, 1] range.

This C++ implementation is released to the public domain. It uses the GPLed FFTW library for computing the FFT.

#include <fftw3.h>
#include <complex>
#include <vector>

void reconstruct_height_map(const float *normal, float *dx, float *dy, int width, int height, float *result)
	typedef std::complex<float> C;
	fftwf_plan plan;

	std::vector<float> nx(width*height), ny(width*height);
	for(int y = 0, i = 0; y < height; ++y)
	for(int x = 0; x < width; ++x, ++i, normal += 3)
		nx[i] = normal[0]/normal[2], ny[i] = normal[1]/normal[2];

	const int half_width = width/2 + 1;
	std::vector<C> Nx(half_width*height), Ny(half_width*height);
	std::vector<C> Dx(half_width*height), Dy(half_width*height);

	plan = fftwf_plan_dft_r2c_2d(height, width, &nx[0], (fftwf_complex*)&Nx[0], FFTW_ESTIMATE);
	fftwf_execute_dft_r2c(plan, &nx[0], (fftwf_complex*)&Nx[0]);
	fftwf_execute_dft_r2c(plan, &ny[0], (fftwf_complex*)&Ny[0]);
	fftwf_execute_dft_r2c(plan, &dx[0], (fftwf_complex*)&Dx[0]);
	fftwf_execute_dft_r2c(plan, &dy[0], (fftwf_complex*)&Dy[0]);

	std::vector<C> F(half_width*height);
	for(int y = 0, i = 0; y < height; ++y)
	for(int x = 0; x < half_width; ++x, ++i)
		float denom = width * height * (norm(Dx[i]) + norm(Dy[i]));
		F[i] = denom > 0 ? - (Dx[i] * Nx[i] + Dy[i] * Ny[i]) / denom : 0;

	plan = fftwf_plan_dft_c2r_2d(height, width, (fftwf_complex*)&F[0], &result[0], FFTW_ESTIMATE);

void reconstruct_height_map1(const float *normal, int width, int height, float *result)
	std::vector<float> dx(width*height), dy(width*height);
	dx[0] = 1, dx[1] = -1;
	dy[0] = 1, dy[width] = -1;
	reconstruct_height_map(normal, &dx[0], &dy[0], width, height, result);

void reconstruct_height_map2(const float *normal, int width, int height, float *result)
	std::vector<float> dx(width*height), dy(width*height);
	dx[width-1] = 1, dx[1] = -1;
	dy[width*(height-1)] = 1, dy[width] = -1;
	reconstruct_height_map(normal, &dx[0], &dy[0], width, height, result);

Note on the sample images

The 8-bit unsigned integer values of the external files were mapped to in-memory floats (and vice versa) by y = 255−1 x − 0.5 for the height maps, and y = 127−1 (x − 128) for the normal maps.

The normal map was deliberately compressed with JPEG to emphasize the precision loss. However, even with lossless compression, the 8-bit quantization of the normal map is enough to produce visible artifacts, especially when wrapping the image.

At last, just for fun, here is a sample where the source image is not a differential of anything at all, yet we are able to construct a height map and a normal map that loosely resemble the original. It does not make any sense, of course.

(a) Lena with B = 1 set and normalized. (b) Reconstructed height map scaled to the [0, 1] interval. (c) Normal map created from the reconstructed height map.

Share on