function Dragon_Curve_adding_corners_gif()
% Programmed in Matlab R2017a.
% Animation of an Dragon Curve unfolding.
% Several versions of the unfolding curve as gif files are created.
%    left 'L' and right 'R' with and without a trails.
%
% 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

greenCol = round([0.1  0.7  0.1]*255)/255; % green
redCol   = round([1    0    0  ]*255)/255; % red

greenFadeWhiteCol = round([0.5  0.9  0.5]*255)/255; % green
redFadeWhiteCol   = round([1    0.5  0.5]*255)/255; % red

greenFadeBlackCol = round([0.1  0.35 0.1]*255)/255; % green
redFadeBlackCol   = round([0.5  0    0  ]*255)/255; % red

% stop acc rot deacc 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)
displacementList = cumsum(speed)/sum(speed); % create position, normalize, 

xLimits = [-0.36,1.19];
yLimits = [-0.36,0.69];

xRange = xLimits(2)-xLimits(1);
yRange = yLimits(2)-yLimits(1);

figHandle = figure(6765301); clf;
set(figHandle,'Units'  ,'pixel');
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)
drawnow;

p1 = 50;  % determined by trial an error
p2 = 300; % determined by trial an error

for mode = {'minimalist','trails'} %
    iFrame    = 0;    % init
    
    x = [0 1]; % initial Line
    y = [0 0]; % initial Line
    
    switch mode{1}
        case {'minimalist'}

            nFramesMove = 1+ (length(displacementList)-1)*maxOrder; % inital Frame + (rotation - redundant frame)*order
            nFramesFadeOut = 5; % fade to white
            nFramesFadeIn = 4; % from white to first image (white not included)
            nFrames= nFramesMove + nFramesFadeOut + nFramesFadeIn;
            
            liWidth = [ones(1,p1)*26 linspace(26,8,p2-p1) linspace(8,2,nFrames-p2)];
            
            displacementIndexList = 1:numel(displacementList); % index "1" (only used in orderr=1) otherwise there would be redundant frames
            % List will be: 1 2 3 ...  end 2 3 ...  end 2 3 ...  end  % 
            % "1" and "end" are full dragon curves; all other numbers are moving
            
        case {'trails'}

            nFramesMove = 1+ (length(displacementList)+2)*maxOrder; % inital Frame + (rotation + intermediate state of the curve)*order
            nFramesFadeOut = 5; % fade to white
            nFramesFadeIn = 4; % from white to first image (white not included)
            nFrames= nFramesMove + nFramesFadeOut + nFramesFadeIn;
            liWidth = [ones(1,p1)*26 linspace(26,8,p2-p1) linspace(8,2,nFrames-p2)]; % LineWidth defined
            displacementIndexList = 0:numel(displacementList)+2; % index "0" (only used in orderr=1) and "+1" for the pure dragon curve without trails
            % List will be: 0 1 2 3 ...  end end+1 1 2 3 ...  end end+1 2 3 ...  end end+1
            % "0" and "end+1" are full dragon curves;
            % "1" and "end" are full dragon curves; but with colorful path plots visible
            % all other numbers are moving
    end
    
    delayTime = ones(1,nFrames)*1/25; % delay per frame

    screenSize = get(groot,'Screensize')-[0 0 5 20]; % [1 1 width height] (minus tolerance for figure borders)
    imageAspectRatio = xRange/yRange;
    MegaPixelTarget = 100*10^6; % Category:Animated GIF files exceeding the 100 MP limit
    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
    scaleReduction = min(...% repeat as often as possible for nice antialiasing
        floor(screenSize(4)/ySize), floor(screenSize(3)/xSize)); 
    if scaleReduction == 0;   error('"MegaPixelTarget" not possible; use smaller target or bigger monitor'); end % check
    
    liWidth = liWidth/sqrt(1000)*sqrt(xSize*scaleReduction); % re-scale if image is not original design size
    
    set(figHandle, 'Position',[1 1 xSize*scaleReduction ySize*scaleReduction]); % big start image for antialiasing later [x y width height]
    set(axesHandle,'Position',[0 0 1 1]); % stretch axis as big as figure, [x y width height]
    xlim(xLimits); ylim(yLimits);
    reducedRGBimageRight  = uint8(ones(ySize,xSize,3,nFrames)); % allocate
    reducedRGBimageLeft = reducedRGBimageRight; % allocate
    for orderr = 1:maxOrder
        
        nPoinitsOld = length(x);
        nPoinitsAdd = nPoinitsOld-1;
        
        addX = zeros(1,nPoinitsAdd); % allocate
        addY = zeros(1,nPoinitsAdd); % allocate
        
        xMean = mean([x(1:end-1); x(2:end)],1); % center on each line segment
        yMean = mean([y(1:end-1); y(2:end)],1); % center on each line segment
        
        xDiff = diff(x); % for scaling to current grid size
        yDiff = diff(y); % for scaling to current grid size
        
        nPoinitsNew = nPoinitsOld+nPoinitsAdd;
        
        oldPos = 1:2:nPoinitsNew;   % index in newX of already existing coners
        addPos = 2:2:nPoinitsNew-1; % index in newX of newly created coners
        
        newX = zeros(1,nPoinitsNew); % allocate
        newY = zeros(1,nPoinitsNew); % allocate
        
        newX(oldPos) = x; % known, non moving base points
        newY(oldPos) = y; % known, non moving base points

        for iDispl = displacementIndexList %
            iFrame = iFrame+1; % next frame
            
            switch mode{1}
                case {'minimalist'}
                    plotTrails = false; %
                    displacement = displacementList(iDispl)*0.5;
                    if iDispl == 1
                        if orderr == 1 % only the very first frame of the animation;
                            delayTime(iFrame) = 0.8; % keep completed order frames longer
                        else % orderr ~= 1 % iDispl=end,orderr=x and iDispl=1,orderr=x+1 are the same image
                            iFrame = iFrame-1; % revert
                            continue % skip
                        end
                    elseif iDispl == displacementIndexList(end)
                        if orderr == maxOrder % only the very last frame of the animation;
                            delayTime(iFrame) = 2.5; % keep max curve on screen for a while
                        else % orderr ~= 1 % iDispl=end,orderr=x and iDispl=1,orderr=x+1 are the same image
                            delayTime(iFrame) = 0.5; % keep completed order frames longer
                        end
                    end
                    
                case {'trails'}
                    % List will be: 0 1 2 3 ...  end end+1 1 2 3 ...  end end+1 2 3 ...  end end+1
                    % "0" and "end+1" are full dragon curves;
                    % "1" and "end" are full dragon curves; but with colorful path plots visible
                    % all other numbers are moving
                    plotTrails = 1; % default
                    if iDispl == 0
                        if orderr == 1 % only the very first frame of the animation
                            displacement = 0;
                            plotTrails = 0; %
                            delayTime(iFrame) = 0.8; % keep completed order frames longer
                        else % orderr ~= 1 % only the very first frame of the animation
                            iFrame = iFrame-1; % revert
                            continue % skip
                        end
                    elseif iDispl == 1
                        delayTime(iFrame) = 0.2; % movement starts, keep frames longer
                        displacement = displacementList(iDispl)*0.5;
                    elseif iDispl == numel(displacementList)
                        delayTime(iFrame) = 0.2; % movement ends, keep frames longer
                        displacement = displacementList(iDispl)*0.5;
                    elseif iDispl == numel(displacementList)+1 % at the end of each modification fade out trails
                        plotTrails   = 0.5; % fade out trails
                        displacement = 0.5;
                    elseif iDispl == numel(displacementList)+2 % at the end of each modification plot pure dragon curve
                        displacement = 0.5;
                        plotTrails   = 0; %
                        if orderr == maxOrder % only the very last frame of the animation;
                            delayTime(iFrame) = 2.5; % keep max curve on screen for a while
                        else % orderr ~= 1 % iDispl=end,orderr=x and iDispl=1,orderr=x+1 are the same image
                            delayTime(iFrame) = 0.3; % keep completed order frames longer
                        end
                    else % normal Frame
                        displacement = displacementList(iDispl)*0.5;
                    end
            end

            addX(1:2:end) = xMean(1:2:end) - displacement*yDiff(1:2:end); % 1st 3rd ... anti-clockwise
            addY(1:2:end) = yMean(1:2:end) + displacement*xDiff(1:2:end); % 1st 3rd ... anti-clockwise
            
            addX(2:2:end) = xMean(2:2:end) + displacement*yDiff(2:2:end); % 2nd 4th ... clockwise
            addY(2:2:end) = yMean(2:2:end) - displacement*xDiff(2:2:end); % 2nd 4th ... clockwise
            
            newX(addPos) = addX; % moving points
            newY(addPos) = addY; % moving points
            
            cla(axesHandle) % fresh frame
            plot(0, 0,'.w'); % dummy, it is needed trust me

            trailsHand = []; % default init
            if plotTrails>0 % not during the "pure dragon curve"-frames  marked with pseudo iDispl "0" and "-1"
                
                trailsXroot = x(1:end-1)+displacement*xDiff;
                trailsYroot = y(1:end-1)+displacement*yDiff;
                if plotTrails == 1
                    colTrailsOdd = redCol;
                    colTrailsEven = greenCol;
                else % plotTrails == 0.5
                    colTrailsOdd = redFadeWhiteCol; % fading out
                    colTrailsEven = greenFadeWhiteCol; % fading out
                    
                end

                signFlip = ones(size(trailsXroot));
                signFlip(2:2:end) = -1; % every second diffent direction
                
                trailsHand = plot(...
                    [trailsXroot;  xMean;  xMean-displacement*yDiff.*signFlip],...
                    [trailsYroot;  yMean;  yMean+displacement*xDiff.*signFlip],...
                    '.-','Color',colTrailsOdd, 'LineWidth',liWidth(iFrame)*0.7, 'MarkerSize',0.7*3.4*liWidth(iFrame));
                
               set(trailsHand(2:2:end),'Color',colTrailsEven); % every second diffent color
            end

            plot(newX,newY,'k.-','LineWidth',liWidth(iFrame),'MarkerSize',3.4*liWidth(iFrame)); % dragon curve
            
            
            dotHandOdd  = []; % default init
            dotHandEven = []; % default init
            if plotTrails>0 % not during the "pure dragon curve"-frames  marked with pseudo iDispl "0" and "-1"
                if plotTrails == 1
                    colDotsOdd  = redCol;
                    colDotsEven = greenCol;
                else % plotTrails == 0.5
                    colDotsOdd  = redFadeBlackCol; % fading out
                    colDotsEven = greenFadeBlackCol; % fading out
                end
                dotHandOdd  = plot(addX(1:2:end),addY(1:2:end),'.','Color',colDotsOdd, 'MarkerSize',liWidth(iFrame)*2); % anti-clockwise displaced points
                dotHandEven = plot(addX(2:2:end),addY(2:2:end),'.','Color',colDotsEven,'MarkerSize',liWidth(iFrame)*2); %      clockwise displaced points
            end
            
            if displacement < 0.25
                text_s = num2str(orderr-1); % order string
                co = interp1([-10 0.1 0.15 10],[0 0 1 1],displacement);
            else
                text_s = num2str(orderr);   % order string
                co = interp1([-10 0.35 0.4 10],[1 1 0 0],displacement);
            end
            text(0.025*xSize*scaleReduction,0.03*ySize*scaleReduction,text_s,... % add order text
                'FontName','Helvetica Narrow','Color',co*[1 1 1],...
                'FontUnits','pixels','Units','pixel',...
                'FontSize',0.09*xSize*scaleReduction,'FontWeight','bold',...
                'HorizontalAlignment','left','VerticalAlignment','baseline');
            
            xlim(xLimits); ylim(yLimits);
            
            %% save animation
            drawnow;
            f = getframe(figHandle);
            reducedRGBimageRight(:,:,:,iFrame) = imReduceSize(f.cdata,scaleReduction); % the size reduction: adds antialiasing

            %% flip top bottom (and red / green)
            objectList = findobj(axesHandle,'Type','line'); %
            valuesY = get(objectList,'YData');
            for i = 1:length(objectList)
                set(objectList(i),'YData',-valuesY{i})
            end
            if ~isempty(dotHandOdd)
                set(dotHandOdd, 'Color',colDotsEven); % flip colors
                set(trailsHand(1:2:end),'Color',colTrailsEven); % every second diffent color
            end
            if ~isempty(dotHandEven)
                set(dotHandEven,'Color',colDotsOdd);  % flip colors
                set(trailsHand(2:2:end),'Color',colTrailsOdd); % every second diffent color
            end
            ylim(fliplr(-yLimits));
            
            drawnow;
            f = getframe(figHandle);
            reducedRGBimageLeft(:,:,:,iFrame) = imReduceSize(f.cdata,scaleReduction); % the size reduction: adds antialiasing
        end
        
        x = newX; % update base structure
        y = newY; % update base structure
    end
    
    whiteImage = ones([ySize,xSize,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
            reducedRGBimageLeft (:,:,:,nFramesMove + iFrame) = whiteImage-(whiteImage-reducedRGBimageLeft (:,:,:,nFramesMove))*brightnes(iFrame); % brightness reduction
            reducedRGBimageRight(:,:,:,nFramesMove + iFrame) = whiteImage-(whiteImage-reducedRGBimageRight(:,:,:,nFramesMove))*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
            reducedRGBimageLeft (:,:,:,nFramesMove+nFramesFadeOut + iFrame) = whiteImage-(whiteImage-reducedRGBimageLeft (:,:,:,1))*brightnes(iFrame); % brightness reduction
            reducedRGBimageRight(:,:,:,nFramesMove+nFramesFadeOut + iFrame) = whiteImage-(whiteImage-reducedRGBimageRight(:,:,:,1))*brightnes(iFrame); % brightness reduction
        end
    end
    
    switch mode{1}
        case {'minimalist'}
            map = gray(16);
        case {'trails'}
            map = createImMap(reducedRGBimageLeft,32,[0 0 0;1 1 1;greenCol;redCol]); % colormap
    end

    for iFrame = 1:nFrames
        imLeft  = rgb2ind(reducedRGBimageLeft (:,:,:,iFrame),map,'nodither'); % rgb to colormap image
        imRight = rgb2ind(reducedRGBimageRight(:,:,:,iFrame),map,'nodither'); % rgb to colormap image
        if iFrame == 1
            imwrite(imLeft, map,fullfile(pathstr, [fname '_' mode{1} '_rectangular_numbered_L.gif']),  'LoopCount',Inf,'DelayTime',delayTime(iFrame)); % individual timings
            imwrite(imRight,map,fullfile(pathstr, [fname '_' mode{1} '_rectangular_numbered_R.gif']), 'LoopCount',Inf,'DelayTime',delayTime(iFrame)); % individual timings
        else
            imwrite(imLeft, map,fullfile(pathstr, [fname '_' mode{1} '_rectangular_numbered_L.gif']),  'WriteMode','append','DelayTime',delayTime(iFrame)); % individual timings
            imwrite(imRight,map,fullfile(pathstr, [fname '_' mode{1} '_rectangular_numbered_R.gif']), 'WriteMode','append','DelayTime',delayTime(iFrame)); % individual timings
        end
    end
    disp([fname '-' mode{1} '_rectangular_numbered_xx.gif  has ' num2str(numel(reducedRGBimageLeft)/3/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 100 MP limit
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 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];