ShapeFromShading.jl

img = generate_surface(shape::AbstractSyntheticShape, albedo::Real = 1.0, illumination_direction::Vector{T} where T <: Real = [0.0, 0.0, 1.0]; img_size::Int=151, noise::Real = 0, relative_noise::Bool = true)

Creates a synthetic image of a given shape with the given illumination properties.

Output

Returns a grayscale image with dimensions img_size in both axes.

Details

Uses the passed properties to calculate the partial differentials of the shape reflectance and uses these to generate a synthetic image of the shape to be used as a consistent test image for SFS testing.

Arguments

The function arguments are described in more detail below.

shape

A AbstractSyntheticShape that specifies the type of shape to be generated. Available shapes are; SynthSphere, Prism, Ripple, Ripple2, Cake, Cake2, SynthVase, Tent, Dem, SynthGaussian and SuperGaussian. Each has a parameter radius, a Real that specifies the radius of the desired shape or its scale within the image. This can acts as either an actual radius (such as for SynthSphere) or rather as a scale factor (such as for Ripple). If left unspecified a default value is used for each.

albedo

A Real that specifies the albedo (amount of light reflected) of the image. If left unspecified a default value of 1.0 is used.

illumination_direction

A Vector{T} where T <: Real that specifies the tilt value to be used by the algorithm. The illumination_direction should be a vector of the form [x,y,z] where x,y,z are in the range [0,1]. If left unspecified a default value of [0.0,0.0,1.0] is used.

img_size

A Int that specifies the dimentions of the final image should be. Defults to 151.

noise

A Real that specifies the standard deviation of noise to be added to the image. Defults to 0.0.

relative_noise

A Bool that specifies if the noise added to the image should be relative to the gradient of the object displayed in the image at that point.

Example

Compute the heightmap for a synthetic image generated by generate_surface.

using Images, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(50), 1.0, [0.5,0.1,0.7])

Reference

1. S. Elhabian, "Hands on Shape from Shading", Computer Vision and Image Processing, 2008.

(scheme::AbstractIntegrator)(p::AbstractArray, q::AbstractArray)

Integrates the gradient field defined by p and q using various methods as described in the following sections.

Output

Returns a heightmap Z which defines the reconstructed surface.

Details

Using the method defined in each type of integrator, each algorithm is run using parameters passed to the struct when creating the integrator. Results can vary wildly depending on the chosen integrator. Available methods are; Frankot, Path, SplitPath, Horn, Durou, Quadratic, TotalVariation, NonConvex1, NonConvex2, AnisotropicDiffusion and MumfordShah. See individual documentation for each method for more detials.

Parameters

Each algorithm has its own set of parameters which are discussed in each the section for each method.

Arguments

p and q:

AbstractArray's which define the gradient field to be integrated.

Example

The following example demonstrates the use of the Frankot integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(SynthSphere(38), img_size = 151)

# Create a Frankot() integrator
frankot = Frankot()

# Calculate the heightmap from the gradients
Z = frankot(p, q)

# Normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
surface(r, r, Z)

Frankot()

Defines the Frankot integrator which contians the Frankot-Chellappa method of integration. The Frankot-Chellappa method is a fast, reliable method for integration surface normals while enforcing integrability using Fourier methods.

Output

Frankot() returns a Frankot integrator which can then be called to run the Frankot-Chellappa method on a gradient field.

Details

Frankot-Chellappa method uses Fourier methods to attempt to solve the Poission equation $\nabla^2z = \partial_up + \partial_vq$. By taking the Fourier transform of both sides we get:

By rearranging the above equation we arrive at an equation for $\hat{z}$;

From which the final surface can be found by taking the inverse Fourier transform of $\hat{z}$.

Due to the way $(\omega_u, \omega_v)$ is defined the algorithm works best when the input dimensions are odd in length. To accommodate this the integrator will pad the edge of the inputs if they are even before running the algorithm. This padding will be removed before returning a value hence output size will be unaffected.

Parameters

Frankot integrator take no parameters.

Example

The following example demonstrates the use of the Frankot integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(SynthSphere(38), img_size = 151)

# Create a Frankot() integrator
frankot = Frankot()

# Calculate the heightmap from the gradients
Z = frankot(p, q)

# Normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
surface(r, r, Z)

Reference

[1] R. T. Frankot and R. Chellappa, "A method for enforcing integrability in shape from shading algorithms," in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 10, no. 4, pp. 439-451, July 1988. doi: 10.1109/34.3909

Path()

Creates a Path() integrator which utilises the average of two path integrals along varying paths. Each path integral reconstructs the surface with accumulating error along the path, hence averaging two different paths can minimise this error, although the method still suffers if the gradient field is not integrable at some points.

Output

Path() returns a Path integrator which can then be called to integrate a gradient field.

Details

Under the assumption that the surface normals are approximately integrable everywhere ($\frac{\partial p}{\partial y}\approx\frac{\partial q}{\partial x}$), then surface can be reconstructed using the path integral defined as:

Which can be broken into two integrals representing the value at each point on the surface as shown below for a path which integrates along the first column then along the row.

The second path used in the algorithm is simply the transpose of the first, integrating along the first row then down the column represented mathematically as:

The algorithm can be written, then discretised as shown below:

It is important to note as mentioned above if there are non-integrable points in the normal field then artefacts can appear in the reconstruction. This is seen in the example below where the otherwise smooth sphere appears "spiky". This can be corrected post reconstruction by smoothing but ideally a different integrator should be used.

Parameters

Path integrator take no parameters.

Example

The following example demonstrates the use of the Path integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(SynthSphere(38), img_size = 151)

# Create a Path() integrator
path = Path()

# Calculate the heightmap from the gradients
Z = path(p, q)

# Normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
surface(r, r, Z)

Reference

[1] D. Forsyth and J. Ponce, Computer vision: a modern approach. Upper Saddle River, N.J: Prentice Hall, 2003, pp. 84-86.

SplitPath()

