% It receives the images containing CA1, CA3, DG and returns the 13
% segmented regions

function [dg_in dg_out_int dg_out_ext dg_out ca1_in ca1_out ca3_in ca3_out ...
    ca1_border_int ca1_border_ext ca3_border_int ca3_border_ext dg_border_int dg_border_ext] = ...
    get_segmented_images(ca1im, ca3im, dgim)

    siz=12;
    
    [ca1_in ca1_out] = binary_segment(ca1im);
    [ca3_in ca3_out] = binary_segment(ca3im);
    [dg_in dg_out] = binary_segment(dgim);
    
    %create and save the border images
    %dg
    dire = 1; norm = 2;
    im = dg_in;
    im(im>0)=1;
    len = size(im,dire);
    wei = repmat((1:len)',1,size(im,norm));
    calc_im = im.*wei;
    res_im = zeros(size(im));  %res_im: contains the mask of the external border
    inn_im = zeros(size(im));  %inn_im: contains the mask of the internal border
    bassolinea=max(calc_im,[],dire);  %coordinates of lower bound
    calc_im(calc_im==0)=size(im,dire);
    altolinea=min(calc_im,[],dire);  %coordinates of upper bound
    for i=siz+1:size(im,norm)-siz
        if bassolinea(i) > 0
            if altolinea(i) > siz && altolinea(i) < len - siz
                res_im(altolinea(i)-siz:altolinea(i)+siz,i-siz:i+siz)=1;
            end
            if bassolinea(i) < len - siz && bassolinea(i) > siz
                res_im(bassolinea(i)-siz:bassolinea(i)+siz,i-siz:i+siz)=1;
            end
        end
    end
    %now, only for dg, calculate the internal border in its own way
    bound = bwboundaries(im,4);
    numPixels = cellfun(@numel,bound);
    [app ix] = max(numPixels);
    sel_bound =bound{ix};
    for i=1:length(sel_bound)
        if sel_bound(i,1) > siz && sel_bound(i,1) < size(im,1) - siz && sel_bound(i,2) > siz && sel_bound(i,2) < size(im,2) - siz %se è nei limiti
            if all(res_im(sel_bound(i,1)-5:sel_bound(i,1)+5,sel_bound(i,2)-5:sel_bound(i,2)+5)==0) %se non è stato messo nel bordo esterno.
              inn_im(sel_bound(i,1)-siz:sel_bound(i,1)+siz,sel_bound(i,2)-siz:sel_bound(i,2)+siz)=1;
            end
        end
    end
    
    dg_border_int = inn_im.*dgim;
    dg_border_ext = res_im.*dgim;
    %ca1
    im = ca1_in;
    im(im>0)=1;
    len = size(im,dire);
    wei = repmat((1:len)',1,size(im,norm));
    calc_im = im.*wei;

    res_im = zeros(size(im));  %res_im: contains the mask of the external border
    inn_im = zeros(size(im));  %inn_im: contains the mask of the internal border
    bassolinea=max(calc_im,[],dire);  %coordinates of lower bound
    calc_im(calc_im==0)=size(im,dire);
    altolinea=min(calc_im,[],dire);  %coordinates of upper bound
    for i=siz+1:size(im,norm)-siz
        if bassolinea(i) > 0
            if altolinea(i) > siz && altolinea(i) < len - siz
                res_im(altolinea(i)-siz:altolinea(i)+siz,i-siz:i+siz)=1;
            end
            if bassolinea(i) < len - siz && bassolinea(i) > siz
                inn_im(bassolinea(i)-siz:bassolinea(i)+siz,i-siz:i+siz)=1;
            end
        end
    end
    ca1_border_int = inn_im.*ca1im;
    ca1_border_ext = res_im.*ca1im;
    %ca3
    dire = 2; norm = 1;
    im = ca3_in;
    im(im>0)=1;
    len = size(im,dire);
    wei = repmat((1:len),size(im,norm),1);
    calc_im = im.*wei;
    res_im = zeros(size(im));  %res_im: contains the external border
    inn_im=zeros(size(im));  %inn_im: contains the internal border
    bassolinea=max(calc_im,[],dire);
    calc_im(calc_im==0)=size(im,dire);
    altolinea=min(calc_im,[],dire);
    for i=siz+1:size(im,norm)-siz
        if bassolinea(i) > 0
            if altolinea(i) > siz && altolinea(i) < len - siz
               res_im(i-siz:i+siz,altolinea(i)-siz:altolinea(i)+siz)=1;
            end
            if bassolinea(i) < len - siz && bassolinea(i) > siz
               inn_im(i-siz:i+siz,bassolinea(i)-siz:bassolinea(i)+siz)=1;
            end
        end
    end
    
    ca3_border_int = inn_im.*ca3im;
    ca3_border_ext = res_im.*ca3im;
    %and, last but not least, divide the outer part of dg
    im=dg_out;
    im(im>0)=1;
    dg_out_int = zeros(size(im)); %internal part output
    dg_out_ext = zeros(size(im)); %external part output
    CC = bwconncomp(im);
    numPixels = cellfun(@numel,CC.PixelIdxList);
    nAreas = sum(numPixels > 1000); %only consider big areas
    mins = zeros(nAreas,1); %minima
    maxs = zeros(nAreas,1); %and maxima position of point in every area
    indexes = zeros(nAreas,1);
    count=0;
    for i=1:numel(numPixels)
        if numPixels(i)>1000
            count=count+1;
            mins(count) = min(CC.PixelIdxList{i});
            maxs(count) = max(CC.PixelIdxList{i});
            indexes(count)=i;
        end
    end
    if nAreas==1
        dg_out_ext(CC.PixelIdxList{indexes(1)}) = dg_out(CC.PixelIdxList{indexes(1)});
    elseif nAreas==2
        [a ix] = min(mins); %ix contains the index of the external part
        dg_out_ext(CC.PixelIdxList{indexes(ix)}) = dg_out(CC.PixelIdxList{indexes(ix)});
        if ix==1, ix=2; else ix=1; end
        dg_out_int(CC.PixelIdxList{indexes(ix)}) = dg_out(CC.PixelIdxList{indexes(ix)});
    elseif nAreas==3
        used = 1:nAreas;
        [a ix1] = min(mins); %ix contains the index of the external part
        dg_out_ext(CC.PixelIdxList{indexes(ix1)}) = dg_out(CC.PixelIdxList{indexes(ix1)});
        used(ix1)=0;
        [a ix2] = max(maxs); %ix contains the index of the external part
        dg_out_ext(CC.PixelIdxList{indexes(ix2)}) = dg_out(CC.PixelIdxList{indexes(ix2)});
        used(ix2)=0;
        [a ix3] = max(used);
        dg_out_int(CC.PixelIdxList{indexes(ix3)}) = dg_out(CC.PixelIdxList{indexes(ix3)});
    else
        used = 1:nAreas;
        [a ix1] = min(mins); %ix contains the index of the external part
        dg_out_ext(CC.PixelIdxList{indexes(ix1)}) = dg_out(CC.PixelIdxList{indexes(ix1)});
        used(ix1)=0;
        [a ix2] = max(maxs); %ix contains the index of the external part
        dg_out_ext(CC.PixelIdxList{indexes(ix2)}) = dg_out(CC.PixelIdxList{indexes(ix2)});
        used(ix2)=0;
        while sum(used)>0
            [a ix3] = max(used);
            dg_out_int(CC.PixelIdxList{indexes(ix3)}) = dg_out(CC.PixelIdxList{indexes(ix3)});
            used(ix3)=0;
        end
    end
	
    
end