function [] = Dragon_Curve_unfolding_gif()
% Programmed in Matlab R2017a.
% Animation of an Dragon Curve unfolding.
% Before anything is drawn the algorithm makes a dry run ('xyLim')
%    to determine proper axis limits to keep the 
%    animation nicely in frame. 
% Then several versions of the unfolding curve as gif files are created.
%    'rectangular' and 'chamfered' with and without a final zoom inwards.
%    all with and without the order number in the lower left.
%
% 2019-05-30 Jahobr - Licensing: CC0 1.0 Universal Public Domain Dedication
 
[pathstr,fname] = fileparts(which(mfilename)); % save files under the same name and at file location
fname = strrep(fname,'_gif','');

maxOrder = 15; % final order of fractal
 
figHandle = figure(345275); clf;
set(figHandle,'Color'  ,'white'); % white background
set(figHandle,'Units'  ,'pixel');
set(figHandle,'MenuBar','none',  'ToolBar','none'); % free real estate for a maximally large image

axesHandle = axes; hold(axesHandle,'on');
axis equal
axis off % invisible axes (no ticks)


% stop acceleration rot de-acceleration stop
accFrames = 8; % frames for acceleration (first frame will be 0 last at full speed, so practicall it is accFrames-2)
speed = [linspace(0,1,accFrames) ones(1,21) linspace(1,0,accFrames)];
speed = speed(1:end-1); % last speed is 0, this does nothing in cumsum; (compensated by +1 frames in center)
angleList = cumsum(speed)/sum(speed) * pi/2; % create position, normalize, scale
 
nFramesZoomOut = 1+ (length(angleList)+1)*maxOrder; % inital Frame + (rotation + intermediate state of the curve)*order
xMinR = -0.1*ones(1,nFramesZoomOut); % initial range
xMaxR = +2.1*ones(1,nFramesZoomOut); % initial range
yMinR = -1.1*ones(1,nFramesZoomOut); % initial range
yMaxR = +1.1*ones(1,nFramesZoomOut); % initial range
 
% LineWidth defined for 900x900px image; for gifs
liWidth = interp1([ 0   0.09  0.54  1]*nFramesZoomOut... during the run
    ,             [24  24     8     2]... defined Linewidth
    ,             1:nFramesZoomOut); % LineWidth sampled for all frames

h_xy = []; % init

