牙齒的根管影像處理

close all;
clear;
clc;

c=0;
for i=22:42
str = int2str(i);
s=['1_1_1_000001' str  '.dcm'];
c=c+1;
if c==1
filename=s;
[X, map] =dicomread(filename);
info = dicominfo(filename);% 00000155
X = dicomread(info);
figure, imshow(X,'DisplayRange',[]);
% I=double(X);
% figure, imshow(I);
% %colormap(map);
% axis off          % Remove axis ticks and numbers

%axis image        % Set aspect ratio to obtain square pixels

k = waitforbuttonpress;
point1 = get(gca,'CurrentPoint');    % button down detected
finalRect = rbbox;                   % return figure units
point2 = get(gca,'CurrentPoint');    % button up detected
point1 = point1(1,1:2);              % extract x and y
point2 = point2(1,1:2);
p1 = min(point1,point2);             % calculate locations
offset = abs(point1-point2);         % and dimensions
x = [p1(1) p1(1)+offset(1) p1(1)+offset(1) p1(1) p1(1)];
y = [p1(2) p1(2) p1(2)+offset(2) p1(2)+offset(2) p1(2)];
hold on
axis manual
plot(x,y)          
I=X(p1(2):p1(2)+offset(2),p1(1):p1(1)+offset(1));
figure, imshow(I,'DisplayRange',[]);
[m n]=size(I);
Z=zeros(m,n);
else
filename=s;
X = dicomread(filename);
info = dicominfo(filename);% 00000155
X = dicomread(info);
I=X(p1(2):p1(2)+offset(2),p1(1):p1(1)+offset(1));
end
Iw = teethbw(I);
f=double(Iw);
Z=Z+f;

end
B = flipud(Z);
figure, surface(B);
-------------------------------------------------------------------


function [FI] = teethbw(I2)
%  close all;
%  clear;clc;

[m n]=size(I2);
%figure, imshow(I2,'DisplayRange',[]);
%I3 = imadjust(I2);
% figure(2), imshow(I3);
% figure(3), imhist(I3);
%BW1 = edge(I3,'prewitt');
BW2 = edge(I2,'canny');
%figure, imshow(BW1);
% figure, imshow(BW2);

% level = graythresh(I3);
% bw = im2bw(I3,level);

bw=BW2;

%figure, imshow(bw);
[labeled,numObjects] = bwlabel(bw,8);
graindata = regionprops(labeled,'basic')
maxArea = max([graindata.Area])
biggestGrain = find([graindata.Area]==maxArea)
for i=1:m
for j=1:n
 if labeled(i,j)==biggestGrain
    bw(i,j)=1;
 else
     bw(i,j)=0;
 end
end
end

% figure, imshow(bw);

  
originalBW=bw;
se = strel('disk',5);      
erodedBW = imclose(originalBW,se);
% figure,imshow(originalBW), figure, imshow(erodedBW)
%bwAreaOpenBW = bwareaopen(erodedBW,50);
% c = [43 185 212];
% r = [38 68 181];

[r c]= find(erodedBW);
% r=ind(:,1);
% c=ind(:,2);
%BWs = bwselect(I2,c,r,4);
% BWs = roipoly(I2,c,r);
% BWs = roifill(I2,c,r);
%BWs=imcontour(I2,3)
%FI=I2;

FI=I2;