clear
clc
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% list directory
nameOfDir = '../imgss/img/'
nameOfOutDir = '../imgss/imgResBestPlaces/'
nameOfMaskDir = '../imgss/imgMask/'

global show;
if (1 == 0)
    show = 0;
    indexes = 25
else
    show = 0;
    indexes = 1:184
end

measures = ones(1,184);


for (index = indexes)
    
    
list=dir(nameOfDir); 
filenames={list(~[list.isdir]).name}

name = filenames{index}

I = imread(strcat(nameOfDir,name));
imageOriginal = rgb2gray(I);

imageOriginal = imresize(imageOriginal, measures(1, index));
imageOriginalO = imageOriginal;





h = fspecial('average', [2,6])
%imageOriginal = medfilt2(imageOriginal, [2,6]);
imageOriginal = imfilter(imageOriginal, h, 'replicate');


BW = edge(imageOriginal,'canny', [0.19 0.61])

se = strel('disk',1);
BWf = imdilate(double(BW), se);
h = fspecial('average', [5,5])
BWfx = imfilter(BWf, h, 'replicate');
BWfx = 1 - BWfx;


for(x = 1:size(BW,2))
    for (y = 1:size(BW,1))
        imageOriginal(y,x) = imageOriginal(y,x) * BWfx(y,x);
    end
end



if (show == 1)
figure
imshow(imageOriginal);
figure
imshow(BWfx);
end

%return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%preprocess mask

%list=dir(nameOfMaskDir); 
%filenames={list(~[list.isdir]).name}

%name = filenames{index}

I = imread(strcat(nameOfMaskDir,name));
if (size(I,3) > 1)
image = rgb2gray(I);
else
    image = I
end

image = imresize(image, measures(1, index));

imMaskDisp = image;

edgeMy = zeros(1, size(I,2));


for(x = 2:size(I,2))
    for (y = 2:size(I,1))
        if (image(y,x) == 255)
            edgeMy(1,x) = y;
            
            break;
        else
            imMaskDisp(y,x) = 0;
        end
        
    end
end


if (show == 1)
figure (555)
imshow(imMaskDisp);
end

b1 = 1;
b2 = size(I,2);
for(y = round(size(I,2)/2):size(I,2))
   if (edgeMy(1,y) == 0)
       b2 = y-1;
       break
   end
end
for(y = round(size(I,2)/2):-1:1)
   if (edgeMy(1,y) == 0)
       b1 = y+1;
       break
   end
end


h = fspecial('average', [10,10])
imSmooth = imfilter(imageOriginal, h, 'replicate');



iter = 0;
diff = 6
maxGrad = 0;
grads = []
body = [];
bodyFin = [];
for(x = b1:5:b2-diff)
    x1 = diff;
    y1 = edgeMy(1,x+diff) - edgeMy(1,x);
    
    %kolmice
    y2 = x1;
    x2 = -y1;
    
    
    
    iter = iter+1;
        
    body(iter, 1) =  round(edgeMy(1,x)+y2*-2);
    body(iter, 2) =  round(x + x2*-2);
    body(iter, 3) =  round(edgeMy(1,x)+y2*2);
    body(iter, 4) =  round(x + x2*2);


    [cx,cy,c] = improfile(imSmooth, [body(iter,2),body(iter,4)],[body(iter,1), body(iter,3)]);
    
    iMax = max(c);
    iMin = min(c);
    
    grad = abs(iMax - iMin);
    
    if (grad > 10)
        grads = [grads grad];
    end
    
end

maxGrad = quantile(grads,0.3);


image = imageOriginalO;
imageA = imageOriginalO;
imageS = imageOriginalO;

for(x = b1:5:b2-diff)
    clear c
    
    x1 = diff;
    y1 = edgeMy(1,x+diff) - edgeMy(1,x);
    
    %kolmice
    y2 = x1;
    x2 = -y1;
    
    
    
    iter = iter+1;
        
    body(iter, 1) =  round(edgeMy(1,x)+y2*-2);
    body(iter, 2) =  round(x + x2*-2);
    body(iter, 3) =  round(edgeMy(1,x)+y2*2);
    body(iter, 4) =  round(x + x2*2);
    
    bodyFin(iter, 1) =  round(edgeMy(1,x));
    bodyFin(iter, 2) =  round(x);
    bodyFin(iter, 3) =  0;
    bodyFin(iter, 4) =  round(edgeMy(1,x)+y2*-2);
    bodyFin(iter, 5) =  round(x + x2*-2);
    bodyFin(iter, 6) =  round(edgeMy(1,x)+y2*2);
    bodyFin(iter, 7) =  round(x + x2*2);


    
    
   
    [cx,cy,c] = improfile(imSmooth, [body(iter,2),body(iter,4)],[body(iter,1), body(iter,3)]);
    
    iMax = max(c);
    iMin = min(c);
     
      
   if ((abs(iMax - iMin) > maxGrad))
        imageA(round(edgeMy(1,x)), round(x)) = 255;
    
        [cx,cy,c] = improfile(imageOriginal, [body(iter,2),body(iter,4)],[body(iter,1), body(iter,3)]);

        [ret, points] = analyzeProfile(c);
        if (ret == 1)
            imageS(round(edgeMy(1,x)), round(x)) = 255;
            bodyFin(iter, 3) =  1;

        end
        
        if (show == 1)
            figure (100)
            imshow(imageS);
        end


    end
end
    
    


itermax = iter
for (iter = 1 + 1:itermax - 1)
   if (bodyFin(iter-1, 3) == 1 &&       bodyFin(iter, 3) == 0 &&        bodyFin(iter+1, 3) == 1)
       
       bodyFin(iter , 3) = 1;
       
   end
