%% Script by http://scalarmotion.wordpress.com/2009/03/15/propeller-image-aliasing/

%% Simulation of image aliasing due to pixel-wise scanning of a propeller. 
%% Assumptions: 
%%      Propeller: 
%%          1. Constant angular velocity 
%%          2. Uniform visual cross-section of the blades (rectangular 
%%          blades) 
%%      Camera: 
%%          1. pixel-wise scanning left-to-right, from top to bottom. 
%%          2. no TV camera like interleaving 
%% Propeller description 
omega = 40; % Angular speed in rotations per second 
length = 200; % Length of the blades in cm 
width = 10; % Cross-section width in cm. (effective projection visible from the camera. assumed constant.) 
numBlades = 4; % number of blade pairs - assume even number of blades 
discWidth = 0; 
%% Camera description 
frameSize = 600*[1 3/4]; % Width & height of the picture frame (in cm) at the plane of the propeller (assuming the propeller is at the center of the image) 
sensorSize = 400*[1 3/4]; % Sensor resolution in pixels 
frameDuration = 1/15; % Time (in seconds) taken to scan all pixels in the sensor 
%% Initializations 
initAngle = 0; % Initial orientation of the propeller. Can be assumed to be anything without loss of generality. 
numPixels = prod(sensorSize); % total number of pixels 
tArrScan = [1:numPixels]'*frameDuration/numPixels; % Time instants at which different pixels are sampled. Starting top-left. 
    % posPixels: position of the pixels. matrix of dimension numPixels x 2. 1 row for each pixel. 2 elements 
    % per pixel for X & Y coordinates. 
posPixels = -0.5+[reshape(repmat([sensorSize(2):-1:1], sensorSize(1), 1), numPixels, 1)-sensorSize(2)/2 mod([0:numPixels-1]', sensorSize(1))+1-sensorSize(1)/2]; 
posPixels = posPixels*frameSize(1)/sensorSize(1); % scaling to translate the pixels on plane of the propeller. 
distPixels = abs(posPixels * [1 i]'); % distance of the pixels from the center of the frame. 
angPixels = angle(posPixels * [i 1]'); % distance of the pixels from the center of the frame. 
%% Camera operation 
img = zeros(numPixels,1); % initialization 
for kk=0:numBlades-1 % do for each balde 
    angleBlade = (2*pi/numBlades*kk)+initAngle+2*pi*mod(tArrScan*omega, 1); % position of the blade when the pixels are being scanned. 
    distPixel2Blade = sum(posPixels .* [-cos(angleBlade) sin(angleBlade)], 2); % distance of the pixels (when they are scanned) from the axis of the blade 
    distPixel2PerpOrigin = sum(posPixels .* [sin(angleBlade) cos(angleBlade)], 2); % distance of the pixels from the axis perpendicular to the blade at origin 
    distPixel2PerpOutEnd = distPixel2PerpOrigin - length; % distance of the pixels from the axis perpendicular to the blade at the ouward end of the blade 
    pixelIsBesideTheBlade = (distPixel2PerpOrigin .* distPixel2PerpOutEnd) <= 0; % Pixel is within the two ends of the blade 
    img = img+ min(abs(distPixel2Blade)<width, pixelIsBesideTheBlade); % add 1 to the pixels that are: (1) within distance width/2 from the line of axis of the blade, and (2) within the two ends of the blade. 
end 
% figure; imagesc(min([reshape(img, sensorSize(1), sensorSize(2))]',1)); colormap(gray); axis equal; % convert the pixels to 2-D image 
%% Create a video to show the results progressively 
tmpimg = zeros(numPixels,1); % initialization 
figure; imagesc([reshape(tmpimg, sensorSize(1), sensorSize(2))]'); colormap(gray); axis equal; 
clf; 
for ii=1:round(numPixels/100):numPixels % only sample 100 frames 
    tmpimg = zeros(numPixels,1); 
    tmpimg(1:ii) = min(1,img(1:ii)); 
    scnrow = ceil(ii/sensorSize(1)); 
    scncol = mod(ii-1,sensorSize(1))+1; 
    scnimg = zeros(numPixels,1); 
    scnimg((scnrow-1)*sensorSize(1)+[1:sensorSize(1)]) = 1; % mark the row 
%     scnimg(scncol:sensorSize(1):end) = 1; % mark the column 
    tmpimg(ii) = 1; 
    bldimg = zeros(numPixels,1); 
    for kk=0:numBlades-1 % do for each balde 
        angleBladeTmp = (2*pi/numBlades*kk)+initAngle+2*pi*mod(tArrScan(ii)*omega, 1); % position of the blade when the pixels are being scanned. 
        distPixel2BladeTmp = posPixels * [-cos(angleBladeTmp) sin(angleBladeTmp)]'; % distance of the pixels (when they are scanned) from the axis of the blade 
        distPixel2PerpOriginTmp = posPixels * [sin(angleBladeTmp) cos(angleBladeTmp)]'; % distance of the pixels from the axis perpendicular to the blade at origin 
        distPixel2PerpOutEndTmp = distPixel2PerpOriginTmp - length; % distance of the pixels from the axis perpendicular to the blade at the ouward end of the blade 
        pixelIsBesideTheBladeTmp = (distPixel2PerpOriginTmp .* distPixel2PerpOutEndTmp) <= 0; % Pixel is within the two ends of the blade 
        bldimg = bldimg+ min(abs(distPixel2BladeTmp)<width, pixelIsBesideTheBladeTmp); % add 1 to the pixels that are: (1) within distance width/2 from the line of axis of the blade, and (2) within the two ends of the blade. 
    end 
    
    pause(1/1000); 
    imagesc(min([reshape(tmpimg+scnimg+bldimg, sensorSize(1), sensorSize(2))]',1)); colormap(gray); axis equal; 
%     imagesc([reshape(max(max(2*tmpimg,scnimg),bldimg), sensorSize(1), sensorSize(2))]'); colormap(gray); axis equal; 
end 
imagesc(min([reshape(img, sensorSize(1), sensorSize(2))]',1)); colormap(gray); axis equal;