How do you shift the data such that the centerline of the jet is at 0?

2 次查看（过去 30 天）
Monique Embury
Monique Embury 2020-1-21
编辑： Adam Danz 2020-2-25
Hi,
I was running a PIV system and didnt realize the camera moved. As I have been working on post processing, I need to shift the plot so it is symmetric about 1040 pixels, 0, and 21 mm. The axis must remian the same and the plot must extend from -1 to 1 still. I pretty much want the axis lined up properly with keeping the same plot.
Any suggections? I have posted my code and a sample of 3 imgaes that need to be shift and 3 images where the camera is where it should be.
shift plots matlab.PNG
correct shift plots matlab.PNG
function [x,y,u_avg,v_avg,CHC_tot] = piv_averages(prefix,suffix,Nstart,Nfinish,interp)
% This program reads in a series of instantaneous PIV vector fields from
% Insight and averages them. The user has the option of excluding
% interpolated vectors, which have CHC > 1. (interp = 0 means do not
% interpolate, while interp = 1 means interpolate).
% Create file name for each image
c_exit=2.776; %speed of sound
x0=1040; %origin of the jet in pixels
y0=53.8; %origin of the jetin pixels
D=923.71; %diameter of jet exit in pixels
v_shift=0;
for i = Nstart:Nfinish
Nstring = int2str(i); % Convert iteration number to a character string
if i < 10
filename_inst = strcat(prefix,'0000',Nstring,suffix);
elseif i < 100
filename_inst = strcat(prefix,'000',Nstring,suffix);
elseif i < 1000
filename_inst = strcat(prefix,'00',Nstring,suffix);
elseif i < 10000
filename_inst = strcat(prefix,'0',Nstring,suffix);
else
filename_inst = strcat(prefix,Nstring,suffix);
end
% Read file name
A_inst = csvread(filename_inst,1,0);
x = A_inst(:,1); % x-position (mm)
y = A_inst(:,2); % y-position (mm)
u = A_inst(:,3); % x-velocity (m/s)
v = A_inst(:,4); % y-velocity (m/s)
chc = A_inst(:,5); % number of good vectors at this location
N = size(x,1); % Length of entire vector array
% Initialize output variables if this is the first file
if i == Nstart
u_tot = zeros(N,1);
v_tot = zeros(N,1);
CHC_tot = zeros(N,1);
end
for j = 1:N
if interp == 0
if chc(j,1) == 1
u_tot(j,1) = u_tot(j,1) + u(j,1);
v_tot(j,1) = v_tot(j,1) + v(j,1);
CHC_tot(j,1) = CHC_tot(j,1) + 1;
end
elseif interp == 1
if chc(j,1) > 0
u_tot(j,1) = u_tot(j,1) + u(j,1);
v_tot(j,1) = v_tot(j,1) + v(j,1);
CHC_tot(j,1) = CHC_tot(j,1) + 1;
end
end
end
end
for j = 1:N
u_avg(j,1) = u_tot(j,1)/CHC_tot(j,1);
v_avg(j,1) = v_tot(j,1)/CHC_tot(j,1);
end
% Set origin to jet exit centerline
x_c = x - (x0);
y_c = y - y0;
% Shift by convective velocity
v = v - v_shift;
% Nondimensionalize variables
x_non = x_c/D; % Nondimensionalize using jet diameter
y_non = y_c/D; % Nondimensionalize using jet diameter
u_non = u_avg/c_exit; % Nondimensionalize using sonic speed
v_non = v_avg/c_exit; % Nondimensionalize using sonic speed
%%%%%%%FOR H/D=2%%%%%%%%%
% 1) Choose a threshold
yThreshold = 400; %accept all vectors that start at or above y = yThreshold;
% 2) identify all quiver arrows that meet the criteria
acceptIdx = y >= yThreshold;
% 3) plot the quiver arrows, but only the ones accepted
% figure(4)
% quiver(x(acceptIdx), y(acceptIdx), u_avg(acceptIdx), v_avg(acceptIdx), 5)
% % add reference line at threshold if you'd like
% line(x,yThreshold);
% Plot nondimensional vector field
figure(2)
quiver(x_non(acceptIdx),y_non(acceptIdx),u_non(acceptIdx),v_non(acceptIdx),4)
axis([-1 1 0 2])
xticks([-1 0 1])
yticks([0 1 2])
set(gca,'YtickLabel',2:-1:0)
set(gca,'XAxisLocation','top','YAxisLocation','left');
xlabel('z/D')
ylabel('x/D')
title('Nondimensional velocity field')
% camroll(90)
% Plot averaged vector field
figure(1)
quiver(x(acceptIdx),y(acceptIdx),u_avg(acceptIdx),v_avg(acceptIdx),4)
axis([0 2000 0 2000]);
set(gca,'YtickLabel',2000:-200:0);
set(gca,'XAxisLocation','top','YAxisLocation','left');
xlabel('z (pixels)')
ylabel('x (pixels)')
title('Dimensional velocity field')
% camroll(90)
% Plot averaged vector field
figure(3)
quiver(x(acceptIdx)/48.11,(y(acceptIdx)/48.11),u_avg(acceptIdx),v_avg(acceptIdx),4)
axis([0 42 0 42])
yticks([0 5 10 15 20 25 30 35 40])
set(gca,'XAxisLocation','top','YAxisLocation','left');
ax = gca;
ax.YTickLabel = flipud(ax.YTickLabel);
xlabel('z (mm)')
ylabel('x (mm)')
title('Dimensional velocity field')
% camroll(90)
% Output averaged data
output_avg = [x y u_avg v_avg CHC_tot];
dlmwrite('average_vels.dat','xyUVC');
dlmwrite('average_vels.dat',output_avg,'-append');
Thank you in advance!
  6 个评论
