This model is much more flexible than commonly used quadratic regression, which lets \(h_{uv}(x_u,x_v)=\gamma _{uv}x_ux_v\) . To model \(P_{uv, 1}, P_{uv, 2}\) , \(N_{uv, 1}, N_{uv, 2}\) 
as flexible one-dimensional non-negative functions, we rely on squaring B-spline expansions [1]} as follows:
\(\begin{aligned}P_{uv, 1}(x_u) = \lbrace \mathbf {s}_u(x_u)^{\text{T}} \theta _{uv, 1}\rbrace ^2 &, \quad P_{uv, 2}(x_v) = \lbrace \mathbf {s}_v(x_v)^{\text{T}} \phi _{uv, 1}\rbrace ^2\\N_{uv, 1}(x_u) = \lbrace \mathbf {s}_u(x_u)^{\text{T}} \theta _{uv, 2}\rbrace ^2 &, \quad N_{uv, 2}(x_v) = \lbrace \mathbf {s}_v(x_v)^{\text{T}} \phi _{uv, 2}\rbrace ^2,\end{aligned}\) 
