Add files via upload

Matlab App v5/function_circular_aperture_cutoff_heisenberg_diffraction.m

@@ -0,0 +1,89 @@
+function [lightfield,del_theta_x,del_theta_y] = function_circular_aperture_cutoff(lightfield, circular_aperture_diameter, wavelength, diffraction_type, low_lim,upp_lim)
+% Function to apply a circular aperture cutoff to a lightfield
+% INPUTS:
+%   - lightfield: Structure containing ray positions and directions
+%   - circular_aperture_diameter: Diameter of the circular aperture
+% OUTPUT:
+%   - new_lightfield: Updated lightfield structure with rays outside the aperture set to NaN
+
+% Start timer for performance measurement
+tic
+
+% Create a logical mask that identifies rays outside the aperture radius
+% Radius is half of the aperture diameter
+mask = sqrt(lightfield.position.x.^2 + lightfield.position.y.^2) > (circular_aperture_diameter / 2);
+
+% % Copy the original lightfield to preserve unchanged rays
+% lightfield = lightfield;
+
+% Apply the mask to all position and direction components
+% Rays outside the aperture are set to NaN in the new lightfield structure
+lightfield.position.x(mask) = nan;
+lightfield.position.y(mask) = nan;
+lightfield.position.z(mask) = nan;
+lightfield.direction.dx(mask) = nan;
+lightfield.direction.dy(mask) = nan;
+lightfield.direction.dz(mask) = nan;
+
+lightfield.position.x = lightfield.position.x(~isnan(lightfield.position.x));
+lightfield.position.y = lightfield.position.y(~isnan(lightfield.position.y));
+lightfield.position.z = lightfield.position.z(~isnan(lightfield.position.z));
+lightfield.direction.dx = lightfield.direction.dx(~isnan(lightfield.direction.dx));
+lightfield.direction.dy = lightfield.direction.dy(~isnan(lightfield.direction.dy));
+lightfield.direction.dz = lightfield.direction.dz(~isnan(lightfield.direction.dz));
+
+if diffraction_type ~= "none"
+    %Calculate distance to each aperture edge
+    aperture_radius = circular_aperture_diameter/2;
+    x_distance = sqrt(aperture_radius^2-lightfield.position.y.^2)-abs(lightfield.position.x); %dimensionalized distance to edge
+    y_distance = sqrt(aperture_radius^2-lightfield.position.x.^2)-abs(lightfield.position.y); %dimensionalized distance to edge
+    del_theta_x = zeros(size(lightfield.position.x)); %Initialize vector for change in theta_x due to diffraction
+    del_theta_y = zeros(size(lightfield.position.y)); %Initialize vector for change in theta_y due to diffraction
+    k = 2*pi/wavelength;
+
+    %Calculate norm distribution characteristics
+    sigma_x = atan(1./(2*x_distance*k));
+    sigma_y = atan(1./(2*y_distance*k));
+
+    %Diffraction model parameters
+    if string(diffraction_type) == 'dynamic limited'
+        %Calculate angle of defelction due to diffraction
+        lower_limit_x = -sigma_x;
+        upper_limit_x = -lower_limit_x;
+        lower_limit_y = -sigma_y;
+        upper_limit_y = -lower_limit_y;
+        %Modeling diffraction in the Monte-Carlo ray-trace environment,
+        %Coffey 1998
+    elseif string(diffraction_type) == "-pi/2 to pi/2 limited"
+        lower_limit_x = -pi/2*ones(size(sigma_x));
+        upper_limit_x = -lower_limit_x;
+        lower_limit_y = -pi/2*ones(size(sigma_y));
+        upper_limit_y = -lower_limit_y;
+    else
+        lower_limit_x = low_lim*ones(size(sigma_x));
+        upper_limit_x = upp_lim*ones(size(sigma_x));
+        lower_limit_y = lower_limit_x;
+        upper_limit_y = upper_limit_x;
+    end %End of type of diffraction
+
+    %Calculate angle due to diffraction
+    parfor counter = 1:numel(lightfield.position.x) % For all rays
+        %Generate values within the truncated range
+        del_theta_x(counter) = truncate_gaussian(0, sigma_x(counter), lower_limit_x(counter), upper_limit_x(counter));
+        del_theta_y(counter) = truncate_gaussian(0, sigma_y(counter), lower_limit_y(counter), upper_limit_y(counter));
+    end
+    %Edge diffraction in Monte Carlo ray tracing, Freniere et al 1999
+    % DOI:10.1117/12.363773
+
+
+
+    %Update new directions in the x and y component by converting from
+    %spherical to cartesian
+
+    lightfield.direction.dx = lightfield.direction.dx + tan(del_theta_x);
+    lightfield.direction.dy = lightfield.direction.dy + tan(del_theta_y);
+    %Display completion message and stop timer
+    fprintf("Lightfield cutoff.");
+end %End of diffraction check
+toc
+end

