defining the center of an image and recognise the contrast disks

so i have these images of a test object. i need to make a code that recognises the contrast disks ( the dark dots around the disk), where the disks are actually present around the whole disk as in the image of the test object ( not the x ray ). i used imcrop and ginput to crop a rectangle in the required area, but that is not good enough because on the parts where the contrast disk is not visible it goes awry also not at all standardised. i need matlab to be able to locate the center of the contrast disks around the bigger disk with minimum user input so the results are repeatable. please help! thank you

[Simon Chan]

Why not use a template to locate the disks? Leeds phantom has a fixed geometry, you can first determine the entire phantom center and its angle of rotation. Finally put the template on it, then you know where are the disks.

[sesilia maidelin]

yes, i think that's what the code below is doing. Thank you!

Answer 1

[DGM]

Consider the example:

A = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/890405/image.png');
B = bwareafilt(A<128,3); % pick the 3 objects with largest area
imshow(B)

% get properties
S = regionprops(B,'area','solidity','centroid','minferetproperties','boundingbox');
C = vertcat(S.Centroid);
ara = vertcat(S.Area);
sol = vertcat(S.Solidity);
mind = vertcat(S.MinFeretDiameter);

% try to discern which objects are which
[~,isring] = max((mind.*ara)./sol);
[~,isbox] = max(sol./ara);

% get relevant geometry information
Bring = S(isring).BoundingBox;
Cring = Bring(1:2) + Bring(3:4)/2; % center of ring
Cbox = C(isbox,:); % center of box
Rring = mind(isring)/2; % approx radius of ring
boxloc = (Cbox-Cring).*[1 1]; % relative position vector of box
baserot = atan2d(boxloc(2),boxloc(1)); % base rotation of entire object

% these are constants extracted from a modified copy
% radii are relative to Rring
Rbox = 0.666*Rring;
Rdot = 0.786*Rring;
rdot = 0.0541*Rring;
dotrot = 30;

