Normal Integration via Inverse Plane Fittingwith Minimum Point-to-Plane Distance1

Xu Cao1 Boxin Shi2,3 Fumio Okura1 Yasuyuki Matsushita1
Osaka University 2 Peking University 3 Peng Cheng Laboratory

Abstract
This paper presents a surface normal integration methodthat solves an inverse problem of local plane fitting. Surface reconstruction from normal maps is essential in photometric shape reconstruction. To this end, we formulate
normal integration in the camera coordinates and jointly
solve for 3D point positions and local plane displacements.
Unlike existing methods that consider the vertical distances
between 3D points, we minimize the sum of squared pointto-plane distances. Our method can deal with both orthographic or perspective normal maps with arbitrary boundaries. Compared to existing normal integration methods,
our method avoids the checkerboard artifact and performs
more robustly against natural boundaries, sharp features,and outliers. We further provide a geometric analysis of thesource of artifacts that appear in previous methods based on
our plane fitting formulation. Experimental results on analytically computed, synthetic, and real-world surfaces showthat our method yields accurate and stable reconstruction
for both orthographic and perspective normal maps1 .

1. Introduction
Reconstructing a surface from its surface normal map,a problem called normal integration [28], is essential forphotometric 3D reconstruction, such as shape from shading [19], photometric stereo [34], and shape from polarization [20]. Among all approaches to normal integration, variational-based methods have been studied for many
years [14, 16, 28], which are based on functional optimization. The objectives are turned into linear systems by either discretizing the functional itself [14] or its optimal
condition [28, 29] in the form of discrete Poisson’s equation. Recently, a non-variational approach based on discrete geometry processing (DGP) [37] was proposed fornormal integration. DGP also solves a linear system to estimate the surface. For all the three linear systems, however,
1 Source code is available at https : / / github . com /hoshino042/NormalIntegration.

Functional [14]Poisson [29]

DGP [37]

Our reconstruction
Figure 1. (Top) Typical problems in solving the linear systems fornormal integration. From left to right: The checkerboard artifact,
Gibbs phenomenon near sharp features, and sensitivity to outliers.
(Bottom) Our results from the same normal maps as above. Thereis no checkerboard artifact, Gibbs phenomenon almost disappears,
and outliers hardly distort the surface.

their least squares approximate solutions suffer from several
problems. For discrete functional, a checkerboard artifact
was reported in [38], as shown in Fig. 1 top-left. Discrete
Poisson’s equation and DGP are sensitive to sharp featuresor outliers, as shown in Fig. 1 top-middle and top-right.
Various remedies have been studied to improve the robustness against these artifacts. First, regularization [13, 14]or smoothness term [38] was used to refine the reconstruction quality. However, introducing a regularization term requires tuning a weighting factor, increasing the effort to obtain a satisfactory result. Second, 2-norm of the residual
vector was replaced by p-norm [6] to increase the robustness; however, it comes with heavier computation. Third,
normal vectors with small elevation angles was treated as
outliers [37], but the heuristic thresholding does not ensure
the robustness to outliers.
Unlike the existing approaches, we focus on the residual vector of the linear systems. Discretizing the functionalor Poisson’s equation by finite difference leads to measuring the difference between the height values of 3D points.
The difference is therefore the vertical distance (i.e., theTangentplane

to outliers and sharp features than these methods.

2. Related Work

Vertical distance
between points

Surfacenormal𝑧

Vertical distanceto the tangent plane(Residual)

Surface

Point-to-plane
distance(Residual)𝑦𝑥
Figure 2. (Left) Vertical distance, i.e., distance along the viewing
direction and (Right) point-to-plane distance, i.e., distance along
the plane normal (Ours).

distance along the viewing direction) between 3D points, as
illustrated in Fig. 2 left. DGP also explicitly computes thevertical distance in its local shaping step2 . These systems
then equate the vertical distance between points to the measured gradient. Geometrically, the residual therefore represents the vertical distance from the estimated point to thetangent plane at its neighboring point.
This paper introduces the perpendicular point-to-plane
distance (i.e., the distance along the plane normal direction) to the problem of normal integration, as illustrated
in Fig. 2 right. We formulate normal integration as aninverse problem of local plane fitting in the camera coordinates. To measure point-to-plane distances, we introduce plane displacements as unknowns, in addition to thepoints’ height values. From a geometric perspective, ourmethod simultaneously moves points on camera rays andplanes along their normal directions to minimize the sumof squared point-to-plane distances. Thanks to the formulation’s flexibility, it is straightforward to apply our method
to orthographic and perspective normal maps with arbitrary
boundaries, or to estimate height values on a different domain, such as the one used by DGP [37]. As a result, ourmethod avoids the checkerboard artifact and performs robustly against sharp features or outliers, as exemplified in
the bottom row of Fig. 1.
Our contributions are summarized as follows• We propose a new method for the problem of normal
integration by jointly optimizing point positions andlocal plane displacements in the camera coordinates.
To our knowledge, this is the first approach performing
normal integration based on point-to-plane distances.
• We show the source of the checkerboard artifact in
the central-difference discretized functional by evolving our plane fitting formulation.
• We point out that discrete Poisson’s equation and DGP
share the same intrinsic principle thus perform similarly, and demonstrate that our method is more robust
2 See Fig. 3 of [37].

