Back to Gallery
  • Remix
  • Share
  • New Entry
on 20 Oct 2024
  • 47
  • 301
  • 0
  • 5
  • 1828
Cite your audio source here (if applicable): Dark Water / Magnus Ludvigsson
drawframe(1);
Write your drawframe function below
% This is an example of how to combine ray tracing and sphere intersection
% with Snell's law for refraction of light.
% Light rays are traced from the location of the light source through the
% surface of the water, and they hit a virtual sphere inside the water.
% No sphere is actually
% drawn; the rays are projected to the surface of the virtual sphere and
% the density of the beams is used to create the caustics. Rays from the
% camera are then drawn to the surface of the water, and these are also
% refracted onto the surface using Snell's law. The density at the
% intersection of the camera-rays with the surface is used to determine (in
% part) the shading of the triangle patch intersecting with teh ray from
% the camera. It is kind of half-way between a standard matlab scene render
% and a true shader in that respect; it is using the patch faces of the
% water surface as the image pixels.
function drawframe(f)
persistent cps ag s po T2 Lt k
if f==1
n=2e5;
N=0:n-1;
t=4*pi*N/(1+5^.5);
p=acos(1-(2*N+1)/n);
po=[cos(t).*sin(p);sin(t).*sin(p);cos(p)]';
po = po + randn(size(po))/3e3;
rng(2,'twister');
po = po./vecnorm(po')';
rd=randn(n,1)/40+1;
% Smoothing around sphere
vn=@(x)(x./vecnorm(x')');
k=convhull(po); % Points on "in" must lie on unit circle
c=@(x)sparse(k(:,x)*[1,1,1],k,1,n,n); % Connectivity
t=c(1)|c(2)|c(3); % Connectivity
ff=spdiags(-sum(t,2)+1,0,t*1.)*35; % Weighting
s=((speye(n)+ff'*ff)\rd); % Solve for s w/regularizer
% Ok... now we want something that can move periodically. So let's convert
% our displacements to phases. Then we can modulate by simply looping over
% 2*pi
ag=linspace(0,2*pi,97);
% Our surface
T2=trisurf(k,p(:,1),p(:,2),p(:,3),'EdgeC','none','AmbientS',1,'DiffuseS',0); % Plot
% Camera & light controls
Lt=[1,-1,1]*1e3; % Light position
cps=[-1,-.4,13];
[x,y,z]=sph2cart(cps(1),cps(2),cps(3));
CP=[x,y,z];
axis equal
camproj p
campos(CP);
camtarget([0,0,0]);
axis off
set(gcf, 'color', 'k');
G=gcf;
G.Position(3:4)=[600,600];
camva(8);
light
end
n=f;
[x,y,z]=sph2cart(cps(1)+ag(n)*0,cps(2),cps(3));
CP=[x,y,z];
o=1500;
s2 = sin((circshift(s, 5)-1)*o+2*ag(n))/200+sin((circshift(s, 11)-1)*o+3*ag(n))/200 + sin((s-1)*o+1*ag(n))/200+1;
p2 = po.*s2;
% Face centers...
fc=squeeze(mean(reshape([p2(k,1);p2(k,2);p2(k,3)],[size(k),3]), 2));
% Face center normals
fcn=fc./vecnorm(fc')';
% Triangulate: this is to get the face normal vector. Eric L. would have a
% better idea of how to do this, I'm sure.
T=triangulation(k,p2(:,1),p2(:,2),p2(:,3));
F=-faceNormal(T);
% Get surface coloration using Snell's law & ray tracing
cf=tot(CP,Lt,F,fc,k);
% Update face surface color
T2.Vertices=p2;
T2.FaceVertexCData = cf;
campos(CP);
drawnow;
frms{n} = getframe(gcf);
end
% Calculate sphere hit locations using Snell's law and ray-sphere
% intersection equations
function [h,I,d] = SNL(F,LF,fc)
% u=.94; % Refraction idx (recip).
u = 0.8;
RD=.6; % Planet radius
% Snell's law, light going into ocean
I=dot(F,LF,2);
t=sqrt(1-min(u.^2*(1-I.^2),1)).*F+u*(LF-I.*F);
t=t./vecnorm(t')';
% Sphere intersection
b=2*sum(t.*fc,2);
c=sum(fc.^2,2)-RD.^2;
D2=(-b-sqrt(b.^2-4*c))/2; % Intersection distance
h=fc+D2.*t; % Intersection point
d=imag(h(:,1))==0&I>0; % Hit index test
h=h(d,:); % hit locations
I=I(d); % ray angle cosines
end
% Calculation of contributions at sphere triangle faces
function v=tot(CP,Lt,F,fc,k)
vn=@(x)(x./vecnorm(x')'); % Vector normalizer
% Raytrace from light to sphere
LF=vn(fc-Lt); % Face to light vectors
[h,I]=SNL(F,LF,fc); % Tray trace to planet surface
% Ray trace from camera to sphere
CF=vn(fc-CP); % Camera -> ocean surface vectors
[C,~,d]=SNL(F,CF,fc);
% Want a color-fade to indicate ambient scattering... get projection in
% direction half way between light and camera
Ltn=CP/norm(CP)+Lt/norm(Lt);
Ltn=Ltn/norm(Ltn);
wt=erf((sum(h.*Ltn,2)-.4)*2)/2+.5; % Taper function between color schemes using erf
% Search for 100 closest points
[Di,D]=knnsearch(real(h),C,'K',99); % 99... save a character
% Make color vectors
g=0*k(:,1)';
v=g;
g(d)=rescale(mean((wt(Di)).^2,2)); % Get fade weighting
h=1-g; % Reverse weighting
v(d)=rescale(sum(I(Di).^2.*exp(-2e3*(D.^2)),2),0.3,1.2); % Caustic intensities
f=sum(F.*CF,2);
v = max(max(v,f'*max(v)/2),sin((acos(sum(F.*CF, 2))))'.^10*.9);
c1=@(x)[x.^8,x.^2.5,x.^0.9];
c2=@(x)[x,x.^1.1,x.^1.8];
c1c=c1(v(:));
c2c=c2(v(:));
h(~d)=1;
v=c1c.*h(:)+c2c.*g(:);
end
Movie
Audio
Remix Tree