Narupon Chattrapiban

This code should be changed a bit.
You might want to change the code in line 71 from

Rbnd = floor(min(yxz)/2)-1;

to

Rbnd = min([xavg,yxz(2)-xavg,yavg,yxz(1)-yavg])-1;

The function can be used as the following.

If you have a 250x250 matrix, M, then the radial profile
originates at (xc,yc) = (125,125) can be obtained from

r = rscan(M,'xavg',xc,'yavg',yc);

The option 'squeezx' and 'squeezy' is used to scan M with
ellipses instead of circles. You can try to change squeezx to
have a different value from squeezy and see what the program
will do. The option 'rot' is used to align the ellipses, 'rlim' is to
change the range of your scan, and 'rstep' is to change the size
of rscan increment. All the variables are described in pixels.
Finally, 'dispflag' and 'dispflagc' can be set to be either 1 or 0.
They are used to either turn the display on or off. Try and see
what they will do.

Anyhow the function rscan here is probably too slow.
You might want to consider the code 'rscan2' below.

For a matrix M, the radial profile 'P' whose coordinates is 'R'
can be obtained from

[R,P] = rscan2(M,rh,N);
figure;plot(R,P);

where 'rh' is the polar coordinates corresponding to M e.g.
x = linspace(-1,1,size(M,2));
y = linspace(-1,1,size(M,1));
[xx,yy] = meshgrid(x,y);
[th,rh]=cart2pol(xx,yy);

and 'N' is the length of R and P you want to have.

try
>> help rscan2

for more examples.

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

function [rAvg,yAvg,rStd,yStd,r,y] = rscan2(Ih,rh,N)

% [rAvg,yAvg,rStd,yStd,r,y] = rscan2(Ih,rh,N)
% Ih -- 2D Matrix
% rh -- polar rho Matrix (e.g. from cart2pol)
% N -- number elements in the output radial profile
%
% Example:
%
% x = linspace(-5,5,500);
% [xx,yy] = meshgrid(x);
% [th,rh]=cart2pol(2*cos(pi/3)*xx-
2*sin(pi/3)*yy,sin(pi/3)*xx+cos(pi/3)*yy);
% figure(1);clf;imagesc(rh)
% Ih = exp(-(rh/2).^2);
% figure(2);imagesc(Ih);
% [rAvg,yAvg] = rscan2(Ih,rh,10000);
% figure(3);plot(rAvg,yAvg);
%
% [th,rh]=cart2pol(xx,yy);
% Ih = exp(-(rh/2).^2);
% figure(2);imagesc(Ih);
% [rAvg,yAvg] = rscan2(Ih,rh,10000);
% figure(3);hold on;plot(rAvg,yAvg,'r--');
%

if nargin < 3, N = numel(rh); end

if N > numel(rh), N = numel(rh); end

dPxl = floor(numel(rh)/N);
nmEL = dPxl*N; % -> prepare for reshape

[r,i] = sort(rh(:));
y = Ih(i);

% v-- strip away element that cannot be reshape
r(nmEL+1:end) = [];
y(nmEL+1:end) = [];

rAvg = mean(reshape(r,dPxl,N),1);
yAvg = mean(reshape(y,dPxl,N),1);
rStd = std(reshape(r,dPxl,N),0);
yStd = std(reshape(y,dPxl,N),0);

if nargout == 0
figurel(10);clf;hold on;box on;
plot(r,y,'g'); plot(rAvg,yAvg,'b--');
end

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

Even though not easy to understand at first.... it works very well for me...excellent

Useful ! Thanks for sharing ! Cheers & Good Luck