Non-variational methods. Early works solve the problem of normal integration by path integration [10, 30, 35].
Another approach is to enforce surface integrability by projecting the gradient field onto the space spanned by a setof basis functions, e.g., Fourier basis functions [9], cosine
functions [11], orthonormal vectors [21], or shapelets [23].
Recently, Xie et al. [37] proposed a discrete geometry processing based normal integration method, which recovers a
discrete polygonal mesh directly from normal maps. Theirfollow-up works dealt with discontinuity preserving [36]and unconnected normal maps [33].
Section 5 demonstrates that the residual vector of DGP’s
linear system has the same geometric meaning as that of discrete Poisson’s equation; therefore, DGP and discrete Poisson’s equation perform similarly.
Variational methods. Pioneered by Ikeuchi [18] and detailed by Horn and Brooks [16], the problem of normal integration has been formulated in a variational manner by
optimizing a functional. Agrawal et al. [1] proposed a general framework to construct a range of solutions to the functional. Instead of minimizing the least squares functional,
different kinds of regularizations were studied to increase
robustness, like L1 norm [6, 27]. Depth measurements werealso used as prior knowledge of the surface [2, 17, 25].
How to better discretize the functional is also an essential issue in practice [3, 28]. Further, boundary conditions
should be carefully handled [8]. A typical method for thediscretization is to vectorize the height map and to approximate the partial derivative by finite difference. Zhu andSmith [38] used 2D Savitzky-Golay filters to approximate
partial derivative. Harker and O’Leary [12, 14] discretized
the functional into a Sylvester equation without vectorizing
the height map. Our method’s critical difference from thevariational methods is that our linear system’s residual vector measures the point-to-plane distance.
Two methods have been studied to reconstruct surfaces
from perspective normal maps. First, Durou et al. [7, 8, 28]
integrated the perspective normal map by taking the logarithm of height values. They then solved the discrete Poisson’s equation similar to the orthographic case, and exponentiate the result to recover the desired height map. Second, Nehab et al. [25], as well as Zhu and Smith [38], derived the relation between the perspective normal map andthe height map based on that the normal vector should be
perpendicular to the surface tangent vectors. To our knowledge, our inverse plane fitting formulation provides a third
alternative to solve the problem of normal integration in perspective case.

3. Problem Statement
We first briefly review the definition of and conventional
methods for the problem of normal integration. A normal map n : Ωn → S 2 ⊂ R3 records the surface orientation at a discrete pixel coordinate u = [u, v]> ∈ Ωn ,where the domain Ωn ⊂ {1, . . . , H} × {1, . . . , W } is defined by the image of dimension H × W . Let p(u) =[x(u), y(u), z(u)]> be the 3D surface point and n(u) =[nx (u), ny (u), nz (u)]> ∈ S 2 be the unit surface normal
vector corresponding to the pixel coordinate u. The problem of normal integration is to reconstruct the height mapz(·) from the normal map n(·).
Variational methods formulate normal integration asheight/shape from gradient, i.e., recovering a height map
from a measured gradient field [p, q]> . Under orthographic
projection,
p=−

nx
,
nz

q=−

ny
.
nz

(1)

We refer the reader to [28] for detailed derivation of [p, q]>
under different types of projections. Variational methods
then find the optimal surface by minimizing the following
functional
ZZ
J (z) =(∂u z − p)2 + (∂v z − q)2 du dv, (2)Ωnto the plane along the plane’s normal direction. We will calld the plane displacement throughout this paper. All points
p ∈ R3 on the plane should satisfy Eq. (4).
Suppose we are given a set of normal vectors {n(u)}over the pixel coordinates. For each normal vector n(u),there is a family of planes parallel to the tangent plane, parameterized by its displacement d(u) asp> n(u) + d(u) = 0.

As the surface is locally approximated by planes, we enforce nearby points of the point p(u) that correspond to
pixel u to lie on a plane asXmin(p> n(u) + d(u))2 ,(6)p∈N (p(u))

where N (p(u)) denotes the set of nearest neighbor points
to p(u). Equation (6) is exactly the same as PlaneSVD [22],a method for estimating normals from known 3D points byfitting planes. In the problem of normal integration, normal
vectors n are given, and the 3D points p are to be estimated,which is why we term our method as “inverse plane fitting”.
To recover the 3D points, we need to jointly consider planeequations at all points asXX(p> n(u) + d(u))2 .
(7)
min
p,d

where ∂u · and ∂v · denote partial derivatives of the functionz : Ωn → R along u and v axes on the image plane. Applying the Euler-Lagrange equation to the functional Eq. (2)leads to the necessary optimal condition∆z = div(p, q) = ∂u p + ∂v q,(3)