% find dot locations
dotangles = linspace(dotrot,180-dotrot,9) + [0; 180] + baserot;
dotangles = flipud(reshape(dotangles.',[],1));
[dotx doty] = pol2cart(dotangles*pi/180,Rdot);
dotxy = [dotx doty] + Cring;

imshow(A); hold on
plot(dotxy(:,1),dotxy(:,2),'o')

%% try to get samples by doing a local average around each point
dotsamples = zeros(numel(dotx),1);
sc = size(A);
[xx yy] = meshgrid(1:sc(2),1:sc(1));
for k = 1:numel(dotx)
    mask = ((xx-dotxy(k,1)).^2+(yy-dotxy(k,2)).^2) <= (rdot*0.5)^2;
    samplesize = nnz(mask);
    dotsamples(k) = sum(A(mask))/samplesize;	
end

dotsamples

dotsamples = 18×1

124.4250 127.0500 129.5641 131.3846 136.0250 132.4872 134.2973 134.8649 137.3250 139.5750

This should work even if the object is rotated or moved/scaled. The thresholding might need to be adjusted if the contrast and levels change.

[sesilia maidelin]

oh wow, thank you so much! this should work great. one more question, i have to calculate the contrast of the roi ( the intensity of the pixels inside the ring ) and one of my main problems was that i couldnt get a circular roi. how do i extract the intensity of the pixels inside these rings?

[DGM]

What object are you talking about when you say "ring"? When I was mentioning "ring" in the code, I was referring to the large ring near the perimeter of the object.

If you mean to find the contrast across the union of all the dots, you can do that like:

kdots = range(dotsamples)./mean(dotsamples)

or using whatever contrast metric you want.

If you want to find the contrast within each dot, you'd have to do that similar to how you'd do the sampling:

dotsampk = zeros(numel(dotx),1);
sc = size(A);
[xx yy] = meshgrid(1:sc(2),1:sc(1));
for k = 1:numel(dotx)
    mask = ((xx-dotxy(k,1)).^2+(yy-dotxy(k,2)).^2) <= (rdot*0.5)^2;
    samp = double(A(mask));
    dotsampk(k) = range(samp)./mean(samp);	
end
dotsampk

Although note that the sample radius is half the estimated dot radius. That's simply to avoid dealing with alignment issues due to limited position accuracy.

[sesilia maidelin]

hey there, sorry i am not clear on something. by dotsamples, do you mean the intensity around the small blue ring or inside it? i was not very coherent when i read your code because of covid. I want to get the contrast between the intensity of the pixels inside the ring and the intensity of the area around said ring. something like this below (sorry for the ugly drawing). Thank you so much for the code, i was able to make good progress with it.

[sesilia maidelin]

hi there, sorry theres a confusion for my question. so i meant the contrast between the area inside the blue ring and a little bit outside of it. so the contrast between that tiny area and the area around it. here's an ugly drawing to explain it. I am going to continue to tinker w your code, but if you could help me with this i would be very thankful!

[DGM]

This finds the contrast between the average value of the inner region of the dot and an annular region around the dot (i.e. its neighboring background). This should reduce the influence of any broader image shading.

First, go ahead and process the image to find the dots. This is the same as before.

A = imread('diskxr_rot2.png');
B = bwareafilt(A<128,3); % pick the 3 objects with largest area
%imshow(~B)

% get properties
S = regionprops(B,'area','solidity','centroid','minferetproperties','boundingbox');
C = vertcat(S.Centroid);
ara = vertcat(S.Area);
sol = vertcat(S.Solidity);
mind = vertcat(S.MinFeretDiameter);

% try to discern which objects are which
[~,isring] = max((mind.*ara)./sol);
[~,isbox] = max(sol./ara);

% get relevant geometry information
Bring = S(isring).BoundingBox;
Cring = Bring(1:2) + Bring(3:4)/2; % center of ring
Cbox = C(isbox,:); % center of box
Rring = mind(isring)/2; % approx radius of ring
boxloc = (Cbox-Cring).*[1 1]; % relative position vector of box
baserot = atan2d(boxloc(2),boxloc(1)); % base rotation of entire object

% these are constants extracted from a modified copy
% radii are relative to Rring
Rbox = 0.666*Rring;
Rdot = 0.786*Rring;
rdot = 0.0541*Rring;
dotrot = 30;

% find dot locations
dotangles = linspace(dotrot,180-dotrot,9) + [0; 180] + baserot;
dotangles = flipud(reshape(dotangles.',[],1));
[dotx doty] = pol2cart(dotangles*pi/180,Rdot);
dotxy = [dotx doty] + Cring;

imshow(A); hold on
plot(dotxy(:,1),dotxy(:,2),'yo')

Then process the dots to find their contrast against their local background.

% get contrast between dots and their neighboring bg
% adjust the mask diameters as you see fit  
% i chose these to allow some misalignment
dotsamples = zeros(numel(dotx),2);
sc = size(A);
[xx yy] = meshgrid(1:sc(2),1:sc(1));
for k = 1:numel(dotx)
    ringmask = ((xx-dotxy(k,1)).^2+(yy-dotxy(k,2)).^2) <= (rdot*2.2)^2 ...
        & ((xx-dotxy(k,1)).^2+(yy-dotxy(k,2)).^2) >= (rdot*1.5)^2;
    mask = ((xx-dotxy(k,1)).^2+(yy-dotxy(k,2)).^2) <= (rdot*0.5)^2;
    
    dotsamples(k,:) = [sum(A(mask))/nnz(mask) sum(A(ringmask))/nnz(ringmask)];
end

kdotsrg = range(dotsamples,2)./mean(dotsamples,2)

kdotsrg = 18×1

0.1219 0.1070 0.0873 0.0695 0.0356 0.0722 0.0543 0.0496 0.0269 0.0074

kdotsst = std(dotsamples,0,2)./mean(dotsamples,2)

kdotsst = 18×1

0.0862 0.0757 0.0617 0.0491 0.0252 0.0510 0.0384 0.0351 0.0190 0.0052

[DGM]

Regarding the question you asked about how the locations are determined:

Once a new image is thresholded, there are only two objects of interest -- the dark ring near the perimeter, and the dark off-center square. The location of the phantom is represented by the center of the ring. The rotation of the phantom is determined by the location of the square with respect to the ring center.

Once the center and base rotation angle are known, the angular displacement of the dots is known. The angle between the square and the first dot is 30d. The angle between each dot is 15d. The radial position of the dots is known relative to the ring radius, so the ring radius will resolve that.

In short, the only information that's needed from the incoming image are the ring center and radius, and the square's location.

Of course, one might ask where those precalculated constants came from. I started by editing a copy of the original image to make it easier to segment.

% extract geometry constants from simplified copy
A = imread('diskxrt.png');
B = bwareafilt(A<128,3);

imshow(B)

S = regionprops(B,'area','solidity','centroid','minferetproperties','boundingbox');
C = vertcat(S.Centroid);
ara = vertcat(S.Area);
sol = vertcat(S.Solidity);
mind = vertcat(S.MinFeretDiameter);

[~,isring] = max((mind.*ara)./sol);
[~,isdot] = max(sol./ara);
isbox = 3;

Bring = S(isring).BoundingBox;
Cring = Bring(1:2) + Bring(3:4)/2; % ring center calculated as center of bbox
%Cring = C(isring,:); % ring center as calculated from image centroid
Cbox = C(isbox,:);
Rring = mind(isring)/2;

boxloc = (Cbox-Cring).*[1 -1];
baserot = atan2d(boxloc(2),boxloc(1));

% these are relative to Rring (as calc via Cring/Bring)
Rbox = norm(Cbox-Cring)/Rring % 0.6474 0.6661

Rbox = 0.6661

Rdot = norm(C(isdot,:)-Cring)/Rring % 0.7670 0.7853

Rdot = 0.7853

rdot = S(isdot).MinFeretDiameter/2/Rring % 0.0541 0.0541

rdot = 0.0541

dotloc = (C(isdot,:)-Cring).*[1 -1];
dotrot = atan2d(dotloc(2),dotloc(1)) - baserot % 30.285 29.499

dotrot = 29.4990

dotangles = linspace(dotrot,180-dotrot,9) + [0; 180] 

dotangles = 2×9

29.4990 44.6242 59.7495 74.8747 90.0000 105.1253 120.2505 135.3758 150.5010 209.4990 224.6242 239.7495 254.8747 270.0000 285.1253 300.2505 315.3758 330.5010

hold on;
plot(C(isring,1),C(isring,2),'yo')
plot(C(isdot,1),C(isdot,2),'*')
plot(C(isbox,1),C(isbox,2),'s')

[sesilia maidelin]

ah i see, thank you so much for your explanation! huge appreciation!