File:AnimatedGears.gif

function animatedGears()
% Source code for drawing gears
% The shape of the gears is not precise, it creates a decent GIF and a SVG.
%
% 2019-05-12 Jahobr

[pathstr,fname] = fileparts(which(mfilename)); % save files under the same name and at file location

RGB.bkgd   = [1   1   1  ]; % white background
RGB.black  = [0   0   0  ]; %
RGB.green  = [0.1 0.7 0.1]; %
RGB.yellow = [1   0.7 0  ]; %
RGB.blue   = [0.2 0.2 1  ]; %
RGB.red    = [1   0.2 0.2]; %
% violetRGB = [0.6 0.2 0.8]; %

RGB = structfun(@(q)round(q*255)/255, RGB, 'UniformOutput',false); % round to values that are nicely uint8 compatible

figHandle = figure(15674459); clf
set(figHandle,'Units','pixel');
set(figHandle,'MenuBar','none',  'ToolBar','none'); % free real estate for a maximally large image
set(figHandle,'Color',RGB.bkgd); % white background
axesHandle = axes;
hold(axesHandle,'on')
axis off % invisible axes (no ticks)
axis equal;

for currentCase = 3:8
    switch currentCase
        case 1 % animatedGears
            saveName = 'animatedGears';
            nFrames = 240;
            
            teeth    = [14, 3*14, 14, 2*14];
            module   = [ 2,    2,  3,    3]; % gear size
            diameter = module.*teeth;
            
            center1 = [0                                       0];
            center2 = [(diameter(1)+diameter(2))/2             0];
            center3 = [center2(1)+(diameter(3)+diameter(4))/2  0];
            
            xLimits = [center1(1)-diameter(1)/2-2*module(1)           center3(1)+diameter(4)/2+module(4)+module(1)]; % use a rim of size "module(1)"
            yLimits = [center3(2)-diameter(4)/2-module(4)-module(1)   center3(2)+diameter(4)/2+module(4)+module(1)]; % use a rim of size "module(1)"
            
            maxMovementOfTheSlowestTooth = 2*pi/teeth(4);
            anglesSlow = linspace(0,maxMovementOfTheSlowestTooth,nFrames+1);
            anglesSlow = anglesSlow(1:end-1); % remove last frame, it would be double
            secondRatio  = teeth(4)/teeth(3);
            anglesMedium = anglesSlow.*secondRatio;
            anglesMedium = anglesMedium + (2*pi/teeth(3)) *0.5; % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
            firstRatio = teeth(2)/teeth(1);
            anglesFast = anglesMedium.*firstRatio;
            anglesFast = anglesFast + (2*pi/teeth(1)) *0.5; % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
            
        case 2 % animated_3_Gear_Row
            saveName = 'animated_3_Gear_Row';
            nFrames  = 80;
            
            teeth    = [14, 2*14, 3*14];
            module   = [1,  1,    1]; % gear size
            diameter = module.*teeth;
            
            center1 = [0                                       0];
            center2 = [(diameter(1)+diameter(2))/2             0];
            center3 = [center2(1)+(diameter(2)+diameter(3))/2  0];

            xLimits = [center1(1)-diameter(1)/2-module(1)*2   center3(1)+diameter(3)/2+module(3)*2]; % use a rim of size "module(1)"
            yLimits = [center3(2)-diameter(3)/2-module(3)*2   center3(2)+diameter(3)/2+module(3)*2];  % use a rim of size "module(1)"
            
            maxMovementOfTheSlowestTooth = 2*pi/teeth(3);
            anglesSlow = linspace(0,maxMovementOfTheSlowestTooth,nFrames+1);
            anglesSlow = anglesSlow(1:end-1); % remove last frame, it would be double
            secondRatio  = teeth(3)/teeth(2);
            anglesMedium = anglesSlow.*secondRatio;
            anglesMedium = anglesMedium + (2*pi/teeth(2)) *0.5; % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
            firstRatio = teeth(2)/teeth(1);
            anglesFast = anglesMedium.*firstRatio;
            anglesFast = anglesFast + (2*pi/teeth(1)) *0.5; % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
            
         case 3 % animated_internal_gear
            saveName = 'animated_internal_gear';
            nFrames  = 80;
            teeth    = [14, 3*14];

            module   = [1,  1]; % gear size
            diameter = module.*teeth;
            
            center1 = [0 (diameter(1)-diameter(2))/2];
            center2 = [0                           0];
            
            xLimits = [center2(1)-diameter(2)/1.5   center2(1)+diameter(2)/1.5]; % 
            yLimits = [center2(2)-diameter(2)/1.5   center2(2)+diameter(2)/1.5]; % 
            
            maxMovementOfTheSlowestTooth = 2*pi/teeth(2);
            anglesSlow = linspace(0,maxMovementOfTheSlowestTooth,nFrames+1);
            anglesSlow = anglesSlow(1:end-1); % remove last frame, it would be double
            ratio  = teeth(2)/teeth(1);
            anglesFast = anglesSlow.*ratio;
            anglesFast = anglesFast + (2*pi/teeth(1)) *0.5; % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
            
         case 4 % animated_rack_and_pinion
            saveName = 'animated_rack_and_pinion';
            nFrames  = 80;
            teeth    = [14, 2*14];

            module   = [1,  1]; % gear size
            diameter = module.*teeth;
            
            center1 = [0              0];
            center2 = [0 -diameter(1)/2];
            
            xLimits = [center1(1)-diameter(1)*1.25       center1(1)+diameter(1)*1.25]; % 
            yLimits = [center1(2)-diameter(1)   center1(2)+diameter(1)/2+module(1)*3]; % 
            
            maxMovementOfTheSlowestTooth = 2*pi/teeth(1);
            anglesSlow = linspace(0,maxMovementOfTheSlowestTooth,nFrames+1);
            anglesSlow = anglesSlow(1:end-1); % remove last frame, it would be double
            sideShift = anglesSlow.*diameter(1)/2;
            sideShift = sideShift + module(1)*3/8; % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
            
        case 5 % animated_two_spur_gears
            saveName = 'animated_two_spur_gears_1_1';
            nFrames  = 80;
            teeth    = [14, 14];
            
        case 6 % animated_two_spur_gears
            saveName = 'animated_two_spur_gears_1_2';
            nFrames  = 80;
            teeth    = [14, 2*14];
            
        case 7 % animated_two_spur_gears
            saveName = 'animated_two_spur_gears_1_3';
            nFrames  = 80;
            teeth    = [14, 3*14];
            
        case 8 % animated_two_spur_gears
            saveName = 'animated_two_spur_gears_1_5';
            nFrames  = 80;
            teeth    = [14, 5*14];
    end
    
    if  contains(saveName,'animated_two_spur_gears') % all versions
        module   = [1,  1]; % gear size
        diameter = module.*teeth;
        
        center1 = [0                           0];
        center2 = [(diameter(1)+diameter(2))/2 0];
        
        xLimits = [center1(1)-diameter(1)/2-module(1)*2   center2(1)+diameter(2)/2+module(2)*2]; % use a rim of size "module(1)"
        yLimits = [center2(2)-diameter(2)/2-module(2)*2   center2(2)+diameter(2)/2+module(2)*2]; % use a rim of size "module(1)"
        
        maxMovementOfTheSlowestTooth = 2*pi/teeth(2);
        anglesSlow = linspace(0,maxMovementOfTheSlowestTooth,nFrames+1);
        anglesSlow = anglesSlow(1:end-1); % remove last frame, it would be double
        ratio  = teeth(2)/teeth(1);
        anglesFast = anglesSlow.*ratio;
        anglesFast = anglesFast + (2*pi/teeth(1)) *0.5; % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
    end
    
    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)
    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
    
    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

    setXYlim(axesHandle,xLimits,yLimits); % set limits and drawnow;
    
    reducedRGBimage = uint8(ones(ySize,xSize,3,nFrames)); % allocate

    liSc = mean([xSize ySize]*scaleReduction)/350; % LineWidth scale; LineWidth is absolut, a bigger images needs thicker lines to keep them in proportion

    for iFrame = 1:nFrames
        cla(axesHandle)
        
        switch currentCase
            case 1 % animatedGears
                drawSpurWheel(center1,teeth(1),module(1),RGB.red   ,liSc,RGB.black,-anglesFast(iFrame));   % left cogwheel  (fast)
                drawSpurWheel(center2,teeth(2),module(2),RGB.blue  ,liSc,RGB.black, anglesMedium(iFrame)); % cogwheel  2 (center)
                drawSpurWheel(center2,teeth(3),module(3),RGB.yellow,liSc,RGB.black, anglesMedium(iFrame)); % cogwheel  3 (center)
                drawSpurWheel(center3,teeth(4),module(4),RGB.green ,liSc,RGB.black,-anglesSlow(iFrame));   % right cogwheel  (slow)
            case 2 % animated_3_Gear_Row
                drawSpurWheel(center1,teeth(1),module(1),RGB.red   ,liSc,RGB.black,-anglesFast(iFrame));   % left cogwheel  (fast)
                drawSpurWheel(center2,teeth(2),module(2),RGB.blue  ,liSc,RGB.black, anglesMedium(iFrame)); % cogwheel  2 (center)
                drawSpurWheel(center3,teeth(3),module(3),RGB.green ,liSc,RGB.black,-anglesSlow(iFrame));   % right cogwheel  (slow)
            case 3 % animated_internal gear
                drawSpurWheel(center1,teeth(1),module(1),RGB.red   ,liSc,RGB.black,-anglesFast(iFrame));   % left cogwheel  (fast)
                drawRingGear(teeth(2),module(2),RGB.green,liSc,RGB.black, -anglesSlow(iFrame));
            case 4 % animated_rack_and_pinion
                drawSpurWheel(center1,teeth(1),module(1),RGB.red   ,liSc,RGB.black,-anglesSlow(iFrame));   % left cogwheel  (fast)
                drawRack(center2,teeth(2),module(2),RGB.green,liSc,RGB.black,-sideShift(iFrame),0);
            case {5,6,7,8} % animated_two_spur_gears
                drawSpurWheel(center1,teeth(1),module(1),RGB.red   ,liSc,RGB.black,-anglesFast(iFrame));   % left cogwheel  (fast)
                drawSpurWheel(center2,teeth(2),module(2),RGB.green ,liSc,RGB.black, anglesSlow(iFrame));   % right cogwheel  (slow)
        end
        
        setXYlim(axesHandle,xLimits,yLimits); % reset limits and drawnow

        f = getframe(figHandle);
        reducedRGBimage(:,:,:,iFrame) = imReduceSize(f.cdata,scaleReduction); % allows subpixel lines       
        
        if iFrame == 1 % SVG
            if ~isempty(which('plot2svg'))
                plot2svg(fullfile(pathstr, [saveName '_Frame1.svg']),figHandle) % by Juerg Schwizer
            else
                disp('plot2svg.m not available; see http://www.zhinst.com/blogs/schwizer/');
            end
        end
    end
    
    switch currentCase % 
        case 1 % animatedGears
            map = createImMap(reducedRGBimage,32,[RGB.bkgd; RGB.black; RGB.red; RGB.blue; RGB.yellow; RGB.green]); % colormap
        case 2 % animated_3_Gear_Row
            map = createImMap(reducedRGBimage,32,[RGB.bkgd; RGB.black; RGB.red; RGB.blue;             RGB.green]); % colormap
        case 3 % animated_internal gear
            map = createImMap(reducedRGBimage,16,[RGB.bkgd; RGB.black; RGB.red;                       RGB.green]); % colormap
        case 4 % animated_rack_and_pinion
            map = createImMap(reducedRGBimage,16,[RGB.bkgd; RGB.black; RGB.red;                       RGB.green]); % colormap
        case {5,6,7,8} % animated_two_spur_gears
            map = createImMap(reducedRGBimage,16,[RGB.bkgd; RGB.black; RGB.red;                       RGB.green]); % colormap
    end
    
    im = uint8(ones(ySize,xSize,1,nFrames)); % allocate
    for iFrame = 1:nFrames
        im(:,:,1,iFrame) = rgb2ind(reducedRGBimage(:,:,:,iFrame),map,'nodither');
    end
    
    imwrite(im,map,fullfile(pathstr, [saveName '.gif']),'DelayTime',1/60,'LoopCount',inf) % save gif
    disp([saveName '.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 100 MP limit
end

function drawSpurWheel(center,toothNumber,module,fillC,linW,linC,startOffset)
% DRAWSPURWHEEL - draw a simple Toothed Wheel
%    center:       [x y]
%    toothNumber:  scalar
%    module:       scalar tooth "size"
%    fillC:        color of filling [r g b]
%    linW:         LineWidth
%    linC:         LineColor
%    startOffset:  start rotation (scalar)[rad]

effectiveRadius = module*toothNumber/2; % effective Radius

outsideRadius =     effectiveRadius+1*  module; %                +---+             +---+
upperRisingRadius = effectiveRadius+0.5*module; %               /     \           /     \
% effective Radius                              %              /       \         /       \
lowerRisingRadius = effectiveRadius-0.5*module; %             I         I       I         I
rootRadius =        effectiveRadius-1.1*module; %     + - - - +         + - - - +         +

angleBetweenTeeth = 2*pi/toothNumber; % angle between 2 teeth
angleOffPoints = (0:angleBetweenTeeth/16:(2*pi));
angleOffPoints = angleOffPoints+startOffset; % apply rotation offset

angleOffPoints(7:16:end)  = angleOffPoints(7:16:end)  + 1/toothNumber^1.2; % hack to create smaller tooth tip
angleOffPoints(11:16:end) = angleOffPoints(11:16:end) - 1/toothNumber^1.2; % hack to create smaller tooth tip

angleOffPoints(8:16:end)  = (angleOffPoints(7:16:end)  + angleOffPoints(9:16:end))/2; % shift the neighbouring tip point in accordingly
angleOffPoints(10:16:end) = (angleOffPoints(11:16:end) + angleOffPoints(9:16:end))/2; % shift the neighbouring tip point in accordingly

angleOffPoints(6:16:end)  = angleOffPoints(6:16:end)  + 1/toothNumber^1.7; % hack to create slender upperRisingRadius
angleOffPoints(12:16:end) = angleOffPoints(12:16:end) - 1/toothNumber^1.7; % hack to create slender upperRisingRadius

radiusOffPoints = angleOffPoints; % allocate with correct site

radiusOffPoints( 1:16:end) = rootRadius;        % center bottom         I
radiusOffPoints( 2:16:end) = rootRadius;        % left bottom           I
radiusOffPoints( 3:16:end) = rootRadius;        % left bottom corner    +
radiusOffPoints( 4:16:end) = lowerRisingRadius; % lower rising bottom      \
radiusOffPoints( 5:16:end) = effectiveRadius;   % rising edge                 \
radiusOffPoints( 6:16:end) = upperRisingRadius; % upper rising edge              \
radiusOffPoints( 7:16:end) = outsideRadius;     % right top  corner                 +
radiusOffPoints( 8:16:end) = outsideRadius;     % right top                         I
radiusOffPoints( 9:16:end) = outsideRadius;     % center top                        I
radiusOffPoints(10:16:end) = outsideRadius;     % left top                          I
radiusOffPoints(11:16:end) = outsideRadius;     % left top  corner                  +
radiusOffPoints(12:16:end) = upperRisingRadius; % upper falling edge             /
radiusOffPoints(13:16:end) = effectiveRadius;   % falling edge                /
radiusOffPoints(14:16:end) = lowerRisingRadius; % lower falling edge       /
radiusOffPoints(15:16:end) = rootRadius;        % right bottom  corner  +
radiusOffPoints(16:16:end) = rootRadius;        % right bottom          I

[X,Y] = pol2cart(angleOffPoints,radiusOffPoints);
X = X+center(1); % center offset
Y = Y+center(2); % center offset
patch(X,Y,fillC,'EdgeColor',linC,'LineWidth',linW)

% % effective Radius
% [X,Y] = pol2cart(angleOffPoints,effectiveRadius);
% X = X+center(1); % center offset
% Y = Y+center(2); % center offset
% plot(X,Y,'-.','Color',linC);

%% shaft
shaftRadius = module*6 /2; % small radius, assuming the effective radius a 6-tooth wheel would have
[X,Y] = pol2cart(angleOffPoints,shaftRadius);
X = X+center(1); % center offset
Y = Y+center(2); % center offset
plot(X,Y,'LineWidth',linW,'Color',linC);
% plot(center(1),center(2),'.','Color',linC)


function drawRingGear(toothNumber,module,fillC,linW,linC,startOffset)
% DRAWRINGGEAR - draw a outer gear
%    center:       [x y]
%    toothNumber:  scalar
%    module:       scalar tooth "size"
%    fillC:        color of filling [r g b]
%    linW:         LineWidth
%    linC:         LineColor
%    startOffset:  start rotation (scalar)[rad]
effectiveRadius = module*toothNumber/2; % effective effectiveRadius

outsideRadius     = effectiveRadius-1*  module; %                +---+             +---+
upperRisingRadius = effectiveRadius-0.5*module; %               /     \           /     \
% effective Radius                              %              /       \         /       \
lowerRisingRadius = effectiveRadius+0.5*module; %             I         I       I         I
rootRadius        = effectiveRadius+1.1*module; %     + - - - +         + - - - +         +

angleBetweenTeeth = 2*pi/toothNumber; % angle between 2 teeth
angleOffPoints = (0:angleBetweenTeeth/16:(2*pi));
angleOffPoints = angleOffPoints+startOffset; % apply rotation offset

%% outerEdge
maxRadius = rootRadius*1.2; % definition of outer line
[Xout,Yout] = pol2cart(angleOffPoints,maxRadius);

%% inner teeth
radiusOffPoints = angleOffPoints; % init

angleOffPoints( 7:16:end) = angleOffPoints(7:16:end)  + 1/toothNumber^1.2; % hack to create smaller tooth tip
angleOffPoints(11:16:end) = angleOffPoints(11:16:end) - 1/toothNumber^1.2; % hack to create smaller tooth tip

angleOffPoints( 8:16:end) = (angleOffPoints(7:16:end)  + angleOffPoints(9:16:end))/2; % shift the neighbouring tip point in accordingly
angleOffPoints(10:16:end) = (angleOffPoints(11:16:end) + angleOffPoints(9:16:end))/2; % shift the neighbouring tip point in accordingly

angleOffPoints( 6:16:end) = angleOffPoints(6:16:end)  + 1/toothNumber^1.7; % hack to create slender tooth
angleOffPoints(12:16:end) = angleOffPoints(12:16:end) - 1/toothNumber^1.7; % hack to create slender tooth

radiusOffPoints( 1:16:end) = rootRadius;        % center bottom         I
radiusOffPoints( 2:16:end) = rootRadius;        % left bottom           I
radiusOffPoints( 3:16:end) = rootRadius;        % left bottom corner    +
radiusOffPoints( 4:16:end) = lowerRisingRadius; % lower rising bottom      \
radiusOffPoints( 5:16:end) = effectiveRadius;   % rising edge                 \
radiusOffPoints( 6:16:end) = upperRisingRadius; % upper rising edge              \
radiusOffPoints( 7:16:end) = outsideRadius;     % right top corner                 +
radiusOffPoints( 8:16:end) = outsideRadius;     % right top                        I
radiusOffPoints( 9:16:end) = outsideRadius;     % center top                       I
radiusOffPoints(10:16:end) = outsideRadius;     % left top                         I
radiusOffPoints(11:16:end) = outsideRadius;     % left top corner                  +
radiusOffPoints(12:16:end) = upperRisingRadius; % upper falling edge             /
radiusOffPoints(13:16:end) = effectiveRadius;   % falling edge                /
radiusOffPoints(14:16:end) = lowerRisingRadius; % lower falling edge       /
radiusOffPoints(15:16:end) = rootRadius;        % right bottom corner   +
radiusOffPoints(16:16:end) = rootRadius;        % right bottom          I

[X,Y] = pol2cart(angleOffPoints,radiusOffPoints);

[Xout,Yout] = poly2cw(Xout,Yout);
[X,   Y   ] = poly2cw(X   ,Y   );
[Xb,Yb] = polybool('subtraction',Xout,Yout, X,Y);
Xb = Xb(~isnan(Xb)); % notNaN
Yb = Yb(~isnan(Yb)); % notNaN
patch(Xb,Yb,fillC,'EdgeColor','none')

plot(X,   Y,   'LineWidth',linW,'Color',linC); % draw teeth outline
plot(Xout,Yout,'LineWidth',linW,'Color',linC); % draw outer circle


function drawRack(center,toothNumber,module,fillC,linW,linC,startOffset,top)
%    center:       [x y]
%    toothNumber:  scalar
%    module:       scalar tooth "size"
%    fillC:        color of filling [r g b]
%    linW:         LineWidth
%    linC:         LineColor
%    startOffset: initial shift
%    top:         1=top 0=bottom
x = (0:toothNumber*4-2)*pi*module/4;
x = x-mean(x)+center(1)+startOffset;

y = zeros(size(x));

y(1:4:end) = y(1:4:end)+1.1*module; %    +###I bottom
y(2:4:end) = y(2:4:end)-1  *module; % +######I tip
y(3:4:end) = y(3:4:end)-1  *module; % +######I tip
y(4:4:end) = y(4:4:end)+1.1*module; %    +###I bottom

x(1:4:end) = x(1:4:end)-0.14*module; % bottom smaller
x(2:4:end) = x(2:4:end)+0.14*module; % tip smaller
x(3:4:end) = x(3:4:end)-0.14*module; % tip smaller
x(4:4:end) = x(4:4:end)+0.14*module; % bottom smaller

x = [x(1)     x  x(end)];
y = [5*module y  5*module];

if ~top
    y = -y; % flip
end
y = y+center(2);
patch(x,y,fillC,'EdgeColor',linC,'LineWidth',linW);


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;


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/scaleReduction ,'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];