Matlab App v5/function_convert_rho_field_to_n_field.m

@@ -0,0 +1,11 @@
+function n_field = function_convert_rho_field_to_n_field(rho,gladstone_dale)
+%Convert rho to n field
+
+%% Load data
+
+load("sample-density.mat");
+rho = rho(:,:,1);
+
+n = gladstone_dale*rho+1;
+end
+

Matlab App v5/function_generate_image_2.m

@@ -0,0 +1,72 @@
+function [rendered_image] = function_generate_image_2(lightfield, x_pixel_number, y_pixel_number, pixel_pitch)
+
+% Remove invalid rays
+invalid = isnan(lightfield.position.x) | isnan(lightfield.position.y);
+lightfield.direction.dx(invalid) = [];
+lightfield.direction.dy(invalid) = [];
+lightfield.position.x(invalid) = [];
+lightfield.position.y(invalid) = [];
+
+tic
+rendered_image = zeros(y_pixel_number, x_pixel_number); %rows by columns
+
+% Compute subpixel locations relative to center of sensor
+x_center = x_pixel_number / 2;
+y_center = y_pixel_number / 2;
+
+% Convert positions to pixel space
+x_pixel_f = (lightfield.position.x / pixel_pitch) + x_center + 1;
+y_pixel_f = (lightfield.position.y / pixel_pitch) + y_center + 1;
+
+% Compute
+alpha = atan(sqrt((lightfield.direction.dx./lightfield.direction.dz).^2+(lightfield.direction.dy./lightfield.direction.dz).^2));
+intensity = cos(alpha).^4;
+
+% Get surrounding pixel indices
+x0 = floor(x_pixel_f);
+y0 = floor(y_pixel_f);
+x1 = x0 + 1;
+y1 = y0 + 1;
+
+% Calculate weights for bilinear interpolation
+wx = x_pixel_f - x0;
+wy = y_pixel_f - y0;
+w00 = (1 - wx) .* (1 - wy); % top-left
+w10 = wx .* (1 - wy);       % top-right
+w01 = (1 - wx) .* wy;       % bottom-left
+w11 = wx .* wy;             % bottom-right
+
+% Filter valid rays that hit the sensor
+valid = x0 >= 1 & x1 <= x_pixel_number & y0 >= 1 & y1 <= y_pixel_number;
+x0 = x0(valid); x1 = x1(valid);
+y0 = y0(valid); y1 = y1(valid);
+w00 = w00(valid); w10 = w10(valid); w01 = w01(valid); w11 = w11(valid);
+intensity = intensity(valid);
+
+%Reshape vector
+pixel_indices = [
+    sub2ind([y_pixel_number, x_pixel_number], y0, x0),  % top-left
+    sub2ind([y_pixel_number, x_pixel_number], y0, x1),  % top-right
+    sub2ind([y_pixel_number, x_pixel_number], y1, x0),  % bottom-left
+    sub2ind([y_pixel_number, x_pixel_number], y1, x1)   % bottom-right
+];
+weights = [
+    intensity .* w00;
+    intensity .* w10;
+    intensity .* w01;
+    intensity .* w11
+];
+
+% Accumulate all contributions using accumarray
+accumulated = accumarray(pixel_indices(:), weights(:), [y_pixel_number * x_pixel_number, 1]);
+rendered_image = reshape(accumulated, [y_pixel_number, x_pixel_number]);
+
+% Flip vertically to match image coordinate convention
+rendered_image = rendered_image(end:-1:1, end:-1:1);
+
+
+fprintf("Image rendered with bilinear interpolation. ");
+toc
+
+figure; imshow(rendered_image(2:end-1,2:end-1), []);
+end