Creates a SplitPath() integrator which utilizes the average of two path integrals along varying paths averaging the value at each step. Each path integral reconstructs the surface with accumlating error along the path, hence averaging two different paths at each step reduces the global error at the cost of local error, although the method still suffers if the gradient field is not integrable at some points it does less so the Path() from which it extends.

Output

SplitPath() returns a SplitPath integrator which can then be called to integrate a gradient field.

Details

Under the assumption that the surface normals are approximately integrable everywhere ($\frac{\partial p}{\partial y}\approx\frac{\partial q}{\partial x}$), then surface can be reconstructed using the path integral defined as:

By expanding on this principle and the discreate summation from Path() we can arrive at the discreate expresion for the value at each point, assuming all values prior to that point have been calculated, as follows:

As with other similar methods (see Horn()) care must be taken with regards to boundaries which can be calculated, to a constant value $z(0,0)$ which is assumed to be the zero point, using:

It is important to note as mentioned above if there are non-integrable points in the normal field then artefacts can appear in the reconstruction. These errors gradually average out but will lead to "streaks" appearing in the reconstruction. This is seen in the example below where the otherwise smooth sphere appears has ripple like structures pointing toward to top right corner. This can be corrected post reconstruction by smoothing but ideally a different integrator should be used. It is also interesting to note the parallels between this method and the Horn and Brooks method, with this method being effectively the forward component of Horn's method. As such this algorithm provided a useful middle ground between direct integration algorithms and iterative algorithms such as the Horn and Brooks method.

Parameters

SplitPath integrator take no parameters.

Example

The following example demonstrates the use of the SplitPath integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(SynthSphere(38), img_size = 151)

# Create a Path() integrator
splitPath = SplitPath()

# Calculate the heightmap from the gradients
Z = splitPath(p, q)

# Normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
surface(r, r, Z)

Reference

[1] D. Forsyth and J. Ponce, Computer vision: a modern approach. Upper Saddle River, N.J: Prentice Hall, 2003, pp. 84-86. [2] B. Horn and M. Brooks, "The variational approach to shape from shading", Computer Vision, Graphics, and Image Processing, vol. 33, no. 2, pp. 174-208, 1986. doi: 10.1016/0734-189x(86)90114-3

Horn(ϵ::Real = 1.0, max_iter::Real = 10000)

Implements the Horn and Brook's method of integrating surface normals. This algorithm offers an iterative solution to the Poisson equation describing the surface providing good reconstructions under most conditions.

Output

Horn() returns a Horn integrator which can then be called to integrate a gradient field.

Details

The Horn and Brook's method attempts to solve the Poisson equation $\nabla^2z = \partial_up + \partial_vq$ by the discretization given below:

Which can be rearranged to give the iterative scheme provided by:

This scheme will always converge to a solution however the rate of convergence may depend upon the initial solution. This implementation will initilize with a zero solution. Neumann boundary conditions are imposed at the edges where the scheme would otherwise go out of bounds.

Parameters

The function parameters are described in more detail below.

Max_iter:

An Int which controls the number of iterations the algorithm will run for.

ϵ:

A Real representing the distance between pixels. This will Control how tall the final reconstruction is relative the array grid.

Example

The following example demonstrates the use of the Horn integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(SynthSphere(38), img_size = 151)

# Create a Horn() integrator
horn = Horn(ϵ = 0.03, max_iter = 10000)

# Calculate the heightmap from the gradients
Z = horn(p, q)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
surface(r, r, Z)

Reference

[1] B. Horn and M. Brooks, "The variational approach to shape from shading", Computer Vision, Graphics, and Image Processing, vol. 33, no. 2, pp. 174-208, 1986. doi: 10.1016/0734-189x(86)90114-3

Durou(ϵ::Real = 1.0, max_iter::Real = 1000)

Implements the Durou and Courteille method of integrating surface normals. This algorithm offers an iterative solution to the Poisson equation describing the surface extending Horn and Brook's method by improving the boundary approximation and providing good reconstructions under most conditions.

Output

Durou() returns a Durou integrator which can then be called to integrate a gradient field.

Details

TheDurou and Courteille's method attempts to solve the Poisson equation $\nabla^2z = \partial_up + \partial_vq$ by the discretization given below:

Which can be rearranged to give the iterative scheme provided by:

This scheme will always converge to a solution however the rate of convergence may depend upon the initial solution. This implementation will initialize with a zero solution. Natural boundary conditions are imposed at the edges using the condition $\partial_uz-p+\partial_v-q=0$. Although faster then the Horn and Brook's method and better at handling boundaries, it can generate a worse solution under some conditions.

Parameters

The function parameters are described in more detail below.

Max_iter:

An Int which controls the number of iterations the algorithm will run for. the range [0,1].

ϵ:

A Real representing the distance between pixels. This will Control how tall the final reconstruction is relative the array grid.

Example

The following example demonstrates the use of the Durou integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(SynthSphere(38), img_size = 151)

# Create a Durou() integrator
durou = Durou(ϵ = 0.03, max_iter = 10000)

# Calculate the heightmap from the gradients
Z = durou(p, q)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
surface(r, r, Z)

Reference

[1] Y. Quéau, J. Durou and J. Aujol, "Normal Integration: A Survey", Journal of Mathematical Imaging and Vision, vol. 60, no. 4, pp. 576-593, 2017. doi: 10.1007/s10851-017-0773-x

Quadratic(z::AbstractArray, λ::AbstractArray = fill(10.0^-6, size(z)), mask::AbstractArray = fill(1.0, size(z)))

Implements the quadratic variational least squared method proposed by Aujol, Durou and Quéau. The algorithm solves the least squares system generated from minimizing a fidelity term $\mathcal{F}(z) = \iint_{(u,v)\in \Omega}\Phi(||\nabla z(u,v) - g(u,v)||)dudv$ and a regularization term $\mathcal{R}(z) = \iint_{(u,v)\in \Omega}\lambda \left[z(u,v)-z^0(u,v)\right]^2dudv$ where $\Phi(s) = s^2$. This method is able to quickly produce good solutions on a smooth surface and can easily handle non-rectangular domains or sub-divided domains and can provide a good starting solution for other algorithms.

