overlap of a scatter plot on a 2D imagesc matrix

Hi,

first of all I apologize if this might be an 'easy' question to answer or something has been discussed already.

I have created a 2D matrix and imaged it with imagesc function, as each pixel of this matrix represents a certain value.

I have to plot or scatter the locus of points corresponding to the middle line encompassed in the nonzero area of this matrix and I have some issues in visualizing it correctly.

I have created x_half_vec and y_half_vec which contain the x coordinate in pixels of the 2D map at half of the non-zero range along the x axis and for each row of the matrix, and its correspective value , respectively.

How can I overlap this line or locus of points on top of this image? I am nor sure x_half_vec and y_half_vec are even necessary

Thank you to anybody will answer me!

Here after the code:

clear all, close all, clc

n0= 1;                                          % Refractive index of material 1
n1= 1.51;                                       % Refractive index of material 2 (n1>n0)

alpha_in=[-30:0.5:65];                          % Source incident angle [deg]
x=size(alpha_in);
x=x(2);

beta=[55:0.1:80];                               % Wedge angle
y=size(beta);
y=y(2);

eta=zeros(y,x);
eta_half=zeros(y,x);
x_half_vec=zeros(1,x);
y_half_vec=zeros(1,y);

% return
n=n1/n0;
alpha_refr=asind((n0/n1).*sind(alpha_in));      % Refracted angle [deg]

r_s=(cosd(alpha_in)-n.*cosd(alpha_refr))./(cosd(alpha_in)+n.*cosd(alpha_refr)); 
r_p=(n.*cosd(alpha_in)-cosd(alpha_refr))./(n.*cosd(alpha_in)+cosd(alpha_refr));         

t_s=1+r_s;         
t_p=1/n.*(1+r_p);

R_s=abs(r_s).^2
R_p=abs(r_p).^2
T_s=(n1.*cosd(alpha_refr))./(n0.*cosd(alpha_in)).*  abs(t_s).^2
T_p=(n1.*cosd(alpha_refr))./(n0.*cosd(alpha_in)).*  abs(t_p).^2

R_eff=0.5*(R_s+R_p);
% return

for i=1:1:y
    for j =1:1:x
        if (alpha_in(j) <= asind(n1*sind(beta(i)-42))) && (alpha_in(j)>= asind(n1*sind(2*beta(i)-138)))
            eta(i,j)=1-R_eff(j);  
            %            disp('True')
            
            [row,col]= find(eta(i,:));               % Non-zero pixel values index
            x_half=round((col(1)+col(end))/2);       % x index in the middle of the range in pixels 
            y_half=eta(i,x_half);                    % correspective y value
            x_half_vec(1,i)=x_half;                  % vector of middle points on x axis in pixels
            y_half_vec(1,j)=y_half;
            
        else
            eta(i,j)=0;
            %             disp('Not satisfied')
        end
    end
end

figure; imagesc(alpha_in,beta,eta); colorbar; 
set(gca,'YDir','normal')
set(gca,'FontSize',16)
caxis ([0.85 1])
xlabel('Incident angle \alpha [deg])','FontSize',18)
ylabel('Wedge angle \beta [deg]','FontSize',18)
title ('Coupling efficiency [%]')

Answer 1

[Mathieu NOE]

hello Giulia

maybe this ?

clear all, close all, clc
n0= 1;                                          % Refractive index of material 1
n1= 1.51;                                       % Refractive index of material 2 (n1>n0)
alpha_in=[-30:0.5:65];                          % Source incident angle [deg]
x=size(alpha_in);
x=x(2);
beta=[55:0.1:80];                               % Wedge angle
y=size(beta);
y=y(2);
eta=zeros(y,x);
eta_half=zeros(y,x);
x_half_vec=zeros(1,x);
y_half_vec=zeros(1,y);
% return
n=n1/n0;
alpha_refr=asind((n0/n1).*sind(alpha_in));      % Refracted angle [deg]
r_s=(cosd(alpha_in)-n.*cosd(alpha_refr))./(cosd(alpha_in)+n.*cosd(alpha_refr)); 
r_p=(n.*cosd(alpha_in)-cosd(alpha_refr))./(n.*cosd(alpha_in)+cosd(alpha_refr));         
t_s=1+r_s;         
t_p=1/n.*(1+r_p);
R_s=abs(r_s).^2;
R_p=abs(r_p).^2;
T_s=(n1.*cosd(alpha_refr))./(n0.*cosd(alpha_in)).*  abs(t_s).^2;
T_p=(n1.*cosd(alpha_refr))./(n0.*cosd(alpha_in)).*  abs(t_p).^2;
R_eff=0.5*(R_s+R_p);
% return
for i=1:1:y
    for j =1:1:x
        if (alpha_in(j) <= asind(n1*sind(beta(i)-42))) && (alpha_in(j)>= asind(n1*sind(2*beta(i)-138)))
           eta(i,j)=1-R_eff(j);  
