A non centered method to compute velocity derivative

I have three vector vx, vy and vz, these vectors contain the values of velocity in every point of my 3d domain and in every time frame. I have a 3d domain and I'm computing derivatives of the velocity in every direction with the Sobel Operators in this way:

s1(:,:,1) = [1 2 1; 0 0 0; -1 -2 -1];

s1(:,:,2) = [2 4 2; 0 0 0; -2 -4 -2];

s1(:,:,3) = [1 2 1; 0 0 0; -1 -2 -1];

s2(:,:,1) = [1 0 -1; 2 0 -2; 1 0 -1];

s2(:,:,2) = [2 0 -2; 4 0 -4; 2 0 -2];

s2(:,:,3) = [1 0 -1; 2 0 -2; 1 0 -1];

s3(:,:,1) = [1 2 1; 2 4 2; 1 2 1];

s3(:,:,3) = -[1 2 1; 2 4 2; 1 2 1];

s1 = s1./32;

s2 = s2./32;

s3 = s3./32;

sx = s2;

sy = s1;

sz = s3;

dvxdx = imfilter(vx,sx,'conv'); Derivative of the x-component of the velocity along direction x

dvxdy = imfilter(vx,sy,'conv'); Derivative of the x-component of the velocity along direction y

dvxdz = imfilter(vx,sz,'conv'); Derivative of the x-component of the velocity along direction z

Since this method is a centered method, it doesn't work very well when I want to compute the derivatives on the boundary of my domain. Does anyone know a not centered method to compute derivatives of the velocity along all the directions given the three vectors vx, vy, vz?

Answer 1

[William Rose]

@Jacopo Rossi,

[edited to crrect typographical errors]

If you use a non-centered-to-the-right method, you will have issues on the right edge. If you use non-centered-to-the-left method, you will have issue on the left edge. You have several options:

1. Apply the technique to the inner region of your volume, far enough away from the edges that you won't have problems. Then the output matrix will be smaller than the original input.
2. Write code to handle the edges as a special case. This can get complicated, and it is not obvious how exactly you should handle the special case. If you don't compute the contributions to the convolution that come from off-edge, it is like assuming the off-edge values are zero, which can give bad results for derivatives.
3. Extend the matrix beyond the edges with an "upside down and backwards" extension. See 1-D example below. Then you filter the extended matrix. Then you keep only the central region of the filtered extended signal. This is what I would do. This is used for derivatives and other filtering of 1-D signals.

Example: Extend a 1D signal by m on each side, so that it can be filtered with a (2m+1)-point-wide centered filter:

t=(0:20)/20; x=cos(2*pi*t)-2*t;     %any waveform or recording
m=3;                     %number of points to add on each end
n=length(x);
%next: compute initial value of extended x
xext=[zeros(1,m),x,zeros(1,m)];
%next: compute the upside-down-and-backwards extensions
for i=1:m
    xext(i)=2*x(1)-x(m+2-i);         %before the beginning
    xext(n+m+i)=2*x(end)-x(end-i);   %tail: "put your behind in your past"-Puumba
end
text=(-m:20+m)/20;       %extended time vector, for plotting
plot(t,x,'-b.',text,xext,'-rx'); legend('original','extended'); grid on

You can extend this idea to 3D.

[William Rose]

@Jacopo Rossi,

I am attaching a script which does what you want.

The first section computes a velocity field, vx(x,y,z), over a volume:

Slices through the velocity field are plotted.

Then the Sobel filters are used to compute the three spatial derivatves. Slices of the derivatives are plotted. The plots (below) show that the derivatives behave badly at the edges of the volume.

Then the velocity field is extended by one sample, all around the sample volume. First, the six faces of the volume are extended by the upside-down-and-backward technique. Then the 12 edges are filled in, by fitting a plane through the three nearest points to each unknown point, and extending the plane to the unknown position. Then the 8 corner points are filled in, by a similar method.

Now we have an extended velocity field. The Sobel filters are convolved across the extended volume to compute the 3 spatial derivatives. The central region of the derivatives are extracted. These correspond to the original sample volume. Slices of the derivatives are plotted (below), which correspond exactly to the earlier slices of derivatives. These plots behave very nicely at the edges. Notice that the vertical axis scales below are different from the vertical axis scales in the previous set of plots.

, in this example, has simple spatial derivatives:

, ,

The plots above indicate that the Sobel operators are working correctly. The vertical scales of the derivative plots above are not correct, because I have not adjusted for the grid spacing. The factors of 1/32 in your operators would have to be altered to get the correct scales for the derivatives in this example.

[Jacopo Rossi]

Hi @William Rose, thank you very much for your detailed solution, I think it's exactly what I need, but I don't find the script you mentioned at the beginning, I'm new to this community and maybe I don't know where to look for it. However I have one doubt, when you talk about the extension of the velocity field by one sample you are referring to a cubic volume (you talked about the 6 faces, 8 corners and 12 edges), right? Is it possible to extend this procedure to a non cubic volume? Because if it is, I think this is perfect for my purpose.

[William Rose]

@Jacopo Rossi,

The code is attached to this message as a file which you should be able to click on and download. I forgot to attach it to my earlier comment.

The example I use in my code is not cubic, but it is rectangular, with dimensions 23 by 21 by 31 before extension, and 25x23x35 after extension.

>> convolutionPastEdge
size(x)=1x21; size(y)=1x23; size(z)=1x33
size(vx)=23x21x33
size(vxExtended)=25x23x35
>> 

[Jacopo Rossi]

Thanks a lot @William Rose, I tested it and it's perfect. Only a question, does the method you used to extend the faces and to fill corners have a particular name?

[William Rose]

@Jacopo Rossi,

I'd call it the Rose method! The face extension method is a generalizaiton of what has been done for 1D signals, in various papers. The methods for the edges and corners are original. You're welcome. Good luck with your work.