Output

Quadratic() returns a Quadratic integrator which can then be called to run the quadratic method on a gradient field.

Details

As mentioned above the quadratic variational least squared method aims to minimize a fidelity term $\mathcal{F}(z) = \iint_{(u,v)\in \Omega}\Phi(||\nabla z(u,v) - g(u,v)||)dudv$ and a regulization term $\mathcal{R}(z) = \iint_{(u,v)\in \Omega}\lambda \left[z(u,v)-z^0(u,v)\right]^2dudv$ where $\Phi(s) = s^2$.

This leads to the minimization problem given by:

By discretizing the problem we can arrive at the functional:

where $\Omega_u^+$ represents the domain where $(u,v)\in\{(u,v)\in\Omega|(u,v)(u+1,v)\in\Omega\}$ etc.. Using this definition the discrete differences can be converted into matrix form, where $p_{u,v}=z_{u+1,v}-z_{u,v}$ is the forward difference in the u direction etc. This data is then stacked into three vectors; $\mathbf{z},\mathbf{p},\mathbf{q}\in\R^{|\Omega|}$. Thus the matrix reresenting each of these is defined as below where m(i) is the mapping of the ith element of this vector to its corresponding point in $(u,v)$ and $D_u^+$ is a $|\Omega|\times|\Omega|$ matrix.

For a 2X2 domain this looks like:

The other three discrete differences matrices are similarly defined from there definitions to be $D_u^-$, $D_v^+$ and $D_v^-$. These can then be used to redefine the minimization problem to be in the form:

Where $\Lambda$ is the $|\Omega|\times|\Omega|$ diagonal matrix containing the values of $\sqrt{\lambda_{u,v}}$. Using the above definitions the negative Laplacian matrix can then be defined as:

Finally the least minimization problem can be represented in the form of a least squares problem of the form $A\mathbf{z}=\mathbf{b}$ where:

This system is then solved using a standard conjugate gradient algorithm where the initialization has only a slight impact on the runtime and no impact on the final solution. The algorithm provides good results on smooth surfaces but struggles in the presence of discontinuities.

Parameters

z:

An AbstractArray which defines the value of $z^0$ the initial solution and prior to be used in the regularization term. Must be provided.

λ:

An AbstractArray the same size as z, defaulting to $10.0^{-6}$ everywhere. This defines theregulization weight at each point. Large values will force the algorithm to keep the solution near to $z^0$ at that position. Can be used to keep the solution near the initial solution or guide the solution to a certain known value at points (i.e. known maxima and minima). This value should be set uniformly small otherwise.

mask:

An AbstractArray the same size as z, which guides the algorithm as to where the valid domain is. Values of 1 will be in the domain $\Omega$ while other values will be ignored and set to $z^0$. This can be used to integrate over sub-domain or to segment the domain into parts. The gen_mask() funtion can be used to generate a mask which will remove non-integrable regions dramatically improving the solution under most condition at the cost of not integrating the entire solution.

Example

The following example demonstrates the use of the Quadratic integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(Prism(75), img_size = 151)

# Create a Quadratic() integrator
quadratic = Quadratic(z=zeros(size(p)))
quadraticMasked = Quadratic(z=zeros(size(p)), mask=gen_mask(p,q,1.0)[:,:,1])

# Calculate the heightmap from the gradients
Z = quadratic(p, q)
Z2 = quadraticMasked(p, q)

# Normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)
Z2 = Z2./maximum(Z2)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
vbox(surface(r, r, Z), surface(r, r, Z2))

Reference

[1] Y. Quéau, J. Durou and J. Aujol, "Variational Methods for Normal Integration", Journal of Mathematical Imaging and Vision, vol. 60, no. 4, pp. 609-632, 2017. doi: 10.1007/s10851-017-0777-6

TotalVariation(z::AbstractArray, α::Real = 1.0, λ::AbstractArray = fill(10.0^-6, size(z)), mask::AbstractArray = fill(1.0, size(z)), max_iter::Int = 100)

Implements the total variational method proposed by Aujol, Durou and Quéau. The algorithm solves the same minimization problem as th Quadraticmethod exept the fidelity terms function $\Phi(s)=||s||_{L_1}$. This method is able to produce good solutions on a smooth and piecewise smooth surface and can easily handle non-rectangular domains or sub-divided domains.

Output

TotalVariation() returns a TotalVariation integrator which can then be called to run the quadratic method on a gradient field.

Details

As discussed above this algorithm used the same fidelity and regularization terms as the Quadratic method exept the $\Phi(s)=s^2$ term is replaced with $\Phi(s)=||s||_{L_1}$. This leads to the minimization problem defined by:

By considering the four posible discreatisations of $∇z(u,v)$ we can generate four posible domains to consider given by; $\Omega^{UV}=\Omega_u^U\cup\Omega_v^V, (U,V)\in\{+,-\}^2$ where $\{+,-\}^2$ refers to all posible combinations of +,-. Using these the following discreate functional can be generated:

Which simplifies to the minimization problem:

This leads to the optimization scheme using an ADMM algorithm defined by:

The z update then can be solved using the linear system defined below, where D_{u,v}^{U,V} and $\Lambda$ are the same at those defined in Quadratic.

Where $\mathbf{P}^{UV^{(k)}}, \mathbf{Q}^{UV^{(k)}}$ are the u and v components of $\mathbf{g}+\mathbf{r}^{UV^{(k)}}-\mathbf{b}^{UV^{(k)}}$. This can be solved using conjugate gradient. Finally the update to $\mathbf{r}^{UV}$ can be computed as:

Parameters

z:

An AbstractArray which defines the value of $z^0$ the initial solution and prior to be used in the regulization term. Must be provided.

α:

