function N =  DetectGrapes( I )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% ICBV091 - Grape detection project
%
%% by Eran Geva and Eran Tomer
%
% notes:
% - the function expects the image's filename, not the image itself.
% - ignore k-means warnings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% read image file and separate into different colors
    function [ RAW R G B ] = readfile( name )
        RAW = imread( strcat( name, '.jpg' ), 'jpeg');
        [ xlen ylen zlen ] = size( RAW );
        R = double( RAW( :, :, 1 ) );
        G = double( RAW( :, :, 2 ) );
        B = double( RAW( :, :, 3 ) );
    end

%% convolution that results with the same size matrix
    function m = cnv( m, ker )
        m = conv2( double( m ), double( ker ), 'same' );
    end

%% finding gradient size and direction
    function [ dx dy gs ga ] = grad( m )
        dy = cnv( m, [ 1 ; -1 ] );
        dx = cnv( m, [ 1 , -1 ] );
        [ ga, gs ] = cart2pol( dx, dy );
    end

%% flattening matrix into a vector
    function vec = flat( m )
        vec = double( reshape( m, numel( m ), 1 ) );
    end

%% clustering for all colors in the image
% done using by matlab's K-Means function
    function m = dokmean( R, G, B, km )
        v = [ flat( R ), flat( G ), flat( B ) ];
        m = reshape( kmeans( v, km, 'emptyaction', 'drop' ), size( R ) );
    end

%% crowded mask -
% finding the density of frequent gradients sizes around a point
    function crowd = remove_crowded( orig, fsize, th_h_colors )
        crowd = ones( xlen, ylen );
        maxcolor = max(max(orig));
        for x1 = 1 : xlen
            for y1 = 1 : ylen
                xfrom = max( 1, x1 - fsize ); xto = min( xlen, x1 + fsize );
                yfrom = max( 1, y1 - fsize ); yto = min( ylen, y1 + fsize );
                boxsize = ( xto - xfrom + 1 ) * ( yto - yfrom + 1 );
                h = histc( flat( orig( xfrom : xto, yfrom : yto ) ), 1 : maxcolor );
                h = sort( h, 'descend' );
                crowd( x1, y1 ) = sum( h( 1 : th_h_colors ) ) / boxsize;
            end
        end
    end

%% angle mask
% finding the dispersing of different gradient angles around a point
    function [ ang1, ang2 ] = top_angles( orig, fsize )
        ang1 = ones( xlen, ylen );
        ang2 = ones( xlen, ylen );
        for x1 = 1 : xlen
            for y1 = 1 : ylen
                xfrom = max( 1, x1 - fsize ); xto = min( xlen, x1 + fsize );
                yfrom = max( 1, y1 - fsize ); yto = min( ylen, y1 + fsize );
                boxsize = ( xto - xfrom + 1 ) * ( yto - yfrom + 1 );
                intr = pi/8;
                h = histc( flat( orig( xfrom : xto, yfrom : yto ) ), -pi - intr/2 : intr :pi + intr/2 );
                h = sort( h, 'descend' ); h( 3 ) = 0;
                ang1( x1, y1 ) = h(1) / boxsize;
                ang2( x1, y1 ) = sum(h(1:2)) / boxsize;
            end
        end
    end

%% creating the plus and cross kernel matrix
    function ppp = makeplus( psize, isx )
        psize = ceil( psize );
        ppp = zeros( psize );
        if ~isx
            ppp( ceil( psize / 2 ), : ) = 1; ppp( :, ceil( psize / 2 ) ) = 1;
        else
            for i = 1 : psize
                ppp( i, i ) = 1; ppp( i, psize - i + 1 ) = 1;
            end
        end
    end

%% cross mask - activating plus or x convulution masking with a threshold
    function mm = cross_filter( m, crossize, thresh, isx )
        if isx; thresh = thresh / s2; crossize = crossize / s2; end
        mm = cnv( m, makeplus( crossize, isx ) ) >= thresh;
    end

%% multicross mask
% repeating plus and x convulution masking a given number of times
    function m = multicross( m, crossize, thresh, times )
        for i = 1 : times
            m = m .* cross_filter( m, crossize, thresh, 1 );
            m = m .* cross_filter( m, crossize, thresh, 0 );
        end
    end

