File:Haus vom Nikolaus all vert.gif

function Haus_vom_Nikolaus()
% The aim is to find all "Eulerian paths" in an house shaped
% undirected graph. This means all edges must be used,
% but no edge can be used twice.
% The start node (1), or root, is at the left foundation.
% The recursive backtracking algorithm finds 44 solution were all edges are
% used. 10 times the algorithm backtracks when not all edges are used
% but the algorithm has no further options.
% All valid solution end in node (2).
%
% by Jahobr 2019-04-04
 
[pathstr,fname] = fileparts(which(mfilename)); % save files under the same name and at file location
cd(pathstr) % move to save folder
 
%     (5)
%   8/   \7
%   /     \
% (4)--6--(3)
%  : \   / :
%  :  \ /  :
%  3   X   5
%  : 2/ \4 :
%  : /   \ :
% (1)--1--(2)
 
% left foundation
possibeConnections(1).neighbours = [2 3 4]; % (nodes or vertices)
possibeConnections(1).edge       = [1 2 3]; % using edge
 
% right foundation
possibeConnections(2).neighbours = [1 3 4]; % (nodes or vertices)
possibeConnections(2).edge       = [1 5 4]; % using edge
 
% right roof-edge
possibeConnections(3).neighbours = [1 2 4 5]; % (nodes or vertices)
possibeConnections(3).edge       = [2 5 6 7]; % using edge
 
% left roof-edge
possibeConnections(4).neighbours = [1 2 3 5]; % (nodes or vertices)
possibeConnections(4).edge       = [3 4 6 8]; % using edge
 
% roof ridge
possibeConnections(5).neighbours = [3 4]; % (nodes or vertices)
possibeConnections(5).edge       = [7 8]; % using edge
 
nodeCoordinates = [-1  1  1 -1  0 ...  % X
    ;              -1 -1  1  1  2];    % Y
 
houseWidth = max(nodeCoordinates(1,:)) - min(nodeCoordinates(1,:));
houseHight = max(nodeCoordinates(2,:)) - min(nodeCoordinates(2,:));
 
nodeCompact = nodeCoordinates*0.85; % smaller set of coordinates; for plots only
nodeCompact(2,3:end) = nodeCompact(2,3:end)+0.05; % for plots only
nodeCompact(end) = nodeCompact(end)+0.1; % for plots only
 
% define root state
currentState.edgeUsed =  false(1,8); % init
currentState.nodeList =  NaN(1,8); % init
currentState.edgeList =  NaN(1,8); % init
currentState.penPos = 1; % current position of the pencil
currentState.recursiveDepth = 0; % current depth in the tree
 
[fullSolutionListNodes, fullSolutionListEdges] = recursiveHausVomNikolaus(currentState,possibeConnections); % call recursive function
nSolutions = size(fullSolutionListNodes,1);
fullSolutionListNodes = [ones(nSolutions,1) fullSolutionListNodes]; % we started at node1; this was not considered so far
isEulerian = all(~isnan(fullSolutionListNodes),2); % is this a valid solution?
disp(fullSolutionListNodes);
disp(fullSolutionListEdges);
fprintf('We found %i eulerian and %i non-Eulerian solutions (%i in total)\n', sum(isEulerian), sum(~isEulerian), nSolutions)
 
 
RGB.white      = [1.0 1.0 1.0];
RGB.black      = [0.0 0.0 0.0];
RGB.brightGrey = [0.7 0.7 0.7];
RGB.darkGrey   = [0.4 0.4 0.4];
RGB.redGrey    = [0.9 0.4 0.4];
RGB.darkRed    = [0.7 0.0 0.0];
 
RGB = structfun(@(q)round(q*255)/255, RGB, 'UniformOutput',false); % round to values that are nicely uint8 compatible
 
nColors = 30; % number of colors on gradient
colMap = copper(round(nColors*1.5));  colMap = colMap(1:nColors,:); % create colormap
 
videoFileNames = {}; % init
 