end


numberOf1 = 0;

for (iter = 1:itermax)
   if (bodyFin(iter, 3) == 1)
        numberOf1 = numberOf1 + 1;
   else
        %check how many
        if (numberOf1 < 4)
            fr = iter - numberOf1 - 1;
            if (fr < 1) fr = 1; end;
            
            for (insIter1 = fr: iter)
                bodyFin(insIter1, 3) = 0;
            end
        end
        
        numberOf1 = 0;
    end
end


for (iter = 1:itermax)  
    if (bodyFin(iter, 1) == 0)
        bodyFin(iter, 3) == 0
    end
end


imageF = imageOriginalO;
for (iter = 1:itermax)  
   if (bodyFin(iter, 3) == 1)
        imageF(bodyFin(iter, 1), bodyFin(iter, 2)) = 255;
   end
end



imageT = imageOriginalO;
for (iter = 1:itermax)  
   if (bodyFin(iter, 3) == 1)
        if (bodyFin(iter, 1) ~= 0)
            [cx,cy,c] = improfile(imageOriginal, [bodyFin(iter,2+3),bodyFin(iter,4+3)],[bodyFin(iter,1+3), bodyFin(iter,3+3)]);
            
            [ret, points] = analyzeProfile(c);
            
            if (points(1,1) ~= 0 && points(2,1) ~= 0)
                imageT(round(cy(points(1,1))), round(cx(points(1,1)))) = 255;
                
                bodyFin(iter, 8) =  (cy(points(1,1)));
                bodyFin(iter, 9) =  (cx(points(1,1)));

                imageT(round(cy(points(2,1))), round(cx(points(2,1)))) = 255;
                                
                bodyFin(iter, 10) =  (cy(points(2,1)));
                bodyFin(iter, 11) =  (cx(points(2,1)));
            else
                bodyFin(iter, 3) = 0;
            end
            
        end
   end
end


iter = 1;
segments = []
totaldists = []
while (iter < itermax)  
    if (bodyFin(iter, 3) == 1)
        dists = [];
        iterFirst = iter;
        while ((bodyFin(iter, 3) == 1) && iter < itermax)  
            LIy = bodyFin(iter, 8);
            LIx = bodyFin(iter, 9);

            MAy = bodyFin(iter, 10);
            MAx = bodyFin(iter, 11);

            if (LIx ~= 0 && LIy ~= 0 && MAx ~= 0 & MAy ~= 0)
                dist = sqrt((LIx - MAx)^2 + (LIy - MAy)^2); 
                dists = [dists, dist];
            end
            
            
            iter = iter + 1
        end
        
        % analyze dists
        meanDist = mean(dists);
        totaldists = [totaldists,dists]
        
         segments = [segments, [iterFirst, iter, meanDist]'];
    else
        iter = iter + 1;
    end
    
end


meanAbs = trimmean(totaldists,40);
for (iter = 1:size(segments,2))
   meanm = segments( 3, iter);
   
   thr = meanAbs/2;
   

   if ((meanm - meanAbs) > thr || ( meanAbs - meanm) > thr)
       for (iterIntr = segments( 1, iter):segments(2, iter))     
            bodyFin(iterIntr, 3) = 0;
       end
   end
end
   


iter = 1;
segments = []
totaldists = []
while (iter < itermax)  
    if (bodyFin(iter, 3) == 1)
        dists = [];
        iterFirst = iter;
%        while ((bodyFin(iter, 3) == 1) && iter < itermax)  
            LIy = bodyFin(iter, 8);
            LIx = bodyFin(iter, 9);

            MAy = bodyFin(iter, 10);
            MAx = bodyFin(iter, 11);

            if (LIx ~= 0 && LIy ~= 0 && MAx ~= 0 & MAy ~= 0)
                dist = sqrt((LIx - MAx)^2 + (LIy - MAy)^2); 
                dists = [dists, dist];
            end
            
            
            iter = iter + 1
 %       end
        
        % analyze dists
        meanDist = mean(dists);
        totaldists = [totaldists,dists]
        
         segments = [segments, [iterFirst, iter, meanDist]'];
    else
        iter = iter + 1;
    end
    
end


%if more segments, and they differs significantly, take only the strongest
%one
meanAbs = trimmean(totaldists,40);
for (iter = 1:size(segments,2))
   meanm = segments( 3, iter);
   
   thr = meanAbs;
   if ((meanm - meanAbs) > thr || ( meanAbs - meanm) > thr/4.0)
       for (iterIntr = segments( 1, iter):segments(2, iter))     
            bodyFin(iterIntr, 3) = 0;
       end
   end
end














measured_dists = [];
imageTT = imageOriginalO;
for (iter = 1:itermax)  
   if (bodyFin(iter, 3) == 1)
       LIy = bodyFin(iter, 8);
       LIx = bodyFin(iter, 9);

       MAy = bodyFin(iter, 10);
       MAx = bodyFin(iter, 11);
       
       if (LIx ~= 0 && LIy ~= 0 && MAx ~= 0 && MAy ~= 0)
            imageTT(round(LIy), round(LIx)) = 255;
            imageTT(round(MAy), round(MAx)) = 255;
            
            measured_dist = sqrt((LIx - MAx)^2 + (LIy - MAy)^2)/measures(1, index); 
            measured_dists = [measured_dists, measured_dist];
       end
   end
end



meanmeasured_dists = mean(measured_dists);

outName = strcat(nameOfOutDir, name, '.txt')
save(outName, 'measured_dists','meanmeasured_dists');



imwrite([imageOriginalO, imageA, imageS, imageF, imageT, imageTT], [nameOfOutDir name '.bmp'], 'bmp')





end



















