# Reconstructing a height map from a normal map

2014-08-20. Last revised on 2015-09-08.

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.

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.

## Through deconvolution

Denote by f(x, y) the height map and n(x, y) = (n_{x}, n_{y}, 1) the normal map that was calculated by some finite difference scheme, that is by convolving the original height map with the kernels d_{x} and d_{y}. 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} ( ǁd_{x} ∗ f − n_{x}ǁ^{2} + ǁd_{y} ∗ f − n_{y}ǁ^{2} )

Taking the Fourier transform F, D_{x}, D_{y}, N_{x}, N_{y} of f, d_{x}, d_{y}, n_{x}, n_{y} respectively, and applying Parseval's theorem:

E = 2^{−1} ( ǁD_{x} F − N_{x}ǁ^{2} + ǁD_{y} F − N_{y}ǁ^{2} )

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

F = (|D_{x}|^{2} + |D_{y}|^{2})^{−1} (conj(D_{x}) N_{x} + conj(D_{y}) N_{y})

For when |D_{x}|^{2} + |D_{y}|^{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.

## Code

Sample use:

The `reconstruct_height_map1`

function assumes one sided differences:

_{x}(x, y) = f(x + 1, y) − f(x, y)

n

_{y}(x, y) = f(x, y + 1) − f(x, y)

The `reconstruct_height_map2`

function assumes central differences:

_{x}(x, y) = f(x + 1, y) − f(x − 1, y)

n

_{y}(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
```
#include
void reconstruct_height_map(const float *normal, float *dx, float *dy, int width, int height, float *result)
{
typedef std::complex C;
fftwf_plan plan;
std::vector 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 Nx(half_width*height), Ny(half_width*height);
std::vector 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]);
fftwf_destroy_plan(plan);
std::vector 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);
fftwf_execute(plan);
fftwf_destroy_plan(plan);
}
void reconstruct_height_map1(const float *normal, int width, int height, float *result)
{
std::vector 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 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.