- /
-
Trapped particles++
on 1 Nov 2024
- 58
- 374
- 0
- 8
- 1973
Cite your audio source here (if applicable): Sound effects from Pixabay.com
drawframe(1);
Write your drawframe function below
% This demo features a simple particle system, where particles are emanated
% from a point source located at the origin and follow trajectories
% described by Brownian motion while being constrained to move inside the
% 3D unit ball. Each particle is assigned its own dynamic properties and
% deforms and changes its color according to its position and lifetime.
% When a particle collides with the unit sphere, a "lightning" that
% connects that particle with the origin is drawn.
% Additionally, the plot in the current axes is captured and
% post-processing effects are applied to it (chromatic aberration, blooming).
function drawframe(f)
persistent O
% Global parameters to fine-tune visual output:
n=75; % number of particles
F=96; % number of frames
col0=[0.2 0.6 1]; % color of particles (at the origin)
col1=[1 0.6 0.2]; % color of particles (on the unit sphere)
step=0.1; % distance travelled by particle at each time step
fov=120; % camera field-of-view
outsize=[320 320]; % size of output frames
% Rename frequently occurring functions with long names to save characters:
tpf=@transformPointsForward;
r=@rescale;
g=@imgaussfilt;
m=@(a,b)mod(a-1,b)+1;
% Do not recalculate animation frames that have been already created!
if size(O,4)==F
imshow(O(:,:,:,f),Border="tight");
return
end
% Start by rendering animation from frame 41 to have a better thumbnail:
f=m(f+40,F);
% Precalculate data for animation:
% Imploding factor (used to make animation loop smoothly):
rho=g([ones(F-20,1);zeros(20,1)],10,FilterSize=[61 1]);
% Calculate coordinates of the particles at each frame by applying
% Brownian motion:
rng twister;
P=randn(F,n,3)*step;
P=imgaussfilt3(P,4,FilterSize=[255 1 1],Padding="circular");
P=rho.*cumsum(P);
% Constrain the particles to lie inside the unit ball:
R=vecnorm(P,2,3).*ones(1,1,3);
out=R>1;
P(out)=P(out)./R(out);
% Generate transformation matrices to make particles wobble:
M3=0.2*rescale(rand(1,1,n),1/2,2).*(eye(3)+(randn(3,3,n)*0.2 .*blkdiag(ones(2),0)) );
M3=imgaussfilt3(M3,[1e-2 1e-2 1],'Padding','circular');
M=affinetform3d;
for i=1:n
M(i)=affinetform3d(blkdiag(M3(:,:,i),1));
end
% Generate the texture image for the particles:
tex=r(fspecial('gaussian',32,8/6))+0.3*r(fspecial('gaussian',32,32/6));
tex=imadjust(tex,[0 0.6]);
% Set figure background to black:
figure(1);
set(gcf,'Pos',[0 0 outsize],'Color',col0/10);
% Remove axes and set camera properties:
cla;
axis equal off;
axis([-1 1 -1 1 -1 1]);
campos([0 0 cscd(fov/2)]);
camup([0 1 0]);
camva(fov);
camproj perspective;
set(gca,YDir="normal",Units="normalized",Position=[0 0 1 1],Clipping="off");
% For each particle that touches the unit sphere,draw a "lightning"
% emanating from the origin and reaching that particle:
rng shuffle twister;
q=0.02*randn(20,n,3);
q([1 end],:,:)=0;
l=linspace(0,1,20)'.*P(f,:,:)+q;
l(:,~out(f,:,1),:)=[];
line(l(:,:,1),l(:,:,2),l(:,:,3),'LineW',1.5,'Color',[col1 0.2]);
% In object space,each particle is represented as one quad with its
% vertices sitting on the corners of the unit square of the xy-plane:
[X0,Y0,Z0]=meshgrid([-1 1],[-1 1],0);
% Get the xyz-coordinates of the unit sphere:
[U0,V0,W0]=sphere(64);
% Render settings for particles:
preset={'EdgeC','none','FaceL','none','FaceA','texturemap','AlphaDataM','none'};
% Loop to draw particles:
for k=1:n
% Take the next matrix and apply the corresponding distortion to the unit square:
q=m(k+f,n);
[X,Y,Z]=tpf(M(q),X0,Y0,Z0);
% Assign color to each particle such that particles that are
% close to the origin are more blue and those that are closer to the
% unit sphere are reddish.
col=interp1([0;1],[col0;col1],min(R(f,k),1).^3 );
% Assign alpha value to each particle such that particles that are
% far from the observer tend to fade away:
alpha=r(P(f,k,3),0.1,1,InputMin=-1,InputMax=1);
% Draw the current particle:
surface(X+P(f,k,1),Y+P(f,k,2),Z+P(f,k,3),'FaceC',col,'AlphaData',tex*alpha,preset{:});
end
% Slightly rotate the unit sphere:
RT=rigidtform3d(5*sin(2*pi*f/F.*[1 1 1]+[0 pi/3 pi/4]),[0 0 0]);
[U,V,W]=tpf(RT,U0,V0,W0);
surface(U,W,V,'FaceC','none','EdgeC',col0,'EdgeA',0.05,'LineW',rand*5);
% Post-Processing part...
A=im2double(frame2im(getframe(gca)));
% Simulate chromatic aberration by zooming each color channels by a
% slightly different factor:
[U,V]=meshgrid(linspace(-1,1,width(A)),linspace(-1,1,height(A)) );
D=-cat(3,U,V)*6;
for c=1:3
A(:,:,c)=imwarp(A(:,:,c),D*(c-1),'linear');
end
% Apply monochromatic Gaussian noise (only for uncompressed video):
% A=max(min(A+randn(size(A,1:2))*3/255,1),0);
% Apply "blooming" to the image:
B=rgb2gray(A);
B=imadjust(B,[0.1 1]);
B=g(B,32,FilterSize=193,Padding="replicate");
B=r(B);
A=A+r(B);
% Restore original frame number:
f=m(f-40,F);
% Display post-processed image and store it in a persistent variable:
imshow(A,Border="tight");
O(:,:,:,f)=A;
end