for mode = {'all','eulerian','non-eulerian'} %
    for ori = {'vert'} %,'hori'
        saveName = [fname '_' mode{1} '_' ori{1}];
        for form = {'gif','mp4','svg'} %
            switch form{1} %
                case {'mp4'}
                    framesPerEdge = 48;
                    FontName = 'Arial';
                    padding = char(8239);
                case {'gif'}
                    framesPerEdge = 24;
                    FontName = 'Arial';
                    padding = char(8239);
                case {'svg'}
                    framesPerEdge = 24;
                    FontName = 'Helvetica';
                    padding = '';
                otherwise
                    error('unknown mode')
            end
            
            switch mode{1}
                case {'all'} % all dead ends
                    solutionListNodes = fullSolutionListNodes; % full list
                    solutionListEdges = fullSolutionListEdges; % full list
                    if strcmp(ori{1},'hori') % horizontal image arrangement
                        nHouseRows    =  5; % nSolutions=54; 5*11=55 (+1 legend)
                        nHouseColumns = 11;
                    else % 'vert' % vertical image arrangement
                        nHouseRows    = 11; % nSolutions=54; 5*11=55 (+1 legend)
                        nHouseColumns =  5;
                    end
                case {'eulerian'}  % only full "valide" path, Eulerian paths
                    solutionListNodes = fullSolutionListNodes(isEulerian,:); % reduced list
                    solutionListEdges = fullSolutionListEdges(isEulerian,:); % reduced list
                    if strcmp(ori{1},'hori') % horizontal image arrangement
                        nHouseRows    = 5; % nSolutions=44; 5*9=45 (+1 legend)
                        nHouseColumns = 9;
                    else % 'vert' % vertical image arrangement
                        nHouseRows    = 9; % nSolutions=44; 5*9=45 (+1 legend)
                        nHouseColumns = 5;
                    end
                case {'non-eulerian'}'  % invalide, non eulerian paths
                    solutionListNodes = fullSolutionListNodes(~isEulerian,:); % reduced list
                    solutionListEdges = fullSolutionListEdges(~isEulerian,:); % reduced list
                    if strcmp(ori{1},'hori')  % not really image arrangement but number of frames
                        nHouseRows    = 2; % nSolutions=10; 2*5=10 (no legend)
                        nHouseColumns = 5;
                    else % 'vert'   % not really image arrangement but number of frames
                        nHouseRows    = 5; % nSolutions=10; 2*5=10 (no legend)
                        nHouseColumns = 2;
                    end
                otherwise
                    error('unknown mode')
            end
            [Y_grid, X_grid] = meshgrid( nHouseRows:-1:1, 1:nHouseColumns ); % all solutions will be plotted in a grid; horizontal for mp4
            maxFrames = framesPerEdge*(8+(2*0.25)); % the maximal number pf printed spline points for full Eulerian paths
            stateList = round(linspace(0.75,9.25,maxFrames));
            colIndexList = ceil(linspace(eps,nColors,maxFrames));
            nSolutions = size(solutionListNodes,1); % update
            
            horzSpacing = houseWidth*1.35;
            vertSpacing = houseHight*1.3;
            X_grid = X_grid.*horzSpacing;
            Y_grid = Y_grid.*vertSpacing;
            
            nFrames = 0;
            for iSol = 1:nSolutions % calculate and store the plot Spline
                
                isNotNan = ~isnan(solutionListNodes(iSol,:)); % is this a valid solution?
                nNodes = sum(isNotNan);
                SolutionCoordinates = nodeCoordinates(:,solutionListNodes(iSol,isNotNan)); % extract the xy-coordinate list of the reached nodes
                
                SolutionCoordInter = interp1(1:nNodes,SolutionCoordinates',1:0.5:nNodes)'; % add 50% intermediate points on the ceter of the edges, to get nice spline curves (creates (2*nNodel)-1 points)
                
                SolutionCoordInter = [SolutionCoordInter(:,1) SolutionCoordInter SolutionCoordInter(:,end)]; % add extra points at start/end to get the spline nice (creates (2*nNodel)+1 points)
                
                SplinePathXY = spcrv(SolutionCoordInter,... % Spline curve by uniform subdivision
                    3, 2000); % order and minimal number of interpolation points
                % The spline has a consistent number of points between
                % [start  all-edge-25%-&-75%-points end] (so 2*nNode-1 sections)
                % The points are a little more crammed at the start
                % If we look between each node:  N1_start  N   N  ...  N  N_End
                %  the ratio of spline points is:      1.25  1   1   1  1.25
                % The total number of spline point is determined by "spcrv"
                nPointsInThisSpline = size(SplinePathXY,2);
                nSplinePointsNeeded = framesPerEdge*(nNodes-1+(2*0.25));
                
                plotValues{iSol} = interp1q((1:nPointsInThisSpline)', SplinePathXY', linspace(1,nPointsInThisSpline,nSplinePointsNeeded)')'; % resample X with known number of points
                nFrames = max(nFrames,nSplinePointsNeeded); % longest spline curve determines the length of the animation
            end
            
            figHandle = figure(192838);
            clf
            set(figHandle,'Units','pixel');
            set(figHandle,'MenuBar','none',  'ToolBar','none'); % free real estate for a maximally large image
            set(figHandle,'Color', RGB.white); % white background
            axesHandle = axes;
            hold(axesHandle,'on')
            axis equal
            axis off % invisible axes (no ticks)
            drawnow;
            
            xLimits = [X_grid(1  )-0.75*houseWidth X_grid(end)+0.75*houseWidth];
            yLimits = [Y_grid(end)-0.60*houseHight Y_grid(1  )+0.80*houseHight];
            xlim(xLimits);
            ylim(yLimits);
            set(axesHandle,'Position',[0 0 1 1]); % stretch axis as big as figure, [x y width height]
            
            xRange = xLimits(2)-xLimits(1);
            yRange = yLimits(2)-yLimits(1);
            
            screenSize = get(groot,'Screensize')-[0 0 5 20]; % [1 1 width height] (minus tolerance for figure borders)
            screenAspectRatio = screenSize(3)/screenSize(4); % width/height
            imageAspectRatio = xRange/yRange;
            switch form{1} %
                case {'mp4','svg'}
                    if imageAspectRatio > screenAspectRatio % width will be the problem
                        xSize = screenSize(3);
                        ySize = xSize/imageAspectRatio; % gif width
                    else % height will be the problem
                        ySize = screenSize(4);
                        xSize = ySize*imageAspectRatio; % gif width
                    end
                    xSize = floor(xSize/16)*16; ySize = floor(ySize/16)*16; % full pixels, multiple-of-16 boundaries (MPEG splits the video into 16x16 squares called macroblocks)
                    scaleReduction = 1;
                case {'gif'}
                    MegaPixelTarget = 51*10^6; % Category:Animated GIF files exceeding the 50 MP limit
                    MegaPixelTarget = MegaPixelTarget*nFrames/maxFrames; % rescale if only shorter paths to match "non-eulerian.gif" and "eulerian.gif"
                    pxPerImage = MegaPixelTarget/nFrames; % pixel per gif frame
                    ySize = sqrt(pxPerImage/imageAspectRatio); % gif height
                    xSize = ySize*imageAspectRatio; % gif width
                    xSize = floor(xSize); ySize = floor(ySize); % full pixels
                    if imageAspectRatio > screenAspectRatio % width will be the problem
                        scaleReduction = floor(screenSize(3)/xSize); % repeat as often as possible
                    else % height will be the problem
                        scaleReduction = floor(screenSize(4)/ySize); % repeat as often as possible
                    end
            end
            linSc = xSize*scaleReduction/nHouseColumns/120;
            
            
            figPos = [1 1 xSize*scaleReduction ySize*scaleReduction]; % big start image for antialiasing later [x y width height]
            set(figHandle, 'Position', figPos);
            if ~all(get(figHandle, 'Position') == figPos);   error('figure Position could not be set');   end % check
            drawnow;
            
            if numel(X_grid) > nSolutions % room for legend
                SolutionCoordinates = nodeCompact(:,[1 4 5 3 4 2 3 1 2]); % one full solution example
                plot(X_grid(end)+SolutionCoordinates(1,:), Y_grid(end)+SolutionCoordinates(2,:),'-','Linewidth',1*linSc,'Color',RGB.darkGrey ); % thin line for legend
            end
            
            for iSol = 1:numel(X_grid)
                xOff = X_grid(iSol);
                yOff = Y_grid(iSol);
                for iNode = 1:5
                    drawCirclePatch(xOff+nodeCompact(1,iNode), yOff+nodeCompact(2,iNode), houseWidth/6, RGB.white, linSc, RGB.darkGrey)
                end
            end
            
            textX = linspace(nodeCoordinates(1,1)-0.1, nodeCoordinates(1,2)+0.1, 9); % prepare text position list
            
            for iSol = 1:nSolutions
                xOff = X_grid(iSol); % plot offset
                yOff = Y_grid(iSol); % plot offset
                
                isNotNan = ~isnan(solutionListNodes(iSol,:)); % is this a valid solution?
                nNodes = sum(isNotNan);
                
                if nNodes == 9 % is this a valide solution?
                    linCol = RGB.brightGrey;
                    txtCol = RGB.black;
                else
                    linCol = RGB.redGrey;
                    txtCol = RGB.darkRed;
                end
                
                if ~strcmp(form{1},'svg') % will not be visible in svg
                    plot(xOff+plotValues{iSol}(1,:), yOff+plotValues{iSol}(2,:),'-','LineWidth',5*linSc,'Color', linCol); % grey background line
                end
                
                for iText = 1:9 % node order as text
                    if ~isnan(solutionListNodes(iSol,iText))
                        txtStr  = num2str(solutionListNodes(iSol,iText));
                    else
                        txtStr  = '-'; % happens in non-Eulerian paths
                    end
                    
                    txtHand(iSol,iText) = text(xOff+textX(iText), yOff-1.43, txtStr,... % node order as text
                        'FontName',FontName,...
                        'Color',txtCol, 'FontUnits','pixel', 'FontSize',floor(20*linSc),...
                        'HorizontalAlignment','center', 'VerticalAlignment','middle');
                end
                
                text(xOff, yOff+1.33, ['(' padding num2str(iSol) padding ')'],... % "House Number"
                    'FontName',FontName,...
                    'Color',txtCol, 'FontUnits','pixel', 'FontSize',floor(18*linSc),...
                    'HorizontalAlignment','center', 'VerticalAlignment','middle');
            end
            
            if numel(X_grid) > nSolutions % room for legend
                for iText = 1:5 % node number as legend in the lower right
                    text(X_grid(end)+nodeCompact(1,iText), Y_grid(end)+nodeCompact(2,iText), num2str(iText),...
                        'FontName',FontName,...
                        'Color',txtCol, 'FontUnits','pixel', 'FontSize',floor(20*linSc),...
                        'HorizontalAlignment','center', 'VerticalAlignment','middle');
                end
            end
            
            %% save animation
            drawnow
            f = getframe(figHandle);
            switch form{1} % Init frame(s)
                case {'mp4'}
                    videoFileNames{end+1} = saveName;
                    vid = VideoWriter(videoFileNames{end},'MPEG-4'); % Prepare the new file
                    vid.FrameRate = 50;
                    vid.Quality = 100; % quality in [0% 100%]
                    open(vid);
                    for iRep = 1:50
                        writeVideo(vid,f.cdata); % save video frames
                    end
                case {'gif'}
                    reducedRGBimage = uint8(ones(ySize,xSize,3,nFrames)); % allocate
                    reducedRGBimage(:,:,:,1) = imReduceSize(f.cdata,scaleReduction); % the size reduction: adds antialiasing
            end
            
            
            currentColIndex = 0;
            for iFrame = 2:nFrames % animate by drawing lines
                if currentColIndex ~= colIndexList(iFrame) % start a new line segment with different color
                    segmentStart = max(1,iFrame-2); % update (-2 for overlap of segments; looks better)
                    for iSol = 1:nSolutions
                        if size(plotValues{iSol},2) >= iFrame % not all path have the same length
                            xOff = X_grid(iSol);
                            yOff = Y_grid(iSol);
                            colHand(iSol) = plot(xOff+plotValues{iSol}(1,segmentStart:iFrame), yOff+plotValues{iSol}(2,segmentStart:iFrame), '-','LineWidth',5*linSc,'Color',colMap(colIndexList(iFrame),:)); % fill the grey template line with color changing line segments
                        end
                    end
                    currentColIndex = colIndexList(iFrame); % update
                else % expand the current line segment
                    for iSol = 1:nSolutions
                        if size(plotValues{iSol},2) >= iFrame % not all path have the same length
                            xOff = X_grid(iSol);
                            yOff = Y_grid(iSol);
                            set(colHand(iSol),'XData',xOff+plotValues{iSol}(1,segmentStart:iFrame),'YData',yOff+plotValues{iSol}(2,segmentStart:iFrame)) % fill the grey template line with color changing line segments
                        end
                    end
                end
                for iSol = 1:nSolutions % make path-text bold
                    if ~isnan(solutionListNodes(iSol,stateList(iFrame)))
                        set(txtHand(iSol,stateList(iFrame)),'FontWeight','bold');
                    end
                end
                
                
                %% save animation
                
                switch form{1} %
                    case {'mp4'}
                        drawnow
                        f = getframe(figHandle);
                        writeVideo(vid,f.cdata); % save video frame
                    case {'gif'}
                        drawnow
                        f = getframe(figHandle);
                        reducedRGBimage(:,:,:,iFrame) = imReduceSize(f.cdata,scaleReduction); % the size reduction: adds antialiasing
                end
                
            end % iFrame = 2:nFrames % animate by drawing lines
            
            
            switch form{1} %
                case {'mp4'}
                    for iRep = 1:50
                        writeVideo(vid,f.cdata); % save video frame
                    end
                    close(vid); disp(['MP4 completed: ' videoFileNames{end}   '.mp4']);
                    
                    
                case {'gif'}
                    if strcmp(mode{1},'eulerian')  % only full "valide" path, Eulerian paths
                        RGBcurrent = rmfield(RGB,{'redGrey','darkRed'}); % remove unused colors
                    else
                        RGBcurrent = RGB;
                    end
                    startMap = cell2mat(struct2cell(RGBcurrent)); % struct2colormap; % list of map colors that are not allowed to be changed
                    
                    map = createImMap(reducedRGBimage,64,startMap); % colormap
                    
                    im = uint8(ones(ySize,xSize,1,nFrames)); % allocate
                    for iFrame = 1:nFrames
                        im(:,:,1,iFrame) = rgb2ind(reducedRGBimage(:,:,:,iFrame),map,'nodither'); % rgb to colormap image
                    end
                    
                    imwrite(im,map,fullfile(pathstr, [saveName '.gif']),'DelayTime',1/25,'LoopCount',inf) % save gif
                    disp([saveName '.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit
                    
                case {'svg'} %
                    drawnow
                    if ~isempty(which('plot2svg'))
                        plot2svg(fullfile(pathstr, [saveName '.svg']),figHandle) % by Juerg Schwizer
                    else
                        disp('plot2svg.m not available; see http://www.zhinst.com/blogs/schwizer/');
                    end
            end
            
        end % form = {'gif','mp4'} %
    end % ori = {'hori','vert'}
end % mode = {'all','eulerian','non-eulerian'}
 
for iVid = 1:numel(videoFileNames)
    disp(['Conversion of ' videoFileNames{iVid} ' from MP4 to WebM (Matlab currently does not support webm)']);
    if exist('ffmpeg.exe','file') == 2 % check if it is a valid path to an existing file
        try
            dosCommand = ['!ffmpeg -i ' videoFileNames{iVid} '.mp4 '...   % Command and source
                '-c:v libvpx-vp9 -crf 30 -b:v 0 -deadline best '... % constant (best) quality; see: https://trac.ffmpeg.org/wiki/Encode/VP9
                '-metadata title="' videoFileNames{iVid} '" '... % see: https://wiki.multimedia.cx/index.php/FFmpeg_Metadata
                '-metadata author="Jahobr" -metadata copyright="CC0 1.0 Universal Public Domain Dedication" '...
                '-metadata license="CC0 1.0 Universal Public Domain Dedication" '...
                videoFileNames{iVid} '.webm'];                            % conversion target
            eval(dosCommand); % run conversion from mp4 to webm
            delete(fullfile(pathstr,[videoFileNames{iVid} '.mp4'])); % delete mp4 version after conversion (to save disc space)
        catch ME
            disp(ME)
        end
    else
        disp(['"ffmpeg.exe" not available in path "' pathstr '" (or other problem); see https://ffmpeg.zeranoe.com/builds/']);
    end
end
 
 
function [solutionListNodes, solutionListEdges] = recursiveHausVomNikolaus(currentState,possibeConnections)
%    5
%   / \
%  4---3
%  : X :
%  1---2
 
currentState.recursiveDepth = currentState.recursiveDepth+1; % current position
solutionListNodes = NaN(0,8); % init empty matrix
solutionListEdges = NaN(0,8); % init empty matrix
newSolutionExplored = false; % default
for possibilityIndex = 1:length(possibeConnections(currentState.penPos).edge) % for all connections
    edgeNum = possibeConnections(currentState.penPos).edge(possibilityIndex); % edge number used in the case
    
    if ~currentState.edgeUsed(edgeNum) % if edge not jet used
        currentStateRecursion = currentState; % prepare to use this edge
        currentStateRecursion.edgeUsed(edgeNum) = true; % use it
        currentStateRecursion.penPos = possibeConnections(currentState.penPos).neighbours(possibilityIndex); % set to next point
        currentStateRecursion.nodeList(currentState.recursiveDepth) = currentStateRecursion.penPos;
        currentStateRecursion.edgeList(currentState.recursiveDepth) = edgeNum;
        
        [solutionListNodesTemp, solutionListEdgesTemp]  = recursiveHausVomNikolaus(currentStateRecursion,possibeConnections);
        newSolutionExplored = true;
        solutionListNodes = [solutionListNodes; solutionListNodesTemp];
        solutionListEdges = [solutionListEdges; solutionListEdgesTemp];
    end
end
if ~newSolutionExplored % no further edges possible; end of recurstion; store valide Solution; -> Backtracking
    % if all(currentState.edgeUsed) % restirct to valide solutions only (all edges used); {not used: create full list; separation is done later}
    solutionListNodes = currentState.nodeList; % store all dead-ends as Solutions
    solutionListEdges = currentState.edgeList; % store all dead-ends as Solutions
    % end
end
 
function im = imReduceSize(im,redSize)
% Input:
%  im:      image, [imRows x imColumns x nChannel x nStack] (unit8)
%                      imRows, imColumns: must be divisible by redSize
%                      nChannel: usually 3 (RGB) or 1 (grey)
%                      nStack:   number of stacked images
%                                usually 1; >1 for animations
%  redSize: 2 = half the size (quarter of pixels)
%           3 = third the size (ninth of pixels)
%           ... and so on
% Output:
%  im:     [imRows/redSize x imColumns/redSize x nChannel x nStack] (unit8)
%
% an alternative is: imNew = imresize(im,1/reduceImage,'bilinear');
%        BUT 'bicubic' & 'bilinear'  produces fuzzy lines
%        IMHO this function produces nicer results as "imresize"
 
[nRow,nCol,nChannel,nStack] = size(im);
 
if redSize==1;  return;  end % nothing to do
if redSize~=round(abs(redSize));             error('"redSize" must be a positive integer');  end
if rem(nRow,redSize)~=0;     error('number of pixel-rows must be a multiple of "redSize"');  end
if rem(nCol,redSize)~=0;  error('number of pixel-columns must be a multiple of "redSize"');  end
 
nRowNew = nRow/redSize;
nColNew = nCol/redSize;
 
im = double(im).^2; % brightness rescaling from "linear to the human eye" to the "physics domain"; see youtube: /watch?v=LKnqECcg6Gw
im = reshape(im, nRow, redSize, nColNew*nChannel*nStack); % packets of width redSize, as columns next to each other
im = sum(im,2); % sum in all rows. Size of result: [nRow, 1, nColNew*nChannel]
im = permute(im, [3,1,2,4]); % move singleton-dimension-2 to dimension-3; transpose image. Size of result: [nColNew*nChannel, nRow, 1]
im = reshape(im, nColNew*nChannel*nStack, redSize, nRowNew); % packets of width redSize, as columns next to each other
im = sum(im,2); % sum in all rows. Size of result: [nColNew*nChannel, 1, nRowNew]
im = permute(im, [3,1,2,4]); % move singleton-dimension-2 to dimension-3; transpose image back. Size of result: [nRowNew, nColNew*nChannel, 1]
im = reshape(im, nRowNew, nColNew, nChannel, nStack); % putting all channels (rgb) back behind each other in the third dimension
im = uint8(sqrt(im./redSize^2)); % mean; re-normalize brightness: "scale linear to the human eye"; back in uint8
 
function drawCirclePatch(x,y,r,fillC,linW,linC)
%    x and y:      coordinates of the center
%    r:            is the radius of the circle
%    fillC:        color of filling [r g b]
%    linW:         LineWidth
%    linC:         LineColor
%    startOffset:  start rotation (scalar)[rad]
angleOffPoints = linspace(0,2.001*pi,50);
xc = x + r*cos(angleOffPoints);
yc = y + r*sin(angleOffPoints);
patch(xc,yc,fillC,'EdgeColor',linC,'LineWidth',linW);
 
function map = createImMap(imRGB,nCol,startMap)
% createImMap creates a color-map including predefined colors.
% "rgb2ind" creates a map but there is no option to predefine some colors,
%         and it does not handle stacked images.
% Input:
%   imRGB:     image, [imRows x imColumns x 3(RGB) x nStack] (unit8)
%   nCol:      total number of colors the map should have, [integer]
%   startMap:  predefined colors; colormap format, [p x 3] (double)
 
imRGB = permute(imRGB,[1 2 4 3]); % step1; make unified column-image (handling possible nStack)
imRGBcolumn = reshape(imRGB,[],1,3,1); % step2; make unified column-image
 
fullMap = double(permute(imRGBcolumn,[1 3 2]))./255; % "column image" to color map
[fullMap,~,imMapColumn] = unique(fullMap,'rows'); % find all unique colors; create indexed colormap-image
% "cmunique" could be used but is buggy and inconvenient because the output changes between "uint8" and "double"
 
nColFul = size(fullMap,1);
nColStart = size(startMap,1);
disp(['Number of colors: ' num2str(nColFul) ' (including ' num2str(nColStart) ' self defined)']);
 
if nCol<=nColStart;  error('Not enough colors');        end
if nCol>nColFul;   warning('More colors than needed');  end
 
isPreDefCol = false(size(imMapColumn)); % init
 
for iCol = 1:nColStart
    diff = sum(abs(fullMap-repmat(startMap(iCol,:),nColFul,1)),2); % difference between a predefined and all colors
    [mDiff,index] = min(diff); % find matching (or most similar) color
    if mDiff>0.05 % color handling is not precise
        warning(['Predefined color ' num2str(iCol) ' does not appear in image'])
        continue
    end
    isThisPreDefCol = imMapColumn==index; % find all pixel with predefined color
    disp([num2str(sum(isThisPreDefCol(:))) ' pixel have predefined color ' num2str(iCol)]);
    isPreDefCol = or(isPreDefCol,isThisPreDefCol); % combine with overall list
end
[~,mapAdditional] = rgb2ind(imRGBcolumn(~isPreDefCol,:,:),nCol-nColStart,'nodither'); % create map of remaining colors
map = [startMap;mapAdditional];