# Create ROIs with a radius and at a 20° angle

6 views (last 30 days)
Latifa Bouguessaa on 3 Jul 2022
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

Image Analyst on 3 Jul 2022
See the FAQ
imshow(grayImage);
hold on;
[rows, columns, numberOfColorChannels] = size(grayImage);
xCenter = columns/2;
yCenter = rows/2;
allTheta = 0 : 20 : 359;
% 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);
Latifa Bouguessaa on 6 Jul 2022
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
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;
% 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')

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by