A non centered method to compute velocity derivative

1 view (last 30 days)
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?

Accepted Answer

William Rose
William Rose on 19 Jan 2022
Edited: William Rose on 20 Jan 2022
[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.
  5 Comments
William Rose
William Rose on 26 Jan 2022
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.

Sign in to comment.

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by