which is in the form of Poisson’s equation. Equations (2)and (3) consider the problem of normal integration on acontinuous domain. As the normal map is measured ona discrete domain, previous methods either discretize thefunctional [14] or Poisson’s equation [28].
We show in Section 5 that these conventional discretization methods minimize a residual by measuring the verticaldistance along viewing directions. On the other hand, ourmethod minimizes the point-to-plane distance along surface
normal directions.

4. Normal Integration via Inverse Plane Fittingp> n + d = 0,(4)

where n ∈ S 2 is the unit normal vector perpendicular tothe plane, and d is the distance from the coordinate originu∈Ωn p∈N (p(u))

To solve this optimization, we introduce two assumptions. First, we will assume that the 2D projections of all
3D points on the image plane are known. This constraint
reduces the degree of freedom of points from 3 to 1, i.e.,
all points can only move along known camera rays. We canthen parameterize a point p as p(z; up ), where up is the
known projection of the 3D point on the image plane, andz is the 3D point’s unknown position on the camera ray.
Second, we will assume that the closeness relationship between 3D points remains the same for their projections on
the image plane. That is,
N (p(u)) = {p(z; up ) | up ∈ N (u)}.

(8)

This assumption holds for continuous surfaces. With thetwo assumptions, the objective Eq. (7) becomes
X X
min(p(z; up )> n(u) + d(u))2 . (9)z,dOur main idea is to approximate the surface locally byplanes perpendicular to normal vectors. We consider theform of a plane equation in a 3D space as(5)u∈Ωn up ∈N (u)There is still one issue left for optimizing Eq. (9),namely, definition of the set of up , which reflects the points
of interest we will solve. We show two design options inthis paper. The first obvious choice is to let {up } = Ωn ,which is the same as conventional variational methods. The
second choice is to place up a half-a-pixel step from u in

CameraCoordinate

𝒑$"

𝒑"#

𝒏(𝒖𝒐𝒐 )

CameraCoordinate𝒏(𝒖𝒐𝒐 )

𝒑##

𝒑""
𝒑#"

𝒑"$
𝑑(𝒖𝒐𝒐 )𝒑#"𝒑""𝑑(𝒖𝒐𝒐 )

Viewing
direction𝒑"#Viewing
direction

𝒖𝒐#
𝒖𝒐𝒐

𝒖$𝒐

𝒖#𝒐

𝒖""
𝒖#"
𝒖𝒐𝒐
𝒖"#

𝒖𝒐$

𝒖##

Pixel Coordinates

Pixel Coordinates

Figure 3. (Left) A local five-point stencil in domain Ωn on theimage plane. We enforce all points projected onto this stencil tolie on a plane with the center normal vector and an unknown displacement. (Right) The domain Ωz used by DGP [37], denoted as
circles. Each pixel in Ωn has four neighbor pixels in Ωz .

Ωn , as illustrated in Fig. 3 right. The second one is the sameas that used by DGP [37], and we denote it as Ωz .
To retrieve the nearest neighbor projections of u, we useball query with radius 1. In the case of {up } = Ωn , thisquery results in a five-point stencil for pixels whose fourneighbor pixels are also in Ωn , as shown in Fig. 3 left. We
use a subscript “oo” to indicate the center pixel u, its corresponding 3D point p, and its height value z. Subscripts“-o”, “+o”, “o+”, “o-” are used for its left / right / upper /lower counterparts in the local five-point stencil. For boundary or corner pixels, N (u) has four or three elements, respectively. The inner summation in Eq. (9) now becomesX(p(z; up )> n(u) + d(u))2
up ∈N (u)(10)
= k[poo

po+

po-

p+o

p+-

p-+

nz 1

p- - ]>

n(u) + d(u)1k22 ,From now on, we will call five-point and four-point planefitting for the case {up } = Ωn and {up } = Ωz , respectively.

4.1. Solution Method
Orthographic case. Under a (scaled) orthographic projection or weak perspective projection, a surface pointp(z; up ) is projected onto the image plane with a scaling>
along xy-axes, i.e., p(z; up ) = [hu>p , z] . The scale factor
h transfers the camera’s xy-coordinates to the pixel coordinates; in the context of finite difference, h is called the step

d

u+o
u>
-o

Stacking equations at all pixels results in a sparse, overdetermined linear system A [ dz ] = b. The number
P of rows of
A or the number of all plane equations is u∈Ωn |N (u)|,which is close to 5|Ωn |; the number of columns of A is
2|Ωn |, as z and d are stacks of |Ωn | height values and plane
displacements, respectively. We can derive a least squares
approximate solution by minimizing the following residual
minkek22z,ds.t.

e = A [ dz ] − b.

(13)

