Create ROIs with a radius and at a 20° angle

Good day,

I hope someone can help me here

I have a DICOM image (Phanton) and I need 18 ROIs (15X15) along a radius (see attached picture)

The radius and center of the phantom are known

Please give me ideas on how to do this

Warm greetings

Latifa

Answer 1

[Image Analyst]

See the FAQ

https://matlab.fandom.com/wiki/FAQ#How_do_I_create_a_circle?

then adapt it like this. Modify my adaptation to use sizes of your choosing.

grayImage = imread('cameraman.tif');
imshow(grayImage);
hold on;
[rows, columns, numberOfColorChannels] = size(grayImage);
xCenter = columns/2;
yCenter = rows/2;
allTheta = 0 : 20 : 359;
radius = 110;
% Show the circle
viscircles([xCenter, yCenter], radius, 'Color','y', 'LineWidth', 2);
boxWidth = 28;
for k = 1 : length(allTheta)
    % Get this theta
    theta = allTheta(k);
    % Get the point on the circle
    x = radius * cosd(theta) + xCenter;
    y = radius * sind(theta) + yCenter;
    % Get the upper left
    xul = x - boxWidth/2;
    yul = y - boxWidth/2;
    pos = [xul, yul, boxWidth, boxWidth];
    rectangle('Position', pos, 'EdgeColor', 'r', 'LineWidth', 2)
    axis square;
    grid on;
end

You can use pos with imcrop() if you want to crop out that little box into its own subimage.

subImage = imcrop(grayImage, pos);

Or you can get the lower right coordinates from

    xul = round(x - boxWidth/2);
    yul = round(y - boxWidth/2);
    xlr = round(xCenter + boxWidth/2);
    ylr = round(yCenter + boxWidth/2);
    subImage = grayImage(yul : yur, xul : xur);

[Image Analyst]

@Latifa Bouguessaa's "Answer" moved here since it's not an answer but really a comment to me:

Hello Sir

Thank you very much for your help.

Now I have one more question. I was able to read in the CT image, segment the phantom in the image and determine the center of the phantom. and with your code I was able to create the ROIs. Now I wanted to sum noise power spectrum from each roi and at the end.

The noise power spectrum for each ROI can be determined like this:

-------------------------------------------------- ----------------------------------------------

Deltax=Spasi(1); %pixel distance in x direction
Deltay=Spasi(2); %pixel distance in y direction
f=(Deltax*Deltay)/(2*n*m); % size of ROI
spe1=mean(x);
spe=abs(fft(spe1)).^2;
nps=f*pe;
[ab]=size(x);
spaki=Spasi(1,1);
sumbu=spaki:spaki:b*spaki;
plot(sumbu, 10*log(nps), 'LineWidth',0.5)

-------------------------------------------------- -------------------------------------------------- -----

How can I now calculate the power noise spectrum for each ROI and then sum them all up?

Thank you again.

Your Latifa

[Image Analyst]

I'd look into pwelch and put that in for each subimage.

To index nps do

nps(k) = f*pe

If my circle code solved your original question, can you please click the "Accept this Answer" link? Thanks in advance. 🙂

[Latifa Bouguessaa]

Hello Sir, thank you for your help. I think I didn't explain the task well. the PNS must be calculated for each ROI. code to calculate the PNS is there. I just don't know how to implement it in the loop. Thank you very much

[Image Analyst]

Make a function to do it, then put it in the loop

grayImage = imread('cameraman.tif');
imshow(grayImage);
hold on;
[rows, columns, numberOfColorChannels] = size(grayImage);
xCenter = columns/2;
yCenter = rows/2;
allTheta = 0 : 20 : 359;
radius = 110;
% Show the circle
viscircles([xCenter, yCenter], radius, 'Color','y', 'LineWidth', 2);
boxWidth = 28;
for k = 1 : length(allTheta)
    % Get this theta
    theta = allTheta(k);
    % Get the point on the circle
    x = radius * cosd(theta) + xCenter;
    y = radius * sind(theta) + yCenter;
    % Get the upper left
    xul = x - boxWidth/2;
    yul = y - boxWidth/2;
    pos = [xul, yul, boxWidth, boxWidth];
    rectangle('Position', pos, 'EdgeColor', 'r', 'LineWidth', 2)
    axis square;
    grid on;
    subImage = imcrop(grayimage, pos);
    pns(k) = ComputePns(subImage); % You write this function.
end

[Latifa Bouguessaa]

Hello Sir.

thank you for your help.

I tried but failed.

now I wanted to sum the gray values of the entire ROIs and display them as one image.

Do you have an idea how to realize that?

Warm greetings

[Image Analyst]

grayImage = imread('cameraman.tif');
imshow(grayImage);
hold on;
[rows, columns, numberOfColorChannels] = size(grayImage);
xCenter = columns/2;
yCenter = rows/2;
allTheta = 0 : 20 : 359;
radius = 110;
% Show the circle
viscircles([xCenter, yCenter], radius, 'Color','y', 'LineWidth', 2);
boxWidth = 28;
for k = 1 : length(allTheta)
    % Get this theta
    theta = allTheta(k);
    % Get the point on the circle
    x = radius * cosd(theta) + xCenter;
    y = radius * sind(theta) + yCenter;
    % Get the upper left
    xul = x - boxWidth/2;
    yul = y - boxWidth/2;
    pos = [xul, yul, boxWidth, boxWidth];
    rectangle('Position', pos, 'EdgeColor', 'r', 'LineWidth', 2)
    axis square;
    grid on;
    subImage = imcrop(grayimage, pos);
    theMeans(k) = sum(subImage(:));
