function Tusi_couple_ellipses()
% source code for Tusi_couple_ellipses
% produces a GIF and a SVG
%
% 2017-04-10 Jahobr (update 14.03.2021)

rBase = 1; % base radius
nFrames = 340;
angleSmall = linspace(0,2*pi,nFrames+1); % define gear position in frames
angleSmall = angleSmall(1:end-1); % remove last frame, it would be double

figHandle = figure(15674455);
whitebg(figHandle,[246 248 249]./250) % color backbground
clf
axesHandle = axes;
hold(axesHandle,'on')
set(figHandle, 'Units','pixel');
set(figHandle, 'Position',[1 1 1000 1000]); % big start image for antialiasing later [x y width height]
set(axesHandle,'Position',[-0.05 -0.05 1.1 1.1]); % stretch axis bigger as figure, easy way to get rid of ticks [x y width height]
xlim([-1.2 1.2]);
ylim([-1.2 1.2]);
axis equal; drawnow;
axis off % invisible axes (no ticks)
set(figHandle, 'Color',[1 1 1]); % white background

reducedRGBimage = uint8(ones(500,500,3,nFrames)); % allocate

for iFrame = 1:nFrames
    
    currentAngle = angleSmall(iFrame);
    
    cla(axesHandle) % fresh frame
    
    
    xs = cos(currentAngle)*0.5*rBase;
    ys = sin(currentAngle)*0.5*rBase;
    
    circle_patch(0,0,rBase,[1 1 1]) % outer circle backbground
    circle_patch(xs,ys,rBase/2,[0.95 0.95 0.95])  % inner circle backbground
    
    plot(rBase*cos(currentAngle),rBase*sin(currentAngle),'.m','MarkerSize',40) % touching point / Instant centre of rotation
    plot(0,0,'+','MarkerSize',15,'LineWidth',4,'Color',[0.5 0.5 0.5]); % center marker
    
    %% draw center circle 
    colD    = [0.5 0.5 0.5]; %
    circle(0,0,rBase/2,colD,4);
    
    
    %% draw inner ellipse
    colG = [0.1 0.7 0.1]; % green
    randiusA = 0.6;
    ellipse(0,0,randiusA,1-randiusA,colG);
    xpG = randiusA*cos(currentAngle);
    ypG =(rBase-randiusA)*sin(currentAngle);
    plot([xs xpG],[ys ypG],':','LineWidth',4,'Color',colG); % radial line

    
    %% draw medium ellipse
    colC    = [0  0.9 0.9]; % cyan    
    randiusA = 0.85;
    h(3) = ellipse(0,0,randiusA,rBase-randiusA,colC,pi/4);
    xpC = xs+sin(currentAngle)*(randiusA-rBase/2);
    ypC = ys+cos(currentAngle)*(randiusA-rBase/2);
    plot([xs xpC],[ys ypC],':','LineWidth',4,'Color',colC); % radial line
    
    %% draw outer ellipse
    colR     = [1 0 0]; % red
    randiusA = 1.35;
    h(1) = ellipse(0,0,randiusA,rBase-randiusA,colR,-pi/4);
    xpR = xs-sin(currentAngle)*(randiusA-rBase/2);
    ypR = ys-cos(currentAngle)*(randiusA-rBase/2);
    plot([xs xpR],[ys ypR],':','LineWidth',4,'Color',colR); % radial line
    
    %% draw verical Line (degenerated ellipse)
    plot([0 0],[rBase,-rBase],'-b','LineWidth',4) % vertical line
    plot([xs 0],[ys rBase*sin(currentAngle)],':b','LineWidth',4); % radial line

    %% final touches
    
    plot(xpG,ypG,'.','MarkerSize',40,'LineWidth',4,'Color',colG); % tracking point
    plot(xpC,ypC,'.','MarkerSize',40,'LineWidth',4,'Color',colC); % tracking point
    plot(xpR,ypR,'.','MarkerSize',40,'LineWidth',4,'Color',colR); % tracking point
    
    circle(0,0,rBase,[0 0 0],5) % outer circle
    circle(xs,ys,rBase/2,[0 0 0],5)  % inner circle
    plot([xs 0],[ys rBase*sin(currentAngle)],'.:b','MarkerSize',40,'LineWidth',4); % tracking point on top of circles
    plot(xs,ys,'.','MarkerSize',40,'LineWidth',4,'Color',colD); % tracking point
    
    

    %% save animation
    drawnow
    pause(0.01)
    saveName = 'Tusi_couple_ellipses';
    f = getframe(figHandle);
    reducedRGBimage(:,:,:,iFrame) = imReduceSize(f.cdata,2); % the size reduction: adds antialiasing
%     
    if iFrame == 72
        plot2svg([saveName '.svg'],figHandle) % by Juerg Schwizer, See http://www.zhinst.com/blogs/schwizer/
    end

end

map = createImMap(reducedRGBimage,128,[0 0 0; 1 1 1; 0.5 0.5 0.5; 0.95 0.95 0.95; colD; colG: colC; colR]); % colormap

im = uint8(ones(500,500,1,nFrames)); % allocate
for iFrame = 1:nFrames
    im(:,:,1,iFrame) = rgb2ind(reducedRGBimage(:,:,:,iFrame),map,'nodither');
end

imwrite(im,map,[saveName '.gif'],'DelayTime',1/50,'LoopCount',inf) % save gif1
disp([saveName '.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 100 MP limit

return
%%

function circle(x,y,r,col,linw)
% x coordinates of the center
% y coordinates of the center
% r is the radius of the circle
angleOffPoints = linspace(0,2*pi,300);
xc = x + r*cos(angleOffPoints);
yc = y + r*sin(angleOffPoints);
plot(xc,yc,'k','LineWidth',linw,'Color',col);

function circle_patch(x,y,r,col)
% x coordinates of the center
% y coordinates of the center
% r is the radius of the circle
angleOffPoints = linspace(0,2*pi,300);
xc = x + r*cos(angleOffPoints);
yc = y + r*sin(angleOffPoints);
patch(xc,yc,col,'EdgeColor','none') % 

function h = ellipse(x,y,a,b,col,theta)
% x coordinates of the center
% y coordinates of the center
% a radius1
% b radius2
if nargin <= 5
    theta = 0;
end
angleOffPoints = linspace(0,2*pi,300);
xe = x + a*cos(angleOffPoints);
ye = y + b*sin(angleOffPoints);

xe_rot =  xe*cos(theta) - ye*sin(theta);
ye_rot =  xe*sin(theta) + ye*cos(theta);
h = plot(xe_rot,ye_rot,'-','LineWidth',4,'Color',col);

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