Similarly, if we work on the domain Ωz , we can derive
the same plane equation as Eq. (12) but with four equations at each pixel. The coefficient matrix A’s dimension
is 4|Ωn | × (|Ωn | + |Ωz |), as there are |Ωn | displacements
and |Ωz | height values.
Perspective case. Suppose a perspective pinhole camera
with intrinsics K projects a surface point p(z; up ) onto theimage plane at up . Its back-projection is a camera ray containing p(u) with an unknown scale z asp(z; up ) = zK−1 ũp ≡ z p̃,p-o ]> n(u) + d(u)1k22 ,where 1 is an all-ones vector.
In the case of {up } = Ωz , N (u) of all pixels has four
elements, no matter whether it is a boundary pixel or not.
Similarly, the inner summation in Eq. (9) becomesX(p(z; up )> n(u) + d(u))2
up ∈N (u)(11)
= k[p++

size. Without loss of generality, we will assume h = 1 inthis paper for notational simplicity. Inserting the point representation into plane equation of Eq. (10) (inner part of thenorm) and moving known terms to the right-hand side, we
obtain the system of linear equations at one pixel as
 >
uoo

  zoo 
nz
1
zo+
u>

o+ 
nz
1
z
o- 
 nx


nz
1   z+o= −  u>(12)o-  [ ny ] .
nz
1
>
z-o

(14)

where ũp = [u, v, 1]> is up in homogeneous coordinates, p̃
is the back-projected point K−1 ũp . In the case of {up } =Ωn , inserting Eq. (14) into the plane equation of Eq. (10)yields five equations at one pixel as
 n> p̃1 oon> p̃o+




zoo
z
1
  zo+o- = 0.
n p̃o1  z+oz-on> p̃+o1dn> p̃-o 1>(15)

In the case of {up } = Ωz , there are four similar equations
at each pixel. Stacking all equations yields a sparse, homogeneous system A [ dz ] = 0. The dimension of A is thesame as that in the corresponding orthographic case. Then
our objective using least squares becomesminkek22z,ds.t.

e = A [ dz ] .

(16)

The non-trivial solution to Eq. (16) is the right-singular vector corresponding to the smallest singular value of A3 .
3 See, for example, A5.3 of [15] for a detailed proof.

Rank, ambiguity, and solver. The system matrix A isrank deficient in the orthographic case and has full-rank in
the perspective case. h In orthographic
case, nullity(A) =i(a) Eq. (12)(b) Eq. (18)(c) Eq. (19)(d) Eq. (20)(e) Eq. (17)11; the null vector is −n◦−1, where n◦−1
is Hadamardzz
inverse of the stacking vector of all nz . The null vector canbe verified in Eq. (12), and implies an offset ambiguity inthe solution. Geometrically understanding, if all points are
moved by 1 unit along the camera rays, the best fit planes
are moved by −nz−1 unit along planes’ normal directions.
We used LSQR [26] to solve Eq. (13) in this paper. Whenthe size of the system matrix A becomes large, solving thenormal equation of Eq. (13) by multigrid methods [4] canbe more efficient.
In perspective case, the solution is up to a scale ambiguity as we are solving a homogeneous system Eq. (16). The
system matrix A can be full-rank because the residual ofplane fitting never achieves 0 unless the surface is a plane.
This intuition is consistent with our observation that A is
only rank-deficient when the surface is a plane.
To efficiently find the solution to Eq. (16), we can compute the eigenvector corresponding to the smallest eigenvalue of A> A, as A’s right singular vectors are A> A’s
eigenvectors. We used ARPACK’s implementation [24] ofimplicitly restarted Lanczos method [32], which is efficientat finding a few eigenvectors of a large sparse matrix.

5. Robustness Analysis and Comparison
This section analyzes three previous orthographic normal integration methods in the language of our plane fitting formulation: Discrete functional [14], discrete Poisson’s equation [29], and DGP [37]. We begin by describing
the source of the checkerboard artifact that appears withcentral-difference-discretized functional. We then show theresidual vectors of the discrete Poisson’s equation and DGP
share the same geometric meaning.
Discrete functional. Central difference is typically used
to discretize the functional Eq. (2), yielding two equations
at one pixel location as(nx
1
2 (z+o − z-o ) = p = − nz
.
(17)
n
y
1
2 (zo+ − zo- ) = q = − nz
Solving the linear system based on Eq. (17) causes thecheckerboard artifact, as shown in Fig. 4 (e). One speculation about the source of the artifact is because Eq. (17)is independent of the center point zoo [38]. We argue that
it is in fact due to the loss of constraints between the left /right and upper / lower point pairs. And as long as we specify their relations, the checkerboard artifact can be avoided
even if there is no constraint using the center point.
To validate our argument, we gradually manipulate our
five-point plane fitting equations Eq. (12) to the discrete

high
high

0low

Figure 4. An evolution from plane fitting to central-difference discretized functional. Each column shows the estimated surface andthe absolute difference map. Central difference implicitly allows alocal offset in each five-point stencil (i.e., fit two parallel planes at
one pixel), introducing the checkerboard artifact.

functional Eq. (17). We first discard the plane constraint
on the center point poo in Eq. (12), yielding four equations
at one pixel as(u − 1)nx + vny + z-o nz + d = 0(left)(u + 1)n + vn + z n + d = 0(right)xy+o z. (18)

un
+
(v
+
1)n
+
z
n
+
d
=
0
(upper)
x
y
o+
z



unx + (v − 1)ny + zo- nz + d = 0(lower)Equation (18) enforces four points on one plane, but now
let us relax the constraint as(left)
(u − 1)nx + vny + z-o nz + d = 0(u + 1)n + vn + z n + d = 0(right)xy+o z, (19)
0unx + (v + 1)ny + zo+ nz + d = 0(upper)unx + (v − 1)ny + zo- nz + d0 = 0(lower)which means we fit two parallel planes for the four points.
In other words, there could be a local offset (i.e., d 6= d0 ) between the left / right and upper / lower point pairs. However,
we can derive the same equations from Eqs. (18) and (19) if
we take the difference between the left and right, upper andlower point 4(nz
2 (z+o − z-o ) = −nx
.
(20)
nz
2 (zo+ − zo- ) = −ny
Dividing both sides by nz , we then arrive at the discrete
functional Eq. (17).
Figure 4 displays the estimated surfaces by the linear systems based on the equations in this evolution.
From Fig. 4 (a) to (b), we discard the constraint on thecenter point. But it does not create the checkerboard artifact, and the reason is that the four neighbor points arestill enforced on the same plane. From Fig. 4 (b) to(c), the checkerboard artifact occurs, implying that relaxing the constraint from Eq. (18) to Eq. (19) introduces the4 This formulation is the same as Zhu and Smith’s method [38] in theorthographic case with central-difference-based discretization.

checkerboard artifact. Further, the checkerboard artifact between Fig. 4 (c) and (d) is indistinguishable, validating that
central difference implicitly allows a local offset between
the left / right and upper / lower point pairs in each five-pointstencil. The checkerboard artifact in Fig. 4 (d) is weaker
than that in (e) because Eq. (20) transfers the vertical distance to the distance along plane normal direction (multiplied by nz ). This phenomenon further validates that measuring distance along normal directions (Eq. (20)) is better
than that along viewing directions (Eq. (17)).
Discrete Poisson. Poisson’s equation (3) is typically discretized at a five-point stencil [1, 8, 28] as
z+o + z-o + zo+ + zo- − 4zoo =qo+ − qop+o − p-o
+
.
2
2
(21)

Quéau et al. [29] point out it is equivalent to consider all
forward / backward differences at one pixel((z+o − zoo = poozo+ − zoo = qoo.
(22)zoo − z-o = poozoo − zo- = qoo
The least squares approximate solutions of two linear systems constructed by stacking Eq. (21) or (22) are the same5 .
We will use Eq. (22) to demonstrate the geometric meaningof the residual of discrete Poisson’s equation.
Consider a center and its right-hand point lying on thetangent plane of the center point(
unx + vny + ẑoo nz + d = 0(center). (23)(u + 1)nx + vny + ẑ+o nz + d = 0(right)We denote by ẑ to stress they are the height values of theintersections of camera rays and the tangent plane, not theones by solving the linear system. Taking the difference
between the two equations yields
ẑ+o − ẑoo = −nx= poo ,nz

(24)

which reveals that the gradient poo is equal to the verticaldistance between two neighbor intersections. We can derive
similar equivalence for other pixel pairs in the five-pointstencil and plug them into Eq. (22) as((
z+o − zoo = ẑ+o − ẑoo
zo+ − zoo = ẑo+ − ẑoo.
zoo − z-o = ẑoo − ẑ-ozoo − zo- = ẑoo − ẑo(25)These equations imply that discrete Poisson’s equation enforces the equality between the vertical distance of estimated neighbor points and neighbor tangent intersections.

natural boundary

1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0

outlier

1.1
1.0
0.9
Poisson 1D
Line fitting

0.8
-1
-0.5
0
0.5
1
-0.5
0
0.5
Figure 5. Poisson equation v.s. line fitting on a 2D semi-circle.
(Left) The boundary point has a large tangent slope (dashed line),deviating estimated points to the tangent along viewing direction.
(Right) We rotated one normal vector as an outlier, whose tangent
line is shown as the dashed line. Poisson’s equation performs moresensitive to the outlier than our line fitting.

The residuals are then the vertical distances from the estimated point to the tangent at its neighbor points.
To demonstrate why this residual is sensitive to natural
boundary or outliers, we compared Poisson’s equation andplane fitting on a 1D function in Fig. 5. Plane fitting degrades to line fitting with the line equation√ xnx +yny +d =0. We used a circle function y = 1 − x2 for x ∈[−1, 1] as the ground truth surface. In Fig. 5 left, we
sampled 9 normal vectors at points evenly distributed in[−0.98, 0.98] and estimate the height values of the 9 points.
The two boundary points are close to the surface’s natural
boundary, so the corresponding tangent has a large slope.
The points estimated by discrete Poisson’s equation deviate largely toward the tangent along viewing direction. The
points estimated by line fitting, on the other hand, are close
to the surface since the point-to-plane distance would notchange largely by moving them on the camera rays.
In Fig. 5 right, we evenly sampled 9 points in [−0.5, 0.5],
i.e., there is no effect of the natural boundary. However,
we rotated one normal vector as an outlier, whose tangent
line is shown as the dashed line. Due to the large verticaldistance introduced by the outlier normal vector, the points
estimated by discrete Poisson’s equation deviate more towards the tangent than line fitting.

Discrete Geometry Processing (DGP). We will showthat the residual of DGP [37] is essentially the same as thatof discrete Poisson’s equation, but on domain Ωz . DGP consists of two steps: (1) local shaping and (2) global blending.
Local shaping defines a local 3D coordinate centered at uoo
on the image plane, and then computes the vertical distancefrom a neighbor pixel to the plane with normal vector noo
passing through uoo (i.e., d = 0). For example, the vertical
5 Equation (21) should be modified for boundary pixels [3].

(local)distance ẑ++of the upper-right pixel u++ is computed as

GT shape &normal map

Poisson [29]Oursfive-point

DGP [37]

Ours
four-point(local)0.5nx + 0.5ny + nz ẑ++
=0
0.5n
+
0.5n
x
y
(local)→ ẑ++
.
=−
nz

(26)
high
high

Global blending then equates the difference between a point
and the average of its four neighbor points to the local vertical distance asz++ −z++ + z+- + z-+ + z- (local)= ẑ++ .
4

5.30 / 3.97

0.86 / 0.69

5.27 / 3.94

0.76 / 0.61

0low
×10−2

(27)

The vertical distance in the local coordinates is the sameas the vertical distance in the camera coordinates, i.e.,(local)ẑ++= ẑ++ − ẑoo . Let zoo be the average of four neighborheight values, then Eq. (27) can be simplified and all four
equations at one pixel location now read((z++ − zoo = ẑ++ − ẑoo
z-+ − zoo = ẑ-+ − ẑoo.
z+- − zoo = ẑ+- − ẑoo
z- - − zoo = ẑ- - − ẑoo
(28)This form is similar to discrete Poisson’s equationin Eq. (25), except that the height values are defined on thedomain Ωz . Therefore the residual of DGP has the same
geometric meaning as that of discrete Poisson’s equation.
Figure 6 compares discrete Poisson’s equation, DGP,
and our methods on two analytically computed normal
maps; see our supplementary material for the strategy
of computing the surface height/normal maps and adding
noise/outliers. As there is a global offset ambiguity in theestimated height maps, we first found an optimal offset between the estimated and the ground truth height maps. We
then computed the Root Mean Squared Error (RMSE) andMean Absolute Error (MAE) as evaluation metrics. Discrete Poisson’s equation and DGP show an indistinguishable error pattern and achieve close RMSE / MAE. Thisphenomenon validates our understanding that the residual
vectors of DGP and discrete Poisson’s equation have thesame geometric meaning. On the other hand, our method
performs more robustly to outliers and sharp features.

high
high

8.32 / 5.21

3.56 / 1.41

8.24 / 5.21

2.79 / 0.64

0low
×10−1

Figure 6. Reconstructed shape and its absolute difference map
from (Top) the normal map with Gaussian noise and 1% outliers ofanisotropic Gaussian surface [14] and (Bottom) a noise-free normal map of sphere with background padded. Discrete Poisson’sequation and DGP shows an indistinguishable error pattern andclose RMSE / MAE (numbers below the error maps). The shapes
from discrete Poisson’s equation or DGP are distorted by outliers,
and Gibbs phenomenon appears near sharp features.

Figure 7 compares state-of-the-art perspective normal integration methods on an analytically computed, synthetic,
and real-world normal map. We computed the ground truth
normal and height map of a S PHERE by a perspective camera model. We used Mitsuba6 to render a normal map andthe corresponding height map of Stanford B UNNY. We thenadded Gaussian noise to and randomly selected 1% normal
vectors as outliers in S PHERE and B UNNY’s normal map.

For the real-world surface, we used the normal map estimated by the state-of-the-art photometric stereo method [5]on the object H ARVEST in DiLiGenT benchmark [31]. See
our supplementary material for the comparison on all objects in DiLiGenT.
To better reflect the linear systems’ features, we only
looked at the least squares approximate solutions of compared methods without any regularization or smoothness
term. As there is a global scale ambiguity in the estimated height maps, we found an optimal scale between
the estimated and ground truth height map before computing RMSE and MAE. We applied the exponentiation [28]to the discrete functional [14] and discrete Poisson’s equation [29]. It is observed that the exponentiation amplifies
the checkerboard artifact of discrete functional and sensitivity to outliers of discrete Poisson’s equation. Spikes occur in both methods. For Zhu and Smith’s method [38], weused two variants as the discretization of partial derivative:
Central difference and the third-order Savitzky-Golay (SG)filter.Zhu and Smith’s method [38] performs robustly to outliers. However, a weaker checkerboard artifact exists in thecentral difference case. Moreover, a higher-order SG filter
does not refine the reconstruction quality compared to central difference. Instead, it introduces severer high-frequency
noise7 . On the other hand, our methods perform stably on

6 Mitsuba Renderer. https://www.mitsuba-renderer.org/
index_old.html, last accessed on November 16, 2020.

7 We empirically found it is the SG smoothing filter that plays a crucialrole in refining the visual quality in Zhu and Smith’s method [38].

6. Comparison of Perspective Methods

GT shape &noisy normal map

Harker et al. [14]discrete functionalQuéau et al. [29]discrete Poisson

Zhu and Smith [38]w/ central difference

Zhu and Smith [38]w/ SG filter

Ours
five-point on Ωn

Ours
four-point on Ωz

high
high

11.56 / 6.07

4.24 / 1.48

0.54 / 0.38

0.35 / 0.20

1.67 / 1.33

0.39 / 0.19

0low
×10−2

high
high

6.02 / 3.65

4.66 / 2.77

3.77 / 1.49

3.77 / 1.57

3.92 / 1.58

3.75 / 1.49

0low
×10−3

high
high

13.38 / 11.24

13.10 / 11.14

12.90 / 11.11

12.91 / 11.11

12.92 / 11.13

12.91 / 11.11

0low
×1

Figure 7. Visual comparison of perspective normal integration methods on normal maps of analytically computed S PHERE, synthetic
B UNNY, and real-world H ARVEST. Every two rows show the ground truth & estimated surface, the normal map as well as the absolute difference map (between ground truth and estimated height map). The two numbers underneath error maps are RMSE and MAE, respectively.
We add Gaussian noise to and randomly select 1% normal vectors as outliers in S PHERE and B UNNY’s normal maps.

all objects. The effect of outliers is hardly observed, andsystematic high-frequency noise never occurs.

7. Conclusion
In this paper we presented, analyzed, and evaluated aninverse plane fitting approach to the problem of normal integration. We formed a linear system that minimizes the sumof squared point-to-plane distances instead of conventionalvertical distances. We showed the source of the checkerboard artifact in the discrete functional-based approach. We
also provided a geometric interpretation of the sensitivity
to outliers or the sharp features of discrete Poisson’s equa-tion and DGP. Our method is straightforward, robust against
outliers and sharp features, hyper-parameter free, yet basedon least squares. We are interested in further studying discontinuity preserving normal integration and depth normal
fusion based on our formulation.

8. Acknowledgments
This work was supported by JSPS KAKENHI Grant
Number JP19H01123, National Natural Science Foundation of China under Grant No. 61872012, 62088102, andBeijing Academy of Artificial Intelligence (BAAI). We
would like to thank Zhuoyu Yang for fruitful discussions.

References[1] Amit Agrawal, Ramesh Raskar, and Rama Chellappa. Whatis the range of surface reconstructions from a gradient
field? In Proc. of European Conference on Computer Vision (ECCV), 2006.
[2] Doris Antensteiner, Svorad Stolc, and Thomas Pock. A review of depth and normal fusion algorithms. Sensors, 2018.
[3] Martin Bahr, Michael Breus, Yvain Quéau, Ali Sharifi
Boroujerdi, and Jean-Denis Durou. Fast and accurate surface normal integration on non-rectangular domains. Computational Visual Media, 2017.
[4] William L Briggs, Van Emden Henson, and Steve F McCormick. A multigrid tutorial. SIAM, 2000.
[5] Guanying Chen, Michael Waechter, Boxin Shi, Kwan-Yee K
Wong, and Yasuyuki Matsushita. What is learned in deep
uncalibrated photometric stereo? In Proc. of European Conference on Computer Vision (ECCV), 2020.
[6] Zhouyu Du, Antonio Robles-Kelly, and Fangfang Lu. Robust surface reconstruction from gradient field using the L1
norm. In Proc. of Biennial Conference of the Australian Pattern Recognition Society on Digital Image Computing Techniques and Applications (DICTA), 2007.
[7] Jean-Denis Durou, Jean-François Aujol, and Frédéric
Courteille. Integrating the normal field of a surface in thepresence of discontinuities. In Workshops of the Computer
Vision and Pattern Recognition (CVPRW), 2009.
[8] Jean-Denis Durou and Frederic Courteille. Integration of anormal field without boundary condition. In Workshops of
the International Conference on Computer Vision (ICCVW),
2007.
[9] Robert T. Frankot and Rama Chellappa. A method for enforcing integrability in shape from shading algorithms. IEEE
Transactions on Pattern Analysis and Machine Intelligence(PAMI), 1988.
[10] Silvano Galliani, Michael Breus, and Yong Chul Ju. Fast androbust surface normal integration by a discrete eikonal equation. In Proc. of British Machine Vision Conference (BMVC),
2012.
[11] Athinodoros S. Georghiades, Peter N. Belhumeur, andDavid J. Kriegman. From few to many: Illumination conemodels for face recognition under variable lighting and pose.
IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), 2001.
[12] Matthew Harker and Paul O’Leary. Least squares surface
reconstruction from measured gradient fields. In Proc. ofComputer Vision and Pattern Recognition (CVPR), 2008.
[13] Matthew Harker and Paul O’Leary. Least squares surface reconstruction from gradients: Direct algebraic methods with
spectral, Tikhonov, and constrained regularization. In Proc.
of Computer Vision and Pattern Recognition (CVPR), 2011.
[14] Matthew Harker and Paul O’leary. Regularized reconstruction of a surface from its measured gradient field. Journal of
Mathematical Imaging and Vision, 2015.
[15] Richard Hartley and Andrew Zisserman. Multiple view geometry in computer vision. Cambridge university press,
2003.

[16] Berthold KP Horn and Michael J Brooks. The variational
approach to shape from shading. Comput. Vis. Graph. Image
Process., 1986.
[17] Itsik Horovitz and Nahum Kiryati. Depth from gradient
fields and control points: Bias correction in photometricstereo. Image and Vision Computing, 2004.
[18] Katsushi Ikeuchi. Constructing a depth map from images.
Technical report, Massachusetts Inst of Tech Cambridge Artificial Intelligence Lab, 1983.
[19] Katsushi Ikeuchi and Berthold KP Horn. Numerical shapefrom shading and occluding boundaries. Artificial intelligence, 1981.
[20] Achuta Kadambi, Vage Taamazyan, Boxin Shi, and Ramesh
Raskar. Polarized 3D: High-quality depth sensing with polarization cues. In Proc. of International Conference on Computer Vision (ICCV), 2015.
[21] Bilge Karacali and Wesley Snyder. Reconstructing discontinuous surfaces from a given gradient field using partial integrability. Computer Vision and Image Understanding, 2003.
[22] Klaas Klasing, Daniel Althoff, Dirk Wollherr, and Martin
Buss. Comparison of surface normal estimation methods forrange sensing applications. In Proc. of International Conference on Robotics and Automation (ICRA)), 2009.
[23] Peter Kovesi. Shapelets correlated with surface normals produce surfaces. In Proc. of International Conference on Computer Vision (ICCV), 2005.
[24] Richard B Lehoucq, Danny C Sorensen, and Chao Yang.
ARPACK users’ guide: solution of large-scale eigenvalue
problems with implicitly restarted Arnoldi methods. SIAM,
1998.
[25] Diego Nehab, Szymon Rusinkiewicz, James Davis, and Ravi
Ramamoorthi. Efficiently combining positions and normals
for precise 3D geometry. ACM Transactions on Graphics(Proc. of ACM SIGGRAPH), 2005.
[26] Christopher C Paige and Michael A Saunders. LSQR:
An algorithm for sparse linear equations and sparse leastsquares. ACM Transactions on Mathematical Software(TOMS), 1982.
[27] Yvain Quéau and Jean-Denis Durou. Edge-preserving integration of a normal field: Weighted least-squares, TV and L1
approaches. In Proc. of International Conference on Scale
Space and Variational Methods in Computer Vision, 2015.
[28] Yvain Quéau, Jean-Denis Durou, and Jean-François Aujol.
Normal integration: a survey. Journal of Mathematical
Imaging and Vision, 2018.
[29] Yvain Quéau, Jean-Denis Durou, and Jean-François Aujol.
Variational methods for normal integration. Journal of Mathematical Imaging and Vision, 2018.
[30] Antonio Robles-Kelly and Edwin R Hancock. A graphspectral method for surface height recovery from needlemaps. In Proc. of Computer Vision and Pattern Recognition(CVPR), 2001.
[31] Boxin Shi, Zhipeng Mo, Zhe Wu, Dinglong Duan, Sai-Kit
Yeung, and Ping Tan. A benchmark dataset and evaluation for non-lambertian and uncalibrated photometric stereo.
IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), 2019.

[32] Danny C Sorensen. Implicit application of polynomial filtersin ak-step arnoldi method. Siam journal on matrix analysis
and applications, 1992.
[33] Miaohui Wang, Wuyuan Xie, and Maolin Cui. Surface reconstruction with unconnected normal maps: An efficientmesh-based approach. In Proc. of ACM International Conference on Multimedia (ACMMM), 2020.
[34] Robert J Woodham. Photometric stereo: A reflectance map
technique for determining surface orientation from image intensity. In Proc. of Image Understanding Systems and Industrial Applications I, 1979.
[35] Zhongquan Wu and Lingxiao Li. A line-integration based
method for depth recovery from surface normals. Computer
Vision, Graphics, and Image Processing, 1988.
[36] Wuyuan Xie, Miaohui Wang, Mingqiang Wei, Jianmin Jiang,
and Jing Qin. Surface reconstruction from normals: A robust
dgp-based discontinuity preservation approach. In Proc. ofComputer Vision and Pattern Recognition (CVPR), 2019.
[37] Wuyuan Xie, Yunbo Zhang, Charlie CL Wang, and Ronald
C-K Chung. Surface-from-gradients: An approach based on
discrete geometry processing. In Proc. of Computer Vision
and Pattern Recognition (CVPR), 2014.
[38] Dizhong Zhu and William A. P. Smith. Least squares surface
reconstruction on arbitrary domains. In Proc. of European
Conference on Computer Vision (ECCV), 2020.