Azimuthally average line plot to 2D image

Hi there,

I am trying to create a 2D image from a line profile.

Effectively, I would like to obtain its radial average and it must be centred in the middle of the image.

I have attached the x and y values of the plot.

I don't know where to begin.

Your help in anyway is appreciated,

Arthur

[Mathieu NOE]

hello

maybe I'm wrong but is the question about doin a kind of circle fit to your data (finding radial average ? )

[Arthur Moya]

@Mathieu NOE, indeed.

Answer 1

[Mathieu NOE]

Hello again

Check this if it works for you.

I had to divide y by 1000 so that x and y are in similar ranges (is one in mm and the other one in microns ?) so circfit would work.

Code below + function.

All the best.

 

a = load('Chip_NPS_fit_x_y.txt');
xs = a(:,1);
ys = a(:,2)/1000;

ind = find(xs>0.15 & xs < 0.2);
[xfit,yfit,Rfit,a] = circfit(xs(ind),ys(ind));

figure(1)
plot(xs,ys,'b',xs(ind),ys(ind),'g')
hold on
rectangle('position',[xfit-Rfit,yfit-Rfit,Rfit*2,Rfit*2],...
    'curvature',[1,1],'linestyle','-','edgecolor','r');
title(sprintf('Best fit: R = %0.3f; Ctr = (%0.3f,%0.3f)',...
    Rfit,xfit,yfit));
plot(xfit,yfit,'g.')
axis equal
hold off

%%%%%%%%%%%%%%%%%%%%%
function   [xc,yc,R,a] = circfit(x,y)
%
%   [xc yx R] = circfit(x,y)
%
%   fits a circle  in x,y plane in a more accurate
%   (less prone to ill condition )
%  procedure than circfit2 but using more memory
%  x,y are column vector where (x(i),y(i)) is a measured point
%
%  result is center point (yc,xc) and radius R
%  an optional output is the vector of coeficient a
% describing the circle's equation
%
%   x^2+y^2+a(1)*x+a(2)*y+a(3)=0
%
%  By:  Izhak bucher 25/oct /1991, 
    x=x(:); y=y(:);
   a=[x y ones(size(x))]\[-(x.^2+y.^2)];
   xc = -.5*a(1);
   yc = -.5*a(2);
   R  =  sqrt((a(1)^2+a(2)^2)/4-a(3));
end

[Arthur Moya]

Thank you so much @Mathieu NOE but I think I may have not explained it as well as I could have. What I am actually looking for is something that looks like the attached image below. Actually, the y units are electrons/pixel while the x units are 1/pixel up to 1/2pixel.

[Mathieu NOE]

ok I got it !

so the idea is to make a surface plot by "rotating" your profile along the y axis

[Mathieu NOE]

here we go :

a = load('Chip_NPS_fit_x_y.txt');
r = a(:,1);
zs = movmean(a(:,2),25); % a bit of smoothing (optional)


theta=linspace(0,2*pi,length(r));

[x,y]=pol2cart(theta, r);
zsr=zs*ones(1,length(r));
figure(1);
s = surf(x,y,zsr);
s.EdgeColor = 'none';
view([0 90]);
axis('equal');
colormap('gray');
colorbar('vert');
axis([-0.3 0.3 -0.3 0.3])

[Arthur Moya]

@Mathieu NOE this is great. I just have 2 more questions in relation to your answer. The result is a surface plot. Would it be possible to convert it into a nxn matrix double. Secondly, the matrix must be a square shape of 2048x2048 pixels. I realised that the surface plot is circular with a maximum radius r. Any advice on how to achieve these two things? Thank you.

Arthur

[Mathieu NOE]

hello

this is it

the output is a matrix of size 2048 x 2048 stored in zq

a = load('Chip_NPS_fit_x_y.txt');
r = a(:,1);
zs = movmean(a(:,2),25); % a bit of smoothing (optional)

figure(1);
plot(r,zs);

 % create the mexican hat 3D plot
 theta=linspace(0,2*pi,length(r));
[x,y]=pol2cart(theta, r);

zsr=zs*ones(1,length(r));
figure(2);
s = surf(x,y,zsr);
s.EdgeColor = 'none';
view([0 90]);
axis('equal');
colormap('gray');
colorbar('vert');

% interpolation on a 2048 x 2048 square grid
% square limit = max(r) / sqrt(2);
ll = max(r) / sqrt(2);
m = linspace(-ll,ll,2048);
[Xq,Yq] = meshgrid(m,m);
zq = griddata(x,y,zsr,Xq,Yq); %the output is a matrix of size 2048 x 2048 stored in zq
figure(3);
s = surf(Xq,Yq,zq);
s.EdgeColor = 'none';
view([0 90]);
axis('equal');
colormap('gray');
colorbar('vert');
xlim([-ll ll]);
ylim([-ll ll]);

[Arthur Moya]

@Mathieu NOE thank you. This works perfectly.

Kind regards,

Arthur

Answer 2

[Image Analyst]

Here is a simple brute force for-loop way. Simply create an image then scan every pixel location determining the radius from the center for that pixel and assigning the signal value there.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format short g;
format compact;
fontSize = 18;
fprintf('Beginning to run %s.m ...\n', mfilename);

data = readmatrix('Chip_NPS_fit_x_y.txt');
r = data(:, 1);
signal = data(:, 2);
nexttile
plot(r, signal, 'b-');
title('Radial Profile', 'FontSize', fontSize);
grid on;
drawnow;

numElements = length(signal)
rows = ceil(2 * numElements / sqrt(2))
middlex = rows / 2;
middley = rows / 2;
grayImage = zeros(rows, rows);
for col = 1 : rows
	fprintf('Assigning column %d.\n', col);
	for row = 1 : rows
		% Get the distance from the center of the image.
		radiusInPixels = sqrt((col - middlex)^2 + (row - middley)^2);
		% Scale it to be in the range 0 - 0.5, like r is.
		% Find the signal value there.
		radius = r(end) * radiusInPixels / numElements;
		[~, index] = min(abs(r - radius));
		grayImage(row, col) = signal(index);
	end
end
nexttile
imshow(grayImage, []);
title('As a 2-D image', 'FontSize', fontSize);
% Maximize window
g = gcf;
g.WindowState = 'maximized';

[Arthur Moya]

@Image Analyst thank you.

This works a charm as well. However, I had told @Mathieu NOE that I needed the final image to be 2048x2048. In your code the final image was 1026x1026.

Kind regards,

Arthur