Matlab App v5/function_generate_image_with_splat.m

@@ -0,0 +1,59 @@
+function [] = function_generate_image_with_splat(lightfield, x_pixel_number, y_pixel_number, pixel_pitch)
+
+% Remove invalid rays
+valid_mask = ~isnan(lightfield.position.x) & ~isnan(lightfield.position.y);
+x = lightfield.position.x(valid_mask);
+y = lightfield.position.y(valid_mask);
+dx = lightfield.direction.dx(valid_mask);
+dy = lightfield.direction.dy(valid_mask);
+
+% Floating-point pixel positions
+x_fp = (x / pixel_pitch) + (x_pixel_number / 2);
+y_fp = (y / pixel_pitch) + (y_pixel_number / 2);
+
+% Compute intensity from direction magnitude
+intensity = max(0, min(1, 1 ./ (1 + sqrt(dx.^2 + dy.^2))));
+
+% Bilinear splatting
+x0 = floor(x_fp);
+y0 = floor(y_fp);
+dx_frac = x_fp - x0;
+dy_frac = y_fp - y0;
+
+% Clip to valid ranges (x0+1 and y0+1 must also be within bounds)
+valid = x0 >= 1 & x0 < x_pixel_number & y0 >= 1 & y0 < y_pixel_number;
+
+x0 = x0(valid);
+y0 = y0(valid);
+dx_frac = dx_frac(valid);
+dy_frac = dy_frac(valid);
+intensity = intensity(valid);
+
+rendered_image = zeros(y_pixel_number, x_pixel_number);
+
+% Accumulate contributions into 4 neighboring pixels
+for i = 1:length(x0)
+    x1 = x0(i) + 1;
+    y1 = y0(i) + 1;
+
+    % Contributions to four neighbors
+    w00 = (1 - dx_frac(i)) * (1 - dy_frac(i));
+    w01 = (1 - dx_frac(i)) * dy_frac(i);
+    w10 = dx_frac(i) * (1 - dy_frac(i));
+    w11 = dx_frac(i) * dy_frac(i);
+
+    rendered_image(y0(i),  x0(i))  = rendered_image(y0(i),  x0(i))  + intensity(i) * w00;
+    rendered_image(y1,     x0(i))  = rendered_image(y1,     x0(i))  + intensity(i) * w01;
+    rendered_image(y0(i),  x1)     = rendered_image(y0(i),  x1)     + intensity(i) * w10;
+    rendered_image(y1,     x1)     = rendered_image(y1,     x1)     + intensity(i) * w11;
+end
+
+% Flip image vertically and horizontally
+rendered_image = flipud(fliplr(rendered_image));
+
+fprintf("Image rendered. ");
+toc
+
+figure; imshow(rendered_image, []);
+
+end

Matlab App v5/function_generate_lightfield.asv