显示 4更早的评论
Monique Embury
Monique Embury 2020-1-21
Well that is not what I wanted to hear. Thank you though. Is there any way to just alter the appearance of the graph so the axis looks correct?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Adam Danz
Adam Danz 2020-1-21
编辑：Adam Danz 2020-1-22
The goal isn't quite clear to me. If you'd just like to scale the quiver plots so that the longest axis fits within -1 to 1, follow this demo.
This first part creates a demo set of inputs that you presumably already have.
% Create quiver plot inputs
[x,y] = meshgrid(0:0.2:2,0:0.2:3);
u = cos(x).*y;
v = sin(x).*y;
% Plot the original data. X values span
% from 0 to 2, y values span from 0 to 3.
clf()
subplot(2,1,1)
quiver(x,y,u,v)
axis equal
This part is what would be applied to your data.
% Determine the range of your data
xRange = range(x(:));
yRange = range(y(:));
scale = max(xRange, yRange);
% Scale data so it's longest axis is within [-1,1]
xS = (x - min(x(:)))/scale*2 - xRange/scale;
yS = (y - min(y(:)))/scale*2 - yRange/scale;
% plot the scaled quiver values
subplot(2,1,2)
quiver(xS,yS,u,v)
axis equal
Or are you trying to find the center point between the two focii?
  18 个评论
显示 16更早的评论
Adam Danz
Adam Danz 2020-1-27
Here's the steps you can take to get you started. See inline comments for details.
% Load the first image of the group of 100 images.
load('average_vels.mat')
% Extract the quiver values
x = output_avg(:,1);
y = output_avg(:,2);
u = output_avg(:,3);
v = output_avg(:,4);
%===============================================
% Produce the quiver plot
clf()
plot(x,y,'k.','MarkerSize',3)
hold on
qh = quiver(x,y,u,v,4,'k');
axis equal
% --Or-- perhaps you'd rather work with the contour plot
% Choose the plot above or this plot below.
% Compute vector length of each quiver
vecLen = sqrt(u.^2 + v.^2);
sortedGrid = sortrows([x,y,vecLen]);
xUnq = unique(x(:));
yUnq = unique(y(:));
vecLenMat = reshape(sortedGrid(:,3),numel(xUnq),numel(yUnq));
clf()
contourf(xUnq,yUnq,vecLenMat,15);
cb = colorbar();
ylabel(cb,'vector magnitude')
axis equal
hold on
plot(x,y,'k.','MarkerSize',3)
% ==============================================
% Here is where you'll manually select the centers
% of each focii with a mouse click.
msgbox('Select the left and then the right focus.') % optional msg to user
LR = ginput(2); % [x,y] coordinates
% You may want to select the closest (x,y) coordinates
% within your grid to the selected points. I'm not sure
% if it's expected or not that the focii are centered
% on a grid point. Let me know if you need help with that.
% ===================================================
% Plot the two selected coordinates
centers = plot(LR(:,1),LR(:,2),'rx');
% Confirm satisfaction with user.
% If they are unsatisfied the code will quit. You can
% change that to make it more useful.
yn = questdlg('Satisfied?',mfilename,'yes','no','yes');
if strcmpi(yn,'no')
delete(centers)
return
end
% Compute the required shift needed to center the center point
% between the two focii (positive means rightward shift is needed)
target = 1040; % where center point should be.
shift = (LR(1,1) + (LR(2,1)-LR(1,1))/2) - target;
% "shift" is the lateral shift needed for all 100 images in this
% groups. Now you can merely loop through those 100 images
% and apply a form of the code below to each image.
%=========================================================
% Shift it
x = x + shift;
% Plot the results if you want
figure()
qh = quiver(x,y,u,v,4,'k');
axis equal
axis tight
xlim([0,inf])
xline(target,'r-') %Show the center point.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Axes Appearance 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by