A Real with defult value of 1.0 which controls the step size. In theory this value should have no impact on final solution but in practice larger values can lead to worse solutions while values which are two small may lead non-convergance in the least square update step, causing the algorithm to hang for long periods of time. Values below 0.25 are not recomended but may work depending on domain size and inputs.

λ:

An AbstractArray the same size as z defulting to $10.0^{-6}$ everywhere. This defines theregulization weight at each point. Large valueas will force the algorithm to keep the solution near to $z^0$ at that position. Can be used to keep the solution near the initial solution or guide the solution to a certian known value at points (i.e. known maxima and minima). This value should be set uniformly small otherwise.

mask:

An AbstractArray the same size as z, which guides the algorithm as to where the valid domain is. Values of 1 will be in the domain $\Omega$ while other values will be ignored and set to $z^0$. This can be used to integrate over sub-domain or to segment the domain into parts. The gen_mask() funtion can be used to generate a mask which will remove non-integrable regions dramatically improving the solution under most condition at the cost of not integrating the entire solution.

Example

The following example demonstrates the use of the TotalVariation integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(Prism(75), img_size = 151)

# Create a TotalVariation() integrator
totalVariation = TotalVariation(z=zeros(size(p)), α=1.0)
totalVariationMasked = TotalVariation(z=zeros(size(p)), α=0.5, mask=gen_mask(p,q,1.0)[:,:,1])

# Calculate the heightmap from the gradients
Z = totalVariation(p, q)
Z2 = totalVariationMasked(p, q)

# Normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)
Z2 = Z2./maximum(Z2)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
vbox(surface(r, r, Z), surface(r, r, Z2))

Reference

[1] Y. Quéau, J. Durou and J. Aujol, "Variational Methods for Normal Integration", Journal of Mathematical Imaging and Vision, vol. 60, no. 4, pp. 609-632, 2017. doi: 10.1007/s10851-017-0777-6

NonConvex1(z::AbstractArray, β::Real = 0.5, λ::AbstractArray = fill(10.0^-6, size(z)), mask::AbstractArray = fill(1.0, size(z)), max_iter::Int = 100)

The first of two non-convex regularization methods proposed by Aujol, Durou and Quéau. The same fidelity and regularization terms to minimize as Quadratic are used, but the convexity of $\Phi$ is sacrificed in order to gain better behaviour around discontinuities and outliers.

Output

NonConvex1() returns a NonConvex1 integrator which can then be called to run the non-convex regularization method method on a gradient field.

Details

The initial minimization problem is the same as in TotalVariation as given below exept $Phi$ is defined using the given function instead of the $L_1$ norm.

The problem can then be discretized to form the following functional where $D_{u,v}^{UV}$ is a $2\times|\Omega|$ matrix formed from stacking the vectors of the finite diferences.

As $f(\mathbf{z})$ is smooth but not convex and $g(\mathbf{z})$ is convex the iPiano algorithm is then used to iterativly solve the minimization of the functional such that:

where $(I+\alpha_1\partial g)^{-1}$ is a proximal operator defined as:

Using the definition of $\Phi$ given above, the derivative of $f(\mathbf{z})$ can be computed as below:

The values of $\alpha_1$ and $\aplha_2$ control the step size of the algorith and are chosen such that $\alpha_2$ is fixed at 0.8 while $\alpha_1$ is chosen using the lazy backtracking method, which uses a Lipschitz constant to caculate a suitably small step size at each step using the below relationship where $\eta>1$ is a constant:

This method, being non-convex, highly relies on a good initial solution and will often only provide a minimal improvment to the solution. A bad initial solution will produce a final solution which does not resemble the surface under most conditions.

Parameters

z:

An AbstractArray which defines the value of $z^0$ the initial solution and prior to be used in the regulization term.

β:

A Real which acts as a hyper-parameter to the function. Large values for β will produce smoother functions but loose discontinuities. Smaller value will preserve discontinuities but lead to staircassing in the solution. Defults to 0.5.

λ:

An AbstractArray the same size as z defulting to $10.0^-6$ everywhere. This defines theregulization weight at each point. Large valueas will force the algorithm to keep the solution near to $z^0$ at that position. Can be used to keep the solution near the initial solution or guide the solution to a certian known value at points (i.e. known maxima and minima). This value should be set uniformly small otherwise.

mask:

An AbstractArray the same size as z, which guides the algorithm as to where the valid domain is. Values of 1 will be in the domain $\Omega$ while other values will be ignored and set to $z^0$. This can be used to integrate over sub-domain or to segment the domain into parts. The gen_mask() funtion can be used to generate a mask which will remove non-integrable regions dramatically improving the solution under most condition at the cost of not integrating the entire solution.

Example

The following example demonstrates the use of the NonConvex1 integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(Prism(75), img_size = 151)

# Create a NonConvex1() integrator
nonConvex1 = NonConvex1(z=zeros(size(p)), β=0.5)
nonConvex1Init = NonConvex1(z=Horn()(p,q), β=0.5)

# Calculate the heightmap from the gradients
Z = nonConvex1(p, q)
Z2 = nonConvex1Init(p, q)

# Normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)
Z2 = Z2./maximum(Z2)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
vbox(surface(r, r, Z), surface(r, r, Z2))

Reference

[1] Y. Quéau, J. Durou and J. Aujol, "Variational Methods for Normal Integration", Journal of Mathematical Imaging and Vision, vol. 60, no. 4, pp. 609-632, 2017. doi: 10.1007/s10851-017-0777-6