%% creating the 2 bagel masks
% the bagel itself and the inside of the bagel
% done to filter points using different regions around it
    function [ inbagel, onbagel ] = make_bagel( boxsize, bagelradius, bagelwidth )
        [ x1 y1 ] = meshgrid( -boxsize:boxsize, -boxsize:boxsize );
        dc = round( sqrt( x1 .^ 2 + y1 .^ 2 ) );
        inbagel = dc < bagelradius - bagelwidth;
        onbagel = abs( dc - bagelradius ) <= bagelwidth;
    end

%% bagel mask
% for each point:
%      count edges inside the bagel around it
%      find the frequent gradient count on the bagel around it
    function [ cin, cpart ] = find_bagels( orig, fsize, bagelradius, bagelwidth, th_h_colors )
        cin = ones( xlen, ylen );
        cpart = ones( xlen, ylen );
        maxcolor = max(max(orig));
        [ inbagel, onbagel ] = make_bagel( fsize, bagelradius, bagelwidth );
        insize = sum(sum(inbagel)); onsize = sum(sum(onbagel));
        orig = [ repmat( 0, fsize, ylen ) ;orig ; repmat( 0, fsize, ylen ) ];
        orig = [ repmat( 0, xlen + fsize * 2, fsize ) orig repmat( 0, xlen + fsize * 2, fsize ) ];
        
        for x1 = 1 : xlen
            for y1 = 1 : ylen
                box = orig( x1 : x1 + fsize*2, y1 : y1 + fsize*2 );
                xsize = min( [ x1+fsize+1 , 2*fsize+1, xlen+fsize-x1 ] );
                ysize = min( [ y1+fsize+1 , 2*fsize+1, ylen+fsize-y1 ] );
                boxratio = fsize * fsize / ( xsize * ysize );
                
                in = box .* inbagel;
                cin( x1, y1 ) = sum(sum(~~in)) / insize * boxratio;
                
                on = box .* onbagel;
                h = histc( flat( on ), 1 : maxcolor );
                h = sort( h, 'descend' );
                cpart( x1, y1 ) = sum( h( 1 : th_h_colors ) ) / onsize * boxratio;
                
            end
        end
    end

%% well... the work itself, within phases.
    function mask = work()
        
        % clustering
        K = dokmean( R, G, B, 10 );
        K = K.^2;
        [ dx dy gs ga ] = grad( K ); gx = double( ~~gs );
        
        % mask initialization
        RGB = R + G + B;
        mask = ones( xlen, ylen );
        
        % color masking
        disp( 'white' ); mask = mask .* ( RGB < 720 );
        disp( 'green' ); mask = mask  & ( G ./ RGB > 0.28 ) & ( G ./ RGB < 0.38 );
        
        % angle mask
        disp( 'angles' ); [ ANG1, ANG2 ] = top_angles( ga, 18 );
        mask = mask .* ( ANG1 > 0.56) .* ( ANG1 < 0.88);
        mask = mask .* ( ANG2 > 0.58) .* ( ANG2 < 0.88);
        
        % crowd mask - simple and histogram based
        disp( 'simple' );SMPL = cnv( gx, ones(21, 21) ) / 441;
        mask = mask .* ( SMPL > 0.14 ) .* ( SMPL < 0.57 );
        disp( 'crowded' ); CROWD = remove_crowded( gs, 10, 2 );
        mask = mask .* ( CROWD > 0.09 ) .* ( CROWD < 0.35 );

        % multicross clean-up
        disp( 'cross' ); mask = multicross( mask, 41, 30, 6 );
        disp( 'cross' ); mask = multicross( mask, 41, 40, 6 );

        % bagels mask - with the average grape size of 19 (9+1+9) pixels
        disp( 'bagels9' ); [ CBGL, CPART ] = find_bagels( gs .* mask, 16, 9, 3, 2 );
        mask = mask .* ( CBGL > 0.02 ) .* ( CBGL < 0.15 );
        mask = mask .* ( CPART > 0.024 ) .* ( CPART < 0.08 );
        
        % multicross clean-up
        disp( 'cross' ); mask = multicross( mask, 41, 30, 6 );
        disp( 'cross' ); mask = multicross( mask, 41, 40, 6 );
    end


s2 = sqrt(2); xlen = 0; ylen = 0;
% if ~exist( 'I', 'var' ); I = '.\labeled\slika1885'; end
[ RAW R G B ] = readfile( I );
N = work();                                                                                                                          %

end







