Help with getting data from image at equidistant positions from the centre - to include the "off circle corners"

Hello, I have an image and I want to get the median of data at each radius value. This is shown by the yellow circles - so each circle (take the median of all pixels on it).

However, in my attemy, I missing all the corners of data that are actually at a larger radius than the image size, but not full circles. Any suggestions how to include this data.

This was my attempt:

         % 1st draw circle of radius R on the flog2 plot
            
            theta = linspace(0,360);    %use n=121 for every 3 degrees (reduce from every degree for speed
            centre = [cx,cy]; %circle centre
            % I= Image
                   
            [X,Y]=ndgrid(1:size(I,1),1:size(I,2));
          
            X=X-cx;Y=Y-cy;%shift coordinate grid
            hold(ax1,"on");
            
            data=[];
            for i=1:cy-1
                radius=i; %circle radius     
                data(i,1)=i;
                y = radius*sind(theta)+centre(2); 
                x = radius*cosd(theta)+centre(1);                 
                
        
                %Get Data on radius
                L=sqrt(X.^2+Y.^2)==radius;
                 data(i,2) = median(I(L));
                 if mod(i,5)==0
                     plot(ax1,x,y,'y'); %draw circle                   
                 end    

                 % Now use interp so use all pixels on circle.   This is
                 % actually what I finally use now               
                 [xx,yy]=pol2cart(theta,radius);   % this is for the interp2 approach
                 x=xx+cx;  y=yy+cy;
                 vals=median(interp2(I,x,y));
                 data(i,3)=vals;
            end
            
            % Now plot data
            hold(ax2,"on");           
            plot(ax2,data(:,1),data(:,3),'y-','LineWidth',3);  

Answer 1

[Matt J]

[M,N]=size(yourImage);

  c=[M+1,N+1]/2; %image center
  
  rho=hypot((1:M)'-c(1), (1:N)-c(2));

G=discretize(rho,0:2:norm([M,N]/2+2) );

Medians=splitapply(@median, yourImage(:),G(:));

[Jason]

Hi Matt, could you give a bit of an explanation of what each line does please|

[Matt J]

The code creates an image G of group labels, which groups pixels together into radial rings. This can then be used in conjunction with splitapply. It may help to plot, to see the effect:

yourImage = double(im2gray(imread('cameraman.tif')));

[M,N,~]=size(yourImage);
[X,Y]=ndgrid((0:M-1)-M/2,(0:N-1)-N/2);
[~,rho] = cart2pol(X,Y);
G=discretize(rho,0:10:norm([M,N]/2+10) );

imshow(G,[])

[Jason]

Thanks Matt, this is very fast! So that Im understanding this, the image G is an image that shows the median value at that location?

I really need the median value at each pixel distance radially, so is this correct?

d=1;
G=discretize(rho,0:d:norm([M,N]/2+d) );

Also, how could I pick off some of the circles and superimpose on the image (just to illustrate a sub set of the circles)?

[Matt J]

So that Im understanding this, the image G is an image that shows the median value at that location?

No, the medians are calculated in this line,

Medians=splitapply(@median, yourImage(:),G(:));

I really need the median value at each pixel distance radially, so is this correct?

I think d=1 is too small if you want the pixel groups to form complete, contiguous rings. The distance between pixels is not uniformly =1 in all directions. At lines radiating at 45 degrees, the distance will be sqrt(2). I would say d=3 would be a safe pick.

[Jason]

OK, so G is the "binning image".

Umm, Im confused now over this method. I really need to identify a central pixel, and then radiate outwards from that and median the circle of pixels in pixel steps.

[Matt J]

Yes, G is the binning image. You have no choice but to bin the pixels into rings that are only approximate. When you have a chessboard, the squares cannot be grouped to form clean, perfect circles. Nor is there necessarily a pixel located at the center of the image. If the image is NxN where N is an even integer, the geometric center lies between the pixels.

[Matt J]

Also, how could I pick off some of the circles and superimpose on the image (just to illustrate a sub set of the circles)?

One way,

yourImage = double(im2gray(imread('cameraman.tif'))); %Example image

[M,N]=size(yourImage);

  c=[M+1,N+1]/2; %image center
  
  rho=hypot((1:M)'-c(1), (1:N)-c(2));

d=3;
G=discretize(rho,0:d:norm([M,N]/2+d) ); %Binning image

Radii=splitapply(@mean, rho(:),G(:)); %Radial distance of each pixel group
Medians=splitapply(@median, yourImage(:),G(:)); %Medians of each pixel group

medianImage=Medians(G); %Median values replicated into rings

region=(60<=rho & rho<=90);
yourImage(region) = medianImage(region); %Replace pixels by median values
imshow(yourImage,[])

plot(Radii,Medians); xlabel Distance; ylabel Median

[Jason]

Thankyou Matt this is awesome.

Coming back to your comment about fractional pixels, does this mean the binning (discretize function) only uses the pixels where it (the circle) passes thru the centre of the pixels.

So for example, the 1st circle corresponding to 1 pixel in x or y from the centre will only have 4 pixels in the calculation rather than the 9 nearest neighbours?

[Matt J]

does this mean the binning (discretize function) only uses the pixels where it (the circle) passes thru the centre of the pixels.

No, the code looks at true, concentric circles of radii,

r = 0:d:norm([M,N]/2+d)

and assigns a pixel to the i-th ring if its center lies in the annulus with inner radius r(i) and outer radius r(i+1).

[Jason]

OK, so If d=1, the 1st cocnentric ring will be 0 to the 1st pixel, then 1st to 2nd pixel?

[Matt J]

I don't know what "1st pixel" and "2nd pixel" means. To see the precise pixel membership, you should look at the binning image G. But, we know that it will depend not only on d but also on whether the dimensions of the image (M,N) are odd or even. For example, with odd dimensions, like M=N=11, the first bin will only be the central pixel,

G=getG(11,11,1)

G = 11×11

8 7 6 6 6 6 6 6 6 7 8 7 6 6 5 5 5 5 5 6 6 7 6 6 5 4 4 4 4 4 5 6 6 6 5 4 3 3 3 3 3 4 5 6 6 5 4 3 2 2 2 3 4 5 6 6 5 4 3 2 1 2 3 4 5 6 6 5 4 3 2 2 2 3 4 5 6 6 5 4 3 3 3 3 3 4 5 6 6 6 5 4 4 4 4 4 5 6 6 7 6 6 5 5 5 5 5 6 6 7

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

whereas for an even-sized image, the first bin will contain the 4 center-most pixels,

G=getG(12,12,1)

G = 12×12

8 8 7 7 6 6 6 6 7 7 8 8 8 7 6 6 5 5 5 5 6 6 7 8 7 6 5 5 4 4 4 4 5 5 6 7 7 6 5 4 3 3 3 3 4 5 6 7 6 5 4 3 3 2 2 3 3 4 5 6 6 5 4 3 2 1 1 2 3 4 5 6 6 5 4 3 2 1 1 2 3 4 5 6 6 5 4 3 3 2 2 3 3 4 5 6 7 6 5 4 3 3 3 3 4 5 6 7 7 6 5 5 4 4 4 4 5 5 6 7

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function G=getG(M,N,d)
 
  c=[M+1,N+1]/2; %image center
  
  rho=hypot((1:M)'-c(1), (1:N)-c(2));

  G=discretize(rho,0:d:norm([M,N]/2+d) );

end

[Jason]

Thats a great example - thanks.

My images are always odd dimesions.

Thanks again for all your help

Jason

[Matt J]

You can also overlay the annuli boundaries to get more visual information about the pixel memberships

plotG(11,11,1)

plotG(12,12,1)

function plotG(M,N,d)
 
  c=[M+1,N+1]/2; %image center
  
  rho=hypot((1:M)'-c(1), (1:N)-c(2));
  
  R=0:d:norm([M,N])/2+d;
  G=discretize(rho, R);
  
  imagesc(G); colormap(1-gray(256)); axis equal;     hold on
  viscircles(repmat(c,numel(R),1), R(:)  );          hold off
  
end

[Jason]

Thats awesome!

Answer 2

[Image Analyst]

See attached demo.

Adapt as needed, for example change mean to median.

[Jason]

Thanks IA, i'll take a look

Answer 3

[Image Analyst]

@Jason, ok try this - it may be a little more intuitive:

grayImage = imread('cameraman.tif');
[rows, columns, numColorChannels] = size(grayImage);
figure('Name', 'Demo by Image Analyst', 'NumberTitle', 'off')
subplot(2, 1, 1);
imshow(grayImage);
axis('on', 'image');
X = 1 : columns;
Y = 1 : rows;
centerX = columns / 2;
centerY = rows / 2;
[x, y] = meshgrid(X, Y);
radii = sqrt((x - centerX) .^ 2 + (y - centerY) .^ 2);
fprintf('There are %d unique radii in this image.\n', numel(unique(radii)))

There are 5924 unique radii in this image.

% Find out how many histogram bins we'll need to be about every pixel.
maxDistance = max(radii(:))

maxDistance = 181.0193

% So the min distance is 0 and the max distance is maxDistance.
% Let's divide this up into numRings "zones".
numRings = 10; % Whatever you want for the thickness of the rings/zones.
zoneBoundaries = linspace(0, maxDistance, numRings);
% Show the lines separating each ring (zone);
centers = repmat([centerX, centerY], numRings, 1);
viscircles(centers, zoneBoundaries);
% For each ring, get a mask and get the median of the pixel values in the mask.
for k = 1 : numRings-1
    innerRadius = zoneBoundaries(k);
    outerRadius = zoneBoundaries(k+1);
    mask = (radii >= innerRadius) & (radii < outerRadius);
    medianValues(k) = median(grayImage(mask));
    fprintf('The median pixel value between distance %.1f and %.1f is %d gray levels.\n',...
        innerRadius, outerRadius, medianValues(k))
end

The median pixel value between distance 0.0 and 20.1 is 17 gray levels.
The median pixel value between distance 20.1 and 40.2 is 48 gray levels.
The median pixel value between distance 40.2 and 60.3 is 68 gray levels.
The median pixel value between distance 60.3 and 80.5 is 104 gray levels.
The median pixel value between distance 80.5 and 100.6 is 155 gray levels.
The median pixel value between distance 100.6 and 120.7 is 154 gray levels.
The median pixel value between distance 120.7 and 140.8 is 147 gray levels.
The median pixel value between distance 140.8 and 160.9 is 148 gray levels.
The median pixel value between distance 160.9 and 181.0 is 150 gray levels.

% Plot medians
subplot(2, 1, 2);
plot(zoneBoundaries(1:end-1), medianValues, 'b.-', 'LineWidth', 2, 'MarkerSize', 25);
grid on;
title('Median Value as a Function of Radius.')
ylabel('Median Value')
xlabel('Radius')

[Jason]

High IA, many thanks for this. Couple of questions.

1: How did you superimpose the viscircles on the image.

2: I really need the numrings such that each median calculation is at exactly 1 pixel from the previous (radially).Would this just be by setting

numRings=maxDistance?

3: I'd also like to use interp for the smaller ring medians since there will be only a small number of pixels on the circle. Not to sure how to do this.

4: Also for the graphic, I wanted to only show a subset of the real calculated circles, so did this - i think its correct

    numRings=maxDistance??? to enable calculations at each pixel radially?
    zoneBoundaries = linspace(0, maxDistance, numRings);
    zoneBoundaries2 = linspace(0, maxDistance, 20);        % Just for graphic
    % Show the lines separating each ring (zone);
    centers = repmat([centerX, centerY], numRings, 1);
    centers2 = repmat([centerX, centerY], 20, 1);          % Just for graphic 
    viscircles(centers2, zoneBoundaries2);

[Image Analyst]

@Jason

1. The function viscircles is the one that overlays circles over an image. There are lots of options so check it out.
2. To have the boundaries be exactly one pixel apart, you'd set the number of rings to the max distance. Keep in mind though that the distances are floating point, factional distances, except for distance purely vertical or horizontal (or a very few other special cases). Even though the pixel centers lie on integer grid coordinates, their distance is not in general going to be an integer.
3. If you're using median, it doesn't really make sense to do an interpolation, does it? Think about it. And if you did, how many points would you interpolate in between? See attached interpolation demos but I don't recommend you interpolate. What is your end goal anyway? Why are you using a non-linear metric like median instead of a linear one like mean? And interpolation would be mixing computations of nonlinear and linear measurements. If you don't have enough pixels in any of the rings, I suggest you just widen the rings. If you wanted, you could smooth the final profile afterwards. I think we need to think about what is the real world physical reason why you want this measurement and if going to subpixel resolution is really needed or would change anything.
4. To do a subset of rings, you'd simply use indexing to pick out every , say, tenth one: viscircles(centers(1:10:end, :), zoneBoundaries(1:10:end));

[Jason]

Thanks I.A