Azimuthally average line plot to 2D image
5 次查看(过去 30 天)
显示 更早的评论
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
2 个评论
Mathieu NOE
2021-3-23
hello
maybe I'm wrong but is the question about doin a kind of circle fit to your data (finding radial average ? )
采纳的回答
Mathieu NOE
2021-3-23
编辑:Image Analyst
2021-3-23
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/559453/image.png)
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
6 个评论
Mathieu NOE
2021-3-24
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]);
更多回答(1 个)
Image Analyst
2021-3-24
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';
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/560323/image.png)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Lighting, Transparency, and Shading 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!