NonConvex2(z::AbstractArray, γ::Real = 1.0, λ::AbstractArray = fill(10.0^-6, size(z, mask::AbstractArray = fill(1.0, size(z)), max_iter::Int = 100)

The second of two non-convex regularization methods proposed by Aujol, Durou and Quéau. The same as NonConvex1 exept $\Phi$ has been replaced with a different non-convex function.

Output

NonConvex2() returns a NonConvex2 integrator which can then be called to run the non-convex regularization method method on a gradient field.

Details

The implimentation is the same as the one in TotalVariation as given below exept $Phi$ is defined as below:

Using the definition of $\Phi$ given above, the derivative of $f(\mathbf{z})$ can be computed as below:

This method, being non-convex, highly relies on a good initial solution and will often only provide a minimal improvment to the solution. A bad initial solution will produce a final solution which does not resemble the surface under most conditions.

Parameters

z:

An AbstractArray which defines the value of $z^0$ the initial solution and prior to be used in the regulization term.

γ:

A Real which acts as a hyper-parameter to the function. Large values for γ will produce smoother functions but loose discontinuities. Smaller value will preserve discontinuities but lead to staircassing in the solution. Defults to 1.0.

λ:

An AbstractArray the same size as z defulting to $10.0^-6$ everywhere. This defines theregulization weight at each point. Large valueas will force the algorithm to keep the solution near to $z^0$ at that position. Can be used to keep the solution near the initial solution or guide the solution to a certian known value at points (i.e. known maxima and minima). This value should be set uniformly small otherwise.

mask:

An AbstractArray the same size as z, which guides the algorithm as to where the valid domain is. Values of 1 will be in the domain $\Omega$ while other values will be ignored and set to $z^0$. This can be used to integrate over sub-domain or to segment the domain into parts. The gen_mask() funtion can be used to generate a mask which will remove non-integrable regions dramatically improving the solution under most condition at the cost of not integrating the entire solution.

Example

The following example demonstrates the use of the NonConvex2 integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(Prism(75), img_size = 151)

# Create a NonConvex2() integrator
nonConvex2 = NonConvex2(z=zeros(size(p)), γ=1.0)
nonConvex2Init = NonConvex2(z=Horn()(p,q), γ=1.0)

# Calculate the heightmap from the gradients
Z = nonConvex2(p, q)
Z2 = nonConvex2Init(p, q)

# Normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)
Z2 = Z2./maximum(Z2)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
vbox(surface(r, r, Z), surface(r, r, Z2))

Reference

[1] Y. Quéau, J. Durou and J. Aujol, "Variational Methods for Normal Integration", Journal of Mathematical Imaging and Vision, vol. 60, no. 4, pp. 609-632, 2017. doi: 10.1007/s10851-017-0777-6

AnisotropicDiffusion(z::AbstractArray, μ::Real = 5.0, ν::Real = 10., λ::AbstractArray = fill(10.0^-6, size(z)), mask::AbstractArray = fill(1.0, size(z)), max_iter::Int = 10)

Defines the anisotropic diffusion method as proposed by Aujol, Durou and Quéau. It utilizes an anisotropic diffusion like process to create a weighted version of the least squares problem solved in the Quadratic method.

Output

AnisotropicDiffusion() returns a AnisotropicDiffusion integrator which can then be called to run the anisotropic diffusion method method on a gradient field.

Details

The minimization problem from Quadratic is modified with the addition weighting term $W(u,v)$ to reach the following minimization problem.

The weighting term can be defined as below where $\mu$ and $\nu$ are parameters which control the impact of gradient of the reconstruction at each iteration and the input gradients.

This lends itself to the creation of the terms $A^{UV}$ and $B^{UV}$ which are $|\Omega|\times|\Omega|$ diagonal matrices, respectively, contianing the values of;

and

where $(U,V)\in\{+,-\}^2$. Using these definititions the original minimization problem can be rewriten as the following iterative scheme:

which is then solved using Cholesky factorization.

Parameters

z:

An AbstractArray which defines the value of $z^0$ the initial solution and prior to be used in the regulization term.

μ:

A Real which acts as a hyper-parameter to the function. Allows for the tuning of discontinuities. Defults to 5.0.

ν:

A Real which acts as a hyper-parameter to the function. Allows for the tuning of discontinuities. Unlike μ is has a minimal impact on the solution and generally does not need adjusting unless an extreme value for μ is used where it can help to balance out the two terms. Defults to 10.0.

λ:

An AbstractArray the same size as z defulting to $10.0^-6$ everywhere. This defines theregulization weight at each point. Large valueas will force the algorithm to keep the solution near to $z^0$ at that position. Can be used to keep the solution near the initial solution or guide the solution to a certian known value at points (i.e. known maxima and minima). This value should be set uniformly small otherwise.

mask:

An AbstractArray the same size as z, which guides the algorithm as to where the valid domain is. Values of 1 will be in the domain $\Omega$ while other values will be ignored and set to $z^0$. This can be used to integrate over sub-domain or to segment the domain into parts. The gen_mask() funtion can be used to generate a mask which will remove non-integrable regions dramatically improving the solution under most condition at the cost of not integrating the entire solution.

Example

The following example demonstrates the use of the AnisotropicDiffusion integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(Prism(75), img_size = 151)

# Create a AnisotropicDiffusion() integrator
anisotropicDiffusion = AnisotropicDiffusion(z=zeros(size(p)), γ=1.0)
anisotropicDiffusionInit = AnisotropicDiffusion(z=Horn()(p,q), γ=1.0)

# Calculate the heightmap from the gradients
Z = anisotropicDiffusion(p, q)
Z2 = anisotropicDiffusionInit(p, q)

# Normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)
Z2 = Z2./maximum(Z2)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
vbox(surface(r, r, Z), surface(r, r, Z2))

Reference

[1] Y. Quéau, J. Durou and J. Aujol, "Variational Methods for Normal Integration", Journal of Mathematical Imaging and Vision, vol. 60, no. 4, pp. 609-632, 2017. doi: 10.1007/s10851-017-0777-6

MumfordShah(z::AbstractArray, μ::Real = 10.0, ϵ::Real = 0.1, λ::AbstractArray = fill(10.0^-6, size(z, mask::AbstractArray = fill(1.0, size(z)), max_iter::Int = 50)

Defines the Mumford-Shah method, the final method proposed by Aujol, Durou and Quéau. It utilizes an adepted version of Mumford and Shah functional to provide a good reconstruction in the presence of discontinuities at the const of a longer runtime then other methods.

Output

MumfordShah() returns a MumfordShah integrator which can then be called to run the Mumford-Shah method method on a gradient field.

Details

This method involves modifying the minimization problem given by the Mumford and Shah functional given below:

Where $K$ is the set of discontinuities. This can then be adapted to our problem by utilizing a Ambrosio-Tortelli approximation to achaive the following functional where $\epsilon\to0$.

This leads to the following discretization of the functional where $\mathbf{w}^{+/-}_{u/v}$ are thevecotr of weights in each forward and backward direction and the matrix $W^{+/-}_{u/v}$ is the diagonal matrix formed from this vector.

This is then solved with a conjugate gradient algorithm and an alternating optimization scheme at each step where the updates are found using the relationships below:

This method can produce good results if the parameters are appropriotly tuned. If the parameters are too large the solution will suffer heavily from staircasing artifacts while setting it too small will result in a smooth solution. Even if the value is chosen correctly the algorithm tends to overfit the final solution to the discontinuities and they will extend into parts of the solution which are actually smooth.

Parameters

z:

An AbstractArray which defines the value of $z^0$ the initial solution and prior to be used in the regulization term.

ϵ:

A Real which acts as a hyper-parameter to the function. Controls how the final solution will converge. Large values will lead to staircasing while small values will over-smooth the surface. Must be relativly small to achieve convergence to a solution. Defults to 0.1.

μ:

A Real which acts as a hyper-parameter to the function. This value controls the smoothness of the final solution. Large values will lead to staircasing while small values will lead to over-smoothed solutions. Defults to 10.0.

λ:

An AbstractArray the same size as z defulting to $10.0^-6$ everywhere. This defines the regulization weight at each point. Large values will force the algorithm to keep the solution near to $z^0$ at that position. Can be used to keep the solution near the initial solution or guide the solution to a certian known value at points (i.e. known maxima and minima). This value should be set uniformly small otherwise.

mask:

An AbstractArray the same size as z, which guides the algorithm as to where the valid domain is. Values of 1 will be in the domain $\Omega$ while other values will be ignored and set to $z^0$. This can be used to integrate over sub-domain or to segment the domain into parts. The gen_mask() funtion can be used to generate a mask which will remove non-integrable regions dramatically improving the solution under most condition at the cost of not integrating the entire solution.

Example

The following example demonstrates the use of the MumfordShah integrator.

using ShapeFromShading, Makie

# Generate synthetic gradients
p, q = synthetic_gradient(Prism(75), img_size = 151)

# Create a MumfordShah() integrator
mumfordShah = MumfordShah(z=zeros(size(p)), γ=1.0)
mumfordShahInit = MumfordShah(z=Horn()(p,q), γ=1.0)

# Calculate the heightmap from the gradients
Z = mumfordShah(p, q)
Z2 = mumfordShahInit(p, q)

# Normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)
Z2 = Z2./maximum(Z2)

# Display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
vbox(surface(r, r, Z), surface(r, r, Z2))

Reference

[1] Y. Quéau, J. Durou and J. Aujol, "Variational Methods for Normal Integration", Journal of Mathematical Imaging and Vision, vol. 60, no. 4, pp. 609-632, 2017. doi: 10.1007/s10851-017-0777-6

DiscreteShape(iterations::Int = 2000, smoothness::Int = 1000, albedo::Real = Inf, illumination_direction::Vector{T} where T <: Real = [Inf, Inf, Inf])

Contins DiscreteShape algorithm which ttempts to produce a heightmap from a grayscale image by minimization of a set of Euler-Lagrange equations. This is done discretely at each point in the image utilizing the second derivates of the surface normals and their Fourier Transforms.

Output

Returns a DiscreteShape functor which can be run on an image to attempt reconstruction of the surface contianed in the image.

Details

The algorithm attempts to minimise the brightness deviation $ϵ$ through the brightness constraint $\epsilon_1$ and smoothness constraint $\epsilon_2$ as defined bellow:

This is minimised using the Euler-Lagrange equations defined as;

and

which become:

which can be further simplified to give:

where $\nabla^2p=p_{xx}+p_{yy}$ and $\nabla^2q=q_{xx}+q_{yy}$ are Laplacians of p and q.

However, for computation we are dealing with a discreate case of these equations which can be defined as below:

where

and

Finally, the algorithm needs to enforce integrability on p and q and retrieve the surface Z. This can be done by taking the Fast Fourier Transform of p and q ()$c_p(\omega_x,\omega_y)$ and $c_q(\omega_x,\omega_y)$) and then using the Inverse Fast Fourier Transform to recover Z and update p and q as per bellow:

where

Arguments

The arguments are described in more detail below.

albedo

A Real that specifies the albedo (amount of light reflected) of the image. If albedo is specified to must the illumination_direction. Defults to Inf which will trigger the algorithm to run estimate_img_properties.

illumination_direction

A Vector{T} where T <: Real that specifies the tilt value to be used by the algorithm. The illumination_direction should be a vector of the form [x,y,z] where x,y,z are int he range [0,1]. If illumination_direction is specified to must the albedo. Defults to [Inf, Inf, Inf] which will trigger the algorithm to run estimate_img_properties.

iterations

An Int that specifies the number of iterations the algorithm is to perform. If left unspecified a default value of 2000 is used.

smoothness

An Int that specifies the strength of the smoothness constraint in the minimised function. Defults to 1000.

Example

Compute the heightmap for a synthetic image generated by generate_surface.

using Images, Makie, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(50), 1.0, [0.2,0,0.9])

#calculate the heightmap (using 500 iterations)
discreteShape = DiscreteShape(iterations = 500, albedo = 1.0, illumination_direction = [0.2,0,0.9])
Z,p,q = discreteShape(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)

Reference

1. S. Elhabian, "Hands on Shape from Shading", Computer Vision and Image Processing, 2008.

DiscreteShapeBound(iterations::Int = 2000, smoothness::Int = 1000, albedo::Real = Inf, illumination_direction::Vector{T} where T <: Real = [Inf, Inf, Inf])

Same as DiscreteShape except it has its initial conditions bound by the image as per bellow where $E$ is the brightness of the image:

Example

Compute the heightmap for a synthetic image generated by generate_surface.

using Images, Makie, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(50), 1.0, [0.2,0,0.9])

#calculate the heightmap (using 500 iterations)
discreteShapeBound = DiscreteShapeBound(iterations = 500, 1.0, [0.2,0,0.9])
Z,p,q = discreteShapeBound(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)

Reference

1. S. Elhabian, "Hands on Shape from Shading", Computer Vision and Image Processing, 2008.

Pentland(slant::Real = Inf, tilt::Real = Inf)

Creates a Pentland fucntor which attempts to produce a heightmap from a grayscale image using Pentland's algorithm. Under the assumptions of Lambertian surface, orthographic projections, the surface being illuminated by distant light sources, the surface is not self-shadowing and the surface has constant albedo, hence it can be ignored. The algorithm employs Tayler series expansion and Fourier transforms to compute a non-iterative solution.

Output

Returns a Pentland functor which can be used to run the Pentland algorithm to reconstruct a surface form a grayscale image.

Details

In Pentlands algorithm the image irradiance is defined as:

This can be reduced using the Taylor series expansion about $p=p_0$ and $p=p_0$ and ignoring the higher order terms becomes:

which for $p_0=q_0=0$ further reduces to:

This gives the following transform identities:

By taking the Fourier transform of both sides if $E(x,y)$ yields the following:

where $F_z$ is the Fourier transform of $Z(x,y)$.

These can be rearranged, and the Inverse Fourier transform used to recover the surface $Z(x,y)$ as per the following:

Arguments

The function arguments are described in more detail below.

slant

A Real that specifies the slant value to be used by the algorithm. The slant should be a value in the range [0,π/2]. If slant is specified to must the tilt. Will defult to Inf which will trigger estimate_img_properties to be run.

tilt

A Real that specifies the tilt value to be used by the algorithm. The tilt should be a value in the range [0,2π]. If tilt is specified to must the slant. Will defult to Inf which will trigger estimate_img_properties to be run.

Example

Compute the heightmap for a synthetic image generated by generate_surface.

using Images, Makie, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(50), 1.0, [0.1,0,0.9])