%            disp('True')
           [row,col]= find(eta(i,:)>eps);           % Non-zero pixel values index
           x_half=round((col(1)+col(end))/2);       % x index in the middle of the range in pixels 
           y_half=eta(i,x_half);                    % correspective y value
           x_half_vec(1,i)=x_half;                  % vector of middle points on x axis in pixels
           y_half_vec(1,j)=y_half;
          
        else
            eta(i,j)=0;
%             disp('Not satisfied')
        end
    end
end
figure; hold on
imagesc(alpha_in,beta,eta); colorbar; 
set(gca,'YDir','normal')
set(gca,'FontSize',16)
caxis ([0.85 1])
xlabel('Incident angle \alpha [deg])','FontSize',18)
ylabel('Wedge angle \beta [deg]','FontSize',18)
title ('Coupling efficiency [%]')
plot(alpha_in(x_half_vec),beta,'r','linewidth',8);
hold off

[Giulia Panusa]

Yes!!!! That's it!

Thank you very much, I was almost there... For some reason I thought I could not assign indeces belonging to vectors of different sizes to one another!

Thank you so much :)

[Image Analyst]

@Giulia Panusa, could you then click the "Accept this answer" link to award @Mathieu NOE his "reputation points". Thanks in advance. 🙂

[Giulia Panusa]

Done ;)

[Mathieu NOE]

Thanks to both of you !

@Image Analyst thanks to take care of my reputation :)

[Giulia Panusa]

Actually, I have another question here..

I need to extract the pixel indices for beta (y axis) which correspond a certain interval range of alpha (x axis).

Now I have a vector (1 x 251), both expressed in pixel position and actual value in degrees for alpha_in.

I need to retrieve values of beta for which I am within a certain range of alpha.

For example: alpha_in >=-2.5 and alpha_in <=2.5 (in pixels they correspond to 56 and 66); what are the corrispective for beta?

[Mathieu NOE]

hello Giulia

this is my suggestion .... I created the lower (black) and upper (magenta) curves - you may remove that from the plot if you wish, it's mostly to show what I'm doing

then if you "zoom" down in your specified alpha range, I could define the new beta range as : beta min is the min values of the "black" segment , and beta max is the max value of the "magenta" segment

this is the demo for a alpha range with an offset of 25 degrees (i.e. alpha_in >=-2.5 + offset and alpha_in <=2.5 +offset) so that the black segment is still visible (not empty) so it proves the method works there too

if you come back to alpha_in >=-2.5 and alpha_in <=2.5 (zero offset) as requested, the black segment is not visible (empty) but it's fairly straitghforward that in that case beta_low = min(beta,[],'all');

Demo for offset = 25 (just for fun and quick test of the code in other alpha range)

beta_min = 56.7000

beta_max = 77.9000

%%%%%%%%%%%%%%%%%%%%%%%%%%%

Demo for offset = 0 (your request)

beta_min = 55

beta_max = 69.8000

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

full code :

clear all, close all, clc
n0= 1;                                          % Refractive index of material 1
n1= 1.51;                                       % Refractive index of material 2 (n1>n0)
alpha_in=[-30:0.5:65];                          % Source incident angle [deg]
x=size(alpha_in);
x=x(2);
beta=[55:0.1:80];                               % Wedge angle
y=size(beta);
y=y(2);
eta=zeros(y,x);
eta_half=zeros(y,x);
x_half_vec=zeros(1,x);
y_half_vec=zeros(1,y);
% return
n=n1/n0;
alpha_refr=asind((n0/n1).*sind(alpha_in));      % Refracted angle [deg]
r_s=(cosd(alpha_in)-n.*cosd(alpha_refr))./(cosd(alpha_in)+n.*cosd(alpha_refr)); 
r_p=(n.*cosd(alpha_in)-cosd(alpha_refr))./(n.*cosd(alpha_in)+cosd(alpha_refr));         
t_s=1+r_s;         
t_p=1/n.*(1+r_p);
R_s=abs(r_s).^2;
R_p=abs(r_p).^2;
T_s=(n1.*cosd(alpha_refr))./(n0.*cosd(alpha_in)).*  abs(t_s).^2;
T_p=(n1.*cosd(alpha_refr))./(n0.*cosd(alpha_in)).*  abs(t_p).^2;
R_eff=0.5*(R_s+R_p);
% return
for i=1:1:y
    for j =1:1:x
        if (alpha_in(j) <= asind(n1*sind(beta(i)-42))) && (alpha_in(j)>= asind(n1*sind(2*beta(i)-138)))
           eta(i,j)=1-R_eff(j);  