@@ -0,0 +1,100 @@
+function lightfield = function_generate_lightfield(number_of_light_rays_per_location, scale, dot_width, dot_height, spacing_x, spacing_y, num_column, num_row, z_location, cone_half_angle,use_gpu)
+
+num_location_per_dot_x = round(dot_width/scale);
+lightsource.location.x = zeros(num_column*num_location_per_dot_x);
+x_values = (scale * (1:num_location_per_dot_x)') + (dot_width + spacing_x) * (0:num_column-1);
+lightsource.location.x = reshape(x_values, [], 1);
+lightsource.location.x = lightsource.location.x-(lightsource.location.x(end)-lightsource.location.x(1))/2; %center around (0,y)
+
+num_location_per_dot_y = round(dot_height/scale);
+lightsource.location.y = zeros(num_row*num_location_per_dot_y);
+y_values = (scale * (1:num_location_per_dot_y)') + (dot_height + spacing_y) * (0:num_row-1);
+lightsource.location.y = reshape(y_values, [], 1);
+lightsource.location.y = lightsource.location.y-(lightsource.location.y(end)-lightsource.location.y(1))/2; %center around (x,0)
+
+%Used to address issue of scale being approximately equal to pixel pitch.
+%Perturb lightsource locations so gridded image doesn't appear.
+
+% lightsource.location.x = lightsource.location.x + 1*scale * (rand(size(lightsource.location.x)) - 0.5);
+% lightsource.location.y = lightsource.location.y + 1*scale * (rand(size(lightsource.location.y)) - 0.5);
+% lightsource.location.x = unifrnd(-1,1,size(lightsource.location.x))*scale+lightsource.location.x;
+% lightsource.location.y = unifrnd(-1,1,size(lightsource.location.y))*scale+lightsource.location.y;
+
+lightsource.location.z = z_location; %Initialize position
+%% Generate uniform lightfield
+
+% Collimated light should have very small divergence/diffusion angle
+cone_half_angle_radians = deg2rad(cone_half_angle);
+tic
+[temp_x_mesh,temp_y_mesh] = meshgrid(lightsource.location.x,lightsource.location.y);
+temp_x_flat = temp_x_mesh(:);
+mu = number_of_light_rays_per_location;
+for counter = 1:size(temp_x_flat,1)
+    num_lightray_array(counter) = min(max(round(normrnd(mu,(2*mu-1)/6)),1),2*mu);
+end
+total_lightray = sum(num_lightray_array);
+
+for i = 1:M
+    for j = 1:N
+        % Calculate the number of rays for the current position
+        num_rays = min(max(round(normrnd(mu, (2*mu - 1) / 6)), 1), 2 * mu);
+
+        % Assign the x, y positions to the arrays
+        x_positions = temp_x_mesh(i,j) * ones(1, num_rays);
+        y_positions = temp_y_mesh(i,j) * ones(1, num_rays);
+
+        % Store the x and y positions in the preallocated arrays
+        lightfield_position_x(current_index : current_index + num_rays - 1) = x_positions;
+        lightfield_position_y(current_index : current_index + num_rays - 1) = y_positions;
+
+        % Update the current index
+        current_index = current_index + num_rays;
+    end
+end
+
+
+% for counter = 1:number_of_light_rays_per_location
+%     lightfield.position.x(:,:,counter) = temp_x_mesh;
+%     lightfield.position.y(:,:,counter) = temp_y_mesh;
+% end
+
+
+
+
+
+
+
+clearvars temp_x_mesh temp_y_mesh;
+lightfield.position.x = lightfield.position.x(:);
+lightfield.position.y = lightfield.position.y(:);
+lightfield.position.z = lightsource.location.z*ones(size(lightfield.position.y));
+
+lightfield.direction.random_cone_angle = (rand(size(lightfield.position.x))-1/2)*cone_half_angle_radians*2;
+lightfield.direction.random_circular_angle = (rand(size(lightfield.position.x))-1/2)*pi*2;
+lightfield.direction.dx = sin(lightfield.direction.random_cone_angle).*cos(lightfield.direction.random_circular_angle);
+lightfield.direction.dy = sin(lightfield.direction.random_cone_angle).*sin(lightfield.direction.random_circular_angle);
+lightfield.direction.dz = (-1)*ones(size(lightfield.direction.dy));
+
+clear lightfield.direction.random_cone_angle lightfield.direction.random_circular_angle
+
+if use_gpu
+    lightfield.position.x = gpuArray(lightfield.position.x);
+    lightfield.position.y = gpuArray(lightfield.position.y);
+    lightfield.position.z = gpuArray(lightfield.position.z);
+    lightfield.direction.dx = gpuArray(lightfield.direction.dx);
+    lightfield.direction.dy = gpuArray(lightfield.direction.dy);
+    lightfield.direction.dz = gpuArray(lightfield.direction.dz);
+end
+
+
+
+% if show_ray_tracing_results
+    % figure;
+    % plot3(lightfield.position.x(1:37:end),lightfield.position.y(1:37:end),lightfield.position.z(1:37:end),'.');
+    % % axis equal;
+    % grid on;
+    % hold on;
+% end
+
+fprintf("Lightfield generated. Number of light rays = %d. ", size(lightfield.position.x,1));
+toc