function[labl, Am, neighbors, maskImg, maskImgB]=slic(im, k, m, seRadius, colopt, imB)


    if ~exist('colopt','var') || isempty(colopt), colopt = 'mean'; end
    if ~exist('mw','var')     || isempty(mw),         mw = 0;      end
    if ~exist('nItr','var')   || isempty(nItr),     nItr = 10;     end
    
    if exist('eim', 'var'), USEDIST = 0; else USEDIST = 1; end
        
    MEANCENTRE = 1;
    MEDIANCENTRE = 2;
    
    if strcmp(colopt, 'mean')
        centre = MEANCENTRE;
    elseif strcmp(colopt, 'median')
        centre = MEDIANCENTRE;        
    else
        error('Invalid colour centre computation option');
    end
    
    [rows, cols, chan] = size(im);
    if chan ~= 3
        error('Image must be colour');
    end
    
    % Convert image to L*a*b* colourspace.  This gives us a colourspace that is
    % nominally perceptually uniform. This allows us to use the euclidean
    % distance between colour coordinates to measure differences between
    % colours.  Note the image becomes double after conversion.  We may want to
    % go to signed shorts to save memory.
    im = rgb2lab(im); 
    imB = rgb2lab(imB);
    [mr, nc, pd] = size(im);

    maskImg = false(mr,nc);
    maskImgB = false(mr,nc);
    
    
    % Apply median filtering to colour components if mw has been supplied
    % and/or non-zero
    if mw
        if length(mw) == 1
            mw(2) = mw(1);  % Use same filtering for L and chrominance
        end
        for n = 1:3
            im(:,:,n) = medfilt2(im(:,:,n), [mw(1) mw(1)]);
            imB(:,:,n) = medfilt2(imB(:,:,n), [mw(1) mw(1)]);
        end
    end
    
    % Nominal spacing between grid elements assuming hexagonal grid
    S = sqrt(rows*cols / (k * sqrt(3)/2));
    
    % Get nodes per row allowing a half column margin at one end that alternates
    % from row to row
    nodeCols = round(cols/S - 0.5);
    % Given an integer number of nodes per row recompute S
    S = cols/(nodeCols + 0.5); 

    % Get number of rows of nodes allowing 0.5 row margin top and bottom
    nodeRows = round(rows/(sqrt(3)/2*S));
    vSpacing = rows/nodeRows;

    % Recompute k
    k = nodeRows * nodeCols;
    
    % Allocate memory and initialise clusters, labels and distances.
    C = zeros(6,k);          % Cluster centre data  1:3 is mean Lab value,
                             % 4:5 is row, col of centre, 6 is No of pixels
    labl = -ones(rows, cols);   % Pixel labels.
    d = inf(rows, cols);     % Pixel distances from cluster centres.
    
    % Initialise clusters on a hexagonal grid
    kk = 1;
    r = vSpacing/2;
    
    for ri = 1:nodeRows
        % Following code alternates the starting column for each row of grid
        % points to obtain a hexagonal pattern. Note S and vSpacing are kept
        % as doubles to prevent errors accumulating across the grid.
        if mod(ri,2), c = S/2; else c = S;  end
        
        for ci = 1:nodeCols
            cc = round(c); rr = round(r);
            C(1:5, kk) = [squeeze(im(rr,cc,:)); cc; rr];
            c = c+S;
            kk = kk+1;
        end
        
        r = r+vSpacing;
    end
    
    % Now perform the clustering.  10 iterations is suggested but I suspect n
    % could be as small as 2 or even 1
    S = round(S);  % We need S to be an integer from now on
    
    for n = 1:nItr
       for kk = 1:k  % for each cluster

           % Get subimage around cluster
           rmin = max(C(5,kk)-S, 1);   rmax = min(C(5,kk)+S, rows); 
           cmin = max(C(4,kk)-S, 1);   cmax = min(C(4,kk)+S, cols); 
           subim = im(rmin:rmax, cmin:cmax, :);  
           assert(numel(subim) > 0)
           
           % Compute distances D between C(:,kk) and subimage
           if USEDIST

               D = dist(C(:, kk), subim, rmin, cmin, S, m);
           else
               D = dist2(C(:, kk), subim, rmin, cmin, S, m, eim, We);
           end

           % If any pixel distance from the cluster centre is less than its
           % previous value update its distance and label
           subd =  d(rmin:rmax, cmin:cmax);
           subl =  labl(rmin:rmax, cmin:cmax);
           updateMask = D < subd;
           subd(updateMask) = D(updateMask);
           subl(updateMask) = kk;
           
           d(rmin:rmax, cmin:cmax) = subd;
           labl(rmin:rmax, cmin:cmax) = subl;           
       end
       
       % Update cluster centres with mean values
       C(:) = 0;
       for r = 1:rows
           for c = 1:cols
              tmp = [im(r,c,1); im(r,c,2); im(r,c,3); c; r; 1];
              C(:, labl(r,c)) = C(:, labl(r,c)) + tmp;
           end
       end
       
       % Divide by number of pixels in each superpixel to get mean values
       for kk = 1:k 
           C(1:5,kk) = round(C(1:5,kk)/C(6,kk)); 
       end
       
       % Note the residual error, E, is not calculated because we are using a
       % fixed number of iterations 
    end
    
    % Cleanup small orphaned regions and 'spurs' on each region using
    % morphological opening on each labeled region.  The cleaned up regions are
    % assigned to the nearest cluster. The regions are renumbered and the
    % adjacency matrix regenerated.  This is needed because the cleanup is
    % likely to change the number of labeled regions.
    if seRadius
        [labl, Am] = mcleanupregions(labl, seRadius);
    else
        labl = makeregionsdistinct(labl);
        [labl, minLabel, maxLabel] = renumberregions(labl);
        Am = regionadjacency(labl);    
    end
    
    
    % Recompute the final superpixel attributes and write information into
    % the Sp struct array.
    N = length(Am);
    disp(N)
    labels = unique(labl);
    nLabels = length(labels);
    neighbors = cell(nLabels,1);
    stdevA = cell(N,1);
    stdevB= cell(N,1);

    imgA_Sp = struct('imgA_L', cell(1,N), 'imgA_A', cell(1,N), 'imgA_B', cell(1,N), ...
                    'std_A_L', cell(1,N), 'std_A_A', cell(1,N), 'std_A_B', cell(1,N),...
                      'r', cell(1,N), 'c', cell(1,N), 'N', cell(1,N));
    imgB_Sp= struct('imgB_L', cell(1,N), 'imgA_A', cell(1,N), 'imgB_B', cell(1, N),...
                   'std_B_L', cell(1,N), 'std_B_A', cell(1,N), 'std_B_B', cell(1,N),...
                    'ro',cell(1,N),'co', cell(1,N), 'N', cell(1, N));              

    [X,Y] = meshgrid(1:cols, 1:rows);
    imgA_L = im(:,:,1);    
    imgA_A = im(:,:,2);    
    imgA_B = im(:,:,3);
    imgB_L=imB(:,:,1);
    imgB_A=imB(:,:,2);
    imgB_B=imB(:,:,3);
 
    for n = 1:N;
       mask = labl==n;
       nm = sum(mask(:));
        if centre == MEANCENTRE     
            imgA_Sp(n).imgA_L = sum(imgA_L(mask))/nm;
            imgA_Sp(n).imgA_A= sum(imgA_A(mask))/nm;
            imgA_Sp(n).imgA_B = sum(imgA_B(mask))/nm;

            imgB_Sp(n).imgB_L = sum(imgB_L(mask))/nm;
            imgB_Sp(n).imgB_A = sum(imgB_A(mask))/nm;
            imgB_Sp(n).imgB_B = sum(imgB_B(mask))/nm;
            
        elseif centre == MEDIANCENTRE
            imgA_Sp(n).imgA_L = median(imgA_L(mask));
            imgA_Sp(n).imgA_A = median(imgA_A(mask));
            imgA_Sp(n).imgA_B = median(imgA_B(mask));

            imgB_Sp(n).imgB_L = median(imgB_L(mask));
            imgB_Sp(n).imgB_A= median(imgB_A(mask));
            imgB_Sp(n).imgB_B = median(imgB_B(mask)); 
        end
        
        imgA_Sp(n).r = sum(Y(mask))/nm;
        imgA_Sp(n).c = sum(X(mask))/nm;
        
        % Compute standard deviations of the colour components of each super
        % pixel. This can be used by code seeking to merge superpixels into
        % image segments.  Note these are calculated relative to the mean colour
        % component irrespective of the centre being calculated from the mean or
        % median colour component values.
        
        %Standard deviation of image A
        imgA_Sp(n).std_A_L = std(imgA_L(mask));
        imgA_Sp(n).std_A_A = std(imgA_A(mask));
        imgA_Sp(n).std_A_B = std(imgA_B(mask));
        imgA_Sp(n).N = nm;  % Record number of pixels in superpixel 
        
        %Standard deviation of image B
        imgB_Sp(n).std_B_L = std(imgB_L(mask)); 
        imgB_Sp(n).std_B_A = std(imgB_A(mask));
        imgB_Sp(n).std_B_B = std(imgB_B(mask));
        imgB_Sp(n).ro = imgA_Sp(n).r;
        imgB_Sp(n).co = imgA_Sp(n).c;
        imgB_Sp(n).N= nm; % Record numbers of pixels in superpixels of image B.

        %Neighnbourhood of superpixels 
        msk = labl == labels(n);
        adjacentPixelMask = imdilate(msk,true(3)) & ~msk;
        neighbors{n} = unique(labl(adjacentPixelMask));

       %Comparison of standard deviation and making of mask image
        if(imgA_Sp(n).std_A_L<=imgB_Sp(n).std_B_L)
             stdevA{n}= 0;
             stdevB{n}=1;
             maskImgB(mask)= true;
             maskImg(mask)= false;
             
        else
               stdevA{n}=1;
               stdevB{n}=0;  
               maskImg(mask) = true;  
               maskImgB(mask)= false;                                               
        end
     
     
    end
    
    