agd.Eikonal.HFM_CUDA.ShapeFromShading

This module solves the shape from shading PDE, with a non-vertical light source. The PDE reads (alpha u_x + beta u_y + gamma) / sqrt(1+u_x^2+u_y^2) = rhs where alpha, beta, gamma are parameters related to the sun direction.

View Source

 1# Copyright 2021 Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
 2# Distributed WITHOUT ANY WARRANTY. Licensed under the Apache License, Version 2.0, see http://www.apache.org/licenses/LICENSE-2.0
 3
 4"""
 5This module solves the shape from shading PDE, with a non-vertical light source. 
 6The PDE reads 
 7(alpha u_x + beta u_y + gamma) / sqrt(1+u_x^2+u_y^2) = rhs
 8where alpha, beta, gamma are parameters related to the sun direction.
 9"""
10
11import os
12import numpy as np
13import cupy as cp
14from . import cupy_module_helper
15from ... import FiniteDifferences as fd
16
17def EvalScheme(cp,u,params,uc=None,mask=None):
18	"""
19	Evaluates the (piecewise) quadratic equation defining the numerical scheme.
20	Inputs :
21	 - uc : plays the role of λ
22	 - params : alpha,beta,gamma,h (grid scale)
23	"""
24	alpha,beta,gamma,h = params
25	if uc is None: uc=u
26	
27	sa,alpha=int(np.sign(alpha)),np.abs(alpha)
28	sb,beta =int(np.sign(beta)), np.abs(beta)
29
30	wx = np.roll(u,-sa,axis=0)
31	wy = np.roll(u,-sb, axis=1)
32	vx = np.minimum(wx,np.roll(u,sa,axis=0))
33	vy = np.minimum(wy,np.roll(u,sb,axis=1))
34
35	residue = (cp*np.sqrt(1+(np.maximum(0,uc-vx)**2+np.maximum(0,uc-vy)**2)/h**2) 
36		+ alpha*(uc-wx)/h+beta*(uc-wy)/h-gamma)
37	return residue if mask is None else np.where(mask,residue,0.)
38
39
40def Solve(rhs,mask,u0,params,niter=300,traits=None):
41	"""
42	Iterative solver for the shape from shading equation.
43	"""
44	float_t = np.float32
45	int_t = np.int32
46	boolatom_t = np.uint8
47	assert len(params)==4
48
49	# Reshape data
50	rhs = cp.asarray(rhs,dtype=float_t)
51	mask = cp.asarray(mask,dtype=boolatom_t)
52	u0 = cp.asarray(u0,dtype=float_t)
53	shape = rhs.shape
54	assert mask.shape==shape and u0.shape==shape
55
56	traits_default = {'side_i':8,'niter_i':8}
57	if traits is None: traits = traits_default
58	else: traits = {**traits_default,**traits}
59
60	shape_i = (traits['side_i'],)*2
61	rhs,mask,u0 = [fd.block_expand(e,shape_i,constant_values=v) 
62		for e,v in [(rhs,np.nan),(mask,False),(u0,0.)] ]
63	
64	# Find active blocks
65	shape_o = rhs.shape[:2]
66	update_o = np.flatnonzero(np.any(mask,axis=(-2,-1))).astype(int_t)
67
68	# Setup the cuda module
69	cuda_path = os.path.join(os.path.dirname(os.path.realpath(__file__)),"cuda")
70	date_modified = cupy_module_helper.getmtime_max(cuda_path)
71	source = cupy_module_helper.traits_header(traits)
72
73	source += ['#include "Kernel_ShapeFromShading.h"',
74	f"// Date cuda code last modified : {date_modified}"]
75	cuoptions = ("-default-device", f"-I {cuda_path}") 
76	source = "\n".join(source)
77	module = cupy_module_helper.GetModule(source,cuoptions)
78
79	def SetCst(name,var,dtype):cupy_module_helper.SetModuleConstant(module,name,var,dtype)
80	SetCst('params',params,float_t)
81	SetCst('shape_o',shape_o,int_t)
82	SetCst('shape_tot', np.array(shape_o)*np.array(shape_i), int_t)
83
84	sfs = module.get_function('JacobiUpdate')
85
86	# Call the kernel
87	rhs,mask,u,update_o = [cp.ascontiguousarray(e) for e in (rhs,mask,u0,update_o)]
88	for i in range(niter):
89		sfs((update_o.size,),shape_i,(u,rhs,mask,update_o))
90
91	return fd.block_squeeze(u,shape)
92
93

