% function object_properties

close all;
clear all;
clc;
addpath('C:\Users\joanna\Documents\IMG_PROC\matlab.contrib');
% path = 'C:\Users\joanna\Documents\AAAData_PrePost_zip\NEW\10000\DICOM\256732_ANONYMO\12710952\7';
path='C:\Users\joanna\Documents\AAAData_PrePost_zip\NEW\10000\DICOM\178281_ANONYMO\12700304\4';
% find all the files in it
av_files = dir(fullfile(path,'*.dcm'));

% put all the files into one series
img=[];
for k=2:20 %size(av_files,1)
    name = av_files(k,1).name;
    X = dicomread(fullfile(path,name));
    img=cat(3,img,X);
end

% figure; imshow(img(:,:,10),[]);
Y = dicominfo(fullfile(path,av_files(1,1).name));

img_temp=img(:,:,14);
figure; imshow(img_temp,[]);

%% Segment all the circular objects
% perform the thresholding procedure
tr = 0.26*max(img(:));
img_tmp=img;
img_tmp(img<tr)=0;
img_tmp(img_tmp~=0)=1;
img_kl=img_tmp;
figure; imshow(img_kl(:,:,14),[]);

%% morphological erosion in 3-D
se(:,:,1)=[0 0 0; 0 1 0;0 0 0];
se(:,:,2)=ones(3);
se(:,:,3)=se(:,:,1);
img_er=imerode(img_kl,se);
figure; imshow(img_er(:,:,14),[]);

%% 3-D holes filling
img_ert=~img_er;
bl=bwlabeln(img_ert,6);
lab=bl(1,1,1);
bl(bl~=lab)=0;
bl(bl~=0)=1;
bl=~bl;
img_fill3=img_er|bl;
figure; subplot(1,2,1); imshow(img_er(:,:,14),[]); subplot(1,2,2); imshow(img_fill3(:,:,14),[]);

%% 2-D holes filling
img_anl=~img_er(:,:,14);
bl=bwlabel(img_anl,4);
lab=bl(1,1);
bl(bl~=lab)=0;
bl(bl~=0)=1;
bl=~bl;
img_fill=img_er(:,:,14)|bl;
figure; subplot(1,2,1); imshow(img_er(:,:,14),[]); subplot(1,2,2); imshow(img_fill,[]);

%% cut the ROI from the image (using boundingbox)
% take into consideration one object only
figure; imshow(img_fill,[]);
[y,x]=ginput(1);
x=round(x);
y=round(y);
bl=bwlabel(img_fill,4);
lab=bl(x,y);
bl(bl~=lab)=0;
bl(bl~=0)=1;
figure; imshow(bl,[]);
rp=regionprops(bl,'BoundingBox')

%% cut the ROI
rp=rp.BoundingBox;
roi=bl(floor(rp(2)):ceil(rp(2)+rp(4)),floor(rp(1)):ceil(rp(1)+rp(3)));
figure; imshow(roi,[]);

%% find the object in image with the biggest area
rp=regionprops(img_fill,'Area'); % we do not obtain the labels
bl=bwlabel(img_fill,4);
rp=regionprops(bl,'Area');
rpa=struct2array(rp)';
f=find(rpa==max(rpa(:)));
bl(bl~=f)=0;
bl(bl~=0)=1;
% figure; imshow(bl,[]);

[B,L] = bwboundaries(bl,'noholes');
figure; imshow(img_temp,[]);

hold on
for k = 1:length(B)
    boundary = B{k};
    plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1)
end

%% estimate the object properties
% Euler number
bl=bwlabel(img_er(:,:,14),8); % CoherentObjects-Holes
rp=regionprops(bl,'EulerNumber'); 
rpe=struct2array(rp)'

bl=bwlabel(img_fill,8);
rp=regionprops(bl,'EulerNumber'); 
rpe=struct2array(rp)'

%% Aspect ratio
% ar=d_min/d_max; a finctio of the largest and the smallest diameter 
% estimate the aspect ration of all the objects in image
[bl,lab]=bwlabel(img_fill,8);
rp=regionprops(bl,'MajorAxisLength'); 
dmax=struct2array(rp)';
rp=regionprops(bl,'MinorAxisLength');
dmin=struct2array(rp)';
ar=dmax./dmin;
% keep the objects with ar value <=2
lab=1:lab';
f=find(ar<=1.5);
lab=lab(f);

bll=zeros(size(bl));
for k=1:length(lab)
    bl_temp=bl;
    bl_temp(bl_temp~=lab(k))=0;
    bl_temp(bl_temp~=0)=1;
    bll=bll|bl_temp;
end
[B,L] = bwboundaries(bll,'noholes');
figure; imshow(img_temp,[]);

hold on
for k = 1:length(B)
    boundary = B{k};
    plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1)
end

%%
% 1. detect the longest object in the sagital aorta image
% path = 'C:\Users\joanna\Documents\AAAData_PrePost_zip\NEW\10000\DICOM\256732_ANONYMO\12710952\7';
% av_files = dir(fullfile(path,'*.dcm'));
% name = av_files(30,1).name;
% img = dicomread(fullfile(path,name));
% figure; imshow(img,[]);
% ................. %

%% Orientation
[bl,lab]=bwlabel(img_fill,8);
rp=regionprops(bl,'Orientation');
or=struct2array(rp);

% 1. detect the aorta in the sagital aorta image
% using the aspect ratio and the orientation features
% .................... %

%% Circularity
% fcirc = 4piArea/perimeter^2
% 3. Find the 20% of the most circular objects in the image, which area is greater than 10pix 
% .................... %

%% Center of gravity
[bl,lab]=bwlabel(img_fill,8);
rp=regionprops(bl,'Centroid');
cg=struct2cell(rp)';
cg=cell2mat(cg);
cg=round(cg);

% plot center of gravity of all the objects in the image
bl=zeros(size(img_fill));
f=sub2ind(size(img_fill),cg(:,2),cg(:,1));
bl(f)=1;
[B,L] = bwboundaries(bl,'noholes');
figure; imshow(img_fill,[]);

hold on
for k = 1:length(B)
    boundary = B{k};
    plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1)
end

%% Convex hull
% find the convexhull of the aorta
[bl,lab]=bwlabel(img_fill,8);
figure; imshow(bl,[]);
[y,x]=ginput(1);
x=round(x);
y=round(y);
lab=bl(x,y);
bl(bl~=lab)=0;
bl(bl~=0)=1;

rp=regionprops(bl,'ConvexHull');
rp=rp.ConvexHull;
% create a convexhull mask
ch = bwconvhull(bl,'objects',8);

% plot the results
[B,L] = bwboundaries(ch,'noholes');
figure; imshow(bl,[]);

hold on
for k = 1:length(B)
    boundary = B{k};
    plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2)
end

%% Convexity
% ratio of the projection area bad the convex hull (mostly limited to 0.99)
% 4. Find all the convex objects
% .................... %

% 5. Try to find in the literature other shape factors usefull for your
% project and implement them
% .................... %