#calculate the heightmap
pentland = Pentland(slant = 0.1106572211739, tilt = 0.0)
Z = pentland(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)

Reference

1. A. Pentland, "Shape Information From Shading: A Theory About Human Perception," [1988 Proceedings] Second International Conference on Computer Vision, Tampa, FL, USA, 1988, pp. 404-413. doi: 10.1109/CCV.1988.590017

Photometric(illumination_direction1::Vector{T} where T <: Real = [Inf, Inf, Inf], illumination_direction2::Vector{T} where T <: Real = [Inf, Inf, Inf], illumination_direction3::Vector{T} where T <: Real = [Inf, Inf, Inf], integration_scheme::AbstractIntegrator = Horn())

Creates a Photometric functor which attempts to produce the gradients and integrate them to retrieve a height map using three-point photometric stereo.

This algorithm solves a system of linear equation formed by varying the lighting conditions on an object while holding the viewing angle constant. Note that the ilumination direction id defined negative of most other algorithms and as such inouts may need to be corrected (see example).

Output

Returns an M by N array (matching dimensions of original image) of Float Z that represents the reconstructed height at the point and the gradients in M by N arrays of Float p and q.

Details

Let ⁠ñ₁, ñ₂ and ñ₃ define the ilumination direction let N be a 3 X 3 matrix whos rows are formed by the row vectors of the illumination direction as per:

