function StarfieldSimulation % Source code that produces a GIF. % Inspired by the "Starfield Simulation" Screensaver for Windows. % % 2017-04-20 Jahobr (update 2019-02-07 Jahobr)  [pathstr,fname] = fileparts(which(mfilename)); % save files under the same name and at file location   uniqueDepth = 200; nFrames = 240; renderDist = 200; % should be smaller or equal then uniqueDepth StarDensity = 0.0005; imageMax = [400, 500]; % pixel (height width) subsampPix = 10; % subsampling per pixel in x and y starRenderLimInPixel = 0.4; % diameter at least a x% of a pixel % if too small it flickers  for currentCase = 1:2     rng(1); % int rand generator to make GIF repeatable; If you dont like the GIF try a differnt value     speed = uniqueDepth/nFrames;     nStars = StarDensity*(2*renderDist * 2*renderDist * uniqueDepth);          x0 = rand(1,nStars)*renderDist*2-renderDist;     y0 = rand(1,nStars)*renderDist*2-renderDist;     sizze = 60+rand(1,nStars)*60;     sizze = [sizze sizze]; % copy space behind each other to loop the gif          depth = rand(1,nStars)*uniqueDepth;     depth = sort(depth); % makes later rendering more efficient     depth = [depth depth+uniqueDepth]; % copy space behind each other to loop the gif          switch currentCase         case 1 % classic             saveName = fname;             x_y =     [x0 x0;y0 y0]; % copy space behind each other to loop the gif             thetaList = zeros(1,nFrames); % rotation         case 2 % rot 90° per gif loop             saveName = [fname '_rot'];             x_y =     [x0 y0;y0 -x0]; % copy space behind each other to loop the gif             thetaList = linspace(0,pi/2,nFrames+1); % rotation             thetaList = thetaList(1:nFrames);     end      rawImage.cdata = uint8(zeros([imageMax 3]));     rawImage.colormap = [];     tic     for iFrame = 1:nFrames                           disp([num2str(iFrame) ' of ' num2str(nFrames)])         image = zeros(imageMax);                  %% calculate perspective                           theta = thetaList(iFrame);         R = [cos(theta) -sin(theta);             sin(theta) cos(theta)];         x_y_rot = R * x_y; % rotate                  plotImCoordX = (x_y_rot(1,:)./depth*imageMax(2))+imageMax(2)/2; % plane of view at depth 1;  focal point at [0 0 0]         plotImCoordY = (x_y_rot(2,:)./depth*imageMax(1))+imageMax(1)/2; % plane of view at depth 1;  focal point at [0 0 0]                  distEuclid = sqrt( x_y_rot(1,:).^2 + x_y_rot(2,:).^2 + depth.^2 );%  euclidean distance focal point at [0 0 0]                  %% Rendering         % It could be done with figure and plot functions but they are not         % designed to display objects on the pixel and sub-pixel scale.         % The result just looks bad.                  % So a own small rendering engine is used:         % check if a particular star rendering is necessary         % lay a grid over the star         % find if grid-points are within the star circle         % if yes, match the grid-point to a pixel and add its portion of         % "brightness" to the image                  startIndex = find(depth>0,1,'first'); % depth was sorted         stopIndex = find(depth<renderDist,1,'last'); % depth was sorted                  for k = startIndex:stopIndex % for all stars that would be in range                          radiusStar = sizze(k)/distEuclid(k)/2;                          if plotImCoordX(k)+radiusStar < 0  % outside to the left                 continue % next star             end             if plotImCoordX(k)-radiusStar > imageMax(2)+1  % outside to the right                 continue % next star             end                          if plotImCoordY(k)+radiusStar < 0  % outside to the bottom                 continue % next star             end             if plotImCoordY(k)-radiusStar > imageMax(1)+1  % outside to the top                 continue % next star             end                          if radiusStar*2>starRenderLimInPixel % diameter at least a defined part of a pixel                 window = ceil(radiusStar *subsampPix); % size of grid to map over stars                                  inCircle = false(2*window,2*window);                                  for ix = 1:window % grid-points in positive x                     for iy = 1:window % grid-points in positive y                                                  distEuclidPixel = sqrt( ((ix-0.5)/subsampPix).^2 + ((iy-0.5)/subsampPix).^2);%  euclidian distance focal point at [0 0 0]                                                  % do a quarter of the circle and flip to -x,y  and -x,-y  and  x,-y (1/8 would be enough actually)                         if  radiusStar>distEuclidPixel % is grid in reach of star                             inCircle(window+ix  ,window+iy  ) = true; % 1 quadrant                             inCircle(window+ix  ,window+1-iy) = true; % 2 quadrant                             inCircle(window+1-ix,window+iy  ) = true; % 3 quadrant                             inCircle(window+1-ix,window+1-iy) = true; % 4 quadrant                                                      end                     end                 end                 xGrid = round((( (1:2*window)-0.5-window) /subsampPix)+plotImCoordX(k)); % create x-grid and round to nearest Image pixel                 yGrid = round((( (1:2*window)-0.5-window) /subsampPix)+plotImCoordY(k)); % create y-grid and round to nearest Image pixel                                  for ix = 1:2*window % pixel grid in x                     if and(xGrid(ix) >= 1, xGrid(ix) <= imageMax(2)) % if in image                         for iy = 1:2*window % pixel grid in y                             if and(yGrid(iy) >= 1, yGrid(iy) <= imageMax(1)) % if in image                                 if inCircle(ix,iy) % if in Circle                                     image(yGrid(iy),xGrid(ix)) = image(yGrid(iy),xGrid(ix))+1/subsampPix^2; % add brightness to pixel                                 end                             end                         end                     end                 end             end         end                  image(image>1) = 1; % fix "over exposure"                  depth = depth-speed; % move                  %% save animation         f = rawImage;         f.cdata(:,:,1) = uint8(image*255);         f.cdata(:,:,2) = uint8(image*255);         f.cdata(:,:,3) = uint8(image*255);                  map = gray(8); % self created map. Otherwise use; [im,map] = rgb2ind(f.cdata,4,'nodither')                  imtemp = rgb2ind(f.cdata,map,'nodither');         im(:,:,1,iFrame) = imtemp;         if iFrame == 1             im(1,1,1,nFrames) = 0; % allocate         end              end          toc     imwrite(im,map,fullfile(pathstr, [saveName '.gif']),'DelayTime',1/30,'LoopCount',inf) %     disp([saveName '.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit end rng('shuffle'); % reset rand