Getting Shape From Shading To Actually Work In Practice
The Mark of Cain approach to recovering surface normals from a single grayscale image sounds elegant on paper. You assume a Lambertian surface with a known lighting direction, then derive a reflectance map and use its isolines to propagate normals across the image. The theory checks out. The practice is where things usually fall apart if you haven't dealt with the edge cases. The algorithm builds a 2D lookup table that maps surface normal to grayscale intensity for a given light direction. That's the reflectance map, R(p, q), where p and q are the partial derivatives of the surface. You find isolines on that map, project them down as directions on the surface, and start marching from boundary pixels inward. Each step gives you a candidate normal, and the shading constraint tells you which one is consistent with the pixel value. I ran into a real problem early on when I was applying this to a scanned stone relief. The surface wasn't smooth. It had fine texture noise that completely broke the propagation. The isolines on the reflectance map are supposed to be continuous, but when your actual image has high-frequency detail, the intensity values jump around and your normal estimates go unstable after just a few pixels from the boundary. What I ended up doing was pre-filtering the image with a bilateral filter instead of a Gaussian. It preserved the edges of the relief while smoothing the micro-texture enough that the propagation held for most of the image. Roughly cut my failure rate from about 40 percent of the surface area down to under 10 percent.
The tricky part nobody mentions is the ambiguity at certain normal orientations. When the surface normal points directly toward or away from the light source, the reflectance map has singularities. At the Specularity point, the gradient of R with respect to p and q goes to zero, which means the isolines converge and become nearly tangent to each other. Your propagation direction becomes ill-defined and you get noise amplification. I found that masking out regions where the predicted intensity would exceed roughly 95 percent of the maximum was a practical way to avoid those zones entirely. You lose some coverage, but the recovered normals in the remaining areas are actually usable. Another thing that trips people up is the assumption that the light direction is known and fixed. In real photography, especially with softbox or diffuse setups, the effective light direction varies across the surface. The algorithm assumes a single distant point light. If your lighting is even slightly off, the entire normal field gets biased. I learned this the hard way when a batch of sculptures looked wrong in 3D. Turns out one light was about 15 degrees higher than I thought, and it shifted all the normals consistently. The fix was to estimate the light direction from known geometric features rather than assuming it. You can calibrate it by finding the region that appears brightest and working backward, or by using a calibration object with a known shape placed in the same lighting. Here's the honest assessment: this method works well on smooth, matte surfaces with controlled lighting. That's about it. On textured surfaces, noisy data, or non-Lambertian materials like polished metal or wet stone, it degrades fast. If you're working with real-world scans where you can't control the environment, you're probably better off combining this with multi-image approaches or switching to photometric stereo if you can get at least three images under different lighting. The single-image case is fundamentally under-constrained anyway, and no amount of clever propagation fixes that.
Implementation-wise, the core loop is straightforward. Compute the reflectance map for your chosen light direction, find isoline directions, set up your boundary conditions, and march inward using a priority queue ordered by distance from the boundary. The bottleneck is usually the interpolation step on the reflectance map, which dominates runtime if you're doing it naively. Precomputing the map on a grid and using bicubic interpolation gets it down to reasonable speeds. A clean Python implementation with NumPy runs in under a minute for a 512x512 image on modest hardware. For anyone actually using this, the biggest win is that it gives you a plausible normal field without needing multiple views or special equipment. The output won't match what a structured light scanner produces, but for artistic applications, quick prototyping, or as an initialization for a more expensive method, it's useful. Just don't expect it to handle anything beyond what the assumptions allow, and spend some time on your preprocessing and light calibration before you trust the results.
Get the Full Details
