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

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.

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!

[darova]

The question is unclear for me. Where is the axis or point you want to shift?

[Monique Embury]

I would like to shift all the data points. In this case shift all data points by 200 pixels to the right. If you only shift the x values you get a plot similar to this:

instead of this the data pots extending to the x axis like this:

Does this make more sense?

[darova]

You want to extrapolate data? Want to fill this part?

[Monique Embury]

yes fill that part.

[darova]

I think it's impossible

[Monique Embury]

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?

Answer 1

[Adam Danz]

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?

[Monique Embury]

I am trying to get the center point of the two focii to be localted at 0 and the plot to extend fully to 1 and -1.

[Adam Danz]

Could you attach the figure or the x,y,u,v variables in a mat file?

[Adam Danz]

In order to play around with the data to determine the best way to find the center between the two focii, I just need the quiver data for a single figure. How are you creating the figures you shared via screenshot? If you attach one of those figure files I could extract the values from that. Or if you're creating the plot via a single call to quiver(), you could save the inputs like this

save('quiverInputs.mat','x','y','u','v') %or whatever variables names they are
quiver(x,y,u,v)

There may be other contextual information that would be useful such as (I'm making this up)

• the two focii should always have the same y coordinate
• the two focii will always be x units apart
• the coordinate of the left focii should always be (5,4)
• etc...

[Monique Embury]

The centerline of the jet SHOULD always be at 1040 pixels. Although this is untrue for this image because the camera shifted wihtout my knowledge during testing. I was able to get MATLAB to spit out the average mat file. I cannot attach it.

"Cannot attach this file because:

• You are limited to 10 daily uploads. If you need to upload additional files, delete one or more files now or wait 24 hours to upload more files."

[Adam Danz]

hmmmm, never saw that one before.

I see you deleted the comment with a bunch of attachements. I also deleted an additional comment of yours in this thread with an attachment. See if that restriction updates in a few minuts.

Alternatively you could upload it somewhere else and provide a download link. But I'll give you a heads up that I probably won't be able to dive into this until later today so no rush.

[Monique Embury]

Okay, I will try to upload them later. Thank you for your time nad help thus far.

[Monique Embury]

https://drive.google.com/file/d/1Lk1zYrgMKFVlWNVuELOkrU85PJ5HN7xF/view?usp=sharing

I shared a Google Drive link since I cannot upload anything still.

[Adam Danz]

Hi Monique, I'm happy to help.

The file is not a mat file. Something must have gone wrong with saving the variables. The file is not a binary mat file. It can be opened in a text editor and merely lists the character, "x,y,U,V,C".

[Monique Embury]

Okay, I tried something else and was able to open it in MATLAB.

https://drive.google.com/open?id=1Lk1zYrgMKFVlWNVuELOkrU85PJ5HN7xF

[Adam Danz]

That worked.

1. Are the two focii at a fixed distance from eachother?
2. Are either of the two focii expected to be at a fixed distance from the left or right edge of the plot?
3. Is the center of the data in the center between the two focii?

I like these kinds of problems because they challenge me but I'm a little stuck at locating these two local minima. The answers to those quesitons above may be helpful.

[Monique Embury]

The center of the data should be located between the two focii. The other two questions only apply in symmetric cases which is rare so I would not use those as a criteria. Also could you show me how to came to that contour plot?

[Adam Danz]

Monique, how many of these images need fixed? I thought of a workaround that would require you to click on the two focii manually and then the code does the rest. But that wouldn't be useful if you have very many images to fix.

The code to reproduce the contour plot is below. The contour levels correspond to the vector magnitudes.

load('average_vels.mat')
x = output_avg(:,1); 
y = output_avg(:,2);
u = output_avg(:,3);
v = output_avg(:,4);

% 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)); 
[cm, ch] = contourf(xUnq,yUnq,vecLenMat,15);
cb = colorbar(); 
ylabel(cb,'vector magnitude')
axis equal

The contour line coordinates can be extracted using this FEX function.

contourTable = getContourLineCoordinates(cm);

The next step would be to isolate the coordinates of the two focii which won't be trivial.

[Monique Embury]

Adam, I have 14 phases to fix this this shift, and then I belive I have another 14, each of those having 100 images. so a total of 2,800 imgaes to fix, but only 28 averages.

[Adam Danz]

If I understood correctly, there are 28 uniquely different camera-shifts and within each erroneous shift there are 100 images that need fixed. Within each of those 100 images, the exact same correction is needed. Also, you don't expect additional corrections needed in the future.

Is that all correct?

If so, I don't think it would be that much work to create an interactive code that loads the first of 100 images and you manually click on the two focii of the first image in that set, then the code corrects for the camera shift for all 100 images. Then that process would repeate for the next set of 100 images and so on.

Does that sound reasonable?

Ideally the code would automatically locate the coordinates of the two focii but I haven't figured out how to do that given these data.

[Monique Embury]

Yes, that sounds correct! That all sounds like what needs to happen.

[Monique Embury]

Thank you again Adam

[Adam Danz]

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.

[Adam Danz]

Glad I could help.

If you have any Qs, let me know.