<p>Rather than fit bilinear patches to your sample points, which are only C0 continuous<span class="math-container">$^\spadesuit$</span> may I instead suggest fitting the surface using Catmull-Rom splines? These will give you C2 continuity. 
Catmull-Rom is an interpolating spline which means it will pass through your sample points, as did the bilinear surface, but will do so smoothly from patch to patch.</p>

<p><span class="math-container">$^\spadesuit$</span> <em>because the first derivatives won't generally match at the boundaries of your patches</em> </p>

<h3>Background/Assumption</h3>

<p>You have a function <span class="math-container">$$F:\Re^2\to\Re^3$$</span> of the form  <span class="math-container">$$F(s,t)=[s,t,f(s,t)]$$</span> where <span class="math-container">$f:\Re^2\to\Re$</span>. The function, <span class="math-container">$f$</span>, can be considered a 'black box' in the sense that we might not be given it explicitly.</p>

<p><span class="math-container">$F$</span> has been regularly sampled as <span class="math-container">$$S_{i,j}=F(x_{scale}.i, y_{scale}.j, f(x_{scale}.i,y_{scale}.j))$$</span> where <span class="math-container">$i,j\in\mathbb{Z}$</span>.  You then want interpolate between these sample points using piecewise bicubic Catmull-Rom patches.</p>

<h3>Constructing patches</h3>

<p>To form a particular bicubic Catmull-Rom patch, <span class="math-container">$C_{i,j}$</span>, which has a given 2x2 set of samples <span class="math-container">$\{S_{i,j},  S_{i+1,j}, S_{i,j+1},  S_{i+1,j+1}\}$</span> at the corners, you need the set of 4x4 samples formed from the 2x2 at the centre and the surrounding 12 samples, i.e.<span class="math-container">$\{S_{i-1,j-1}$</span> <em>through</em> <span class="math-container">$S_{i+2,j+2}\}$</span>,  as shown below.
(My apologies that this is 'hand placed' so is far from ideal)
<a href="https://i.sstatic.net/M48A3.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/M48A3.png" alt="enter image description here"></a></p>

<p>If we just consider the behaviour of the surface at sample, <span class="math-container">$S_{i,j}$</span>, then the partial derivative in the <span class="math-container">$i$</span> direction is given by 
<span class="math-container">$$\frac{\partial C_{ij}}{\partial i}=\frac{S_{i+1,j}-S_{i-1,j}}{2}$$</span>
and similarly, in the <span class="math-container">$j$</span> direction by
<span class="math-container">$$\frac{\partial C_{ij}}{\partial j}=\frac{S_{i,j+1}-S_{i,j-1}}{2}$$</span></p>

<p>Remembering these are vectors, you can use the cross product to generate a surface normal at any of the <span class="math-container">$S_{i,j}$</span> points<span class="math-container">$^\clubsuit$</span>. (Though you may need an extra set of boundary samples)</p>

<blockquote>
  <p><em>Ali: You mentioned that both partial derivatives are vectors(in i,j directions). could you write each vector components plz</em></p>
</blockquote>

<p>Certainly: If we say <span class="math-container">$S_{i,j}=[x_{i,j},y_{i,j},z_{i,j}]^T$</span> then <span class="math-container">$$\frac{\partial C_{ij}}{\partial i}=\frac{1}{2}.[(x_{i+1,j}-x_{i-1,j}),(y_{i+1,j}-y_{i-1,j}),(z_{i+1,j}-z_{i-1,j})]^T$$</span>
<span class="math-container">$$\frac{\partial C_{ij}}{\partial j}=\frac{1}{2}.[(x_{i,j+1}-x_{i,j-1}),(y_{i,j+1}-y_{i,j-1}),(z_{i,j+1}-z_{i,j-1})]^T$$</span></p>

<p>If you've sampled the surface 'finely' enough you could probably choose to render the surface by just using 2 or 4 triangles per original patch. </p>