%            disp('True')
           [row,col]= find(eta(i,:)>eps);           % Non-zero pixel values index
           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           % lower and upper limit curves  processing
           x1(1,i) = col(1);
           y1(1,i) = row(1);
           x2(1,i) = col(end);
           y2(1,i) = row(end);
           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           x_half=round((col(1)+col(end))/2);       % x index in the middle of the range in pixels 
           y_half=eta(i,x_half);                    % correspective y value
           x_half_vec(1,i)=x_half;                  % vector of middle points on x axis in pixels
           y_half_vec(1,j)=y_half;
          
        else
            eta(i,j)=0;
%             disp('Not satisfied')
        end
    end
end
figure; hold on
imagesc(alpha_in,beta,eta); colorbar; 
set(gca,'YDir','normal')
set(gca,'FontSize',16)
caxis ([0.85 1])
xlabel('Incident angle \alpha [deg])','FontSize',18)
ylabel('Wedge angle \beta [deg]','FontSize',18)
title ('Coupling efficiency [%]')
plot(alpha_in(x_half_vec),beta,'r','linewidth',12);
alpha_in1 = alpha_in(x1);
alpha_in2 = alpha_in(x2);
plot(alpha_in1,beta,'m','linewidth',12);
plot(alpha_in2,beta,'k','linewidth',12);
hold off

%% zoom within range alpha : alpha_in >=-2.5 & alpha_in <=2.5
offset = 0; % 0 = your request , 25 just for quick test the code in other alpha ranges
ind_al = find(alpha_in >=-2.5 +offset & alpha_in <=2.5  +offset);
alpha2 = alpha_in(ind_al);
eta2 = eta(:,ind_al);

% lower and upper limit curves  processing
ind_al12 = find(alpha_in1 >=-2.5 +offset & alpha_in1 <=2.5 +offset);
ind_al22 = find(alpha_in2 >=-2.5 +offset & alpha_in2 <=2.5 +offset);

alpha_in12 = alpha_in1(ind_al12);
beta12 =  beta(ind_al12);
alpha_in22 = alpha_in2(ind_al22);
beta22 =  beta(ind_al22);

figure; hold on
imagesc(alpha2,beta,eta2); colorbar; 
set(gca,'YDir','normal')
set(gca,'FontSize',16)
caxis ([0.85 1])
xlabel('Incident angle \alpha [deg])','FontSize',18)
ylabel('Wedge angle \beta [deg]','FontSize',18)
title ('Coupling efficiency [%]')
plot(alpha_in12,beta12,'m','linewidth',12);
plot(alpha_in22,beta22,'k','linewidth',12);
hold off

%% beta range
if isempty(beta22)
    beta_min = min(beta,[],'all');
else
    beta_min = min(beta22);
end

beta_max = max(beta12);

[Giulia Panusa]

Dear Mathieu,

sorry for this late response, and thank you so much for your help. Your feedback helped me to get what I wanted. I have modified a bit the code you suggested and did two functions in order to retrieve the beta intervals. Just out of curiosity I post them here:

Function1:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% n0= 1; % Refractive index of material 1

% n1= 1.51; % Refractive index of material 2 (n1>n0)

% delta_x=0.5;

% alpha_in=[-50:delta_x:65]; % Source incident angle [deg]

% x=size(alpha_in);

% x=x(2);

%

% beta_angle=[55:0.1:80]; % Wedge angle

% y=size(beta_angle);