for mode = {'xyLim','zoom','rectangular','chamfered'} %
    
    x_y = [0 0;0 1];  % inital line [x1 x2... ; y1 y2...]
    x_y_plot  = x_y;  % inital line
    iFrame    = 0;    % init
    
    switch mode{1}
        case {'xyLim','rectangular','chamfered'}
            nFramesZoomIn = 0; % no zoom back in
            nFramesFadeOut = 5; % fade to white 
            nFramesFadeIn = 4; % from white to first image (white not included)
        case 'zoom'
            nFramesZoomIn = 800-nFramesZoomOut; % zoom back in
            nFramesFadeOut = 0; % done automatically
            nFramesFadeIn = 4;
    end
    nFrames = nFramesZoomOut+nFramesZoomIn+nFramesFadeOut+nFramesFadeIn;
    
    delayTime = ones(1,nFrames)*1/25; % delay per frame
    delayTime(1:2) = 0.4; % keep first two frames longer
    delayTime(nFramesZoomOut) = 2.5; % keep max curve on screen for a while
    
    MegaPixelTarget = 100*10^6; % Category:Animated GIF files exceeding the 100 MP limit
    xySize = floor(sqrt(MegaPixelTarget/nFrames)); %  gif size in pixel
    screenSize = get(groot,'Screensize')-[0 0 5 20]; % [1 1 width height] (minus tolerance for figure borders)
    scaleReduction = min(...% reduction for nice antialiasing (for big numvers you need a 4K monitor or a virtural combination of several monitors using "Nvidia Surround" to fit the figure)
        floor(screenSize(4)/xySize), floor(screenSize(3)/xySize)); 
    if scaleReduction == 0;   error('"MegaPixelTarget" not possible; use smaller target or bigger monitor'); end % check
    
    xyFigS = xySize*scaleReduction; % full Figure size; limited by your monitor
    liWidth = liWidth/sqrt(900)*sqrt(xyFigS); % re-scale if image is not 900x900px
    figPosition = [1 1 xyFigS xyFigS]; % big start image for antialiasing later [x y width height]
    set(figHandle, 'Position',figPosition); % big start image; reduction allows for subpixel LineWidth
    if ~all(get(figHandle, 'Position')==figPosition)
        error('figure Position could not be set')
    end
    reducedRGBimage1 = ones([xySize,xySize,3,nFrames],'uint8'); % allocate
    reducedRGBimage2 = ones([xySize,xySize,3,nFrames],'uint8'); % allocate
 
    for orderr = 1:maxOrder % create zoom out; generating the dragon curve
        
        rot_x_y = x_y(:,end); % point of rotation
        
        x_y_new = fliplr(x_y_plot); 
        x_y_new = x_y_new - rot_x_y*ones(1,size(x_y_new,2)); % shift to zero at point of rotation
        
        for currA = [-2*pi angleList 2.5*pi]
            
            if and(currA == -2*pi, orderr ~= 1) % only the very first frame of the animation
                continue
            end
            
            iFrame = iFrame+1; % next frame
            [x_y_newRot] = rotateCordiantes(x_y_new,currA); % rotate
            x_y_newRot = x_y_newRot + rot_x_y*ones(1,size(x_y_newRot,2)); % shift back
            
            if currA == 2.5*pi % end of every rotation
                
                x_y_new = fliplr(x_y); % use unmodified curve
                x_y_new = x_y_new - rot_x_y*ones(1,size(x_y_new,2)); % shift to zero at point of rotation
                [x_y_newRot] = rotateCordiantes(x_y_new,currA); % rotate
                x_y_newRot = x_y_newRot + rot_x_y*ones(1,size(x_y_newRot,2)); % shift back
                x_y = [x_y  x_y_newRot(:,2:end)]; % extend curve
                x_y = round(x_y); % eliminate any rounding errors by the rotation (not really a problem)
                
                switch mode{1}
                    case {'xyLim','rectangular'}
                        x_y_plot = x_y; % no modification
                    case {'chamfered','zoom'}
                        np = size(x_y,2); 
                        x_y_plot = interp1(1:np,x_y',1:0.25:np)'; % interpolate point at: 1 1.25 1.5 1.75 2 ..
                        x_y_plot = x_y_plot(:,[1 (2:2:end) end]); % use only: first 1.25 1.75 1.25 1.75 ... last
                end
            end
            
            switch mode{1}
                
                case 'xyLim'
                    xMinR(iFrame:end) = min([xMinR(iFrame) x_y(1,:) x_y_newRot(1,:)]); % record extreme range
                    xMaxR(iFrame:end) = max([xMaxR(iFrame) x_y(1,:) x_y_newRot(1,:)]); % record extreme range
                    yMinR(iFrame:end) = min([yMinR(iFrame) x_y(2,:) x_y_newRot(2,:)]); % record extreme range
                    yMaxR(iFrame:end) = max([yMaxR(iFrame) x_y(2,:) x_y_newRot(2,:)]); % record extreme range
                    
                case {'rectangular','chamfered','zoom'}
                    cla(axesHandle)
                    setXYlim(axesHandle, [xMinR(iFrame) xMinR(iFrame)+maxR(iFrame)], [yMinR(iFrame) yMinR(iFrame)+maxR(iFrame)]) % use pre-created zoom
                    
                    
                    if currA == -2*pi % very first frame
                        h_xy = plot(x_y_plot(1,:), x_y_plot(2,:), 'k.-','LineWidth',liWidth(iFrame),'MarkerSize',3.4*liWidth(iFrame)); % just line, no rotation point marker
                    elseif currA == 2.5*pi % end of every rotation
                        h_xy = plot(x_y_plot(1,:), x_y_plot(2,:), 'k.-','LineWidth',liWidth(iFrame),'MarkerSize',3.4*liWidth(iFrame)); % just line, no rotation point marker
                    else
                        h_xy = plot(x_y_plot(1,:), x_y_plot(2,:), 'b.-','LineWidth',liWidth(iFrame),'MarkerSize',3.4*liWidth(iFrame)); % non rotating original line
                        plot(x_y_newRot(1,:), x_y_newRot(2,:), 'k.-','LineWidth',liWidth(iFrame),'MarkerSize',3.4*liWidth(iFrame)); %  rotating section
                        plot(rot_x_y(1),rot_x_y(2),'.r','MarkerSize',liWidth(iFrame)*6); % rotation point marker
                    end
                    
                    drawnow; %
                    f1 = getframe(figHandle);
                    reducedRGBimage1(:,:,:,iFrame) = imReduceSize(f1.cdata,scaleReduction); % the size reduction: adds antialiasing
                    
                    if currA < pi*0.25
                        text_s = num2str(orderr-1); % order string
                        co = interp1([-100 pi*0.1 pi*0.15 100],[0 0 1 1],currA);
                    else
                        text_s = num2str(orderr);   % order string
                        co = interp1([-100 pi*0.35 pi*0.4 100],[1 1 0 0],currA);
                    end
                    tHand = text(0.07*xyFigS,0.066*xyFigS,text_s,... % add order text
                        'FontName','Helvetica Narrow','Color',co*[1 1 1],...
                        'FontUnits','pixels','Units','pixel',...
                        'FontSize',0.12*xyFigS,'FontWeight','bold',...
                        'HorizontalAlignment','left','VerticalAlignment','baseline');

                    drawnow; %
                    f2 = getframe(figHandle);
                    reducedRGBimage2(:,:,:,iFrame) = imReduceSize(f2.cdata,scaleReduction); % the size reduction: adds antialiasing
            end
        end
    end
 
    if strcmpi(mode{1},'zoom')
        aFrames = 20; % frames for acceleration
        cFrames = 20; % frames for constant speed
        dFrames = 70; % frames for deceleration
        RangeChangeSpeed = [linspace(0,1,aFrames) ones(1,cFrames) linspace(1,0,dFrames) zeros(1,nFramesZoomIn-aFrames-cFrames-dFrames)];
        Range_Change  = cumsum(RangeChangeSpeed)/sum(RangeChangeSpeed); % create position, normalize, scale
        maxR(end); % start Range
        endRange = 15; % end Range
        rangeList = maxR(end)-Range_Change*(maxR(end)-endRange);
        
        liWidth(end); % start Width with final value of building it up
        endLiSizeGif = 10/sqrt(900)*sqrt(xyFigS); % end Width % re-scale if image is not 900x900px; 
        liWidthZoom = liWidth(end)-(Range_Change.^6)*(liWidth(end)-endLiSizeGif); % create new linewith range
        
        aFrames = 20; % frames for acceleration
        cFrames = 10; % frames for constant speed
        dFrames = 80; % frames for deceleration
        finMovSpeed = 0.05; % by trial and error
        PosChangeSpeed = [linspace(0,1,aFrames) ones(1,cFrames) linspace(1,finMovSpeed,dFrames) ones(1,nFramesZoomIn-aFrames-cFrames-dFrames)*finMovSpeed];
        Pos_Change  = cumsum(PosChangeSpeed)/sum(PosChangeSpeed); % create position, normalize, scale
        
        xPos = xMinR(end)+ Pos_Change*(-100-xMinR(end)); % stretch from generic shape to start and end point
        yPos = yMinR(end)+ Pos_Change*(50-yMinR(end));   % stretch from generic shape to start and end point
        co = [zeros(1,nFramesZoomIn-5) linspace(0,1,5)]; % color Black to the end, then fade to white
        
        figure(43542); clf; hold on; grid on
        title('limits zoom in')
        plot(1:nFramesZoomIn,xPos,'r', 1:nFramesZoomIn,xPos+rangeList,'r', 1:nFramesZoomIn,yPos,'b', 1:nFramesZoomIn,yPos+rangeList,'b');
        legend({'x limits','x limits','y limits','y limits'})
        
        figure(figHandle); % back to main figure
        for iFrame = 1:nFramesZoomIn %
            
            set(h_xy,'LineWidth',liWidthZoom(iFrame),'MarkerSize',3.4*liWidthZoom(iFrame),'Color',co(iFrame)*[1 1 1]); % set linewidth
            set(tHand,'Color',[1 1 1]); % hide text
            setXYlim(axesHandle, [xPos(iFrame) xPos(iFrame)+rangeList(iFrame)], [yPos(iFrame) yPos(iFrame)+rangeList(iFrame)]) % use pre-created zoom
            drawnow
            f1 = getframe(figHandle);
            reducedRGBimage1(:,:,:,nFramesZoomOut + iFrame) = imReduceSize(f1.cdata,scaleReduction); % the size reduction: adds antialiasing
            
            if iFrame<=23 % at the start of the zoom in -> fade out the number
                col = interp1([0 17 23],[0 0 1],iFrame);
                set(tHand,'Color',col*[1 1 1]); %
            else
                set(tHand,'String',''); % just get rid of the number from now on
            end
            drawnow
            f2 = getframe(figHandle);
            reducedRGBimage2(:,:,:,nFramesZoomOut + iFrame) = imReduceSize(f2.cdata,scaleReduction); % the size reduction: adds antialiasing
        end
    end
    
    switch mode{1} %
        case {'rectangular','chamfered','zoom'}
            % The gif is looping automatically. To avoid a sudden change
            % from the last frame to the next first frame some smooth
            % transition frames are added at the end
            whiteImage = ones([xySize,xySize,3,1],'uint8')*255; % allocate
            if nFramesFadeOut>0
                brightnes = linspace(1,0,nFramesFadeOut+1);
                brightnes = brightnes(2:end); % first "1" not needed
                for iFrame = 1:nFramesFadeOut % fade to white
                    reducedRGBimage1(:,:,:,nFramesZoomOut+nFramesZoomIn + iFrame) = whiteImage-(whiteImage-reducedRGBimage1(:,:,:,nFramesZoomOut+nFramesZoomIn))*brightnes(iFrame); % brightness reduction
                    reducedRGBimage2(:,:,:,nFramesZoomOut+nFramesZoomIn + iFrame) = whiteImage-(whiteImage-reducedRGBimage2(:,:,:,nFramesZoomOut+nFramesZoomIn))*brightnes(iFrame); % brightness reduction
                end
            end
            if nFramesFadeIn>0
                brightnes = linspace(0,1,nFramesFadeIn+2);
                brightnes = brightnes(2:end-1); % first "0" and last "1" not needed
                for iFrame = 1:nFramesFadeIn % from white to frame1
                    reducedRGBimage1(:,:,:,nFramesZoomOut+nFramesZoomIn+nFramesFadeOut + iFrame) = whiteImage-(whiteImage-reducedRGBimage1(:,:,:,1))*brightnes(iFrame); % brightness reduction
                    reducedRGBimage2(:,:,:,nFramesZoomOut+nFramesZoomIn+nFramesFadeOut + iFrame) = whiteImage-(whiteImage-reducedRGBimage2(:,:,:,1))*brightnes(iFrame); % brightness reduction
                end
            end
    end
    
    switch mode{1}
        case 'xyLim'
            % getting a smooth zoom is the most difficult part
            % the image is defined as square, but the full height and width
            % is not always required. Extra space is assigned at the side,
            % which requires it next.
            
            constPadding = 0.015; % 1.5% padding around the actually used area
            
            paddingFactorLowLeft  = interp1([1 [0.15 0.2 0.8 0.95 1]*nFramesZoomOut], [1 1 1.04 1.04 1 1]+constPadding, 1:nFramesZoomOut);
            xMinR = xMinR.*paddingFactorLowLeft; % start and end fit nicely with constPadding but during movement some padding is required around the dragon curve
            yMinR = yMinR.*paddingFactorLowLeft; % start and end fit nicely with constPadding but during movement some padding is required around the dragon curve

            paddingFactorTopRight = interp1([1 [0.15 0.2 0.8 0.95 1]*nFramesZoomOut], [1 1 1.02 1.02 1 1]+constPadding, 1:nFramesZoomOut);
            xMaxR = xMaxR.*paddingFactorTopRight; % start and end fit nicely with constPadding but during movement some padding is required around the dragon curve
            yMaxR = yMaxR.*paddingFactorTopRight; % start and end fit nicely with constPadding but during movement some padding is required around the dragon curve
            
            xRange = xMaxR-xMinR;  yRange = yMaxR-yMinR;

            maxR = max(xRange,yRange);
            assignableX = (maxR-xRange); % assignable space in x-direction
            assignableY = (maxR-yRange); % assignable space in y-direction
            
            figure(43543); clf; hold on; grid on
            title('required limits')
            plot([1:nFramesZoomOut nFramesZoomOut:-1:1],[xMinR fliplr(xMaxR)],'r', [1:nFramesZoomOut nFramesZoomOut:-1:1],[yMinR fliplr(yMaxR)],'b');
            plot(assignableX,'m'); plot(assignableY,'k');
            legend({'x limits','y limits','assignableX','assignableY'})
 
            for iFrame = 1:nFramesZoomOut
                xRange = xMaxR-xMinR;  yRange = yMaxR-yMinR;
                maxR = max(xRange,yRange);
                assignableX = (maxR-xRange); % assignable space in x-direction
                assignableY = (maxR-yRange); % assignable space in y-direction
                
                if assignableX(iFrame)>0 % x has wiggle room
                    needAtRight = find(xMaxR >= xMaxR(iFrame)+ assignableX(iFrame),1,'first'); % when will space be required next?
                    needAtLeft  = find(xMinR <= xMinR(iFrame)- assignableX(iFrame),1,'first'); % when will space be required next?
                    
                    if and(isempty(needAtRight),isempty(needAtLeft))
                        nVal = xMinR(iFrame)-assignableX(iFrame)/2;
                        xMinR(iFrame:end) = min(xMinR(iFrame:end),nVal); % assign space
                        nVal = xMaxR(iFrame)+assignableX(iFrame)/2;
                        xMaxR(iFrame:end) = max(xMaxR(iFrame:end),nVal); % assign space
                    elseif or(isempty(needAtLeft), needAtLeft>needAtRight)
                        nVal = xMaxR(iFrame)+assignableX(iFrame);
                        xMaxR(iFrame:end) = max(xMaxR(iFrame:end),nVal); % assign space
                    elseif or(isempty(needAtRight), needAtLeft<needAtRight)
                        nVal = xMinR(iFrame)-assignableX(iFrame);
                        xMinR(iFrame:end) = min(xMinR(iFrame:end),nVal); % assign space
                    end
                end
                if assignableY(iFrame)>0 % y has wiggle room
                    needAtRight = find(yMaxR >= yMaxR(iFrame)+ assignableY(iFrame),1,'first'); % when will space be required next?
                    needAtLeft  = find(yMinR <= yMinR(iFrame)- assignableY(iFrame),1,'first'); % when will space be required next?
                    
                    if and(isempty(needAtRight),isempty(needAtLeft))
                        nVal = yMinR(iFrame)-assignableY(iFrame)/2;
                        yMinR(iFrame:end) = min(yMinR(iFrame:end),nVal); % assign space
                        nVal = yMaxR(iFrame)+assignableY(iFrame)/2;
                        yMaxR(iFrame:end) = max(yMaxR(iFrame:end),nVal); % assign space
                    elseif or(isempty(needAtRight), needAtLeft<needAtRight)
                        nVal = yMinR(iFrame)-assignableY(iFrame);
                        yMinR(iFrame:end) = min(yMinR(iFrame:end),nVal); % assign space
                    elseif or(isempty(needAtLeft), needAtLeft>needAtRight)
                        nVal = yMaxR(iFrame)+assignableY(iFrame);
                        yMaxR(iFrame:end) = max(yMaxR(iFrame:end),nVal); % assign space
                    end
                end
            end
            
            % the zomm [min max] covers now the neccessary area but is very "stuttering"
            figure(43544); clf; hold on; grid on
            title('limits with assigned space/ filtered')
            plot([1:nFramesZoomOut nFramesZoomOut:-1:1],[xMinR fliplr(xMaxR)],'r', [1:nFramesZoomOut nFramesZoomOut:-1:1],[yMinR fliplr(yMaxR)],'b');
            % create smooth zoom with moving average filter & Zero-phase distortion
            wind = 30; % window size
            xMinR = (filtfilt(ones(1,wind)/wind,1,(xMinR)));
            xMaxR = (filtfilt(ones(1,wind)/wind,1,(xMaxR)));
            yMinR = (filtfilt(ones(1,wind)/wind,1,(yMinR)));
            yMaxR = (filtfilt(ones(1,wind)/wind,1,(yMaxR)));
 
            % Filter form end to start (right to left). This creates space anticipatory rather then lagging.
            % As a result the dragon_curve does not reach out of the frame 
            wind = 24;  b=ones(1,wind)/wind;  a=1; % filter parameters
            xMinR = fliplr(filter(b,a,fliplr(xMinR),filtic(b,a,ones(1,wind)*xMinR(end),ones(1,wind)*xMinR(end))));
            xMaxR = fliplr(filter(b,a,fliplr(xMaxR),filtic(b,a,ones(1,wind)*xMaxR(end),ones(1,wind)*xMaxR(end))));
            yMinR = fliplr(filter(b,a,fliplr(yMinR),filtic(b,a,ones(1,wind)*yMinR(end),ones(1,wind)*yMinR(end))));
            yMaxR = fliplr(filter(b,a,fliplr(yMaxR),filtic(b,a,ones(1,wind)*yMaxR(end),ones(1,wind)*yMaxR(end))));
            xRange = xMaxR-xMinR;  yRange = yMaxR-yMinR;
            maxR = max(xRange,yRange);
            
            plot([1:nFramesZoomOut nFramesZoomOut:-1:1],[xMinR fliplr(xMaxR)],'r:', [1:nFramesZoomOut nFramesZoomOut:-1:1],[yMinR fliplr(yMaxR)],'b:');
            legend({'x limits','y limits','x filtered','y filtered'})
            figure(figHandle); % back to main figure
            
        otherwise %  {'rectangular','chamfered','zoom'}
            % create animation
            colormapImage = permute(reducedRGBimage1(:,:,:,:),[1 2 4 3]); % step1; make unified image
            colormapImage = reshape(colormapImage,[xySize,xySize*(nFrames) 3 1]); % step2; make unified image
            map = createImMap(colormapImage,32,[0 0 0;1 1 1;1 0 0;0 0 1]); % colormap

            for iFrame = 1:nFrames
                im1 = rgb2ind(reducedRGBimage1(:,:,:,iFrame),map,'nodither');
                im2 = rgb2ind(reducedRGBimage2(:,:,:,iFrame),map,'nodither');
                
                if iFrame == 1
                    imwrite(im1,map,fullfile(pathstr, [fname '_' mode{1} '_unnum.gif']),   'LoopCount',Inf,'DelayTime',delayTime(iFrame)); % individual timings
                    imwrite(im2,map,fullfile(pathstr, [fname '_' mode{1} '_numbered.gif']),'LoopCount',Inf,'DelayTime',delayTime(iFrame)); % individual timings
                else
                    imwrite(im1,map,fullfile(pathstr, [fname '_' mode{1} '_unnum.gif']),   'WriteMode','append','DelayTime',delayTime(iFrame)); % individual timings
                    imwrite(im2,map,fullfile(pathstr, [fname '_' mode{1} '_numbered.gif']),'WriteMode','append','DelayTime',delayTime(iFrame)); % individual timings
                end
            end
            disp([fname '.gif  has ' num2str(numel(reducedRGBimage1)/3/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 100 MP limit        
    end
end


function [xy] = rotateCordiantes(xy,anglee)
% [x1 x2 x3 ... ; y1 y2 y3 ...] coordinates to rotate
% anglee angle of rotation in [rad]
rotM = [cos(anglee) -sin(anglee); sin(anglee) cos(anglee)];
xy = rotM*xy;


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 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];


function setXYlim(axesHandle,xLimits,yLimits)
% set limits; practically the axis overhangs the figure all around, to
% hide rendering error at line-ends.
% Input:
%   axesHandle:        
%   xLimits, yLimits:  [min max]
overh = 0.05; % 5% overhang all around; 10% bigger in x and y
xlim([+xLimits(1)*(1+overh)-xLimits(2)*overh  -xLimits(1)*overh+xLimits(2)*(1+overh)])
ylim([+yLimits(1)*(1+overh)-yLimits(2)*overh  -yLimits(1)*overh+yLimits(2)*(1+overh)])
set(axesHandle,'Position',[-overh -overh  1+2*overh 1+2*overh]); % stretch axis as bigger as figure, [x y width height]
drawnow;