{'fraction_non_alphanumeric': 0.04154290256453165, 'fraction_numerical': 0.010709952725599298, 'mean_word_length': 4.04516673701984, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
What's In A Patch, II: Visualizing generic surfaces

We continue the development of a linear algebraic framework for the shape-from-shading problem, exploiting the manner in which tensors arise when scalar (e.g. image) and vector (e.g. surface normal) fields are differentiated multiple times. In this paper we apply that framework to develop Taylor expansions of the normal field and build a boot-strapping algorithm to find these polynomial surface solutions (under any light source) consistent with a given patch to arbitrary order. A generic constraint on the image derivatives restricts these solutions to a 2-D subspace, plus an unknown rotation matrix. The parameters for the subspace and rotation matrix encapsulate the ambiguity in the shading problem.


Introduction
We have been developing a linear-algebraic approach to the shape-fromshading problem [7], and this paper is a continuation of that effort. 1 The key intuition underlying the research effort is that, if one were to 'drill down' in derivatives for the surface, then this should correspond to analagous derivatives for the image (Fig. 1). We earlier considered which image (derivatives) are most likely given certain shape (normal) derivatives, for similiar levels of differentiation. We now combine the derivatives in the Taylor sense, to calculate a representation of the full space of possible surface patches that could correspond to a given image patch, up to some order of differentiation. The algorithmic formulation allows us to 'bootstrap' the next-order structure from the previously calculated structure, similar to the use of recursion in series solutions to ordinary differential equations. It is important because, as shown in Fig. 3, the differential order of the image and of the surface must somehow be coupled. For convenience, we work in ambient Euclidean space; for a related abstract treatment using covariant derivatives, see [9].
In general there are infinite surfaces that can correspond to a given image patch, given the ill-posedness of the problem. The contribution here is that we are able to 'control' it with linear algebraic machinery [7] plus a generic lighting constraint [5]. Importantly, the notion of control that emerges is algorithmic, and indicates that, while the generic assumptions can structure solutions, even when they are linked to a particular (differential) level ambiguities remain. In particular, for a smooth Lambertian image patch (with no additional information regarding the light source or boundary), we develop a characterization of the set of possible underlying surface patches. Unlike previous approaches, we do not restrict the underlying surface patch to be a polynomial of a fixed degree, but rather allow the surface patch to account for all intensity variation in the image patch. Working with Taylor approximations, we let the order be dependent on each image patch and introduce an algorithm that realizes the set of generic surface patches corresponding to that image patch. In effect this provides a visualization of the ambiguity in the problem well beyond bas-relief [2]; see Fig. 1.
An illustration of the family of solutions is shown in Fig. 2, for a Lambertian image patch of a cylinder. This is, of course, a very special object Figure 1: By Taylor's theorem, increased derivatives at a point provide additional information within a neighborhood. We apply this idea to the shape inference problem, by noting that derivatives of the image D j I can be related to derivatives on the surface D j N , although the map is not one-to-one. In this paper we characterize the ambiguity in possible surface reconstructions even with a generic lighting assumption. The results further illustrate the rich interactions between linear algebra and differential geometry. in which the surface normal variation is restricted to the radial direction and the curvature forms are low rank; see discussion in [7]. But it is also an important object in shape-from-shading research. Algorithmically it has been invoked to motivate a mean-curvature prior [1] and used to estimate human 'reflectance functions' [13]. Importantly, even in this special case, there is enormous variation in the perceived shape [14,11,10] and it plays a key role in light-source identification algorithms (e.g., [12]). Nevertheless, even among these many possibilities, neither the the bas-relief nor the generic lighting parameters, are invoked, so the variation is truly impressive. These concepts are developed in the course of this paper; we end with the algorithm that was used to compute these surfaces.

Background
A background review on shape-from-shading algorithms is provided in [17] and in the companion paper [7]. Here we concentrate on the few papers that are explicitly based on a patch model. For this approach the idea is to solve for local patches individually and then "stitch" them together [3,15,8].
In [4], the authors formulate the problem as solving a (large) system of polynomials using modern homotopy solvers over a triangulation. This is feasible for small images, involves (up to) quartic interactions, and leads to exact recovery of all possible solutions. For general patches, one can model the associated pixel values with various degrees of underlying surface complexity, represented as a Taylor polynomial in either heights or normals. [15] assume the image patch derives from a second-order surface and, therefore seek a quadratic solution; in [8] the image patch is modeled from a third-order surface. Assuming a local surface patch is exactly modeled by a quadratic, there are in general only four solutions to the local image formation model, i.e., the coefficients of the quadratic and the light source. If the image patch is large enough (i.e., number of pixels) relative to the number of coefficients of the Taylor polynomial, then the local patch can be determined up to a four-fold ambiguity [15,8]. Clearly over-fitting can be a problem if the image patch is taken to be too large; e.g., errors will arise in fitting a quadratic surface to an image patch that arose from a quartic surface. In general, fixing the underlying surface complexity while considering successively larger image patches creates overfitting, whereas fixing the constraints (image patch size) while increasing the degree of the Taylor polynomial leads to increasing am- Figure 2: A set of possible solutions to a cylindrical image patch. The rows depict variation in the generic solution space; the columns correspond to the bas-relief ambiguity (e.g. changing the attitude of the tangent plane and light source). The richness of these possible solutions, given that we will show they all are 'generic,' is well beyond what is normally expected [13]. The Taylor approximating polynomialsĪ andN have degree 5. Generic parameters (c 1 , c 2 ) chosen by c 1 = c 2 varying linearly from 0 to 1.  biguity. See Figure 3. The variability in possible solutions is also biologically relevant [16], where researchers have attempted to identify the response to an image for a given surface. The variability that we describe must be kept in mind when attempting to assess the uniqueness (or lack thereof) in such neural responses.

Conditions on Generic Normal Fields Associated with a Lambertian Image Patch
Branching off from the companion paper [7], we now want to determine the set of possible underlying generic surface patches corresponding to a given image patch. We assume orthogonal projection, so the surface S can be thought of as height field over the image plane. For more background on the setup and relevant notation, see [7]. The image patch I(x, y) centered at the origin is modeled by a Taylor approximationĪ(x, y): x 2 D 2 p I(e 1 , e 1 ) + xy D 2 p I(e 1 , e 2 ) + y 2 D 2 p I(e 2 , e 2 ) + third and higher order terms The e i represent the standard basis vectors. As in [7], we use the notation D j p I(·, ·, . . . , ·) to represent the j th derivative of the function I(x, y) at the point p. It is a multilinear j−form that requires j vector inputs to return a scalar value in R (See Appendix A in [7]). Thus, D j p I : R 2 j → R for each j. Our goal is to understand the under-determined map fromĪ(x, y) to a Taylor approximation of the surface normal field,N (x, y).N (x, y) is a multivariate polynomial with coefficients D j p N (·, ·, . . . , ·) : R 2 j → R 3 . Relating the two Taylor approximations can be understood by relating the coefficients. Thus, we seek an algorithm that takes the known values {D j p I} n j=1 as inputs and and outputs the {D j p N } n j=1 . Building on the earlier results, this one to many map carries the ambiguity when going from the Taylor approximation of the image to the Taylor approximation of the surface. For notational simplicity, we will drop the p subscript for the rest of the analysis.
A subtlety arises because the normal vector is unit length. First, represent the surface as a height function over the image plane: F (x, y) = {x, y, f (x, y)}. Using subscripts to denote partial differentation, the associated normal field is

Image intensity gives a projection of the normal field
A Lambertian image intensity is given by I(x, y) = αL · N (x, y). Assuming constant albedo, we set α = 1. We apply derivative operators j times to both sides to obtain: To appreciate D j N , think of e.g., D 2 N : this requires two R 2 vectors {v 1 , v 2 } as input and outputs a vector in R 3 . This output is the "change in the change in the normal field" as we differentiate first in the v 1 direction and then in the v 2 direction. Note that, due to the symmetry of the derivative, the order of the inputs does not matter.
The square root term in the denominator creates difficulties in relating the {D j I} n j=1 to the parameters {f x , f y , f xx , . . .} (or other surface parameters). According to the above equation 3, we see that the relationships between {D j I} n j=1 and {D j N } n j=1 is a projection along the (unknown) light source vector L. We will need two more linearly independent projections in order to uniquely define the remainder of {D j N } n j=1 and thus recover the normal fieldN (x, y). As we now show, one of these additional projections will be set by the unit length condition ofN (x, y). The final projection can be freely set and represents the ambiguity in the problem.

Normalization constraints yield another projection of the normal field
The normalization constraint can be expanded as a system of linear constraints in the Taylor series. By enforcing this linear system of constraints, we can ensure an approximately unit length normal field (up to error O(x n+1 )) in the following manner.

8
The normalization constraint is: Here, we write ·, · as the standard dot product in R 3 and we write N 0 = N (p) as the normal vector at the center of our patch. The above equation can be differentiated in an arbitrary image vector direction u and evaluated at p: Differentiate in another direction v: We do not choose {u, v} before we differentiate; we could keep these directions unknown and general. That is, we consider the D j N as a (1, j) tensor -a linear machine seeking j vectors and outputting a vector in R 3 . We create D j N , D k N as a new (0, j + k) tensor in the following way: · Construct D j N ⊗ D k N as a (2, j + k) tensor. The two contravariant parts correspond to the R 3 vectors D j N , D k N once j + k inputs have been chosen.
· Lower an index associated with the unique contravariant part (in R 3 ) of D k N to get a (1, j + k + 1) tensor.
· Contract the two indices associated with the R 3 parts (we now have one covariant and one contravariant) to perform the dot product. Thus, where C is the contraction operator and lowers the appropriate index. Examining the Equations 5, 6, a pattern emerges: if {D j N } j<m were known, then we could calculate D m N , N 0 . This key point will allow us to solve for D m N , N 0 for each m inductively and thereby gain knowledge of the projection of the D j N coefficients onto the central normal N 0 . Continuing to take derivatives and rearranging, we get the following proposition.
Proposition. The constraint N , N = 1 can be Taylor approximated up to order k by enforcing series of linear constraints, for every j ≤ k.
Here Π is the set of combinations of a objects chosen from j objects. These combinations arise from the application of the product rule multiple times. We let V represent the space R 2 j of j 2D image vector inputs into the tensor D j N and then π a (V ) represents a subset of a inputs out of the j possible ones. Thus, π a (V ) C represents the remaining j − a inputs. Note that since the order of the inputs doesn't matter (but whether they get fed to D a N or D j−a N does), we use combinations.
In conclusion, if {D j N } j<m were known, we could acquire the projections of D m N onto the vector N 0 . As D m N is fully defined when its projection onto three linearly independent vectors is known, it remains to search for one more projection. Unfortunately, there is no other information in the shape from shading problem that allows us to directly set a third projection of the D j N coefficients. This inherent ambiguity in the problem is due to the normal field being a higher dimensional entity (taking values on S 2 ) than the intensity function (taking values on R).

Using generic lighting to obtain a third projection
We now introduce a device that will allow us to calculate a third projection of D j N . Let G(x, y) : Ω → R be any smooth function and let b be any direction in R 3 not in the span of {L, N 0 }. Suppose we choose to set b T D j N = D j G, ∀j ∈ {1, . . . , n}. Provided we ensure L T D j N = D j I, ∀j ∈ {1, . . . , n} and equation 8 holds for each j, we will construct a Taylor seriesN (x, y) that is approximately unit length and approximately matches the image. (For a discussion of these Taylor remainder errors, please see the Appendix.) Integrating this normal field would provide a surface patch matching the original image patch for any G; however, some choices of G may be better (i.e. more robust and probable) than others.
A natural choice for G comes from the generic framework and notation developed in [5,6], which we now develop for our case. Define Here, we have unfolded (see Appendix in [7]) the various tensors D j N into 3 × 2 j matrices and then appended them together. Let m = n i=1 2 i . Then, β is a linear map from R m to R 3 and Y is a linear map from R m to R . Now, consider the following rendering function: Following [5,6], we call L the generic variable, β the scene parameters, and Y the observations. L is an unknown vector in R 3 . Following [5,6], we assume a Gaussian noise model on the observations Y : whereŶ is the ideal rendered observation and T ∼ N (0, Σ) with Σ = diag(σ 2 ) for some σ ∈ R + . For the noise model, we have Appling Bayes' theorem and integrating over the generic variable L yields the posterior distribution: = k (fidelity) (prior probability) (genericity) (16) where P β (β), P L (L 0 ) are prior distributions on the surface and light source parameters, L 0 is the light source that can best account for the observations given a chosen β and A is a matrix with the following elements: We now seek the solutions that maximize the posterior probability P (β|Y ). From equation 15, we maximize by choosing β and L 0 so that ||Y −f (L 0 , β)|| = 0 while at the same time setting det(A) = 0.
Where, ||Y − f (L 0 , β)|| = 0, the condition that det(A) = 0 is equivalent to the constraint that β is a low rank 3 × m matrix. Under this condition, β is determined (up to two constants c 1 , c 2 ) by its projection onto two linearly independent vectors. (We ignore the rank 1 case, as it's infinitesimally unlikely compared to the rank 2 case.) β's projection onto two linearly independent vectors can already be obtained, as the components of β are each D j N . Thus, if we restrictN (x, y) to the generic solutions, we can solve for it unambiguously up to the unknowns {c 1 , c 2 , L 0 , N 0 }. We now show this.

Representing unknown lighting and tangent plane orientation via change of basis
As the projections of D j N are onto the vectors {L, N 0 }, we will work in a basis defined by those vectors. Define l t to be the unit length projection of L 0 onto the tangent plane perpendicular to N 0 . That is, Note that P is an unknown orthogonal matrix, since we don't know either the normal or the direction of the light source. However, rather than computing the Taylor surfaceN (x, y) in the standard R 3 basis, we will instead compute the modified Taylor surface P TN (x, y). This is merely considering the output ofN (x, y) in a different frame. In this fashion, we solve for a family of surfaces that will all match the Taylor image polynomial. To obtain a single member of that family, we choose an element Q ∈ SO 3 (R) and Figure 4: A factorization to separate the unknown portion of the surface from its geometry. We precompose P T with the unfolded D j N tensors of size 3 × 2 ⊗j . This is equivalent to orthogonally changing the basis of the output space.
multiply to get Q (P TN (x, y)). This is equivalent to choosing a normal and light source for the scene. Thus our new goal is to solve for P TN (x, y) by solving for the coefficients of the Taylor series P T D j N , 1 ≤ j ≤ n; see Fig. 4. We do this in an inductive manner, as we will need the {P T D j N } j<k in order to solve for P T D k N .

Algorithm
In this section, we will put the pieces described in the above sections together in order to create an inductive algorithm that can solve for all generic surfaces corresponding to a single image patch. To do this efficiently, we use an unfolding of the tensors D j I and D j N .
We recall that D j N is a (1, j) tensor, that is, a multilinear map from R 2 j to R 3 . It has a matrix representation -a rank 1 unfolding -whose dimensions are is 3 × 2 j . To calculate the action of this tensor on our inputs {v 1 , v 2 , . . . , v j }, v i ∈ R 2 , we apply its matrix representation to the Kronecker product of the inputs w = v 1 ⊗ v 2 . . . ⊗ v j . We seek these matrix representations. Let r j i stand for row i of P T D j N .  Figure 2, we show many solution surfaces to a given image. Across columns, we change the rotation matrix P , which amounts to changing the light source and central normal (e.g. bas-relief ambiguity). The Taylor approximating polynomialsĪ andN have degree 5. Across rows, we change the values of the generic parameters c 1 , c 2 , chosen by c 1 = c 2 varying linearly from −1 to 1. However, all solutions are considered 'generic" according to Freeman's definition [5]. 14 Figure 6: Similar to Figure 2, we show many solution surfaces to a given image. Across columns, we change the rotation matrix P , which amounts to changing the light source and central normal (e.g. bas-relief ambiguity). The Taylor approximating polynomialsĪ andN have degree 5. Across rows, we change the values of the generic parameters c 1 , c 2 , chosen by c 1 = c 2 varying linearly from −1 to 1. However, all solutions are considered 'generic" according to Freeman's definition [5]. 15 Suppose {P T D j N } j<k were known and the linear combination constants {c 1 , c 2 } were chosen. We will use an inductive algorithm to define the matrices P T D j N consecutively in j.
By equation 8, we can calculate r j+1 1 by noting that an orthogonal transformation does not change inner products. Thus, the RHS of equation 8 can be calculated and r j+1 1 can be set equal to it. Next, we define r 2 = l T t D j+1 N in the following manner: Thus, given r j 1 from the normalization constraints and D j I from the image information, we can find the next row r j 2 uniquely. Note that we are assuming that l t exists, which it will at every regular point. It remains to define r j 3 as the linear combination of the previous two rows r j 3 = c 1 r j 1 + c 2 r j 2 . Now, we have defined P T D j N by its rows {r j 1 , r j 2 , r j 3 } and we continue on to P T D j+1 N .
It remains to define the base case: P T D 1 N . From equation 5, we know its first row r 1 1 must be the 0 vector. From equation 25, we find that r 1 2 is just the weighted brightness gradient r 1 2 = D 1 I √ 1−I 2 . Finally, due to the generic assumption, we set r 1 3 = c 2 r 1 2 . Proposition. Following the algorithm described above (and summarized in Algorithm 1) yields coefficients {D j N } m 1 defining a multivariate Taylor poly-nomialN (x, y) that is generic according to [6], matches exactly the image Taylor approximationĪ(x, y), and is unit length (up to error O(h m )).
See the Appendix for precise details regarding the unit length error.

Discussion
There have been many attempts to relate local image derivatives to local surface derivatives in Lambertian shading, but due to the 1 + f 2 x + f 2 y in the denominator, the derivatives become more and more complex and the analysis soon becomes intractable. There have also been approaches towards representing the surface in a different way (principal directions basis, covariant derivatives, stereographic projections) but all tend to gain complexity as more derivatives are considered. In this work, we have described a method of representation that does not increase in analytic complexity as more derivatives are considered -this allows us to calculate all generic Taylor expansions (of any order) of a surface for a given image.
Essentially, we project each set of Taylor coefficients D j N onto an unknown plane (the "visible attribute plane") defined by some unknown rotation matrix P . This is a plane in R 3 that is spanned by the light source l t and the central normal N 0 . We know the projections as the D j N , N 0 is determined by the unit normal constraint. The D j N , l t is a function of the image and D j N , N 0 . We do not know this visible attribute plane but we can calculate these projections onto it from the image. Thus, any two normal fields with the same Taylor coefficients D j N will have equivalent projections onto this plane at each differential level j. The ambiguity stems from the fact that we cannot know the heights of these tensors above the plane (the projection of D j N along b).
Two final remarks: First, it is possible to make the problem more well posed (less ambiguous) if information from the boundaries is included. Sec-  ond, our algorithm can be modified to generate all Taylor polynomial normal fieldsN (x, y) (instead of only the generic ones) for a given image patch by modifying the above algorithm: simply select r j 3 arbitrarily for each j. That is, rather than selecting r j 3 as a linear combination of r j 1 and r j 2 , r j 3 can be set equal to any v ∈ R 2 j (satisfying the appropriate symmetry according to differentiation, e.g. N xy = N yx ). Thus, we see that the generic assumption significantly reduces the ambiguity of our Taylor expansion: we either must choose 2 generic constants or m j=1 j + 1 other values. In summary, image patches are related to surface patches by some rendering function. Lambertian models with orthographic projection are the simplest such functions, and amount to projection of surface normals against a light source vector. This suggests a linear algebraic approach. We here take this one step further, for the inverse problem, by asking which surface patches could be consistent with a given image patch up to some number of derivatives. This allows identification of stucture at different levels, although these 'Taylor tensors' can only be determined up to unknown generic constants c 1 , c 2 and unknown pose P ∈ SO 3 . Nevertheless, it leads to an algorithm for calculating every possible consistent normal field, even with a generic lighting assumption. Perhaps not surprisingly, such visualizations go well beyond the bas-relief ambiguity, even with generic constraints.