<p>Alternatively, you can recursively subdivide your quad mesh to get a finer set of triangles. To double the number of points in the 'i' direction, simply insert a new point between each <span class="math-container">$S_{i,j}$</span> and <span class="math-container">$S_{i+1,j}$</span> that equals ...
<span class="math-container">$$\frac{9(S_{i,j}+S_{i+1,j}) - (S_{i-1,j}+S_{i+2,j})}{16}$$</span>
 Once you have the finer grid in 'i' simply repeat with 'j' with the appropriate equivalent set of points. You can repeat this ad nauseam to get whatever triangle resolution you desire ;-)</p>

<p><span class="math-container">$^\clubsuit$</span> (Copying some text from the comments below): I've made a couple of valid assumptions here.  Noting that if a given first (partial) derivative is non-zero (which is guaranteed in your case because of the way you construct the surface), then it also is locally tangential to that surface. Further, we know your two partial derivates aren't parallel, so the cross product will construct a valid normal. 
FWIW For some arrangements of control points, i.e. in more general use cases, you can get so-called 'degenerate' surfaces where you can't naively apply the above rules, but you're safe here. </p>

CC BY-SA 4.0
2019-11-06T23:41:44.917Z
Thanks Simon. You mentioned that both partial derivatives are vectors(in i,j directions). could you write each vector components plz.
ali
CC BY-SA 4.0
2019-11-07T22:43:07.147Z
Thanks for the update. My math is not very good. I thought derivative is the rate of change of "function value" with respect to one of its arguments(unless it is a vector-valued function and you are assuming that the point S(x,y,z) is the function value); the gradient is a vector but not the partial derivatives, is that correct?
ali
CC BY-SA 4.0
2019-11-08T10:56:07.927Z
*"I thought derivative is the rate of change of "function value" with respect to one of its arguments"*   That is indeed correct.  But for a surface, _if_ a given first (partial) derivative is non-zero (which is guaranteed in your case because of the way you construct the surface), then it also is locally tangential to the surface. Further, we know *your* two partial derivates _aren't_ parallel, so the cross product will construct a valid normal. (FWIW For some arrangements of control points, in more general use cases, you can get so-called 'degenerate' surfaces, but you're safe here)
Simon F
CC BY-SA 4.0
2019-11-08T14:24:22.910Z
"but for a surface, if a given first (partial) derivative is non-zero then it also is locally tangential...". I am not making a fuss :) but your argument already assumes partial derivative to be a vector(otherwise it won't be tangential) whereas I think it is scalar.
ali
CC BY-SA 4.0
2019-11-08T17:22:00.090Z
But you are building a surface in 3D.  I'll add another note
Simon F
CC BY-SA 4.0
2019-11-08T19:45:41.470Z
Perfect. It's all clear now. I had missed the capital F function.
ali
CC BY-SA 4.0
2019-11-06T14:59:45.760Z
2019-11-08T17:44:55.533Z
2019-11-08T17:44:55.533Z
Simon F
<p>I am trying to calculate the derivatives of a sampled data set by following an example in <a href="https://books.google.co.uk/books/about/Data_Visualization.html?id=GPwYCwAAQBAJ&amp;printsec=frontcover&amp;source=kp_read_button&amp;redir_esc=y#v=onepage&amp;q=%20gllightfv&amp;f=false" rel="nofollow noreferrer">this book</a>. The data is a set of 2d quads(green area in below) and their vertices(4 each). The data is in the following example is sampled from <span class="math-container">$z=f(x,y)=e^{-(x^2+y^2)}$</span> but, in general, the underlying function might not be known.</p>

<p><a href="https://i.sstatic.net/ySAtb.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ySAtb.png" alt="enter image description here"></a> </p>

<p>The book states that the derivatives of any points are calculated by:</p>

<p><a href="https://i.sstatic.net/iHqWJ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/iHqWJ.png" alt="enter image description here"></a></p>

<p><a href="https://i.sstatic.net/KdZ65.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/KdZ65.png" alt="enter image description here"></a> are the quadrilateral basis functions:</p>

<p><a href="https://i.sstatic.net/opYtc.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/opYtc.png" alt="enter image description here"></a></p>

<p>and</p>

<p><a href="https://i.sstatic.net/xrwjd.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/xrwjd.png" alt="enter image description here"></a></p>

<p>I can calculate df/dx and df/dy. </p>

<p>But where I struggle is calculating df/dz. I assume</p>

<p><a href="https://i.sstatic.net/YKxwY.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/YKxwY.png" alt="enter image description here"></a> </p>

<p>and therefore df/dz = 0. But this cannot be correct since df/dz should be 1(by making the above analytic function implicit and differentiating it for z).</p>

CC BY-SA 4.0
2019-10-28T08:47:30.970Z
Can I ask some questions for clarification?So you have a function, f(x,y), that you have sampled on a regular grid of (x,y) locations, to create a set of *bilinear* patches?  Are those corner points, therefore, of the form (xi, yi, f(xi, yi))? That is Z is pointing up the screen?
Simon F
CC BY-SA 4.0
2019-10-28T13:36:43.940Z
Exactly as you said
ali
CC BY-SA 4.0
2019-10-28T14:08:33.923Z
(Sigh. I wish the comments supported equations.)  I don't understand why you are trying to manufacture a "df/dz".  You have a surface, which to reduce confusion we'll call G(x,y) = (x, y, f(x,y)), and you need to compute on that surface the two partial derivatives (vectors),  dG/dx and dG/dy. These are (1, 0, df/dx) and (0, 1, df/dy) respectively,   where, unless I'm mistaken, df/dx is -2*x*f(x,y) etc  If you just want your piecewise surface to look smooth say using Gouraud shading, evaluate the tangents at the sample points, create your normals (via cross products) and shade.
Simon F
CC BY-SA 4.0
2019-10-28T18:52:56.707Z
Thanks Simon. What you said is correct. But you don't always have the analytic formula and all you have is sampled data. The second formula above(taken from the book) I assume is supposed to calculate partial deliveries for with respect to each x,y,z. It does but not for z as expected to be 1
ali
CC BY-SA 4.0
2019-10-29T09:00:51.280Z
Ok. Would it be ok to edit your original question to make that clearer?  As for constructing derivatives etc, assuming your surface is meant to be at least C1 continuous, if it were up to me, I'd use a different underlying piecewise interpolation.
Simon F
CC BY-SA 4.0
2019-10-30T15:56:44.480Z
@SimonF comments do support tex equations. For example: $\frac{\alpha_1}{\gamma}$
joojaa
CC BY-SA 4.0
2019-11-05T11:30:41.127Z
@SimonF, thanks. You mentioned: " I'd use a different underlying piecewise interpolation". So would you mind explaining what interpolation you use in this case.
ali
CC BY-SA 4.0
2019-11-05T14:15:53.900Z
I'll try to write something tomorrow.
Simon F
<sampling><interpolation>
computing derivatives of sampled data
ali
# computing derivatives of sampled data

(Asked by: username_0 on 2019-10-26)

I am trying to calculate the derivatives of a sampled data set by following an example in [this book](https://books.google.co.uk/books/about/Data_Visualization.html?id=GPwYCwAAQBAJ&printsec=frontcover&source=kp_read_button&redir_esc=y#v=onepage&q=%20gllightfv&f=false). The data is a set of 2d quads(green area in below) and their vertices(4 each). The data is in the following example is sampled from $z=f(x,y)=e^{-(x^2+y^2)}$ but, in general, the underlying function might not be known.

[![enter image description here](https://i.sstatic.net/ySAtb.png)](https://i.sstatic.net/ySAtb.png)

The book states that the derivatives of any points are calculated by:

[![enter image description here](https://i.sstatic.net/iHqWJ.png)](https://i.sstatic.net/iHqWJ.png)

[![enter image description here](https://i.sstatic.net/KdZ65.png)](https://i.sstatic.net/KdZ65.png) are the quadrilateral basis functions:

[![enter image description here](https://i.sstatic.net/opYtc.png)](https://i.sstatic.net/opYtc.png)

and

[![enter image description here](https://i.sstatic.net/xrwjd.png)](https://i.sstatic.net/xrwjd.png)

I can calculate df/dx and df/dy.

But where I struggle is calculating df/dz. I assume

[![enter image description here](https://i.sstatic.net/YKxwY.png)](https://i.sstatic.net/YKxwY.png)

and therefore df/dz = 0. But this cannot be correct since df/dz should be 1(by making the above analytic function implicit and differentiating it for z).

**Tags:** sampling, interpolation

> **Comments:**
> **username_1**: Can I ask some questions for clarification?So you have a function, f(x,y), that you have sampled on a regular grid of (x,y) locations, to create a set of *bilinear* patches?  Are those corner points, therefore, of the form (xi, yi, f(xi, yi))? That is Z is pointing up the screen?  
> **username_0**: Exactly as you said  
> **username_1**: (Sigh. I wish the comments supported equations.)  I don't understand why you are trying to manufacture a "df/dz".  You have a surface, which to reduce confusion we'll call G(x,y) = (x, y, f(x,y)), and you need to compute on that surface the two partial derivatives (vectors),  dG/dx and dG/dy. These are (1, 0, df/dx) and (0, 1, df/dy) respectively,   where, unless I'm mistaken, df/dx is -2*x*f(x,y) etc  If you just want your piecewise surface to look smooth say using Gouraud shading, evaluate the tangents at the sample points, create your normals (via cross products) and shade.  
> **username_0**: Thanks username_1. What you said is correct. But you don't always have the analytic formula and all you have is sampled data. The second formula above(taken from the book) I assume is supposed to calculate partial deliveries for with respect to each x,y,z. It does but not for z as expected to be 1  
> **username_1**: Ok. Would it be ok to edit your original question to make that clearer?  As for constructing derivatives etc, assuming your surface is meant to be at least C1 continuous, if it were up to me, I'd use a different underlying piecewise interpolation.  
> **username_2**: @username_1 comments do support tex equations. For example: $\frac{\alpha_1}{\gamma}$  
> **username_0**: @username_1, thanks. You mentioned: " I'd use a different underlying piecewise interpolation". So would you mind explaining what interpolation you use in this case.  
> **username_1**: I'll try to write something tomorrow.  

---

## 1 answers

### Answer by username_1 on 2019-11-06. Score: 1 | Accepted answer

Rather than fit bilinear patches to your sample points, which are only C0 continuous$^\spadesuit$ may I instead suggest fitting the surface using Catmull-Rom splines? These will give you C2 continuity.
Catmull-Rom is an interpolating spline which means it will pass through your sample points, as did the bilinear surface, but will do so smoothly from patch to patch.

$^\spadesuit$ *because the first derivatives won't generally match at the boundaries of your patches*

### Background/Assumption

You have a function $$F:\Re^2\to\Re^3$$ of the form $$F(s,t)=[s,t,f(s,t)]$$ where $f:\Re^2\to\Re$. The function, $f$, can be considered a 'black box' in the sense that we might not be given it explicitly.

$F$ has been regularly sampled as $$S\_{i,j}=F(x\_{scale}.i, y\_{scale}.j, f(x\_{scale}.i,y\_{scale}.j))$$ where $i,j\in\mathbb{Z}$. You then want interpolate between these sample points using piecewise bicubic Catmull-Rom patches.

### Constructing patches

To form a particular bicubic Catmull-Rom patch, $C\_{i,j}$, which has a given 2x2 set of samples $\{S\_{i,j}, S\_{i+1,j}, S\_{i,j+1}, S\_{i+1,j+1}\}$ at the corners, you need the set of 4x4 samples formed from the 2x2 at the centre and the surrounding 12 samples, i.e.$\{S\_{i-1,j-1}$ *through* $S\_{i+2,j+2}\}$, as shown below.
(My apologies that this is 'hand placed' so is far from ideal)
[![enter image description here](https://i.sstatic.net/M48A3.png)](https://i.sstatic.net/M48A3.png)

If we just consider the behaviour of the surface at sample, $S\_{i,j}$, then the partial derivative in the $i$ direction is given by
$$\frac{\partial C\_{ij}}{\partial i}=\frac{S\_{i+1,j}-S\_{i-1,j}}{2}$$
and similarly, in the $j$ direction by
$$\frac{\partial C\_{ij}}{\partial j}=\frac{S\_{i,j+1}-S\_{i,j-1}}{2}$$

Remembering these are vectors, you can use the cross product to generate a surface normal at any of the $S\_{i,j}$ points$^\clubsuit$. (Though you may need an extra set of boundary samples)

> *Ali: You mentioned that both partial derivatives are vectors(in i,j directions). could you write each vector components plz*

Certainly: If we say $S\_{i,j}=[x\_{i,j},y\_{i,j},z\_{i,j}]^T$ then $$\frac{\partial C\_{ij}}{\partial i}=\frac{1}{2}.[(x\_{i+1,j}-x\_{i-1,j}),(y\_{i+1,j}-y\_{i-1,j}),(z\_{i+1,j}-z\_{i-1,j})]^T$$
$$\frac{\partial C\_{ij}}{\partial j}=\frac{1}{2}.[(x\_{i,j+1}-x\_{i,j-1}),(y\_{i,j+1}-y\_{i,j-1}),(z\_{i,j+1}-z\_{i,j-1})]^T$$

If you've sampled the surface 'finely' enough you could probably choose to render the surface by just using 2 or 4 triangles per original patch.

Alternatively, you can recursively subdivide your quad mesh to get a finer set of triangles. To double the number of points in the 'i' direction, simply insert a new point between each $S\_{i,j}$ and $S\_{i+1,j}$ that equals ...
$$\frac{9(S\_{i,j}+S\_{i+1,j}) - (S\_{i-1,j}+S\_{i+2,j})}{16}$$
Once you have the finer grid in 'i' simply repeat with 'j' with the appropriate equivalent set of points. You can repeat this ad nauseam to get whatever triangle resolution you desire ;-)

$^\clubsuit$ (Copying some text from the comments below): I've made a couple of valid assumptions here. Noting that if a given first (partial) derivative is non-zero (which is guaranteed in your case because of the way you construct the surface), then it also is locally tangential to that surface. Further, we know your two partial derivates aren't parallel, so the cross product will construct a valid normal.
FWIW For some arrangements of control points, i.e. in more general use cases, you can get so-called 'degenerate' surfaces where you can't naively apply the above rules, but you're safe here.

> **Comments:**
> **username_0**: Thanks username_1. You mentioned that both partial derivatives are vectors(in i,j directions). could you write each vector components plz.  
> **username_0**: Thanks for the update. My math is not very good. I thought derivative is the rate of change of "function value" with respect to one of its arguments(unless it is a vector-valued function and you are assuming that the point S(x,y,z) is the function value); the gradient is a vector but not the partial derivatives, is that correct?  
> **username_1**: *"I thought derivative is the rate of change of "function value" with respect to one of its arguments"*   That is indeed correct.  But for a surface, _if_ a given first (partial) derivative is non-zero (which is guaranteed in your case because of the way you construct the surface), then it also is locally tangential to the surface. Further, we know *your* two partial derivates _aren't_ parallel, so the cross product will construct a valid normal. (FWIW For some arrangements of control points, in more general use cases, you can get so-called 'degenerate' surfaces, but you're safe here)  
> **username_0**: "but for a surface, if a given first (partial) derivative is non-zero then it also is locally tangential...". I am not making a fuss :) but your argument already assumes partial derivative to be a vector(otherwise it won't be tangential) whereas I think it is scalar.  
> **username_1**: But you are building a surface in 3D.  I'll add another note  
> **username_0**: Perfect. It's all clear now. I had missed the capital F function.  

---