def EvalScheme(cp, u, params, uc=None, mask=None): View Source

18def EvalScheme(cp,u,params,uc=None,mask=None):
19	"""
20	Evaluates the (piecewise) quadratic equation defining the numerical scheme.
21	Inputs :
22	 - uc : plays the role of λ
23	 - params : alpha,beta,gamma,h (grid scale)
24	"""
25	alpha,beta,gamma,h = params
26	if uc is None: uc=u
27	
28	sa,alpha=int(np.sign(alpha)),np.abs(alpha)
29	sb,beta =int(np.sign(beta)), np.abs(beta)
30
31	wx = np.roll(u,-sa,axis=0)
32	wy = np.roll(u,-sb, axis=1)
33	vx = np.minimum(wx,np.roll(u,sa,axis=0))
34	vy = np.minimum(wy,np.roll(u,sb,axis=1))
35
36	residue = (cp*np.sqrt(1+(np.maximum(0,uc-vx)**2+np.maximum(0,uc-vy)**2)/h**2) 
37		+ alpha*(uc-wx)/h+beta*(uc-wy)/h-gamma)
38	return residue if mask is None else np.where(mask,residue,0.)

Evaluates the (piecewise) quadratic equation defining the numerical scheme. Inputs :

uc : plays the role of λ
params : alpha,beta,gamma,h (grid scale)

def Solve(rhs, mask, u0, params, niter=300, traits=None): View Source

41def Solve(rhs,mask,u0,params,niter=300,traits=None):
42	"""
43	Iterative solver for the shape from shading equation.
44	"""
45	float_t = np.float32
46	int_t = np.int32
47	boolatom_t = np.uint8
48	assert len(params)==4
49
50	# Reshape data
51	rhs = cp.asarray(rhs,dtype=float_t)
52	mask = cp.asarray(mask,dtype=boolatom_t)
53	u0 = cp.asarray(u0,dtype=float_t)
54	shape = rhs.shape
55	assert mask.shape==shape and u0.shape==shape
56
57	traits_default = {'side_i':8,'niter_i':8}
58	if traits is None: traits = traits_default
59	else: traits = {**traits_default,**traits}
60
61	shape_i = (traits['side_i'],)*2
62	rhs,mask,u0 = [fd.block_expand(e,shape_i,constant_values=v) 
63		for e,v in [(rhs,np.nan),(mask,False),(u0,0.)] ]
64	
65	# Find active blocks
66	shape_o = rhs.shape[:2]
67	update_o = np.flatnonzero(np.any(mask,axis=(-2,-1))).astype(int_t)
68
69	# Setup the cuda module
70	cuda_path = os.path.join(os.path.dirname(os.path.realpath(__file__)),"cuda")
71	date_modified = cupy_module_helper.getmtime_max(cuda_path)
72	source = cupy_module_helper.traits_header(traits)
73
74	source += ['#include "Kernel_ShapeFromShading.h"',
75	f"// Date cuda code last modified : {date_modified}"]
76	cuoptions = ("-default-device", f"-I {cuda_path}") 
77	source = "\n".join(source)
78	module = cupy_module_helper.GetModule(source,cuoptions)
79
80	def SetCst(name,var,dtype):cupy_module_helper.SetModuleConstant(module,name,var,dtype)
81	SetCst('params',params,float_t)
82	SetCst('shape_o',shape_o,int_t)
83	SetCst('shape_tot', np.array(shape_o)*np.array(shape_i), int_t)
84
85	sfs = module.get_function('JacobiUpdate')
86
87	# Call the kernel
88	rhs,mask,u,update_o = [cp.ascontiguousarray(e) for e in (rhs,mask,u0,update_o)]
89	for i in range(niter):
90		sfs((update_o.size,),shape_i,(u,rhs,mask,update_o))
91
92	return fd.block_squeeze(u,shape)

Iterative solver for the shape from shading equation.