I've made a program that converts a normal map texture (DirectX DDS supported directly) into a height map image (PNG)

Here is a high-level overview of my algorithm:

• Loading the Normal Map: The normal map image is loaded and converted to a floating-point representation in the range [-1, 1].
  
  
• Preprocessing the Normal Map: The normal map is tiled into a larger map to handle boundary conditions more effectively during the processing.
  
  
• Defining Subregions: The larger normal map is divided into smaller subregions to allow for parallel processing.
  
  
• Computing Height Maps for Subregions: Each subregion's height map is computed using Poisson integration, which involves solving a large sparse linear system.
  
  What is the Poisson Equation?
  The Poisson equation is a fundamental equation in physics and engineering, often used to describe the potential field caused by a given charge or mass density distribution. In the context of image processing, it's used to reconstruct a surface from gradient information.
  
  In our case, the Poisson equation helps to convert the normal map into a height map by solving for the height values at each pixel based on the gradients.
  
  The Equation
  The Poisson equation in two dimensions can be written as:
  
  ∂²h/∂x² + ∂²h/∂y² = f(x, y)
  Here:
  
  • h is the height function we want to find.
  • f(x, y) is a function representing the known gradients.
    
  
  Discretizing the Poisson Equation
  To solve this on a computer, we convert the continuous equation into a discrete form using finite differences. This involves approximating the second derivatives with differences between neighboring pixel values.
  For a pixel at position (i, j):
  
  (h(i+1, j) - 2h(i, j) + h(i-1, j))/Δx² + (h(i, j+1) - 2h(i, j) + h(i, j-1))/Δy² = f(i, j)

  Assuming uniform grid spacing (Δx = Δy = 1 for simplicity):
  
  h(i+1, j) + h(i-1, j) + h(i, j+1) + h(i, j-1) - 4*h(i, j) = f(i, j)

  This forms a system of linear equations where each pixel’s height value depends on its neighbors.
  
  Matrix Representation
  
  We represent this system as a matrix equation:
  
  A * h = b
  
  Here:
  
  • A is a sparse matrix encoding the relationships between neighboring pixels.
  • h is the vector of unknown height values.
  • b is the vector derived from the gradient information f(x, y).
    
  
  Solving the System
  To find the height map, we solve the linear system:
  
  A * h = b
  using numerical methods. This gives us the height values at each pixel, which we then reshape into the height map.
  
  
• Combining Subregion Height Maps: The height maps from all subregions are combined into a single height map using a weighted blending technique. I use cosine smoothing in this version which allows me to minimize the boundary effect between subregions. Previously, I've used gaussian bluring and 50% overlap, but this method gives better results - you can see the difference between non-FFT and FFT version (less repeated boundary artifacts)
  
  
  
• Enhancing High-Frequency Details: High-frequency details are enhanced using multi-scale unsharp masking to improve the visual quality of the height map.
  
  
• Clipping Histogram Levels: The histogram levels of the height map are adjusted to improve contrast.
  
  
• S-Curve Histogram Adjustment: An S-curve adjustment is applied to smooth the height map
  
  
• Removing Spikes: Spikes in the height map are removed using a combination of median filtering and Gaussian blurring.
  
  
• Saving the Height Map: The final height map is saved as an image file.
  

Adobe (around 40 seconds)

This code (20 min calculation)

This code FFT version (30 seconds)