% y=y(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [coupling_efficiency,x_half_vec,y_half_vec,half_angle] = test_beta_2(n0,n1,n,x,y,alpha_in,beta_angle,eta) % inputs to run this function

eta=zeros(y,x);

eta_half=zeros(y,x);

x_half_vec=zeros(1,x);

y_half_vec=zeros(1,y);

n=n1/n0;

alpha_refr=asind((n0/n1).*sind(alpha_in)); % Refracted angle [deg]

r_s=(cosd(alpha_in)-n.*cosd(alpha_refr))./(cosd(alpha_in)+n.*cosd(alpha_refr));

r_p=(n.*cosd(alpha_in)-cosd(alpha_refr))./(n.*cosd(alpha_in)+cosd(alpha_refr));

R_s=abs(r_s).^2

R_p=abs(r_p).^2

R_eff=0.5*(R_s+R_p);

% eta=zeros(y,x);

eta_half=zeros(y,x);

x_half_vec=zeros(1,x);

y_half_vec=zeros(1,y);

for i=1:1:y

for j =1:1:x

if (alpha_in(j) <= asind(n1*sind(beta_angle(i)-42))) && (alpha_in(j)>= asind(n1*sind(2*beta_angle(i)-138)))

eta(i,j)=1-R_eff(j);

% [row,col]= find(eta(i,:)); % Non-zero pixel values index

[row,col]=find(eta(i,:)>eps);

x_half=round((col(1)+col(end))/2); % x index in the middle of the range in pixels

y_half=eta(i,x_half); % correspective y value

x_half_vec(1,i)=x_half; % vector of middle points on x axis in pixels

y_half_vec(1,j)=y_half;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% lower and upper limit curves processing

x1(1,i) = col(1);

y1(1,i) = row(1);

x2(1,i) = col(end);

y2(1,i) = row(end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

else

eta(i,j)=0;

end

end

end

%% Display coupling efficiency map

figure; hold on

imagesc(alpha_in,beta_angle,eta); colorbar;

set(gca,'YDir','normal')

set(gca,'FontName','Palatino Linotype','FontSize',14,'LineWidth',1); colorbar

caxis ([0.85 1])

xlabel('Incident angle \alpha [deg])','FontSize',18)

ylabel('Wedge angle \beta [deg]','FontSize',18)

title ('Coupling efficiency [%]')

half_angle=alpha_in(x_half_vec);

plot(half_angle,beta_angle,'r','linewidth',3); % Plot middle curve

ylim([55 80])

alpha_in1 = alpha_in(x1); % Find correspective value of alpha in the lower index of the range we are studying

alpha_in2 = alpha_in(x2); % Find correspective value of alpha in the upper index of the range we are studying

plot(alpha_in1,beta_angle,'m','linewidth',3); % Plot of upper limit

plot(alpha_in2,beta_angle,'k','linewidth',3); % Plot of lower limit

ylim([55 80])

hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

coupling_efficiency=eta;

end

Function 2:

% This script computes the intervals of beta

function [xx1_vector,xx2_vector,yy1_vector,yy2_vector] = intervals_definition(alpha_in,half_angle,beta_angle,coupling_efficiency)

step_field=5;

offset=[0:step_field:40];

s_o=size(offset);

s_o=s_o(2);

beta_min_vec=zeros(1,s_o);

beta_max_vec=zeros(1,s_o);

for i=1:1:s_o

ind_al = find(alpha_in >=-2.5 +offset(i) & alpha_in <=2.5 +offset(i)); % It finds the indices for alpha encompassed within a certain range

alpha2 = alpha_in(ind_al); % Actual value of alpha in that range

eta2 = coupling_efficiency(:,ind_al); % Coupling efficiency for alpha within that range (submatrix of eta)

ind_al_beta_half = find(half_angle >=-2.5 +offset(i) & half_angle <=2.5 +offset(i)); % Finds indices of beta in correspondence of half angle within a certain range

alpha_beta=half_angle(ind_al_beta_half); % Values of alpha in that range (needs to be consistent with alpha2)

beta_half=beta_angle(ind_al_beta_half); % Values of beta in that range

beta_min=min(beta_half); % Minimum value of beta

beta_max=max(beta_half); % Maximum value of beta

beta1=find(beta_angle==beta_min); % beta range min in pixels

beta2=find(beta_angle==beta_max); % beta range max in pixels

beta_range=beta_half;

xx1=find(alpha_in==offset(i)-2.5)

xx2=find(alpha_in==offset(i)+2.5)

beta_min=round(beta_min);

yy1=find (beta_angle==beta_min)

beta_max=round(beta_max);

yy2=find (beta_angle==beta_max)

xx1_vector(1,i)=xx1;

xx2_vector(1,i)=xx2;

yy1_vector(1,i)=yy1;

yy2_vector(1,i)=yy2;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Thank you again for your support!

Best,

Giulia.

[Mathieu NOE]

hello Giulia

thanks for the update - glad you have a solution that works

Good luck for the future !

[Giulia Panusa]

You too! :)