10um linescan with 100x objective (modified from https://github.com/omsai/microscope-optical-input)

raster.ipynb

raster.R

## Calculate power distribution of laser scanning photo-activation device
##
## Laser beam symbols used are from the beam power formula:
## [1] en.wikipedia.org/wiki/Gaussian_beam#Power_through_an_aperture
##
## References:
## [2] radiantzemax.com/kb-en/Goto50125.aspx
## [4] en.wikipedia.org/wiki/Beam_diameter#D4.CF.83_or_second_moment_width
## [5] en.wikipedia.org/wiki/Normal_distribution
##
## Theory contributors in alphabetical order:
## Browne, M. (Andor Technology)
## Karunarathne, A. (Washington University in St. Louis)
## Magidson, V. (Wadsworth Institute, NYSDOH)
##
## Code Author:
## Pariksheet Nanda (Andor Technology)

library(spatstat) # calulate image convolution
library(ggplot2)  # display complex plots
library(reshape2) # reshaping matrix as data.frame for ggplot

## Create the round laser beam profile at the specimen plane

# Gaussian laser beam power (micro-Watts) at back focal plane, Symbol: Ps
# Measured by user's PD-300 laser meter using single pixel FRAP
Ps <- 4.8 / 40

To <- 0.84 # Percent transmission through above objective at 445nm 
Tm <- 0.8 # Percent transmission through media, etc (estimated)
# Gaussian laser beam power (micro-Watts), Symbol: Po
Po <- Ps * To * Tm

# FWHM Beam diameter / waist size (microns) at objective, Symbol: Wo
# Measured at objective specimen plane using FRAPPA slide
Wo_FWHM <- 0.8

# Convert FWHM to 1/e^2 [2]
Wo <- Wo_FWHM * 0.8493218 * 2

camera_pixels <- 512 # e2v CCD97 pixel
fov_width <- 76.6 # microns
# Pixel calibration (microns / pixel)
cal <- fov_width / camera_pixels

# Beam axis pixel length on camera (pixels), Symbol: px_len
# Wo for a single mode beam is 4 sigma wide [4]
px_len_2sig <- Wo / cal
# Extend the beam diameter out to 6 sigma to contain 99% of Po
px_len <- (6 / 4.) * px_len_2sig

px_fit <- floor(px_len) / px_len
px_to_edge <- floor(3 / px_fit)
edge <- px_fit * px_to_edge
steps <- floor(px_len / px_fit)
#x <- seq(-edge, edge, length=steps)  # keep for debugging 1D construction

# Map 1D arrays radially
x <- seq(-edge, edge, length=steps)
mx <- matrix(rep(x, length(x)), length(x), length(x)) # square matrix of x
my <- t(mx)
z <- sqrt(mx^2 + my^2)

# Create 2D gaussian beam profile
norm <- function(x, mean=0, sigma=1) {
    # Returns single value from Gaussian / Normal distribution curve [5]
    (1/(sigma * sqrt(2 * pi))) * exp(-0.5 * ((x - mean) / sigma) ^ 2)
}

beam_profile <- norm(z)

# norm() gives a 1D probability distribution, the area of which is 1.  The 2D
# beam profile thus needs to be re-normalized back to 1, for it to be
# multiplied by the Po scalar, obj_mag^2, and transmission losses
beam_profile <- beam_profile / sum(beam_profile)
beam_profile <- beam_profile * Po

## Setup ROI

w_microns <- 10
w <- round(w_microns / cal) # pixels
h <- 1 # pixels

length <- ncol(beam_profile)  # length in pixels of beam_profile square
x <- w + length - length %% 2
y <- h + length - length %% 2
roi <- matrix(1, x, y)

## Apply the laser beam over the ROI

roi <- convolve.im(as.im(roi), as.im(beam_profile))

## Show results

# Area calculation
w_um <- w * cal
h_um <- h * cal
x_um <- x * cal
y_um <- y * cal
actual_area <- x_um * y_um
drawn_area <- w_um * h_um

# Power calculation
peak_pixel_power <- max(roi) # uW
actual_power <- sum(roi) # uW
drawn_power <- sum(roi[length/2:w+length/2,length/2:h+length/2])
density_drawn_power <- drawn_power / (drawn_area / 1000**2) # uW/mm^2

# Energy calculation
dwell_time <- 80  # micro-seconds
repeats <- 1
peak_pixel_energy <- max(roi) * dwell_time * repeats / 1000 # uJ
actual_energy <- actual_power * dwell_time * repeats / 1000 # uJ
drawn_energy <- drawn_power * dwell_time * repeats / 1000 # uJ
density_drawn_energy <- drawn_energy / (drawn_area / 1000**2) # uJ/mm^2

sprintf('Laser Beam at Back Focal Plane-')
sprintf('Measured Power: %3.0f uW', Ps)
sprintf('')
sprintf('Laser Beam at Objective Specimen Plane-')
sprintf('Calculated Power: %3.0f uW', Po)
sprintf('Measured Width: %1.1f um, in FWHM', Wo_FWHM)
sprintf('')
sprintf('Beam Profile Pixelation at Back Focal Plane-')
sprintf('Calibration: %1.3f um/pixel', cal)
sprintf('Size: %d x %d pixels, covering 3 standard deviations',
        nrow(beam_profile), ncol(beam_profile))

m <- as.matrix(beam_profile)
mt <- melt(t(m))
names(mt)[names(mt)=="Var1"] <- "y"
names(mt)[names(mt)=="Var2"] <- "x"
mt$x <- (mt$x - mean(mt$x)) * cal
mt$y <- (mt$y - mean(mt$y)) * cal
ggplot(mt, aes(x=x, y=y, fill=value)) +
  scale_fill_distiller(palette="Spectral") +
  geom_raster(hjust=1, vjust=1) +
  coord_fixed(expand=FALSE) +
  labs(list(
    title=expression(paste('Beam profile power of laser (', mu, 'W)')),
    x=expression(paste(mu, 'm')),
    y=expression(paste(mu, 'm'))
  )) +
  geom_blank()

# 2D plot of laser at specimen plane
m <- as.matrix(roi)
mt <- melt(t(m))
names(mt)[names(mt)=="Var1"] <- "y"
names(mt)[names(mt)=="Var2"] <- "x"
mt$x <- (mt$x - mean(mt$x)) * cal
mt$y <- (mt$y - mean(mt$y)) * cal
ggplot(mt, aes(x=x, y=y, fill=value)) +
  scale_fill_distiller(palette="Spectral") +
  geom_raster(hjust=1,vjust=1) +
  coord_fixed(expand=FALSE) +
  labs(list(
    title=expression(paste('Laser power at specimen plane (', mu, 'W)')),
    x=expression(paste('Specimen plane (', mu, 'm)')),
    y=expression(paste('Specimen plane (', mu, 'm)'))
    )) +
  geom_blank() +
  annotate("rect", color="black", alpha=0,
           xmin=-w_um/2, xmax=w_um/2,
           ymin=-h_um/2, ymax=h_um/2)

# 1D cross section of laser at specimen plane
m <- t(as.matrix(roi))
y <- m[nrow(m)/2,]
x <- 1:length(y)
x <- (x - mean(x)) * cal
line <- data.frame(x, y)
ggplot(line, aes(x=x, y=y)) +
  geom_line(color="red") +
  geom_vline(xintercept=c(-w_um/2,+w_um/2)) +
  labs(list(
    title="Cross section of laser spread at specimen plane",
    x=expression(paste('ROI length (', mu, 'm)')),
    y=expression(paste('Laser power (', mu, 'W)'))
  ))