end

but theMeans is a vector of 20 elements long. How do you plan on making an image out of it?

[Latifa Bouguessaa]

thanks SIR for the help

I have to calculate the noise power spectrum for each ROI

[Latifa Bouguessaa]

Hello Sir

Thanks for the quick solutions

I have now written my code in such a way that I calculate the NSP using the formula (see appendix) at the end. But I'm not entirely satisfied with the NPS curve. I do not know why.

thank you again

My Code:

clear all

clc

% ------------------------------------------------------------

% ------------------------------------------------------------

% Open CT image

Citra = dicomread('Testdaten_0180.dcm');%Please edit this file name

info = dicominfo('Testdaten_0180.dcm'); %Please edit this file name

Potong=info.RescaleIntercept;

Miring=info.RescaleSlope;

FV=info.ReconstructionDiameter;

Spasi=info.PixelSpacing;

imshow(Citra,'DisplayRange',[]);

% ------------------------------------------------------------

% ------------------------------------------------------------

% Transform CT data to HU

[a b]=size(Citra);

D1=int16(Citra);

D2=zeros(a,b);

D2=int16(D2);

for k=1:a

for l=1:b

D2(k,l)=D1(k,l)*Miring+Potong;

end

end

figure(1)

imshow(D2,'DisplayRange',[]);

% ------------------------------------------------------------

% ------------------------------------------------------------

% Automatic segmentation

B=zeros(a,b);

for i=1:a

for j=1:b

if D2(i,j) < (800-1024)

B(i,j)=0;

else

B(i,j)=1;

end

end

end

Obyek=sum(sum(B));

if Obyek > 1

BR = bwmorph(B,'remove');

BL=bwlabel(BR);

NM=max(max(BL));

for put=1:NM

CB=zeros(a,b);

for i=1:a

for j=1:b

if BL(i,j) == put

CB(i,j)=1;

end

end

end

DB = imfill(CB,'holes');

LA(put)=bwarea(DB);

end

terbesar=max(LA);

for iter=1:NM

if (LA(iter)==terbesar)

urutan =iter;

end

end

CB=zeros(a,b);

for i=1:a

for j=1:b

if BL(i,j) == urutan

CB(i,j)=1;

end

end

end

DB = imfill(CB,'holes');

else

DB=zeros(a,b);

end

figure(2)

imshow(D2,'DisplayRange',[]);

DBP = bwmorph(DB,'remove');

[n m]=size(DBP);

p=0;

for i=1:n

for j=1:m

if DBP(i,j)==1

p=p+1;

x(p)=j;

y(p)=i;

end

end

end

hold on

scatter(x,y,'w', '.')

% ------------------------------------------------------------

% ------------------------------------------------------------

% Determination center of the phantom image

x0=round(a/2); y0=round(b/2);

N=0;

xpos=0; ypos=0;

for i=1:a

for j=1:b

if DB(i,j) > 0

N=N+1;

xpos=xpos+i;

ypos=ypos+j;

end

end

end

yb=round(xpos/N);

xb=round(ypos/N);

plot(xb,yb,'w.', 'MarkerSize', 30)

[rows, columns, numberOfColorChannels] = size(D2);

xCenter = xb;

yCenter = yb;

allTheta = 0 : 10 : 359;

radius = 100;

% Show the circle

viscircles([xCenter, yCenter], radius, 'Color','y', 'LineWidth', 2);

boxWidth = 19;

co = 1;

for k = 1 : length(allTheta)

% Get this theta

theta = allTheta(k);

% Get the point on the circle

x = radius * cosd(theta) + xCenter;

y = radius * sind(theta) + yCenter;

% Get the upper left

xul = x - boxWidth/2;

yul = y - boxWidth/2;

pos = [xul, yul, boxWidth, boxWidth];

rectangle('Position', pos, 'EdgeColor', 'r', 'LineWidth', 2)

axis square;

grid on;

subImage(:,:,co) = imcrop(D2, pos);

co=co+1;

end

figure(4)

for i=1:36

subplot(6,6,i)

imagesc(subImage(:,:,i))

end

%Rauschleistungsspektrum berechnen

Deltax=Spasi(1);

Deltay=Spasi(2);

f=(Deltax*Deltay)/(2*size(subImage,1)*size(subImage,2));

for i=1:size(subImage,3)

spe1=mean(subImage(:,:,i));

spe=abs(fft(spe1)).^2;

nps(i,:)=f.*spe;

end

nps = mean(nps);

%nps=movmean(nps,3);

[a b]=size(subImage(:,:,1));

figure (5)

spaki=Spasi(1,1);

SpatialFrequency=(1/((2*boxWidth+1)*spaki));

[f g]=size(nps);

sumbu=0:SpatialFrequency:(g-1)*SpatialFrequency;

plot(sumbu,nps ,'LineWidth',1.2)

xlabel('Spatial Frqequen (mm)^-1')

ylabel('NPS')

set(gcf, 'color', 'w')