agd.ExportedCode.Notebooks_NonDiv.ShapeFromShading

 1# Code automatically exported from notebook ShapeFromShading.ipynb in directory Notebooks_NonDiv
 2# Do not modify
 3import numpy as np
 4import matplotlib.pyplot as plt
 5
 6def EvalScheme(cp,u,uc,params):
 7    """
 8    Evaluates the (piecewise) quadratic equation defining the numerical scheme.
 9    Inputs :
10     - uc : plays the role of λ
11    """
12    alpha,beta,gamma,h = params
13    wx = np.roll(u,-1,axis=0)
14    wy = np.roll(u,-1,axis=1)
15    vx = np.minimum(wx,np.roll(u,1,axis=0))
16    vy = np.minimum(wy,np.roll(u,1,axis=1))
17
18    return (cp*np.sqrt(1+(np.maximum(0,uc-vx)**2+np.maximum(0,uc-vy)**2)/h**2) +
19            alpha*(uc-wx)/h+beta*(uc-wy)/h-gamma)
20
21def LocalSolve(cp,vx,vy,wx,wy,params):
22    """
23    Solve the (piecewise) quadratic equation defining the numerical scheme.
24    Output: solution λ.
25    """
26    alpha,beta,gamma,h = params
27    # Trying with two active positive parts 
28    
29    # Quadratic equation coefficients.
30    # a lambda^2 - 2 b lambda + c =0
31    a = (2.*cp**2 - (alpha+beta)**2)
32    b = cp**2 *(vx+vy) - (alpha+beta)*(alpha*wx+beta*wy+h*gamma)
33    c = cp**2*(h**2+vx**2+vy**2)-(gamma*h+alpha*wx+beta*wy)**2
34    
35    delta = b**2 - a*c
36    good = np.logical_and(delta>=0,a!=0)
37    u = 0*cp;
38    # TODO : Is that the correct root selection ?
39    u[good] = (b[good]+np.sqrt(delta[good]))/a[good] 
40
41    vmax = np.maximum(vx,vy)
42    good = np.logical_and(good,u>=vmax)
43    
44    # Trying with one active positive part
45    # TODO : restrict computations to not good points to save cpu time ?
46    
47    vmin = np.minimum(vx,vy)
48    a = (cp**2 - (alpha+beta)**2)
49    b = cp**2 *vmin - (alpha+beta)*(alpha*wx+beta*wy+h*gamma)
50    c = cp**2*(h**2+vmin**2)-(gamma*h+alpha*wx+beta*wy)**2
51
52    delta = b**2 - a*c
53    ggood = np.logical_and(np.logical_and(delta>=0,a!=0), 1-good)
54    u[ggood] = (b[ggood] +np.sqrt(delta[ggood]))/a[ggood]
55    
56    good = np.logical_or(good,np.logical_and(ggood,u>=vmin))
57    
58    # No active positive part
59    # equation becomes linear, a lambda - b = 0
60    a = alpha+beta+0.*cp
61    b = alpha*wx+beta*wy +gamma*h - cp*h
62    bad = np.logical_not(good)
63    u[bad]=b[bad]/a[bad]
64    return u
65    
66def JacobiIteration(u,Omega,c,params):
67    """
68    One Jacobi iteration, returning the pointwise solution λ to the numerical scheme.
69    """
70    wx = np.roll(u,-1,axis=0)
71    wy = np.roll(u,-1,axis=1)
72    vx = np.minimum(wx,np.roll(u,1,axis=0))
73    vy = np.minimum(wy,np.roll(u,1,axis=1))
74    
75#    sol=LocalSolve(c,vx,vy,wx,wy,params)
76    sol = u+LocalSolve(c,vx-u,vy-u,wx-u,wy-u,params)
77    u[Omega] = sol[Omega]
78
79def OneBump(x,y):
80    bump = 0.5-3.*((x-0.5)**2+(y-0.5)**2)
81    return np.maximum(bump, np.zeros_like(x))
82
83def ThreeBumps(x,y):
84    bump1 = 0.3-3*((x-0.4)**2+(y-0.5)**2)
85    bump2 = 0.25-3*((x-0.65)**2+(y-0.6)**2)
86    bump3 = 0.25-3*((x-0.6)**2+(y-0.35)**2)
87    return np.maximum.reduce([bump1,bump2,bump3,np.zeros_like(bump1)])
88
89def Volcano(x,y):
90    r = np.sqrt((x-0.5)**2+(y-0.5)**2)
91    volcano = 0.05+1.5*(1+x)*(r**2-6*r**4)
92    return np.maximum(volcano, np.zeros_like(x))
93
94def GenerateRHS(height,params):
95    α,β,γ,h = params
96    hx,hy = np.gradient(height,h)
97    Intensity = (α*hx+β*hy+γ)/np.sqrt(1+hx**2+hy**2)
98    Omega = height>0
99    return Intensity,Omega