Create ROIs with a radius and at a 20° angle

10 views (last 30 days)
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

Accepted Answer

Image Analyst
Image Analyst on 3 Jul 2022
See the FAQ
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);
  8 Comments
Latifa Bouguessaa
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
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')

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by