File:Ellipses Astroid Tusi couple a.gif

function Ellipses_Astroid_Tusi_couple()
% source code for drawing Ellipses_Astroid_Tusi_couple
% produces a GIF and a SVG
%
% 2017-05-25 Jahobr (reworked 2018-06-10)

p = 1;
nFrames = 320;
xySize = 400; % size in pixel
scaleReduction = 2; % the size reduction: adds antialiasing

angleSmall = linspace(0,2*pi,nFrames+1); % define gear position in frames
angleSmall = angleSmall(1:end-1); % remove last frame, it would be double
angleSmall = circshift(angleSmall,round(nFrames*0.89)); % shift a view frames to get good visibility of all lines in frame 1
[pathstr,fname] = fileparts(which(mfilename)); % save files under the same name and at file location

figHandle = figure(15674455);
clf
axesHandle = axes;
hold(axesHandle,'on')
set(figHandle, 'Units','pixel');
set(figHandle, 'position',[1 1 xySize*scaleReduction xySize*scaleReduction]); % big start image for nicer 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]
set(figHandle,'GraphicsSmoothing','on') % requires at least version 2014b
xlim(p*[-1.15 1.15]);
ylim(p*[-1.15 1.15]);
axis equal; drawnow;

nE = 7; % number ellipses

colE = [... color of ellipses
    1.0  0.5  0  ;... % 1 orange
    0.8  0.8  0  ;... % 2 yellow
    0.6  1.0  0  ;... % 3 yellow green
    0    0.7  0  ;... % 4 green
    0    0.8  1.0;... % 5 cyan
    0.0  0    1.0;... % 6 blue
    0.4  0    1.0];...% 7 violet
colE   = round(colE.*255)./255; % round to uint8 position

fullCircCol    = round([0.0 0.0 0.0].*255)./255; % color of full circle
halfCircCol    = round([0.2 0.2 0.2].*255)./255; % color of 0.5  circle
quarterCircCol = round([0.4 0.4 0.4].*255)./255; % color of 0.25 circle
astroidCol     = round([1.0 0.0 0.0].*255)./255; % color astroid
touchCol       = round([1.0 0.0 1.0].*255)./255; % touching point

%% draw ellipses

ellOffset = linspace(0,p,nE);
for iE = 1:nE
    ellipse(0,0,p-ellOffset(iE),ellOffset(iE),colE(iE,:),5,'-');
end
for iE = nE:-1:1 % Cosmetic enhancement: a set of thinner ellipses to improve line crossings
    ellipse(0,0,p-ellOffset(iE),ellOffset(iE),colE(iE,:),4,'-');
end
for iE = 1:nE % Cosmetic enhancement: a set of thinner ellipses to improve line crossings
    ellipse(0,0,p-ellOffset(iE),ellOffset(iE),colE(iE,:),3,'-');
end
t = linspace(0,2*pi,300);
x = p/4*(3*cos(t) + cos(3*t));
y = p/4*(3*sin(t) - sin(3*t));
plot(x,y,'-','color',astroidCol,'linewidth',5); % plot Astroid

for iE = nE:-1:1 % Cosmetic enhancement: a set of thinner ellipses to improve line crossings and overlap with Astroid
    ellipse(0,0,p-ellOffset(iE),ellOffset(iE),colE(iE,:),2,'-');
end
for iE = 1:nE % Cosmetic enhancement: a set of thinner ellipses to improve line crossings and overlap with Astroid
    ellipse(0,0,p-ellOffset(iE),ellOffset(iE),colE(iE,:),1,'-');
end

ellipse(0,0, p,p, fullCircCol,5,'-'); % outer circle
reducedRGBimage = uint8(ones(xySize,xySize,3,nFrames)); % allocate
h_del = []; % init

for iFrame = 1:nFrames
    delete(h_del) % fresh frame    
    
    currentAngle = angleSmall(iFrame);
    
    h_del(1) = plot( p*cos(currentAngle), p*sin(currentAngle),'.','color',touchCol,'MarkerSize',40); % touching point / Instant centre of rotation
    h_del(2) = ellipse(cos(currentAngle)*0.5,sin(currentAngle)*0.5, p/2,p/2, halfCircCol,5,'-');  % half circle
    h_del(3) = plot([p*cos(currentAngle) 0],[0 p*sin(currentAngle)],'-','color',halfCircCol,'linewidth',7); % trammel bar
    h_del(4) = ellipse(cos(currentAngle)*0.75,sin(currentAngle)*0.75, p/4,p/4, quarterCircCol,5,'-');  % quarter circle
    
    xDot = linspace(p*cos(currentAngle), 0, nE);
    yDot = linspace(0, p*sin(currentAngle), nE);
    
    for iE = 1:nE
        h_del(4+iE) = plot(xDot(iE),yDot(iE),'o', 'LineWidth',1.5,...
            'MarkerSize',14,...
            'MarkerEdgeColor',halfCircCol,...
            'MarkerFaceColor',colE(iE,:));
    end
    
    x = p/4*(3*cos(currentAngle) + cos(3*currentAngle));
    y = p/4*(3*sin(currentAngle) - sin(3*currentAngle));
    
    h_del(5+iE) = plot([x cos(currentAngle)*0.75],[y sin(currentAngle)*0.75],'-','color',quarterCircCol,'linewidth',7);  % radius quarter circle
    h_del(6+iE) = plot(   cos(currentAngle)*0.75 ,   sin(currentAngle)*0.75 ,'.','color',quarterCircCol,'MarkerSize',36);  % center quarter circle
    h_del(7+iE) = plot(x,y,'o', 'LineWidth',1.5,...
        'MarkerSize',14,...
        'MarkerEdgeColor',quarterCircCol,...
        'MarkerFaceColor',astroidCol); % plot point on Astroid

    %% save animation
    drawnow
    % pause(0.01)
    f = getframe(figHandle);
    reducedRGBimage(:,:,:,iFrame) = imReduceSize(f.cdata,scaleReduction); % the size reduction: adds antialiasing
    
    if iFrame == 1 % svg; create colormap; allocate
        if ~isempty(which('plot2svg'))
            plot2svg(fullfile(pathstr, [fname '_a.svg']),figHandle) % by Juerg Schwizer, See http://www.zhinst.com/blogs/schwizer/
        else
            disp('plot2svg.m not available; see http://www.zhinst.com/blogs/schwizer/');
        end
    end
    
end

map = createImMap(reducedRGBimage,64,[1 1 1;colE;fullCircCol;halfCircCol;quarterCircCol;astroidCol;touchCol]); % create color map; enforce self defined colores
im = uint8(ones(xySize,xySize,1,nFrames)); % allocate
for iFrame = 1:nFrames
    im(:,:,1,iFrame) = rgb2ind(reducedRGBimage(:,:,:,iFrame),map,'nodither');
end

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


function h = ellipse(x,y,a,b,col,linw,lSty)
% x coordinates of the center
% y coordinates of the center
% a radius1
% b radius2
angleOffPoints = linspace(0,2.001*pi,300);
xe = x + a*cos(angleOffPoints);
ye = y + b*sin(angleOffPoints);
h = plot(xe,ye,lSty,'linewidth',linw,'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:
%  imNew:  unit8([imRows/redSize x imColumns/redSize x nChannel x nStack])
%
% 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 colores; 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 colores
    [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];