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