%test - dicom image process =

close all; 

clear all ; 

clc ;

m=dicominfo('bmode.dcm');

X=dicomread(im);

a_1=rgb2gray(X(29:568,235:790,:,1)); %extraction of first image // the range to extract jst the image

a_2=rgb2gray(X(29:568,235:790,:,2)); %extraction of the second image

figure, 

subplot(1,2,1);imshow(X(:,:,:,1));

subplot(1,2,2);imshow(X(:,:,:,2));

impixelinfo;

%figure, 

%i=[1:36];

%montage(X(29:568,235:790,:,i));

%detect feature points 

imagepoints1=detectSURFFeatures(a_1,'MetricThreshold',1000);

imagepoints2=detectSURFFeatures(a_2,'MetricThreshold',1000);

%extract feature discryptors

feature_1=extractFeatures(a_1,imagepoints1);

feature_2=extractFeatures(a_2,imagepoints2);

%plot the extracted features

figure;

subplot(1,2,1);

imshow(a_1);

hold on; 

plot(selectStrongest(imagepoints1,500));

subplot(1,2,2);

imshow(a_2);

hold on; 

plot(selectStrongest(imagepoints2,500));

title('extracted features on the image');

%MATCH EXTRACTED FEATURES

%indexPairs_1=matchFeatures(feature_1,feature_2,'maxRatio',0.9); 

%ommitted as too many outliers are produced

indexPairs_1=matchFeatures(feature_1,feature_2,'Unique',true);

matchedPoints1 = imagepoints1(indexPairs_1(:, 1));

matchedPoints2 = imagepoints2(indexPairs_1(:, 2));

figure; 

showMatchedFeatures(a_1,a_2,matchedPoints1,matchedPoints2);

title('features matched in the original image');

%outliers delete

c=[];

for z=1:422 

    if sumsqr(matchedPoints2(z).Location-matchedPoints1(z).Location)>10

       c=[c,z];

    end

end

matchedPoints1(c)=[];

matchedPoints2(c)=[];

figure; 

showMatchedFeatures(a_1,a_2,matchedPoints1,matchedPoints2);

title('features matched after deletion of outliers');

%create point cloud