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