and I₁, I₂ and I₃ be the intensity values at a point (x,y) and let Ĩ=[I₁,I₂,I₃]′ be the column vector formed by these values.

The surface normal at the point (x,y) can now be directly calculated using:

where n is the surface normal.

Arguments

The function arguments are described in more detail below.

illumination_direction1

A Vector{T} where T <: Real that specifies the first illumination vector to be used by the algorithm. The illumination_direction should be a vector of the form [x,y,z] where x,y,z are int he range [0,1]. Defults to [Inf,Inf,Inf] which triggers estimate_img_properties to be called.

illumination_direction2

A Vector{T} where T <: Real that specifies the second illumination vector to be used by the algorithm. The illumination_direction should be a vector of the form [x,y,z] where x,y,z are int he range [0,1]. Defults to [Inf,Inf,Inf] which triggers estimate_img_properties to be called.

illumination_direction3

A Vector{T} where T <: Real that specifies the third illumination vector to be used by the algorithm. The illumination_direction should be a vector of the form [x,y,z] where x,y,z are int he range [0,1]. Defults to [Inf,Inf,Inf] which triggers estimate_img_properties to be called.

integration_scheme

A keywork argument of type AbstractIntegrator which specifies which integration scheme is to be used to reconstruct the surface from the calculated normals. Defutls to Horn().

using Images, Makie, ShapeFromShading

#generate synthetic images
n₁ = [0,0,1.0]
n₂ = [1.0,0.2,1.0]
n₃ = [0.2,0.9,1.0]
img1, img2, img3 = generate_photometric(SynthSphere(50), 1.0, n₁, n₂, n₃)

#calculate the heightmap
photometric = Photometric(-n₁, -n₂, -n₃, Horn())
Z,p,q = photometric(img1, img2, img3)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)

Reference

1. R. Woodham, "Photometric Method For Determining Surface Orientation From Multiple Images", Optical Engineering, vol. 19, no. 1, 1980. doi:10.1117/12.7972479

Shah(slant::Real = Inf, tilt::Real = Inf, iterations::Int = 200)

Creates a Shah fucntor which attempts to produce a heightmap from a grayscale image using Shah's algorithm.

