function outputFile = make_rotating_4d_cube_video(outputFile, ... durationSeconds, haunterImageFile) %MAKE_ROTATING_4D_CUBE_VIDEO Render a faceted, eldritch tesseract. % % outputFile = MAKE_ROTATING_4D_CUBE_VIDEO() % creates "Z:\cube\rotating_4d_shining_cube.mp4". % % outputFile = MAKE_ROTATING_4D_CUBE_VIDEO(outputFile, durationSeconds, % haunterImageFile) lets you choose the destination, duration and monster % artwork. The defaults are: % Resolution: 720 x 576 pixels % Frame rate: 25 frames per second % Duration: 10 seconds (250 frames) % Rotation: one complete 360-degree rotation in the X-W plane % Appearance: dark translucent facets with shifting internal colour, % red striations, travelling glints and luminous edges % Haunter: one centred replacement frame every 1/13 second % (13 flashes per second, scheduled without timing drift), % growing from 12% size to full-frame by its final flash % % By default, "haunter_of_the_dark.png" is read from the same folder as % this M-file. The supplied project ZIP contains the required artwork. % % The last stored frame deliberately precedes the first by one angular % step. This avoids duplicating the first frame and makes the result loop % seamlessly. No toolboxes are required. % if nargin < 1 || isempty(outputFile) outputFile = fullfile('Z:\cube', 'rotating_4d_shining_cube.mp4'); end if nargin < 2 || isempty(durationSeconds) durationSeconds = 10; end if nargin < 3 || isempty(haunterImageFile) programFolder = fileparts(mfilename('fullpath')); haunterImageFile = fullfile(programFolder, ... 'haunter_of_the_dark.png'); end if ~(ischar(outputFile) || (isstring(outputFile) && isscalar(outputFile))) error('outputFile must be a character vector or string scalar.'); end outputFile = char(outputFile); if ~(isnumeric(durationSeconds) && isscalar(durationSeconds) && ... isfinite(durationSeconds) && durationSeconds > 0) error('durationSeconds must be a positive finite scalar.'); end if ~(ischar(haunterImageFile) || ... (isstring(haunterImageFile) && isscalar(haunterImageFile))) error('haunterImageFile must be a character vector or string scalar.'); end haunterImageFile = char(haunterImageFile); if ~isfile(haunterImageFile) error(['Haunter artwork not found: %s\nKeep ', ... 'haunter_of_the_dark.png beside this M-file, or pass its ', ... 'full path as the third argument.'], haunterImageFile); end frameRate = 25; frameWidth = 720; frameHeight = 576; supersample = 2; numberOfFrames = round(durationSeconds * frameRate); flashRate = 13; initialHaunterScale = 0.12; if numberOfFrames < 2 error('The duration must produce at least two video frames.'); end % Map the mathematically exact 13 Hz event times to the nearest available % 25 fps video frames. The resulting gaps are one or two frames long, % but the schedule never drifts and gives exactly 13 flashes per complete % second. Beginning at frame one also preserves the cadence across a loop. actualDuration = numberOfFrames / frameRate; flashTimes = (0:(ceil(actualDuration * flashRate) - 1)) / flashRate; flashTimes = flashTimes(flashTimes < actualDuration); flashFrameNumbers = unique(round(flashTimes * frameRate) + 1); flashFrameNumbers = flashFrameNumbers( ... flashFrameNumbers <= numberOfFrames); flashFrameMask = false(1, numberOfFrames); flashFrameMask(flashFrameNumbers) = true; flashOrdinalForFrame = zeros(1, numberOfFrames); flashOrdinalForFrame(flashFrameNumbers) = 1:numel(flashFrameNumbers); haunterFrame = prepareHaunterFrame(haunterImageFile, ... frameWidth, frameHeight); outputFolder = fileparts(outputFile); if ~isempty(outputFolder) && ~isfolder(outputFolder) mkdir(outputFolder); end [vertices, edges, faces] = makeTesseract(); % A fixed initial orientation exposes all four dimensions clearly. initialRotation = ... planeRotation(4, 1, 2, deg2rad(14)) * ... planeRotation(4, 2, 3, deg2rad(-19)) * ... planeRotation(4, 3, 4, deg2rad(17)); % A fixed 3-D camera orientation makes the perspective projection easy % to read while the actual animation remains a true 4-D X-W rotation. cameraRotation = ... planeRotation(3, 1, 2, deg2rad(22)) * ... planeRotation(3, 2, 3, deg2rad(-18)) * ... planeRotation(3, 1, 3, deg2rad(11)); % Measure the entire animation before rendering, so the drawing scale is % constant and the tesseract does not appear to breathe in and out. maximumExtent = 0; for frameNumber = 1:numberOfFrames angle = 2 * pi * (frameNumber - 1) / numberOfFrames; rotation4D = planeRotation(4, 1, 4, angle); rotatedVertices = rotation4D * initialRotation * vertices; [projected, ~] = project4Dto2D(rotatedVertices, cameraRotation); maximumExtent = max(maximumExtent, max(abs(projected(:)))); end drawingScale = 0.41 * min(frameWidth, frameHeight) / maximumExtent; writer = VideoWriter(outputFile, 'MPEG-4'); writer.FrameRate = frameRate; writer.Quality = 95; open(writer); writerCleanup = onCleanup(@() safelyCloseWriter(writer)); fprintf('Creating %s\n', outputFile); fprintf('%d frames at %d fps, %d x %d pixels\n', ... numberOfFrames, frameRate, frameWidth, frameHeight); fprintf(['Haunter inserts: %d single frames at %g Hz ', ... '(photosensitivity warning)\n'], numel(flashFrameNumbers), flashRate); fprintf('Haunter growth: %.0f%% to full-frame over successive flashes\n', ... 100 * initialHaunterScale); for frameNumber = 1:numberOfFrames angle = 2 * pi * (frameNumber - 1) / numberOfFrames; if flashFrameMask(frameNumber) % Each successive one-frame apparition is larger. Scaling by % flash ordinal, rather than absolute frame number, guarantees % that the final apparition is exactly full-frame even when the % 13 Hz cadence does not place a flash on the movie's last frame. flashOrdinal = flashOrdinalForFrame(frameNumber); if numel(flashFrameNumbers) > 1 growthProgress = (flashOrdinal - 1) / ... (numel(flashFrameNumbers) - 1); else growthProgress = 1; end haunterScale = initialHaunterScale + ... (1 - initialHaunterScale) * growthProgress; frame = makeScaledHaunterFrame(haunterFrame, haunterScale, ... frameWidth, frameHeight); else rotation4D = planeRotation(4, 1, 4, angle); rotatedVertices = rotation4D * initialRotation * vertices; [projected, depth, vertices3D] = ... project4Dto2D(rotatedVertices, cameraRotation); frame = renderShiningTesseract(projected, vertices3D, depth, ... edges, faces, angle, drawingScale, frameWidth, frameHeight, ... supersample); end writeVideo(writer, frame); if mod(frameNumber, frameRate) == 0 || frameNumber == numberOfFrames fprintf('Rendered %d of %d frames\n', frameNumber, numberOfFrames); end end close(writer); clear writerCleanup outputFile = char(java.io.File(outputFile).getCanonicalPath()); fprintf('Finished: %s\n', outputFile); end function frame = makeScaledHaunterFrame(fullFrame, scale, ... frameWidth, frameHeight) % Resize the prepared artwork and centre it on black without toolboxes. if scale >= 1 - eps frame = fullFrame; return end scaledWidth = max(1, min(frameWidth, round(frameWidth * scale))); scaledHeight = max(1, min(frameHeight, round(frameHeight * scale))); [queryX, queryY] = meshgrid(linspace(1, frameWidth, scaledWidth), ... linspace(1, frameHeight, scaledHeight)); scaledImage = zeros(scaledHeight, scaledWidth, 3); source = double(fullFrame) / 255; for channel = 1:3 scaledImage(:, :, channel) = interp2(source(:, :, channel), ... queryX, queryY, 'linear'); end scaledImage = uint8(round(255 * min(1, max(0, scaledImage)))); frame = zeros(frameHeight, frameWidth, 3, 'uint8'); firstColumn = floor((frameWidth - scaledWidth) / 2) + 1; firstRow = floor((frameHeight - scaledHeight) / 2) + 1; lastColumn = firstColumn + scaledWidth - 1; lastRow = firstRow + scaledHeight - 1; frame(firstRow:lastRow, firstColumn:lastColumn, :) = scaledImage; end function frame = prepareHaunterFrame(imageFile, frameWidth, frameHeight) % Load, centre-crop and resize the supplied artwork using base MATLAB only. [source, ~, alpha] = imread(imageFile); if size(source, 3) == 1 source = repmat(source, 1, 1, 3); elseif size(source, 3) > 3 source = source(:, :, 1:3); end source = normaliseImage(source); if ~isempty(alpha) alpha = normaliseImage(alpha); source = source .* repmat(alpha, 1, 1, 3); end [sourceHeight, sourceWidth, ~] = size(source); targetAspect = frameWidth / frameHeight; sourceAspect = sourceWidth / sourceHeight; if sourceAspect > targetAspect croppedWidth = max(1, round(sourceHeight * targetAspect)); firstColumn = floor((sourceWidth - croppedWidth) / 2) + 1; source = source(:, firstColumn:(firstColumn + croppedWidth - 1), :); elseif sourceAspect < targetAspect croppedHeight = max(1, round(sourceWidth / targetAspect)); firstRow = floor((sourceHeight - croppedHeight) / 2) + 1; source = source(firstRow:(firstRow + croppedHeight - 1), :, :); end [sourceHeight, sourceWidth, ~] = size(source); [queryX, queryY] = meshgrid(linspace(1, sourceWidth, frameWidth), ... linspace(1, sourceHeight, frameHeight)); resized = zeros(frameHeight, frameWidth, 3); for channel = 1:3 resized(:, :, channel) = interp2(source(:, :, channel), ... queryX, queryY, 'linear'); end % Lift shadow detail slightly so the creature registers during its single % 1/25-second appearance while retaining the near-black atmosphere. resized = min(1, max(0, resized)) .^ 0.72; frame = uint8(round(255 * resized)); end function image = normaliseImage(image) % Convert integer, logical or floating-point image data to doubles in [0, 1]. if isinteger(image) image = double(image) / double(intmax(class(image))); elseif islogical(image) image = double(image); else image = double(image); if max(image(:)) > 1 image = image / 255; end end image = min(1, max(0, image)); end function [vertices, edges, faces] = makeTesseract() % Sixteen vertices at every combination of -1 and +1 in four dimensions. vertices = zeros(4, 16); for vertexIndex = 0:15 coordinates = bitget(uint16(vertexIndex), 1:4); vertices(:, vertexIndex + 1) = 2 * double(coordinates(:)) - 1; end % Thirty-two edges. Column three records the dimension in which the two % endpoints differ, allowing each family of parallel edges its own hue. edges = zeros(32, 3); edgeIndex = 0; for vertexIndex = 0:15 for dimension = 1:4 otherVertex = bitxor(uint16(vertexIndex), ... bitshift(uint16(1), dimension - 1)); if vertexIndex < double(otherVertex) edgeIndex = edgeIndex + 1; edges(edgeIndex, :) = [vertexIndex + 1, ... double(otherVertex) + 1, dimension]; end end end % Twenty-four square faces. The final two columns record the pair of % varying dimensions, used to give each family of facets its own hue. faces = zeros(24, 6); faceIndex = 0; for dimensionA = 1:3 for dimensionB = dimensionA + 1:4 fixedDimensions = setdiff(1:4, [dimensionA, dimensionB]); for fixedState = 0:3 baseVertex = uint16(0); for fixedBit = 1:2 if bitget(uint16(fixedState), fixedBit) baseVertex = bitset(baseVertex, ... fixedDimensions(fixedBit), 1); end end cornerA = bitset(baseVertex, dimensionA, 1); cornerB = bitset(baseVertex, dimensionB, 1); oppositeCorner = bitset(cornerA, dimensionB, 1); faceIndex = faceIndex + 1; faces(faceIndex, :) = [double([baseVertex, cornerA, ... oppositeCorner, cornerB]) + 1, dimensionA, dimensionB]; end end end end function rotation = planeRotation(numberOfDimensions, axisA, axisB, angle) % Rotation in the plane defined by axisA and axisB. rotation = eye(numberOfDimensions); cosine = cos(angle); sine = sin(angle); rotation(axisA, axisA) = cosine; rotation(axisA, axisB) = -sine; rotation(axisB, axisA) = sine; rotation(axisB, axisB) = cosine; end function [projected2D, depth, vertices3D] = ... project4Dto2D(vertices4D, cameraRotation) % Perspective projection from 4-D to 3-D, followed by 3-D to 2-D. distance4D = 4.5; perspective4D = distance4D ./ (distance4D - vertices4D(4, :)); vertices3D = vertices4D(1:3, :) .* perspective4D; vertices3D = cameraRotation * vertices3D; depth = vertices3D(3, :); distance3D = 5.5; perspective3D = distance3D ./ (distance3D - depth); projected2D = vertices3D(1:2, :) .* perspective3D; end function frame = renderShiningTesseract(points, vertices3D, depth, edges, ... faces, phase, drawingScale, frameWidth, frameHeight, supersample) % Render an obsidian-like solid with changing, apparently internal light. largeWidth = frameWidth * supersample; largeHeight = frameHeight * supersample; canvas = zeros(largeHeight, largeWidth, 3, 'uint8'); x = supersample * (frameWidth / 2 + drawingScale * points(1, :)); y = supersample * (frameHeight / 2 - drawingScale * points(2, :)); minimumDepth = min(depth); depthSpan = max(depth) - minimumDepth; if depthSpan < eps depthSpan = 1; end % The six face families are deliberately sombre. Their highlights shift % between bruised violet, deep blue, green, blood-red and cold white. facetColours = [ ... 75, 10, 100; ... 8, 48, 78; ... 105, 16, 42; ... 8, 68, 56; ... 42, 16, 110; ... 12, 72, 102]; faceDepth = mean(depth(faces(:, 1:4)), 2); [~, faceOrder] = sort(faceDepth, 'ascend'); lightDirection = [cos(phase); ... 0.65 * sin(2 * phase + 0.7); 0.85]; lightDirection = lightDirection / norm(lightDirection); viewDirection = [0; 0; 1]; halfVector = lightDirection + viewDirection; halfVector = halfVector / norm(halfVector); for orderIndex = 1:size(faces, 1) face = faces(faceOrder(orderIndex), :); vertexIndices = face(1:4); family = faceFamily(face(5), face(6)); sideA = vertices3D(:, vertexIndices(2)) - ... vertices3D(:, vertexIndices(1)); sideB = vertices3D(:, vertexIndices(4)) - ... vertices3D(:, vertexIndices(1)); normal = cross(sideA, sideB); normalLength = norm(normal); if normalLength < eps continue end normal = normal / normalLength; diffuse = abs(dot(normal, lightDirection)); specular = abs(dot(normal, halfVector)) ^ 30; internalPulse = 0.5 + 0.5 * sin(3 * phase + ... 1.13 * family + 0.31 * orderIndex); depthBrightness = (faceDepth(faceOrder(orderIndex)) - ... minimumDepth) / depthSpan; baseColour = facetColours(family, :); colour = baseColour * (0.28 + 0.55 * diffuse + ... 0.22 * internalPulse + 0.18 * depthBrightness) + ... [205, 225, 255] * (0.75 * specular); colour = uint8(round(min(255, colour))); opacity = min(0.42, 0.10 + 0.12 * diffuse + ... 0.13 * internalPulse + 0.12 * specular); canvas = fillPolygonAlpha(canvas, x(vertexIndices), ... y(vertexIndices), colour, opacity); % A narrow highlight slides over each facet. Different whole-number % rates keep the pattern complex while preserving the seamless loop. glintRate = 1 + mod(family + orderIndex, 3); glintPosition = 0.5 + 0.42 * sin(glintRate * phase + ... 0.83 * family + 0.19 * orderIndex); glintStart = (1 - glintPosition) * ... [x(vertexIndices(1)); y(vertexIndices(1))] + ... glintPosition * [x(vertexIndices(4)); y(vertexIndices(4))]; glintEnd = (1 - glintPosition) * ... [x(vertexIndices(2)); y(vertexIndices(2))] + ... glintPosition * [x(vertexIndices(3)); y(vertexIndices(3))]; glintStrength = 0.18 + 0.58 * specular + 0.20 * internalPulse; glintColour = uint8(round(min(255, glintStrength * ... [190, 215, 255] + 0.45 * double(colour)))); canvas = drawLine(canvas, glintStart(1), glintStart(2), ... glintEnd(1), glintEnd(2), glintColour, 1); % The Trapezohedron is described as almost black with red striations. striationPosition = 0.5 + 0.34 * sin(2 * phase + ... 0.57 * family + 0.11 * orderIndex); stripeStart = (1 - striationPosition) * ... [x(vertexIndices(1)); y(vertexIndices(1))] + ... striationPosition * [x(vertexIndices(2)); y(vertexIndices(2))]; stripeEnd = (1 - striationPosition) * ... [x(vertexIndices(4)); y(vertexIndices(4))] + ... striationPosition * [x(vertexIndices(3)); y(vertexIndices(3))]; stripeColour = uint8(round((0.18 + 0.26 * internalPulse) * ... [255, 42, 68])); canvas = drawLine(canvas, stripeStart(1), stripeStart(2), ... stripeEnd(1), stripeEnd(2), stripeColour, 1); end edgeDepth = (depth(edges(:, 1)) + depth(edges(:, 2))) / 2; [~, drawOrder] = sort(edgeDepth, 'ascend'); % X, Y and Z use spectral colours; W remains pale gold so the truly % fourth-dimensional connections can still be followed. dimensionColours = [ ... 68, 202, 255; ... 58, 238, 186; ... 186, 112, 255; ... 255, 205, 92]; for orderIndex = 1:size(edges, 1) edge = edges(drawOrder(orderIndex), :); averageDepth = (depth(edge(1)) + depth(edge(2))) / 2; brightness = 0.38 + 0.62 * ... (averageDepth - minimumDepth) / depthSpan; brightness = brightness * (0.88 + 0.12 * ... sin(2 * phase + 0.41 * orderIndex)); colour = uint8(round(brightness * dimensionColours(edge(3), :))); bloomColour = uint8(round(0.25 * double(colour))); canvas = drawLine(canvas, x(edge(1)), y(edge(1)), ... x(edge(2)), y(edge(2)), bloomColour, 6); canvas = drawLine(canvas, x(edge(1)), y(edge(1)), ... x(edge(2)), y(edge(2)), colour, 2); end % Draw the sixteen vertices after the edges so their intersections remain % crisp. Nearer vertices are brighter than distant ones. [~, vertexOrder] = sort(depth, 'ascend'); for orderIndex = 1:numel(vertexOrder) vertex = vertexOrder(orderIndex); brightness = 0.50 + 0.50 * ... (depth(vertex) - minimumDepth) / depthSpan; pulse = 0.88 + 0.12 * sin(3.1 * phase + 0.67 * vertex); colour = uint8(round(255 * brightness * pulse * ... [0.88, 0.95, 1.00])); haloColour = uint8(round(0.24 * double(colour))); canvas = drawDisc(canvas, round(x(vertex)), round(y(vertex)), ... 9, haloColour); canvas = drawDisc(canvas, round(x(vertex)), round(y(vertex)), ... 4, colour); end % Exact 2x box downsampling, requiring no Image Processing Toolbox. frameSum = uint16(canvas(1:2:end, 1:2:end, :)) + ... uint16(canvas(2:2:end, 1:2:end, :)) + ... uint16(canvas(1:2:end, 2:2:end, :)) + ... uint16(canvas(2:2:end, 2:2:end, :)); frame = uint8((double(frameSum) + 2) / 4); end function family = faceFamily(dimensionA, dimensionB) % Map the six possible pairs of varying dimensions to palette rows. pair = sort([dimensionA, dimensionB]); pairs = [1, 2; 1, 3; 1, 4; 2, 3; 2, 4; 3, 4]; family = find(all(pairs == pair, 2), 1); end function image = fillPolygonAlpha(image, polygonX, polygonY, colour, opacity) % Scan-line polygon fill with alpha blending; no toolbox functions required. [imageHeight, imageWidth, ~] = size(image); minimumY = max(1, ceil(min(polygonY))); maximumY = min(imageHeight, floor(max(polygonY))); numberOfVertices = numel(polygonX); planeSize = imageHeight * imageWidth; for row = minimumY:maximumY scanY = row + 0.5; intersections = zeros(1, numberOfVertices); count = 0; for vertex = 1:numberOfVertices nextVertex = mod(vertex, numberOfVertices) + 1; y1 = polygonY(vertex); y2 = polygonY(nextVertex); crosses = (y1 <= scanY && y2 > scanY) || ... (y2 <= scanY && y1 > scanY); if crosses count = count + 1; intersections(count) = polygonX(vertex) + ... (scanY - y1) * ... (polygonX(nextVertex) - polygonX(vertex)) / (y2 - y1); end end if count < 2 continue end intersections = sort(intersections(1:count)); for pairIndex = 1:2:(count - 1) firstColumn = max(1, ceil(intersections(pairIndex))); lastColumn = min(imageWidth, ... floor(intersections(pairIndex + 1))); if firstColumn > lastColumn continue end indices = sub2ind([imageHeight, imageWidth], ... row * ones(1, lastColumn - firstColumn + 1), ... firstColumn:lastColumn); for channel = 1:3 channelIndices = indices + (channel - 1) * planeSize; existing = double(image(channelIndices)); blended = (1 - opacity) * existing + ... opacity * double(colour(channel)); image(channelIndices) = uint8(round(blended)); end end end end function image = drawLine(image, x1, y1, x2, y2, colour, radius) % Draw a thick line into a uint8 RGB image without toolbox functions. numberOfSteps = max(2, ceil(max(abs(x2 - x1), abs(y2 - y1))) + 1); x = round(linspace(x1, x2, numberOfSteps)); y = round(linspace(y1, y2, numberOfSteps)); [imageHeight, imageWidth, ~] = size(image); for offsetY = -radius:radius for offsetX = -radius:radius if offsetX^2 + offsetY^2 > radius^2 continue end xx = x + offsetX; yy = y + offsetY; valid = xx >= 1 & xx <= imageWidth & yy >= 1 & yy <= imageHeight; if ~any(valid) continue end indices = sub2ind([imageHeight, imageWidth], yy(valid), xx(valid)); image = colourPixels(image, indices, colour); end end end function image = drawDisc(image, centreX, centreY, radius, colour) % Draw a filled circular vertex marker. [imageHeight, imageWidth, ~] = size(image); [offsetX, offsetY] = meshgrid(-radius:radius, -radius:radius); inside = offsetX.^2 + offsetY.^2 <= radius^2; x = centreX + offsetX(inside); y = centreY + offsetY(inside); valid = x >= 1 & x <= imageWidth & y >= 1 & y <= imageHeight; indices = sub2ind([imageHeight, imageWidth], y(valid), x(valid)); image = colourPixels(image, indices, colour); end function image = colourPixels(image, indices, colour) % Retain the brightest channel value where primitives overlap. planeSize = size(image, 1) * size(image, 2); for channel = 1:3 channelIndices = indices + (channel - 1) * planeSize; image(channelIndices) = max(image(channelIndices), colour(channel)); end end function safelyCloseWriter(writer) % Ensure a partially rendered video is closed if an error interrupts rendering. try close(writer); catch % The writer was already closed. end end