Under the assumption that the albedo is constant, and the surface is Lambertian, the algorithm employs discrete approximations of p and q using finite differences in order to linearize the reflectance map in terms of Z by taking the first term of its Taylor series expansion.

Output

Returns an M by N array (matching dimensions of original image) of Float Z that represents the reconstructed height at the point and the gradients in M by N arrays of Float 'p' and 'q'.

Details

Given the reflectance map for the surface defined as bellow:

where $i_x=\frac{I(1)}{I(3)}=\frac{\cos\tau\sin\sigma}{\cos\sigma}=\cos\tau \tan\sigma$ and $i_y=\frac{I(2)}{I(3)}=\frac{\sin\tau\sin\sigma}{\cos\sigma} =\sin\tau\tan\sigma$.

and p and q are discretely approximated as:

Shah linearized the function $f=E-R=0$ in terms of $Z$ in the vicinity of $Z^{k-1}$ by crating a system of linear equations which can be solved iteratively using the Jacobi iterative scheme, simplifying the Taylor series expansion to the first order to get the following:

which by letting $Z^n(x,y)=Z(x,y)$ gives:

where,

which as $Z^0(x,y)=0$, allows the algorithm to iteratively solve for $Z(x,y)$.

Arguments

The function arguments are described in more detail below.

iterations

An Int that specifies the number of iterations the algorithm is to perform. If left unspecified a default value of 200 is used.

slant

A Real that specifies the slant value to be used by the algorithm. The slant should be a value in the range [0,π/2]. If slant is specified to must the tilt. Will defult to Inf which will trigger estimate_img_properties to be run.

tilt

A Real that specifies the tilt value to be used by the algorithm. The tilt should be a value in the range [0,2π]. If tilt is specified to must the slant. Will defult to Inf which will trigger estimate_img_properties to be run.

Example

Compute the heightmap for a synthetic image generated by generate_surface.

using Images, Makie, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(50), 1.0, [0.1,0,0.9])

#calculate the heightmap
shah = Shah(slant = 0.1106572211739, tilt = 0.0)
Z = shah(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)

Reference

1. T. Ping-Sing and M. Shah, "Shape from shading using linear approximation", Image and Vision Computing, vol. 12, no. 8, pp. 487-498, 1994. doi:10.1016/0262-8856(94)90002-7

Deterministic(p::AbstractArray, q::AbstractArray, integration_factor::Real = 10, smoothness::Real = 20, β::Real = 1.0, integration_scheme::AbstractIntegrator = Horn())

Impliments a deterministic method for reconstructing the surface using a gradient descent method. This version makes use of a Optim.jl solver.

Output

Returns a Deterministic functor which can be called on an image to reconstruct the surface normals from the image.

Details

Under the assumption that the object is illuminated from directly above ([0,0,1]) and that the surface is Lambertian, the eikonal equation, shown below, can be used to describe the surface in an image.

As the surface must be integrable, the integrability condition; $p_y=q_x$ must hold on the surface. Using this a discreate itegrability function can be derived in the form:

A discreate functional for the energy can now be derived from these in the form:

Finally a smoothness term can be introduced as below. This both aids in the quality of the reconstruction and reduces the number of minima in the functional helping to ensure a good solution is found.

The problem is then minimized using gradient decsent method to attampt to find a solution minimizing the energy function. Gradient descent is terminated when $|\nabla\epsilon|<\beta\sqrt{2|\Omega|}$.

Parameters

p

An AbstractArray which contains initial solution for p, the x gradient of the surface. Better inital solution will produce better results.

q

An AbstractArray which contains initial solution for q, the y gradient of the surface. Better inital solution will produce better results.

integration_factor

A Real which controls the impact of the integration term on the energy function. Can be used to tune the final solution. Defults to 10.

smoothness

A Real which controls the impact of the smoothness term on the energy function. Can be used to tune the final solution. Defults to 20.

β

A Real which controls the value of $|\nabla\epsilon|$ which will terminate the algorithm. If algorithm terminates on first iteration value will need reducing and if algorithm fails to converge the value will need increasing. Defults to 1.0.

integration_scheme

A AbstractIntegrator which will be used to integrate the calculated gradients. Defults to Horn().

Example

The following example demonstrates the use of the Deterministic method.

using ShapeFromShading, Makie, Images

#generate synthetic image and initial solution
img = generate_surface(SynthSphere(38))
pin, qin = synthetic_gradient(SynthSphere(38))

#calculate the heightmap
deterministic = Deterministic(p=pin,q=qin, integration_factor = 10, smoothness = 50, β = 1.0, integration_scheme = Horn())
Z,p,q = deterministic(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(abs.(Z))

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
surface(r, r, Z)

Reference

[1] A. Crouzil, X. Descombes and J. Durou, "A multiresolution approach for shape from shading coupling deterministic and stochastic optimization," in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 25, no. 11, pp. 1416-1421, Nov. 2003. doi: 10.1109/TPAMI.2003.1240116

DeterministicCustom(p::AbstractArray, q::AbstractArray, integration_factor::Real = 10, smoothness::Real = 20, β::Real = 1.0, integration_scheme::AbstractIntegrator = Horn())

Same as Deterministic exept impliments a custom gradient descent and line search algorithm.

Example

The following example demonstrates the use of the DeterministicCustom method.

using ShapeFromShading, Makie, Images

#generate synthetic image and initial solution
img = generate_surface(SynthSphere(38))
pin, qin = synthetic_gradient(SynthSphere(38))

#calculate the heightmap
deterministicCustom = DeterministicCustom(p=pin,q=qin, integration_factor = 10, smoothness = 50, β = 1.0, integration_scheme = Horn())
Z,p,q = deterministicCustom(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(abs.(Z))

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:4
surface(r, r, Z)

Reference

[1] A. Crouzil, X. Descombes and J. Durou, "A multiresolution approach for shape from shading coupling deterministic and stochastic optimization," in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 25, no. 11, pp. 1416-1421, Nov. 2003. doi: 10.1109/TPAMI.2003.1240116