Holography and Light Scattering in Python¶
Release:  3.0.0 

HoloPy
is a python based tool for working with digital
holograms and light scattering. HoloPy can be used to analyze holograms in two complementary ways:
 Backward propagation of light from a digital hologram to reconstruct 3D volumes.
 This approach requires no prior knowledge about the scatterer
 Forward propagation of light from a scattering calculation of a predetermined scatterer.
 Comparison to a measured hologram with Bayesian inference allows precise measurement of scatterer properties and position.
HoloPy provides a powerful and userfriendly python interface to fast scattering and optical propagation theories implemented in Fortran and C code. It also provides a set of flexible objects that make it easy to describe and analyze data from complex experiments or simulations.
HoloPy started as a project in the Manoharan Lab at Harvard University. If you use HoloPy, you may wish to cite one or more of the sources listed in References and credits. We also encourage you to sign up for our User Mailing List to keep up to date on releases, answer questions, and benefit from other users’ questions.
User Guide¶
Skip to the Loading Data tutorial if you already have HoloPy installed and want to get started quickly.
Getting Started¶
Installation¶
As of version 3.0, HoloPy supports only Python 3. We recommend using the anaconda distribution of Python, which makes it easy to install the required dependencies. In the future (hopefully by the time you read this documentation), HoloPy will be available on condaforge, so that you can install it with:
conda install c condaforge holopy
in a shell, terminal, or command prompt. Once you have HoloPy installed, open an IPython console or Jupyter Notebook and run:
import holopy
If this line works, skip to Using HoloPy before diving into the tutorials.
Windows Support¶
To find a HoloPy conda package for Windows, use
anaconda search holopy
You may need to manually install any missing dependencies (eg. emcee) from condaforge like
conda install c condaforge emcee
You can also build HoloPy from source by following the instructions for Installing HoloPy for Developers.
Dependencies¶
HoloPy’s hard dependencies can be found in requirements.txt. Optional dependencies for certain calculations include:
Using HoloPy¶
You will probably be most comfortable using HoloPy in Jupyter (resembles Mathematica) or Spyder (resembles Matlab) interfaces. One perennially tricky issue concerns matplotlib backends. HoloPy is designed to be used with an interactive backend. In the console, try running:
from holopy import test_disp
test_disp()
You should see a window pop up with an image, and you should be able to change the square to a circle or diamond by using the left/right arrow keys. If you can, then you’re all set! Check out our Loading Data tutorial to start using HoloPy. If you don’t see an image, or if the arrow keys don’t do anything, you can try setting your backend with one of the following:
%matplotlib tk
%matplotlib qt
%matplotlib gtk
%matplotlib gtk3
Note that these commands will only work in an IPython console or Jupyter
Notebook. If the one that you tried gave an ImportError
, you should restart
your kernel and try another. Note that there can only be one matplotlib backend
per ipython kernel, so you have the best chance of success if you restart your
kernel and immediately enter the %matplotlib
command before doing anything
else. Sometimes a backend will be chosen for you (that cannot be changed later)
as soon as you plot something, for example by running test_disp()
or
show()
. Trying to set to one of the above backends that is not installed
on your system will result in an error, but will also prevent you from setting a different
backend until you restart your kernel.
An additional option in Spyder is to change the backend through the menu: Tools > Preferences > IPython console > Graphics. It will not take effect until you restart your kernel, but it will then remember your backend for future sessions, which can be convenient.
An additional option in jupyter is to use %matplotlib
nbagg
to use inline interactive plots.
Loading Data¶
HoloPy can work with any kind of image data, but we use it for digital holograms, so our tutorials will focus mostly on hologram data.
Loading and viewing a hologram¶
We include a couple of example holograms with HoloPy. Lets start by loading and viewing one of them
import holopy as hp
from holopy.core.io import get_example_data_path
imagepath = get_example_data_path('image01.jpg')
raw_holo = hp.load_image(imagepath, spacing = 0.0851)
hp.show(raw_holo)
(Source code, png, hires.png, pdf)
The first few lines just specify where to look for an image. The most important line actually loads the image so that you can work with it:
raw_holo = hp.load_image(imagepath, spacing = 0.0851)
HoloPy can import any image format that can be handled by Pillow.
The spacing argument tells holopy about the scale of your image. Here, we had
previously measured that each pixel is a square with side length 0.0851 microns.
In general, you should specify spacing
as the distance between adjacent pixel centres.
You can also load an image without specifying a spacing value if you just want
to look at it, but most holopy calculations will give erroneous results on such an image.
The final line simply displays the loaded image on your screen
with the builtin HoloPy show()
function. If you don’t see anything, you may need to set your matplotlib backend.
Refer to Using HoloPy for instructions.
Correcting Noisy Images¶
The raw hologram has some nonuniform illumination and an artifact in the
upper right hand corner from dust somewhere in the optics. These types of
things can be removed if you are able to take a background image with the same optical setup but
without the object of interest. Dividing the raw hologram by the background using bg_correct()
can usually improve the image a lot.
from holopy.core.process import bg_correct
bgpath = get_example_data_path('bg01.jpg')
bg = hp.load_image(bgpath, spacing = 0.0851)
holo = bg_correct(raw_holo, bg)
hp.show(holo)
(Source code, png, hires.png, pdf)
Often, it is beneficial to record multiple background images. In this case,
we want an average image to pass into bg_correct()
as our background.
bgpath = get_example_data_path(['bg01.jpg', 'bg02.jpg', 'bg03.jpg'])
bg = hp.core.io.load_average(bgpath, refimg = raw_holo)
holo = bg_correct(raw_holo, bg)
hp.show(holo)
Here, we have used load_average()
to construct an average of the three background
images specified in bgpath
. The refimg
argument allows us to specify a reference
image that is used to provide spacing and other metadata to the new, averaged image.
If you are worried about stray light in your optical train, you should
also capture a darkfield image of your sample, recorded with no laser illumination.
A darkfield image is specified as an optional third argument to bg_correct()
.
dfpath = get_example_data_path('df01.jpg')
df = hp.load_image(dfpath, spacing = 0.0851)
holo = bg_correct(raw_holo, bg, df)
hp.show(holo)
Some convenient tools for manipulating image data are included within HoloPy. See the HoloPy Tools page for details.
Telling HoloPy about your Experimental Setup¶
Recorded holograms are a product of the specific experimental setup that produced them. The image only makes sense when considered with information about the experimental conditions in mind. When you load an image, you have the option to specify some of this information in the form of metadata that is associated with the image. In fact, we already saw an example of this when we specified image spacing earlier. The sample in our image was immersed in water, which has a refractive index of 1.33. It was illuminated by a red laser with wavelength of 660 nm and polarization in the xdirection. We can tell HoloPy all of this information when loading the image:
raw_holo = hp.load_image(imagepath, spacing=0.0851, medium_index=1.33, illum_wavelen=0.660, illum_polarization=(1,0))
You can then view these metadata values as attributes of raw_holo
, as in raw_holo.medium_index
.
However, you must use a special function update_metadata()
to edit them. If we forgot to
specify metadata when loading the image, we can use update_metadata()
to add it later:
holo = hp.core.update_metadata(holo, medium_index=1.33, illum_wavelen=0.660, illum_polarization=(1,0))
Note
Spacing and wavelength must both be written in the same units  microns in the example above. Holopy has no builtin length scale and requires only that you be consistent. For example, we could have specified both parameters in terms of nanometers or meters instead.
HoloPy images are stored as xarray DataArray objects. xarray is a powerful tool that makes it easy to keep track of metadata and extra image dimensions, distinguishing between a time slice and a volume slice, for example. While you do not need any knowledge of xarray to use HoloPy, some familiarity will make certain tasks easier. This is especially true if you want to directly manipulate data before or after applying HoloPy’s builtin functions.
Saving and Reloading Holograms¶
Once you have backgrounddivided a hologram and associated it with metadata, you might want to save it so that you can skip those steps next time you are working with the same image:
hp.save('outfilename', holo)
saves your processed image to a compact HDF5 file. In fact, you can use save()
on any holopy object. To reload your same hologram with metadata you would write:
holo = hp.load('outfilename')
If you would like to save your hologram to an image format for easy visualization, use:
hp.save_image('outfilename', holo)
Additional options of save_image()
allow you to control how image intensity is scaled.
Images saved as .tif (the default) will still contain metadata, which will be retrieved if
you reload with load()
, but not load_image()
Note
Although HoloPy stores metadata even when writing to .tif image files, it is still recommended that
holograms be saved in HDF5 format using save()
. Floating point intensity values are rounded
to 8bit integers when using save_image()
, resulting in information loss.
NonSquare Pixels¶
The holograms above make use of several default assumptions. When you load an image like:
raw_holo = hp.load_image(imagepath, spacing = 0.0851)
you are making HoloPy assume a square array of evenly spaced grid points. If your pixels are not square, you can provide pixel spacing values in each direction:
raw_holo = hp.load_image(imagepath, spacing = (0.0851, 0.0426))
Most displays will default to displaying square pixels but if you
use HoloPy’s builtin show()
function to display the image, your hologram will display
with pixels of the correct aspect ratio.
Reconstructing Data (Numerical Propagation)¶
A hologram contains information about the electric field amplitude and phase at the detector plane. Shining light back through a hologram allows reconstruction of the electric field at points upstream of the detector plane. HoloPy performs this function mathematically by numerically propagating a hologram (or electric field) to another position in space. This allows you to reconstruct 3D sample volumes from 2D images. The light source is assumed to be collimated.
Example Reconstruction¶
import numpy as np
import holopy as hp
from holopy.core.io import get_example_data_path, load_average
from holopy.core.process import bg_correct
imagepath = get_example_data_path('image01.jpg')
raw_holo = hp.load_image(imagepath, spacing = 0.0851, medium_index = 1.33, illum_wavelen = 0.66, )
bgpath = get_example_data_path(['bg01.jpg','bg02.jpg','bg03.jpg'])
bg = load_average(bgpath, refimg = raw_holo)
holo = bg_correct(raw_holo, bg)
zstack = np.linspace(0, 20, 11)
rec_vol = hp.propagate(holo, zstack)
hp.show(rec_vol)
(Source code, png, hires.png, pdf)
We’ll examine each section of code in turn. The first block:
import numpy as np
import holopy as hp
from holopy.core.io import get_example_data_path, load_average
from holopy.core.process import bg_correct
loads the relevant modules from HoloPy and NumPy. The second block:
imagepath = get_example_data_path('image01.jpg')
raw_holo = hp.load_image(imagepath, spacing = 0.0851, medium_index = 1.33, illum_wavelen = 0.66)
bgpath = get_example_data_path(['bg01.jpg','bg02.jpg','bg03.jpg'])
bg = load_average(bgpath, refimg = raw_holo)
holo = bg_correct(raw_holo, bg)
reads in a hologram and divides it by a corresponding background image. If this is unfamiliar to you, please review the Loading Data tutorial.
Next, we use numpy’s linspace to define a set of distances between the image plane and the reconstruction plane at 2micron intervals to
propagate our image to. You can also propagate to a single distance
or to a set of distances obtained in some other fashion.
The actual propagation is accomplished with propagate()
:
zstack = np.linspace(0, 20, 11)
rec_vol = hp.propagate(holo, zstack)
Here, HoloPy has projected the hologram image through space to each of the distances contained in zstack
by using the metadata that we
specified when loading the image. If we forgot to load optical metadata with the image,
we can explicitly indicate the parameters for propagation to obtain an identical object:
rec_vol = hp.propagate(holo, zstack, illum_wavelen = 0.660, medium_index = 1.33)
Visualizing Reconstructions¶
You can display the reconstruction with show()
:
hp.show(rec_vol)
Pressing the left and right arrow keys steps through volumes slices  propagation to different zplanes. (Don’t use the down arrow key; it will mess up the stepping due to a peculiarity of Matplotlib. If this happens, close your plot window and show it again. Sorry.). If the left and right arrow keys don’t do anything, you might need to set your matplotlib backend. Refer to Using HoloPy for instructions.
Reconstructions are actually comprised of complex numbers. show()
defaults to showing you the amplitude of the image. You can get different, and
sometimes better, contrast by viewing the phase angle or imaginary part of the
reconstruction:
hp.show(rec_vol.imag)
hp.show(np.angle(rec_vol))
These phase sensitive visualizations will change contrast as you step through because you hit different places in the phase period. Such a reconstruction will work better if you use steps that are an integer number of wavelengths in medium:
med_wavelen = holo.illum_wavelen / holo.medium_index
rec_vol = hp.propagate(holo, zstack*med_wavelen)
hp.show(rec_vol.imag)
Cascaded Free Space Propagation¶
HoloPy calculates reconstructions by performing a convolution of the hologram with the reference wave over the distance to be propagated. By default, HoloPy calculates a single transfer function to perform the convolution over the specified distance. However, a better reconstruction can sometimes be obtained by iteratively propagating the hologram over short distances. This cascaded free space propagation is particularly useful when the reconstructions have fine features or when propagating over large distances. For further details, refer to Kreis 2002.
To implement cascaded free space propagation in HoloPy, simply pass a cfsp
argument
into propagate()
indicating how many times the hologram should be iteratively
propagated. For example, to propagate in three steps over each distance, we write:
rec_vol = hp.propagate(holo, zstack, cfsp = 3)
Reconstructing Point Source Holograms¶
Holograms are typically reconstructed optically by shining light back through them. This corresponds mathematically to propagating the field stored in the hologram to some different plane. The propagation performed here assumes that the hologram was recorded using a point source (diverging spherical wave) as the light source. This is also known as lensfree holography. Note that this is different than propagation calculations where a collimated light source (plane wave) is used. For recontructions using a plane wave see Reconstructing Data (Numerical Propagation).
This pointsource propagation calculation is an implementation of the algorithm that appears in Jericho and Kreuzer 2010. Curently, only square input images and propagation through media with a refractive index of 1 are supported.
Example Reconstruction¶
import holopy as hp
import numpy as np
from holopy.core.io import get_example_data_path
from holopy.propagation import ps_propagate
from scipy.ndimage.measurements import center_of_mass
imagepath = get_example_data_path('ps_image01.jpg')
bgpath = get_example_data_path('ps_bg01.jpg')
L = 0.0407 # distance from light source to screen/camera
cam_spacing = 12e6 # linear size of camera pixels
mag = 9.0 # magnification
npix_out = 1020 # linear size of output image (pixels)
zstack = np.arange(1.08e3, 1.18e3, 0.01e3) # distances from camera to reconstruct
holo = hp.load_image(imagepath, spacing=cam_spacing, illum_wavelen=406e9, medium_index=1) # load hologram
bg = hp.load_image(bgpath, spacing=cam_spacing) # load background image
holo = hp.core.process.bg_correct(holo, bg+1, bg) # subtract background (not divide)
beam_c = center_of_mass(bg.values.squeeze()) # get beam center
out_schema = hp.core.detector_grid(shape=npix_out, spacing=cam_spacing/mag) # set output shape
recons = ps_propagate(holo, zstack, L, beam_c, out_schema) # do propagation
hp.show(abs(recons[:,350:550,450:650])) # display result
(Source code, png, hires.png, pdf)
We’ll examine each bsection of code in turn. The first block:
import holopy as hp
import numpy as np
from holopy.core.io import get_example_data_path
from holopy.propagation import ps_propagate
from scipy.ndimage.measurements import center_of_mass
loads the relevant modules. The second block:
imagepath = get_example_data_path('ps_image01.jpg')
bgpath = get_example_data_path('ps_bg01.jpg')
L = 0.0407 # distance from light source to screen/camera
cam_spacing = 12e6 # linear size of camera pixels
mag = 9.0 # magnification
npix_out = 1020 # linear size of output image (pixels)
zstack = np.arange(1.08e3, 1.18e3, 0.01e3) # distances from camera to reconstruct
defines all parameters used for the reconstruction. Numpy’s linspace was used to define a set of distances at 10micron intervals to propagate our image to. You can also propagate to a single distance or to a set of distances obtained in some other fashion. The third block:
holo = hp.load_image(imagepath, spacing=cam_spacing, illum_wavelen=406e9, medium_index=1) # load hologram
bg = hp.load_image(bgpath, spacing=cam_spacing) # load background image
holo = hp.core.process.bg_correct(holo, bg+1, bg) # subtract background (not divide)
beam_c = center_of_mass(bg.values.squeeze()) # get beam center
out_schema = hp.core.detector_grid(shape=npix_out, spacing=cam_spacing/mag) # set output shape
reads in a hologram and subtracts the corresponding background image. If this is unfamiliar to you, please review the Loading Data tutorial. The third block also finds the center of the reference beam and sets the size and pixel spacing of the output images.
Finally, the actual propagation is accomplished with
ps_propagate()
and a cropped region of the result is
displayed. See the Reconstructing Data (Numerical Propagation) page for details on
visualizing the reconstruction results.
recons = ps_propagate(holo, zstack, L, beam_c, out_schema) # do propagation
hp.show(abs(recons[:,350:550,450:650])) # display result
Magnification and Output Image Size¶
Unlike the case where a collimated beam is used as the illumination
and the pixel spacing in the reconstruction is the same as in the
original hologram, for lensfree reconstructions the pixel spacing
in the reconstruction can be chosen arbitrarily. In order to magnify
the reconstruction the spacing in the reconstruction plane should be
smaller than spacing in the original hologram. In the code above, the
magnification of the reconstruction can be set using the variable
mag
, or when calling ps_propagate()
directly the desired
pixel spacing in the reconstruction is specified through the
spacing of out_schema
. Note that the output spacing will not be
the spacing of out_schema
exactly, but should be within a few
percent of it. We recommend calling get_spacing()
on recons
to get the actual spacing used.
Note that the total physical size of the plane that is reconstructed
remains the same when different output pixel spacings are used. This
means that reconstructions with large output spacings will only have
a small number of pixels, and reconstructions with small output
spacings will have a large number of pixels. If the linear size (in
pixels) of the total reconstruction plane is smaller than
npix_out
, the entire reconstruction plane will be returned.
However, if the linear size of total reconstruction plane is
larger than npix_out
, only the center region of the
reconstruction plane with linear size npix_out
is returned.
In the current version of the code, the amount of memory needed to
perform a reconstruction scales with mag
^{2}. Presumably this
limitation can be overcome by implementing the steps described in the
Convolution section of the Appendix of
Jericho and Kreuzer 2010.
Scattering Calculations¶
Optical physicists and astronomers have worked out how to compute the scattering of light from many kinds of objects. HoloPy provides an easy interface for computing scattered fields, intensities, scattering matrices, crosssections, and holograms generated by microscopic objects.
A Simple Example¶
Let’s start by calculating an inline hologram generated by a plane wave scattering from a microsphere.
import holopy as hp
from holopy.scattering import calc_holo, Sphere
sphere = Sphere(n = 1.59, r = .5, center = (4, 4, 5))
medium_index = 1.33
illum_wavelen = 0.660
illum_polarization = (1,0)
detector = hp.detector_grid(shape = 100, spacing = .1)
holo = calc_holo(detector, sphere, medium_index, illum_wavelen, illum_polarization)
hp.show(holo)
We’ll examine each section of code in turn. The first few lines :
import holopy as hp
from holopy.scattering import calc_holo, Sphere
load the relevant modules from HoloPy that we’ll need for doing our calculation. The next line describes the scatterer we would like to model:
sphere = Sphere(n = 1.59, r = .5, center = (4, 4, 5))
We will be scattering light off a Scatterer
object,
specifically a Sphere
. A Scatterer
object
contains information about the geometry (position, size, shape) and optical
properties (refractive index) of the object that is scattering light. We’ve
defined a spherical scatterer with radius 0.5 microns and index of refraction
1.59. This refractive index is approximately that of polystyrene. Next, we need
to describe how we are illuminating our sphere, and how that light will be
detected:
medium_index = 1.33
illum_wavelen = 0.66
illum_polarization = (1,0)
detector = hp.detector_grid(shape = 100, spacing = .1)
We are going to be using red light (wavelength = 660 nm in vacuum) polarized in the xdirection to illuminate a sphere immersed in water (refractive index = 1.33). Refer to Units and Coordinate System if you’re confused about how the wavelength and polarization are specified.
The scattered light will be collected at a detector, which is frequently a
digital camera mounted onto a microscope. We defined our detector as a 100 x
100 pixel array, with each square pixel of side length .1 microns. The
shape
argument tells HoloPy how many pixels are in the detector and affects
computation time. The spacing
argument tells HoloPy how far apart each
pixel is. Both parameters affect the absolute size of the detector.
After getting everything ready, the actual scattering calculation is straightforward:
holo = calc_holo(detector, sphere, medium_index, illum_wavelen, illum_polarization)
hp.show(holo)
Congratulations! You just calculated the inline hologram generated at the detector plane by interference between the scattered field and the reference wave. For an inline hologram, the reference wave is simply the part of the field that is not scattered or absorbed by the particle.
You might have noticed that our scattering calculation requires much of the same
metadata we specified when loading an image. If we have an experimental image
from the system we would like to model, we can use that as an argument in
calc_holo()
instead of our detector
object created from
detector_grid()
. HoloPy will calculate a hologram image with pixels at
the same positions as the experimental image, and so we don’t need to worry
about making a detector_grid()
with the correct shape
and spacing
arguments.
from holopy.core.io import get_example_data_path
imagepath = get_example_data_path('image0002.h5')
exp_img = hp.load(imagepath)
holo = calc_holo(exp_img, sphere)
Note that we didn’t need to explicitly specify illumination information when
calling calc_holo()
, since our image contained saved metadata and HoloPy
used its values. Passing an image to a scattering function is particularly
useful when comparing simulated data to experimental results, since we can
easily recreate our experimental conditions exactly.
So far all of the images we have calculated are holograms, or the interference
pattern that results from the superposition of a scattered wave with a reference
wave. Holopy can also be used to examine scattered fields on their own. Simply
replace calc_holo()
with calc_field()
to look at scattered
electric fields (complex) or calc_intensity()
to look at field
amplitudes, which is the typical measurement in a light scattering experiment.
More Complex Scatterers¶
Coated Spheres¶
HoloPy can also calculate holograms from coated (or multilayered) spheres. Constructing a coated sphere differs only in specifying a list of refractive indices and outer radii corresponding to the layers (starting from the core and working outwards).
coated_sphere = Sphere(center=(2.5, 5, 5), n=(1.59, 1.42), r=(0.3, 0.6))
holo = calc_holo(exp_img, coated_sphere)
hp.show(holo)
If you prefer thinking in terms of the thickness of subsequent layers, instead
of their distance from the center, you can use LayeredSphere
to achieve
the same result:
from holopy.scattering import LayeredSphere
coated_sphere = LayeredSphere(center=(2.5, 5, 5), n=(1.59, 1.42), t=(0.3, 0.3))
Collection of Spheres¶
If we want to calculate a hologram from a collection of spheres, we must
first define the spheres individually, and then combine them into a
Spheres
object:
from holopy.scattering import Spheres
s1 = Sphere(center=(5, 5, 5), n = 1.59, r = .5)
s2 = Sphere(center=(4, 4, 5), n = 1.59, r = .5)
collection = Spheres([s1, s2])
holo = calc_holo(exp_img, collection)
hp.show(holo)
Adding more spheres to the cluster is as simple as defining more
sphere objects and passing a longer list of spheres to the
Spheres
constructor.
Nonspherical Objects¶
To define a nonspherical scatterer, use Spheroid
or Cylinder
objects. These axisymmetric scatterers are defined by two dimensions, and can describe scatterers that are elongated or squashed along one direction.
By default, these objects are aligned with the zaxis, but they can be rotated into any orientation by passing a set of Euler angles to the rotation
argument when defining the scatterer. See Rotations of Scatterers for information on how these angles are defined.
As an example, here is a hologram produced by a cylinder aligned with the vertical axis (xaxis
according to the HoloPy Coordinate System).
Note that the hologram image is elongated in the horizontal direction since the sides of the cylinder scatter light more than the ends.
import numpy as np
from holopy.scattering import Cylinder
c = Cylinder(center=(5, 5, 7), n = 1.59, d=0.5, h=2, rotation=(0,np.pi/2, 0))
holo = calc_holo(exp_img, c)
hp.show(holo)
Customizing Scattering Calculations¶
While the examples above will be sufficient for most purposes, there are a few additional options that are useful in certain scenarios.
Scattering Theories in HoloPy¶
HoloPy contains a number of scattering theories to model the scattering from different kinds of scatterers. By default, scattering from single spheres is calculated using Mie theory, which is the exact solution to Maxwell’s equations for the scattered field from a spherical particle, originally derived by Gustav Mie and (independently) by Ludvig Lorenz in the early 1900s.
A scatterer composed of multiple spheres can exhibit multiple scattering
and coupling of the nearfields of neighbouring particles. Mie theory doesn’t include
these effects, so Spheres
objects are by default calculated using the
SCSMFO package from Daniel Mackowski.
This calculation uses Tmatrix methods to give the exact solution to Maxwell’s equation
for the scattering from an arbitrary arrangement of nonoverlapping spheres.
Sometimes you might want to calculate scattering from multiple spheres
using Mie theory if you are worried about computation time or if you are
using multilayered spheres (HoloPy’s implementation of the multisphere theory
can’t currently handle coated spheres). You can specify Mie theory manually when
calling the calc_holo()
function:
from holopy.scattering import Mie
holo = calc_holo(exp_img, collection, theory = Mie)
Similarly, HoloPy calculates scattering from cylindrical or spheroidal particles by using Tmatrix code from Michael Mishchenko, but these scatterer types are not compatible with Mie theory.
Holopy can also access a discrete dipole approximation (DDA) theory to model arbitrary nonspherical objects. See the Scattering from Arbitrary Structures with DDA tutorial for more details. It is fairly easy to add your own scattering theory to HoloPy. See Adding a new scattering theory for details. If you think your new scattering theory may be useful for other users, please consider submitting a pull request.
Detector Types in HoloPy¶
The detector_grid()
function we saw earlier creates holograms that
display nicely and are easily compared to experimental images. However, they can
be computationally expensive, as they require calculations of the electric field
at many points. If you only need to calculate values at a few points, or if your
points of interest are not arranged in a 2D grid, you can use
detector_points()
, which accepts either a dictionary of coordinates or
indvidual coordinate dimensions:
x = [0, 1, 0, 1, 2]
y = [0, 0, 1, 1, 1]
z = 1
coord_dict = {'x': x, 'y': y, 'z': z}
detector = hp.detector_points(x = x, y = y, z = z)
detector = hp.detector_points(coord_dict)
The coordinates for detector_points()
can be specified in terms of either
Cartesian or spherical coordinates. If spherical coordinates are used, the
center
value of your scatterer is ignored and the coordinates are
interpreted as being relative to the scatterer.
Static light scattering calculations¶
Scattering Matrices¶
In a static light scattering measurement you record the scattered intensity at a number of locations. A common experimental setup contains multiple detectors at a constant radial distance from a sample (or a single detector on a goniometer arm that can swing to multiple angles.) In this kind of experiment you are usually assuming that the detector is far enough away from the particles that the farfield approximation is valid, and you are usually not interested in the exact distance of the detector from the particles. So, it’s most convenient to work with amplitude scattering matrices that are angledependent. (See [Bohren1983] for further mathematical description.)
import numpy as np
from holopy.scattering import calc_scat_matrix
detector = hp.detector_points(theta = np.linspace(0, np.pi, 100), phi = 0)
distant_sphere = Sphere(r=0.5, n=1.59)
matr = calc_scat_matrix(detector, distant_sphere, medium_index, illum_wavelen)
Here we omit specifying the location (center) of the scatterer. This is
only valid when you’re calculating a farfield quantity. Similarly, note
that our detector, defined from a detector_points()
function,
includes information about direction but not distance. It is typical
to look at scattering matrices on a semilog plot. You can make one as follows:
import matplotlib.pyplot as plt
plt.figure()
plt.semilogy(np.linspace(0, np.pi, 100), abs(matr[:,0,0])**2)
plt.semilogy(np.linspace(0, np.pi, 100), abs(matr[:,1,1])**2)
plt.show()
You are usually interested in the intensities of the scattered fields, which are proportional to the modulus squared of the amplitude scattering matrix. The diagonal elements give the intensities for the incident light and the scattered light both polarized parallel and perpendicular to the scattering plane, respectively.
Scattering CrossSections¶
The scattering cross section provides a measure of how much light from an
incident beam is scattered by a particular scatterer. Similar to calculating
scattering matrices, we can omit the position of the scatterer for calculation
of cross sections. Since cross sections integrates over all angles, we can also
omit the detector
argument entirely:
from holopy.scattering import calc_cross_sections
x_sec = calc_cross_sections(distant_sphere, medium_index, illum_wavelen, illum_polarization)
x_sec returns an array containing four elements. The first element is the scattering cross section, specified in terms of the same units as wavelength and particle size. The second and third elements are the absorption and extinction cross sections, respectively. The final element is the average value of the cosine of the scattering angle.
Scattering from Arbitrary Structures with DDA¶
The discrete dipole approximation (DDA) lets us calculate scattering
from any arbitrary object by representing it as a closely packed array
of point dipoles. In HoloPy you can make use of the DDA by specifying
a general Scatterer
with an indicator function (or set of
functions for a composite scatterer containing multiple media).
HoloPy uses ADDA to do the actual DDA calculations, so you will need to install ADDA and be able to run:
adda
at a terminal for HoloPy DDA calculations to succeed. To install ADDA, first download or clone the code from GitHub. In a terminal window, go to the directory ’adda/src’ and compile using one of three options:
make seq
or:
make
or:
make OpenCL
make seq
will not take advantage of any parallel processing of the cores
on your computer. make
uses mpi for parallel processing. make OpenCL
uses
OpenCL for parallel processing. If the make does not work due to missing packages,
you will have to download the specified packages and install them.
Next, you must modify your path in your .bashrc or /bash_profile (for mac). Add the line:
export PATH=$PATH:userpath/adda/src/seq
or:
export PATH=$PATH:userpath/adda/src/mpi
or:
export PATH=$PATH:userpath/adda/src/OpenCL
where you should use the path that matches the make you chose above.
A lot of the code associated with DDA is fairly new so be careful; there are probably bugs. If you find any, please report them.
Defining the geometry of the scatterer¶
To calculate the scattering pattern for an arbitrary object, you first need an indicator function which outputs ‘True’ if a test coordinate lies within your scatterer, and ‘False’ if it doesn’t. The indicator function is an argument of the constructor of your scatterer.
For example, if you wanted to define a dumbbell consisting of the union of two overlapping spheres you could do so like this:
import holopy as hp
from holopy.scattering import Scatterer, Sphere, calc_holo
import numpy as np
s1 = Sphere(r = .5, center = (0, .4, 0))
s2 = Sphere(r = .5, center = (0, .4, 0))
detector = hp.detector_grid(100, .1)
dumbbell = Scatterer(lambda point: np.logical_or(s1.contains(point), s2.contains(point)),
1.59, (5, 5, 5))
holo = calc_holo(detector, dumbbell, medium_index=1.33, illum_wavelen=.66, illum_polarization=(1, 0))
Here we take advantage of the fact that Spheres can tell us if a point
lies inside them. We use s1
and s2
as purely geometrical
constructs, so we do not give them indices of refraction, instead
specifying n when defining dumbbell
.
HoloPy contains convenient wrappers for many builtin ADDA constructions.
The dumbbell defined explicitly above could also have been defined with the HoloPy Bisphere
class instead.
Similar classes exist to define an Ellipsoid
, Cylinder
, or Capsule
.
Mutiple Materials: A Janus Sphere¶
You can also provide a set of indicators and indices to define a scatterer containing multiple materials. As an example, lets look at a janus sphere consisting of a plastic sphere with a high index coating on the top half:
from holopy.scattering.scatterer import Indicators
import numpy as np
s1 = Sphere(r = .5, center = (0, 0, 0))
s2 = Sphere(r = .51, center = (0, 0, 0))
def cap(point):
return(np.logical_and(np.logical_and(point[...,2] > 0, s2.contains(point)),
np.logical_not(s1.contains(point))))
indicators = Indicators([s1.contains, cap],
[[.51, .51], [.51, .51], [.51, .51]])
janus = Scatterer(indicators, (1.34, 2.0), (5, 5, 5))
holo = calc_holo(detector, janus, medium_index=1.33, illum_wavelen=.66, illum_polarization=(1, 0))
We had to manually set up the bounds of the indicator functions here because the automatic bounds determination routine gets confused by the cap that does not contain the origin.
We also provide a JanusSphere
scatterer which is very
similar to the scatterer defined above, but can also take a rotation
angle to specify other orientations:
from holopy.scattering import JanusSphere
janus = JanusSphere(n = [1.34, 2.0], r = [.5, .51], rotation = (np.pi/2, 0),
center = (5, 5, 5))
Bayesian inference of Parameter Values¶
Scattering calculations can inform us about the hologram produced by a specific scatterer, but they can’t tell us anything about what type of scatterer produced an experimentally measured hologram. For this inverse scattering problem, we turn to a Bayesian inference approach. We calculate the holograms produced by many similar scatterers, and we evaluate which ones are closest to our measured hologram. We can then use known information about the scatterers to determine the parameters of the scatterer (for example, its size and refractive index) that are most likely to have produced the observed hologram.
In the following example, we will infer the size, refractive index, and position of a spherical scatterer. The approach and formalism is described in more detail in [Dimiduk2016]. For more information on Bayesian inference in general, see [Gregory2005].
Here is the full example. We’ll go through it stepbystep afterward:
import holopy as hp
import numpy as np
from holopy.core.io import get_example_data_path, load_average
from holopy.core.process import bg_correct, subimage, normalize
from holopy.scattering import Sphere, calc_holo
from holopy.inference import prior, AlphaModel, tempered_sample
# load an image
imagepath = get_example_data_path('image01.jpg')
raw_holo = hp.load_image(imagepath, spacing = 0.0851, medium_index = 1.33, illum_wavelen = 0.66, illum_polarization = (1,0))
bgpath = get_example_data_path(['bg01.jpg','bg02.jpg','bg03.jpg'])
bg = load_average(bgpath, refimg = raw_holo)
data_holo = bg_correct(raw_holo, bg)
# process the image
data_holo = subimage(data_holo, [250,250], 200)
data_holo = normalize(data_holo)
# Set up the prior
s = Sphere(n=prior.Gaussian(1.5, .1), r=prior.BoundedGaussian(.5, .05, 0, np.inf),
center=prior.make_center_priors(data_holo))
# Set up the noise model
noise_sd = data_holo.std()
model = AlphaModel(s, noise_sd=noise_sd, alpha=1)
result = tempered_sample(model, data_holo)
result.values()
hp.save('examplesampling.h5', result)
The first few lines import the code needed to compute holograms and do parameter estimation.
import holopy as hp
import numpy as np
from holopy.core.io import get_example_data_path, load_average
from holopy.core.process import bg_correct, subimage, normalize
from holopy.scattering import Sphere, calc_holo
from holopy.inference import prior, AlphaModel, tempered_sample
Preparing Data¶
Next, we load a hologram from a file using the same steps as those in Loading Data
# load an image
imagepath = get_example_data_path('image01.jpg')
raw_holo = hp.load_image(imagepath, spacing = 0.0851, medium_index = 1.33, illum_wavelen = 0.66, illum_polarization = (1,0))
bgpath = get_example_data_path(['bg01.jpg','bg02.jpg','bg03.jpg'])
bg = load_average(bgpath, refimg = raw_holo)
data_holo = bg_correct(raw_holo, bg)
You will notice that the hologram data is localized to a region near the center
of the image. We only want to compare calculated holograms to this region, so we
will crop our image with subimage()
. We also need to normalize the data
so that its mean is 1, since calculations return a normalized result. Since our
image is background divided, its mean is already very close to 1, but it is good
to get in the habit of normalizing anyway.
# process the image
data_holo = subimage(data_holo, [250,250], 200)
data_holo = normalize(data_holo)
Note
It is often useful to test an unfamiliar technique on data for which you know the expected outcome.
Instead of actual data, you could use a hologram calculated from calc_holo()
, and modulated
by random noise with add_noise()
.
Defining a Probability Model¶
Priors¶
We know that the hologram was produced by a spherical scatterer, so we want to
define a Sphere
object like we did in the Scattering Calculations
tutorial. However, in this case we don’t know what parameters to specify for the
sphere (since that is what we’re trying to find out). Instead, we write down a
probabilistic statement of our prior information about the sphere. In
statistics, we call this a prior. For the case we are investigating here, we
would probably have some best guess and uncertainty about the size and index of
the particle, obtained from the supplier or from prior work with the particle.
We will guess the radius to be 0.5 micrometers (with 50 nm error) and refractive
index to be 1.5 (with 0.1 error). We also need to provide a prior for the
position of the sphere. We can use a hough()
transform to get a pretty
good guess of where the particle is in x and y, but it is difficult to determine
where it is in z.
Note
One trick to get a better estimate of z position is to numerically propagate
the hologram backwards in space with propagate()
, and look for where
the interference fringes vanish.
Let’s turn our information about priors into code by defining our scatterer:
s = Sphere(n=prior.Gaussian(1.5, .1), r=prior.BoundedGaussian(.5, .05, 0, np.inf),
center=prior.make_center_priors(data_holo))
We use a Gaussian distribution as the prior because it is the most conservative
distribution we can use if all we know is some expected value and some
uncertainty on that expected value. For the radius we also know that it must be
nonnegative, so we can bound the Gaussian at zero. The
make_center_priors()
function automates generating priors for a sphere
center using center_find()
(based on a hough transform). It assigns
Gaussian priors for x and y, and picks a large uniform prior for z to represent
our ignorance about how far the particle is from the imaging plane. In this case
the center prior will be:
[Gaussian(mu=11.4215, sd=0.0851),
Gaussian(mu=9.0945, sd=0.0851),
Uniform(lower_bound=0, upper_bound=170.2)]
Likelihood¶
Next we need to define a model that tells HoloPy how probable it is that we would see the data we observed given some hypothetical scatterer position, size and index. In the language of statistics, this is referred to as a likelihood. In order to compute a likelihood, you need some estimate of how noisy your data is (so that you can figure out how likely it is that the differences between your model and data could be explained by noise). Here we use the standard deviation of the data, which is an overestimate of the true noise, since it also includes variation due to our signal.
noise_sd = data_holo.std()
model = AlphaModel(s, noise_sd=noise_sd, alpha=1)
Note
alpha
is a model parameter that scales the scattered beam intensity
relative to the reference beam. It is often less than 1 for reasons that are
poorly understood. If you aren’t sure what value it should take in your
system, you can allow alpha
to vary by giving it a prior, as we did with
the sphere parameters.
Sampling the Posterior¶
Finally, we can sample the posterior probability for this model. Essentially, a
set of proposed scatterers are randomly generated according to the priors we
specified. Each of these scatterers is then evaluated in terms of how well it
matches the experimental hologram data_holo
. A Monte Carlo algorithm
iteratively produces and tests sets of scatterers to find the scatterer
parameters that best reproduce the target hologram. We end up with a
distribution of values for each parameter (the posterior) that represents our
updated knowledge about the scatterer when accounting for the expected
experimental hologram. To do the actual sampling, we use
tempered_sample()
(ignoring any RuntimeWarnings about invalid values):
result = tempered_sample(model, data_holo)
The above line of code may take a long time to run (it takes 1015 mins on our 8core machines). If you just want to quickly see what the results look like, try:
result = tempered_sample(model, data_holo, nwalkers=10, samples=100, max_pixels=100)
This code should run very quickly, but its results cannot be trusted for any
actual data. Nevertheless, it can give you an idea of what format the results
will take. In our last line of code, we have adjusted three parameters to make
the code run faster: nwalkers
describes the number of scatterers produced in
each generation. samples
describes how many generations of scatterers to
produce. Together, they define how many scatterering calculations must be
performed. For the values chosen in the fast code, a Monte Carlo steady state
will not yet have been achieved, so the resulting posterior distribution is not
very meaningful. max_pixels
describes the maximum number of pixels compared
between the experimental holgoram and the test holograms. It turns out that
holograms contain a lot of redundant information owing to their symmetry, so a
subset of pixels can be analyzed without loss of accuracy. However, 100 pixels
is probably too few to capture all of the relevant information in the hologram.
You can get a quick look at our obtained values with:
result.values()
result.values()
gives you the maximum a posteriori probability (MAP) value as
well as onesigma credibility intervals (or you can request any other sigma with
an argument to the function). You can also look only at central measures:
result.MAP
result.mean
result.median
Since calculation of useful results takes a long time, you will usually want to save them to an hdf5 file:
hp.save('examplesampling.h5', result)
References¶
[Gregory2005]  Gregory, P. (2005) Bayesian Logical Data Analysis. Cambridge University Press 
Fitting Models to Data¶
In addition to Bayesian inference, HoloPy can also do simpler leastsquares fits to determine the scatterer parameters that best match an experimentally measured hologram. The main advantage of this technique is that it can be much faster. The drawback is that good intial guesses of each parameter are required to obtain accurate results.
Note
The HoloPy fitting methods have been superseded by the Bayesian inference techniques described in the Bayesian inference of Parameter Values tutorial. We strongly recommend that approach unless you have a good reason that fitting is preferable in your particular situation.
A Simple Fit¶
We start by loading and processing data just as we did for the parameter inference in the previous tutorial.
import holopy as hp
import numpy as np
from holopy.core.io import get_example_data_path, load_average
from holopy.core.process import bg_correct, subimage, normalize
from holopy.scattering import Sphere, calc_holo
# load an image
imagepath = get_example_data_path('image01.jpg')
raw_holo = hp.load_image(imagepath, spacing = 0.0851, medium_index = 1.33, illum_wavelen = 0.66, illum_polarization = (1,0))
bgpath = get_example_data_path(['bg01.jpg','bg02.jpg','bg03.jpg'])
bg = load_average(bgpath, refimg = raw_holo)
data_holo = bg_correct(raw_holo, bg)
# process the image
data_holo = subimage(data_holo, [250,250], 200)
data_holo = normalize(data_holo)
Define a Model¶
The model specification is a little bit different from the inference case.
First, we define a parameterized scatterer including initial guesses and absolute bounds
using the Parameter
class. Note that the bounds here are not uncertainty values as in
the inference case, but instead represent the full allowed range of a parameter (like the Uniform
prior).
The center
coordinates must be specified as (x, y, and z, in that order).
Here, we will keep particle radius and refractive index fixed. Fitting works best when there are only a few uncertain parameters.
You can find guesses for x and y coordinates with center_find()
, and guess z with propagate()
.
In this image (uncropped version), the particle’s center is near (24, 22, 15), with coordinates in microns.
from holopy.fitting import fit, Model
from holopy.fitting import Parameter as par
par_s = Sphere(center = (par(guess = 24, limit = [15,30]),
par(22, [15, 30]), par(15, [10, 20])), r = .5, n = 1.58)
Then this parametrized scatterer, along with a desired scattering calculation, is used to define a model:
model = Model(par_s, calc_holo, alpha = par(.6, [.1, 1]))
alpha
is an additional fitting parameter first introduced in [Lee2007] (see References and credits for additional details).
To see how well the guess in your model lines up with the hologram you are fitting to, use :
guess_holo = calc_holo(data_holo, par_s, scaling=model.alpha)
Run the Fit¶
Once you have all of that set up, running the fit is almost trivially simple:
result = fit(model, data_holo)
We can see just the fit results with result.scatterer.center
.
The initial guess of the sphere’s position (24, 22, 15)
was corrected by the fitter to (24.17,21.84,16.42). Notice that we have achieved subpixel position resolution!
From the fit,
result.scatterer
gives the scatterer that best matches the hologram,
result.alpha
is the alpha for the best fit. result.chisq
and
result.rsq
are statistical measures of the the goodness of the fit.
You can also compute a hologram of the final fit result to compare to the data with:
result_holo = calc_holo(data_holo, result.scatterer, scaling=result.alpha)
Finally, we save the result with:
hp.save('result.h5', result)
Speeding up Fits with Random Subset Fitting¶
A hologram usually contains far more information than is needed to determine the number of parameters you are interested in. Because of this, you can often get a significantly faster fit with no little or no loss in accuracy by fitting to only a random fraction of the pixels in a hologram.
result = fit(model, data_holo, random_subset=.01)
You will want to do some testing to make sure that you still get acceptable answers with your data, but our investigations have shown that you can frequently use random fractions of .1 or .01 with little effect on your results and gain a speedup of 10x or greater.
Advanced Parameter Specification¶
Complex Index of Refraction¶
You can specify a complex index with:
from holopy.fitting import ComplexParameter
Sphere(n = ComplexParameter(real = par(1.58, step = 0.01), imag = 1e4))
This will fit to the real part of index of refraction while holding
the imaginary part fixed. You can fit to it as well by specifying
imag = par(1e4)
instead of imag = 1e4
. In a case like this
where we are providing a small imaginary part for numerical stability,
you would not want to fit to it. However fitting to an imaginary index
component could be useful for a metal particle. Setting the key word argument step = 0.01
specifies the the step size used in calculating
the numerical derivatives of this parameter. Specifying a small step
size is often necessary when fitting for an index of refraction.
Tying Parameters¶
You may desire to fit holograms with tied parameters, in which several physical quantities that could be varied independently are constrained to have the same (but nonconstant) value. A common example involves fitting a model to a multiparticle hologram in which all of the particles are constrained to have the same refractive index, but the index is determined by the fitter. This may be done by defining a Parameter and using it in multiple places :
from holopy.scattering import Spheres
n1 = par(1.59)
sc = Spheres([Sphere(n = n1, r = par(0.5e6), \
center = [10., 10., 20.]), \
Sphere(n = n1, r = par(0.5e6), center = [9., 11., 21.])])
Developer’s Guide¶
Installing HoloPy for Developers¶
If you are going to hack on holopy, you probably only want to compile the scattering extensions.
For Mac and Linux:
Download or clone the latest version of HoloPy from Git Hub at https://github.com/manoharanlab/holopy.
Let’s say you downloaded or cloned HoloPy to
/home/me/holopy
. Then open a terminal, cd
to /home/me/holopy
and run:
python setup.py build_ext inplace
This puts the extensions inside the source tree, so that you can work
directly from /home/me/holopy
. You will need to add
/home/me/holopy
to your python_path for python to find the
module when you import it.
Note for Mac users: gfortran may put its library in a place python can’t find it. If you get errors including something like can't find /usr/local/libgfortran.3.dynlib
you can symlink them in from your install. You can do this by running:
sudo ln s /usr/local/gfortran/lib/libgfortran.3.dynlib /usr/local/lib
sudo ln s /usr/local/gfortran/lib/libquadmath.3.dynlib /usr/local/lib
For Windows: Installation on Windows is still a work in progress, but we have been able to get HoloPy working on Windows 10 with an AMD64 architecture (64bit) processor.
Install Anaconda with Python 3.6 and make sure it is working.
Install the C compiler. It’s included in Visual Studio 2015 Community. Make sure it is working with a C helloworld.
From now on, make sure any command prompt window invokes the right environment conditions for compiling with VC. To do this, make sure
C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\vcvarsall.bat
is added to the system path variable. This batch detects your architecture, then runs another batch that sets the path include the directory with the correct version of the VC compiler.Install cython and made sure it works.
Install Intel’s Fortran compiler. A good place to start is the trial version of Parallel Studio XE. Make sure it is working with a Fortran helloworld.
Install mingw32make, which does not come with Anaconda by default.
Download or clone HoloPy from https://github.com/manoharanlab/holopy.
Open the command prompt included in Intel’s Parallel Studio. Run
holopy/setup.py
. It is necessay to use Intel’s Parallel Studio command prompt to avoid compiling errors.Install the following dependencies that don’t come with Anaconda:
conda install xarray dask netCDF4 bottleneck conda install c astropy emcee=2.2.1
Open an iPython console where holopy is installed and try
import holopy
.
If the above procedure doesn’t work, or you find something else that does, please let us know so that we can improve these instructions.
How HoloPy Stores Data¶
Images in HoloPy are stored in the format of xarray DataArrays. Spatial
information is tracked in the DataArray’s dims
and coords
fields
according to the HoloPy Coordinate System. Additional dimensions are
sometimes specified to account for different zslices, times, or field
components, for example. Optical parameters like refractive index and
illumination wavelength are stored in the DataArray’s attrs
field.
The detector_grid()
function simply creates a 2D image composed entirely
of zeros. In contrast, the detector_points()
function creates a DataArray
with a single dimension named ‘point’. Spatial coordinates (in either Cartesian
or spherical form) track this dimension, so that each data value in the array
has its own set of coordinates unrelated to its neighbours. This type of
onedimensional organization is sometimes used for 2D images as well. Inference
and fitting methods typically use only a subset of points in an image (see
Speeding up Fits with Random Subset Fitting), and so it makes sense for them to keep track of lists of
location coordinates instead of a grid. Furthermore, HoloPy’s scattering
functions accept coordinates in the form of a 3xN array of coordinates. In both
of these cases, the 2D image is flattened into a 1D DataArray like that created
by detector_points()
. In this case the single dimension is ‘flat’ instead
of ‘point’. HoloPy treats arrays with these two named dimensions identically,
except that the ‘flat’ dimension can be unstacked to restore a 2D image or 3D
volume.
HoloPy’s use of DataArrays sometimes assigns smaller DataArrays in attrs
,
which can lead to problems when saving data to a file. When saving a DataArray
to file, HoloPy converts any DataArrays in attrs
to numpy arrays, and keeps
track of their dimension names separately. HoloPy’s save_image()
writes a
yaml dump of attrs
(along with spacing information) to the
imagedescription
field of .tif file metadata.
Bayesian inference of Parameter Values returns a lot of information, which is stored in the form of a SamplingResult
object.
This object stores the model and EmceeStrategy
that were used in the inference calculation as attributes.
An additional attribute named dataset
is an xarray Dataset
that contains both the data used in the inference calculation, as well as the raw output.
The parameter values at each step of the sampling chain and the calculated logprobabilities at each step are stored here under the samples
and lnprobs
namespaces.
Adding a new scattering theory¶
Adding a new scattering theory is relatively straightforward. You just need to define a new scattering theory class and implement one or two methods to compute the raw scattering values:
class YourTheory(ScatteringTheory):
def _raw_fields(self, positions, scatterer, medium_wavevec, medium_index, illum_polarization):
# Your code here
def _raw_scat_matrs(self, scatterer, pos, medium_wavevec, medium_index):
# Your code here
def _raw_cross_sections(self, scatterer, medium_wavevec, medium_index, illum_polarization):
# Your code here
You can get away with just defining one of _raw_scat_matrs or _raw_fields if you just want holograms, fields, or intensities. If you want scattering matrices you will need to implement _raw_scat_matrs, and if you want cross sections, you will need to implement _raw_cross_sections. We seperate out _raw_fields from _raw_scat_matrs because we want to provide a faster fields implementation for mie and multisphere (and you might want to for your theory).
You can look at the Mie theory in HoloPy for an example of calling Fortran functions to compute scattering (C functions will look similar from the python side) or DDA for an an example of calling out to an external command line tool by generating files and reading output files.
Adding a new inference model¶
To perform inference, you need a noise model. You can make a new noise model by inheriting from NoiseModel
. This class has all the machinery to compute likelihoods of observing data given some set of parameters and assuming Gaussian noise.
To implement a new model, you just need to implement one function: _forward
.
This function receives a dictionary of parameter values and a data shape schema (defined by detector_grid()
, for example) and needs to return simulated data of shape specified. See the _forward
function in AlphaModel
for an example of how to do this.
If you want to use some other noise model, you may need to override _lnlike and define the probablity given your uncertainty. You can reference _lnlike in NoiseModel
.
Running Tests¶
HoloPy comes with a suite of tests that ensure everything has been
built correctly and that it’s able to perform all of the calculations
it is designed to do. To run these tests, navigate to the root of the
package (e.g. /home/me/holopy
) and run:
python run_nose.py
HoloPy Tools¶
Holopy contains a number of tools to help you with common tasks when analyzing holograms. This page provides a summary of the tools available, while full descriptions can be found in the relevant code reference.
General Image Processing Tools¶
The tools described here are frequently used when analyzing holgrams. They are available from the holopy.core.process
namespace.
The normalize()
function divides an image by its average, returning an
image with a mean pixel value of 1. Note that this is the same normalization
convention used by HoloPy when calculating holograms with calc_holo()
.
Cropping an image introduces difficulties in keeping track of the relative
coordinates of features within an image and maintaining metadata. By using the
subimage()
function, the image origin is maintained in the cropped image,
so coordinate locations of features (such as a scatterer) remain unchanged.
Since holograms of particles usually take the form of concentric rings, the
location of a scatterer can usually be found by locating the apparent center(s)
of the image. Use center_find()
to locate one or more centers in an
image.
You can remove isolated dead pixels with zero intensity (e.g. for a background
division) by using zero_filter()
. This function replaces the dead pixel
with the average of its neighbours, and fails if adjacent pixels have zero
intensity.
The add_noise()
function allows you to add Gaussiancorrelated random
noise to a calculated image so that it more closely resembles experimental data.
To find gradient values at all points in an image, use image_gradient()
.
To simply remove a planar intensity gradient from an image, use
detrend()
. Note that this gives a mean pixel value of zero.
Frequency space analysis provides a powerful tool for working with images. Use
fft()
and ifft()
to perform fourier transforms and inverse fourier
transforms, respectively. These make use of scipy.fftpack
functions, but are
wrapped to correctly interpret HoloPy objects. HoloPy also includes a Hough
transform (hough()
) to help identify lines and other features in your
images.
Math Tools¶
HoloPy contains implementations of a few mathematical functions related to
scattering calculations. These functions are available from the
holopy.core.math
namespace.
To find the distance between two points, use cartesian_distance()
.
To rotate a set of points by arbitrary angles about the three coordinate axes,
use rotate_points()
. You can also calculate a rotation matrix with
rotation_matrix()
to save and use later.
To convert spherical coordinates into Cartesian coordinates, use
to_cartesian()
. To convert Cartesian coordinates into spherical
coordinates, use to_spherical()
.
When comparing data to a model, the chisquared and rsquared values provide
measures of goodnessoffit. You can access these through chisq()
and
rsq()
.
Concepts¶
Units¶
HoloPy does not enforce any particular set of units. As long as you are consistent, you can use any set of units, for example pixels, meters, or microns. So if you specify the wavelength of your red imaging laser as 658 then all other units (x, y, z position coordinates, particle radii, etc.) must also be specified in nanometers.
Coordinate System¶
For image data (data points arrayed in a regular grid in a single plane), HoloPy defaults to placing the origin, (0,0), at the top left corner as shown below. The xaxis runs vertically down, the yaxis runs horizontally to the right, and the zaxis points out of the screen, toward you. This corresponds to the way that images are treated by most computer software.
In sample space, we choose the z axis so that distances to objects from the camera/focal plane are positive (have positive z coordinates). The price we pay for this choice is that the propagation direction of the illumination light is then negative. In the image above, light travels from a source located in front of the screen, through a scatterer, and onto a detector behind the screen.
More complex detector geometries will define their own origin, or ask you to define one.
Rotations of Scatterers¶
Certain scattering calculations in HoloPy require specifying the orientation of a scatterer (such as a Janus sphere) relative to the HoloPy coordinate system. We do this in the most general way possible by specifying three Euler angles and a reference orientation. Rotating a scatterer initially in the reference orientation through the three Euler angles \(\alpha\), \(\beta\), and \(\gamma\) (in the active transformation picture) yields the desired orientation. The reference orientation is specified by the definition of the scatterer.
The Euler rotations are performed in the following way:
 Rotate the scatterer an angle \(\alpha\) about the HoloPy \(z\) axis.
 Rotate the scatterer an angle \(\beta\) about the HoloPy \(y\) axis.
 Rotate the scatterer an angle \(\gamma\) about the HoloPy \(z\) axis.
The sense of rotation is as follows: each angle is a rotation in the clockwise direction about the specified axis, viewed along the positive direction of the axis from the origin. This is the usual sense of how rotations are typically defined in math:
holopy package¶
Module contents¶
HoloPy is a set of tools for working with digital holograms and light scattering. It contains tools for working loading data and associating it with expermiental metadata, reconstructing holograms, calculating light scattering, fitting scattering models to data, and visualizing images and calculations.
Subpackages¶
holopy.core package¶
Subpackages¶
holopy.core.io package¶
Common entry point for holopy io. Dispatches to the correct load/save functions.

default_extension
(inf, defext='.h5')¶

get_example_data
(name)¶

get_example_data_path
(name)¶

load
(inf, lazy=False)¶ Load data or results
Parameters: inf (string) – String specifying an hdf5 file containing holopy data Returns: obj – The array object contained in the file Return type: xarray.DataArray

load_average
(filepath, refimg=None, spacing=None, medium_index=None, illum_wavelen=None, illum_polarization=None, normals=None, channel=None, image_glob='*.tif')¶ Average a set of images (usually as a background)
Parameters:  filepath (string or list(string)) – Directory or list of filenames or filepaths. If filename is a directory, it will average all images matching image_glob.
 refimg (xarray.DataArray) – reference image to provide spacing and metadata for the new image.
 spacing (float) – Spacing between pixels in the images. Used preferentially over refimg value if both are provided.
 medium_index (float) – Refractive index of the medium in the images. Used preferentially over refimg value if both are provided.
 illum_wavelen (float) – Wavelength of illumination in the images. Used preferentially over refimg value if both are provided.
 illum_polarization (listlike) – Polarization of illumination in the images. Used preferentially over refimg value if both are provided.
 normals (listlike) – Orientation of detector. Used preferentially over refimg value if both are provided.
 image_glob (string) – Glob used to select images (if images is a directory)
Returns: averaged_image – Image which is an average of images
Return type: xarray.DataArray

load_image
(inf, spacing=None, medium_index=None, illum_wavelen=None, illum_polarization=None, normals=None, channel=None, name=None)¶ Load data or results
Parameters:  inf (single or list of basestring or files) – File to load. If the file is a yaml file, all other arguments are ignored. If inf is a list of image files or filenames they are all loaded as a a timeseries hologram
 channel (int (optional)) – number of channel to load for a color image (in general 0=red, 1=green, 2=blue)
Returns: obj
Return type: The object loaded,
holopy.core.marray.Image
, or as loaded from yaml

pack_attrs
(a, do_spacing=False)¶

save
(outf, obj)¶ Save a holopy object
Will save objects as yaml text containing all information about the object unless outf is a filename with an image extension, in which case it will save an image, truncating metadata.
Parameters:  outf (basestring or file) – Location to save the object
 obj (
holopy.core.holopy_object.HoloPyObject
) – The object to save
Notes
Marray objects are actually saved as a custom yaml file consisting of a yaml header and a numpy .npy binary array. This is done because yaml’s saving of binary array is very slow for large arrays. HoloPy can read these ‘yaml’ files, but any other yaml implementation will get confused.

save_image
(filename, im, scaling='auto', depth=8)¶ Save an ndarray or image as a tiff.
Parameters:  im (ndarray or
holopy.image.Image
) – image to save.  filename (basestring) – filename in which to save image. If im is an image the function should default to the image’s name field if no filename is specified
 scaling ('auto', None, or (NoneInt, NoneInt)) – How the image should be scaled for saving. Ignored for float output. It defaults to auto, use the full range of the output format. Other options are None, meaning no scaling, or a pair of integers specifying the values which should be set to the maximum and minimum values of the image format.
 depth (8, 16 or 'float') – What type of image to save. Options other than 8bit may not be supported for many image types. You probably don’t want to save 8bit images without some kind of scaling.
 im (ndarray or

unpack_attrs
(a)¶
holopy.core.process package¶
Routines for image processing. Useful for preprocessing raw holograms prior to extracting final data or postprocessing reconstructions.
The centerfinder module is a group of functions for locating the centers of holographic ring patterns. The module can find the center of a singlesphere holographic pattern, a dimer holographic pattern, or the centers of multiple (wellseparated: clearly separate ring patterns with separate centers) single spheres or dimers. The intended use is for determining an initial parameter guess for hologram fitting.
We thank the Grier Group at NYU for suggesting the use of the Hough transform. For their independent implementation of a Houghbased holographic feature detection algorithm, see: http://physics.nyu.edu/grierlab/software/circletransform.pro For a case study and further reading, see: F. C. Cheong, B. Sun, R. Dreyfus, J. AmatoGrill, K. Xiao, L. Dixon & D. G. Grier, Flow visualization and flow cytometry with holographic video microscopy, Optics Express 17, 1307113079 (2009).

center_find
(image, centers=1, threshold=0.5, blursize=3.0)¶ Finds the coordinates of the center of a holographic pattern. The coordinates returned are in pixels (row number, column number). Intended for finding the center of single particle or dimer holograms which basically show concentric circles. The optional threshold parameter (between 0 and 1) gives a bound on what magnitude of gradients to include in the calculation. For example, threshold=.75 means ignore any gradients that are less than 75% of the maximum gradient in the image. The optional blursize parameter sets the size of a Gaussian filter that is applied to the image. This step improves accuracy when small features in the image have large gradients (e.g. dust particles on the camera). Without blurring, these features may be incorrectly identified as the hologram center. For best results, blursize should be set to the radius of features to be ignored, but smaller than the distance between hologram fringes. To skip blurring, set blursize to 0.
Parameters:  image (ndarray) – image to find the center(s) in
 centers (int) – number of centers to find
 threshold (float (optional)) – fraction of the maximum gradient below which all other gradients will be ignored (range 0.99)
 blursize (float (optional)) – radius (in pixels) of the Gaussian filter that is applied prior to Hough transform
Returns: res – row(s) and column(s) of center(s)
Return type: ndarray
Notes
When threshold is close to 1, the code will run quickly but may lack accuracy. When threshold is set to 0, the gradient at all pixels will contribute to finding the centers and the code will take a little bit longer.

hough
(col_deriv, row_deriv, centers=1, threshold=0.25)¶ Following the approach of a Hough transform, finds the pixel which the most gradients point towards or away from. Uses only gradients with magnitudes greater than threshold*maximum gradient. Once the pixel is found, uses a brightnessweighted average around that pixel to refine the center location to return. After the first center is found, the sourrounding area is blocked out and another brightest pixel is searched for if more centers are required.
Parameters:  col_deriv (numpy.ndarray) – ycomponent of image intensity gradient
 row_deriv (numpy.ndarray) – xcomponent of image intensity gradient
 centers (int) – number of centers to find
 threshold (float (optional)) – fraction of the maximum gradient below which all other gradients will be ignored (range 0.99)
Returns: res – row and column of center or centers
Return type: ndarray

image_gradient
(image)¶ Uses the Sobel operator as a numerical approximation of a derivative to find the x and y components of the image’s intensity gradient at each pixel.
Parameters: image (ndarray) – image to find the gradient of Returns:  gradx (ndarray) – xcomponents of intensity gradient
 grady (ndarray) – ycomponents of intensity gradient
Handles Fourier transforms of HoloPy images by using scipy’s fftpack. Tries to correctly interpret dimensions from xarray.

fft
(a, overwrite=False, shift=True)¶ More convenient Fast Fourier Transform
An easier to use fft function, it will pick the correct fft to do based on the shape of the Marray, and do the fftshift for you. This is intended for working with images, and thus for dimensions greater than 2 does slicewise transforms of each “image” in a multidimensional stack
Parameters:  a (ndarray) – The array to transform
 overwrite (bool) – Allow this function to overwrite the Marry you pass in. This may improve performance slightly. Default is not to overwrite
 shift (bool) – Whether to preform an fftshift on the Marry to give low frequences near the center as you probably expect. Default is to do the fftshift.
Returns: fta – The fourier transform of a
Return type: ndarray

ft_coord
(c)¶

ft_coords
(cs)¶

get_spacing
(c)¶

ifft
(a, overwrite=False, shift=True)¶ More convenient Inverse Fast Fourier Transform
An easier to use ifft function, it will pick the correct ifft to do based on the shape of the Marry, and do the fftshift for you. This is indendended for working with images, and thus for dimensions greater than 2 does slicewise transforms of each “image” in a multidimensional stack
Parameters:  a (ndarray) – The array to transform
 overwrite (bool) – Allow this function to overwrite the Marry you pass in. This may improve performance slightly. Default is not to overwrite
 shift (bool) – Whether to preform an fftshift on the Marry to give low frequences near the center as you probably expect. Default is to do the fftshift.
Returns: ifta – The inverse fourier transform of a
Return type: ndarray

ift_coord
(c)¶

ift_coords
(cs)¶

transform_metadata
(a, inverse)¶
Image enhancement through background subtraction, contrast adjustment, or detrending

add_noise
(image, noise_mean=0.1, smoothing=0.01, poisson_lambda=1000)¶ Add simulated noise to images. Intended for use with exact calculated images to make them look more like noisy ‘real’ measurements.
Real image noise usually has correlation, so we smooth the raw random variable. The noise_mean can be controlled independently of the poisson_lambda that controls the shape of the distribution. In general, you can stick with our default of a large poisson_lambda (ie for imaging conditions not near the shot noise limit).
Defaults are set to give noise vaguely similar to what we tend to see in our holographic imaging.
Parameters:  image (ndarray or Image) – The image to add noise to.
 smoothing (float) – Fraction of the image size to smooth by. Should in general be << 1
 poisson_lambda (float) – Used to compute the shape of the noise distribution. You can generally leave this at its default value unless you are simulating shot noise limited imaging.
Returns: noisy_image – A copy of the input image with noise added.
Return type: ndarray

bg_correct
(raw, bg, df=None)¶ Correct for noisy images by dividing by a background. The calculation used is (rawdf)/(bgdf).
Parameters:  raw (xarray.DataArray) – Image to be background divided.
 bg (xarray.DataArray) – background image recorded with the same optical setup.
 df (xarray.DataArray) – dark field image recorded without illumination.
Returns: corrected_image – A copy of the background divided input image.
Return type: xarray.DataArray

detrend
(image)¶ Remove linear trends from an image.
Performs a 2 axis linear detrend using scipy.signal.detrend
Parameters: image (ndarray) – Image to process Returns: image – Image with linear trends removed Return type: ndarray

normalize
(image)¶ Normalize an image (NumPy array) by dividing by the pixel average. This gives the image a mean value of 1.
Parameters: image (ndarray) – The array to normalize Returns: normalized_image – The normalized image Return type: ndarray

simulate_noise
(shape, mean=0.1, smoothing=0.01, poisson_lambda=1000)¶ Create an array of correlated noise. The noise_mean can be controlled independently of the poisson_lambda that controls the shape of the distribution. In general, you can stick with our default of a large poisson_lambda (ie for imaging conditions not near the shot noise limit).
Defaults are set to give noise vaguely similar to what we tend to see in our holographic imaging.
Parameters:  shape (int or array_like of ints) – shape of noise array
 smoothing (float) – Fraction of the image size to smooth by. Should in general be << 1
 poisson_lambda (float) – Used to compute the shape of the noise distribution. You can generally leave this at its default value unless you are simulating shot noise limited imaging.
Returns: noisy_image – A copy of the input image with noise added.
Return type: ndarray

subimage
(arr, center, shape)¶ Pick out a region of an image or other array
Parameters:  arr (numpy.ndarray) – The array to subimage
 center (tuple of ints or floats) – The desired center of the region, should have the same number of elements as the arr has dimensions. Floats will be rounded
 shape (int or tuple of ints) – Desired shape of the region. If a single int is given the region will be that dimension in along every axis. Shape should be even
Returns: sub – Subset of shape shape centered at center. For marrays, marray.origin will be set such that the upper left corner of the output has coordinates relative to the input.
Return type: numpy.ndarray or
RegularGrid
marray object

zero_filter
(image)¶ Search for and interpolate pixels equal to 0. This is to avoid NaN’s when a hologram is divided by a BG with 0’s.
Parameters: image (ndarray) – Image to process Returns: image – Image where pixels = 0 are instead given values equal to average of neighbors. dtype is the same as the input image Return type: ndimage
Module contents¶
Loading, saving, and basic processing of data.
holopy.core contains code to load images and holopy yamls into
marray
objects. It also contains the machinery for saving
all HoloPy objects as holopy yaml files. Finally, it provides some
basic mathematical operations, mostly as higher level wrappers around
numpy or scipy routines.
Main use cases are
 Image or other data file + metadata =>
Image
or otherMarray
object  Raw
Image
+ processing => processedImage
object  Any
HoloPyObject
from calculations or processing => achival yaml text or text/binary result
holopy.core.holopy_object module¶
Root class for all of holopy. This class provides serialization to and from yaml text file for all holopy objects.
yaml files are structured text files designed to be easy for humans to read and write but also easy for computers to read. HoloPy uses them to store information about experimental conditions and to describe analysis procedures.

class
HoloPyObject
¶ Bases:
holopy.core.holopy_object.Serializable
Ancestor class for all HoloPy classes.
HoloPy object’s purpose is to provide the machinery for saving to and loading from HoloPy yaml files

classmethod
from_yaml
(loader, node)¶

like_me
(filter_none=True, **kwargs)¶

classmethod
to_yaml
(dumper, data)¶

classmethod

class
Serializable
¶ Bases:
yaml.YAMLObject
Base class for any object that wants a nice clean yaml output

classmethod
to_yaml
(dumper, data)¶

classmethod

class
SerializableMetaclass
(name, bases, kwds)¶ Bases:
yaml.YAMLObjectMetaclass

ordered_dump
(dumper, tag, data)¶
holopy.core.math module¶

cartesian_distance
(p1, p2=[0, 0, 0])¶ Return the Cartesian distance between points p1 and p2.
Parameters: p2 (p1,) – Coordinates of point 1 and point 2 in Ndimensions Returns: dist – Cartesian distance between points p1 and p2 Return type: float64

chisq
(fit, data)¶ Calculate perpoint value of chisquared comparing a bestfit model and data.
Parameters:  fit (array_like) – Values of bestfit model to be compared to data
 data (array_like) – Data values to be compared to model
Returns: chisq – Chisquared per point
Return type: Notes
chisquared is defined as
\[\chi^2 = \frac{1}{N}\sum_{\textrm{points}} (\textrm{fit}  \textrm{data})^2\]where \(N\) is the number of data points.

rotate_points
(points, theta, phi, psi)¶ Rotate an array of points about Euler angles in a z, y, z convention.
Parameters:  points (arraylike (n,3)) – Set of points to be rotated
 phi, psi (theta,) – Euler rotation angles in z, y, z convention. These are not the same as angles in spherical coordinates.
Returns: rotated_points – Points rotated by Euler angles
Return type:

rotation_matrix
(alpha, beta, gamma, radians=True)¶ Return a 3D rotation matrix
Parameters:  beta, gamma (alpha,) – Euler rotation angles in z, y, z convention
 radians (boolean) – Default True; treats input angles in radians
Returns: rot – Rotation matrix. To rotate a column vector x, use np.dot(rot, x.)
Return type: array(3,3)
Notes
The Euler angles rotate a vector (in the active picture) by alpha clockwise about the fixed lab z axis, beta clockwise about the lab y axis, and by gamma about the lab z axis. Clockwise is defined as viewed from the origin, looking in the positive direction along an axis.
This breaks compatability with previous conventions, which were adopted for compatability with the passive picture used by SCSMFO.

rsq
(fit, data)¶ Calculate correlation coeffiction Rsquared comparing a bestfit model and data.
Parameters:  fit (array_like) – Values of bestfit model to be compared to data
 data (array_like) – Data values to be compared to model
Returns: rsq – Correlation coefficient Rsquared.
Return type: Notes
Rsquared is defined as
\[R^2 = 1  \frac{\sum_{\textrm{points}} (\textrm{data}  \textrm{fit})^2}{\sum_{\textrm{points}} (\textrm{data}  \bar{\textrm{data}})^2}\]where \(\bar{\textrm{data}}\) is the mean value of the data. If the model perfectly describes the data, \(R^2 = 1\).

to_cartesian
(r, theta, phi)¶ Returns Cartesian coordinates of a point given in spherical polar coordinates.
Parameters: theta, phi (r,) – Spherical polar coordinates of point. Returns: cartesian_coords – Dictionary of Cartesian coordinates of point with keys ‘x’, ‘y’, and ‘z’. Return type: dict

to_spherical
(x, y, z)¶ Return the spherical polar coordinates of a point in Cartesian coordinates.
Parameters: y, z (x,) – Cartesian coordinates of point Returns: spherical_coords – Dictionary of spherical polar coordinates (r, theta, phi) of point with keys ‘r’, ‘theta’, ‘phi’. Return type: dict Notes
theta is the polar angle measured from the z axis with range (0, pi). phi is the azimuthal angle with range (0, 2 pi).
holopy.core.metadata module¶
Classes for defining metadata about experimental or calculated results.

copy_metadata
(old, new, do_coords=True)¶

data_grid
(arr, spacing=None, medium_index=None, illum_wavelen=None, illum_polarization=None, normals=None, name=None, z=0)¶ Create a set of detector points along with other experimental metadata.
Returns: Return type: DataArray object Notes
Use of higherlevel detector_grid() and detector_points() functions is recommended.

default_norms
(coords, n)¶

detector_grid
(shape, spacing, normals=None, name=None)¶ Create a rectangular grid of pixels to represent a detector on which scattering calculations are to be performed.
Parameters:  shape (int or listlike (2)) – If int, detector is a square grid of shape x shape points. If array_like, detector has shape[0] rows and shape[1] columns.
 spacing (int or listlike (2)) – If int, distance between square detector pixels. If array_like, spacing[0] between adjacent rows and spacing[1] between adjacent columns.
 normals (listlike or None) – If listlike, detector orientation.
 name (string) –
Returns: grid – DataArray of zeros with coordinates calculated according to shape and spacing
Return type: DataArray object
Notes
Typically used to define a set of points to represent the pixels of a digital camera in scattering calculations.

detector_points
(coords={}, x=None, y=None, z=None, r=None, theta=None, phi=None, normals='auto', name=None)¶ Returns a onedimensional set of detector coordinates at which scattering calculations are to be done.
Parameters:  coords (dict, optional) – Dictionary of detector coordinates. Default: empty dictionary. Typical usage should not pass this argument, giving other parameters (Cartesian x, y, and z or polar r, theta, and phi coordinates) instead.
 y (x,) – Cartesian x and y coordinates of detectors.
 z (int or array_like, optional) – Cartesian z coordinates of detectors. If not specified, assume z = 0.
 r (int or array_like, optional) – Spherical polar radial coordinates of detectors. If not specified, assume r = infinity (farfield).
 theta (int or array_like, optional) – Spherical polar coordinates (polar angle from z axis) of detectors.
 phi (int or array_like, optional) – Spherical polar azimuthal coodinates of detectors.
 normals (string, optional) – Default behavior: normal in +z direction for Cartesian coordinates, r direction for polar coordinates. Nondefault behavior not currently implemented.
 name (string) –
Returns: grid – DataArray of zeros with calculated coordinates.
Return type: DataArray object
Notes
Specify either the Cartesian or the polar coordinates of your detector. This may be helpful for modeling static light scattering calculations. Use detector_grid() to specify coordinates of a grid of pixels (e.g., a digital camera.)

flat
(a)¶

from_flat
(a)¶

get_extents
(im)¶

get_spacing
(im)¶

get_values
(a)¶

make_coords
(shape, spacing, z=0)¶

primdim
(a)¶

sphere_coords
(a, origin=(0, 0, 0), wavevec=1)¶

to_vector
(c)¶

update_metadata
(a, medium_index=None, illum_wavelen=None, illum_polarization=None, normals=None)¶ Returns a copy of an image with updated metadata in its ‘attrs’ field.
Parameters:  a (xarray.DataArray) – image to update.
 medium_index (float) – Updated refractive index of the medium in the image.
 illum_wavelen (float) – Updated wavelength of illumination in the image.
 illum_polarization (listlike) – Updated polarization of illumination in the image.
 normals (listlike) – Updated detector orientation of the image.
Returns: b – copy of input image with updated metadata. The ‘normals’ field is not allowed to be empty.
Return type: xarray.DataArray
holopy.core.utils module¶
Misc utility functions to make coding more convenient

dict_without
(d, keys)¶ Exclude a list of keys from a dictionary.
Silently ignores any key in keys that is not in the dict (this is intended to be used to make sure a dict does not contain specific keys) :param d: The dictionary to operate on :type d: dict :param keys: The keys to exclude :type keys: list(string) :param returns: A copy of dict without any of the specified keys :type returns: d2

ensure_array
(x)¶

ensure_listlike
(x)¶

is_none
(o)¶ Check if something is None.
This can’t be done with a simple is check anymore because numpy decided that array is None should do an element wise comparison.
Parameters: o (object) – Anything you want to see if is None

mkdir_p
(path)¶ Equivalent to mkdir p at the shell, this function makes a directory and its parents as needed, silently doing nothing if it exists.

repeat_sing_dims
(indict, keys='all')¶
holopy.fitting package¶
Module contents¶
Fit models of scattering to data
Make precision measurements of a scattering system by fitting a model of it to data
The fitting module is used to:
holopy.fitting.fit module¶
Routines for fitting a hologram to an exact solution

class
FitResult
(parameters, scatterer, chisq, rsq, converged, time, model, minimizer, minimization_details)¶ Bases:
holopy.core.holopy_object.HoloPyObject
The results of a fit.
You should not make objects of this class directly, they will be given to you by
fit()
Parameters:  parameters (array(float)) – The fitted values for each parameter
 scatterer (
Scatterer
) – The best fit scatterer  chisq (float) – The \(\chi^2\) goodness of fit
 rsq (float) – The \(R^2\) goodness of fit
 converged (bool) – Did the minimizer converge
 time (float) – Time in seconds the fit took
 minimizer (
Minimizer
) – Them minimizer used in the fit  minimization_details (object) – Additional information returned by the minimizer about the minimization

alpha
¶

fitted_holo
(schema)¶

classmethod
from_summary
(summary, scatterer_cls)¶ Build a FitResult from a summary.
The returned FitResult will be incomplete because summaries do not contain all of the minimizer information that full results do
Parameters: summary (dict) – A dict as from cls.summary containing information about a fit.

next_model
()¶ Construct a model to fit the next frame in a time series

niter
()¶

summary
()¶ Put just the essential components of a fit result in a dictionary

summary_misc
= ['rsq', 'chisq', 'time', 'converged', 'niter']¶

fit
(model, data, minimizer=<class 'holopy.fitting.minimizer.Nmpfit'>, random_subset=None)¶ fit a model to some data
Parameters:  model (
Model
object) – A model describing the scattering system which leads to your data and the parameters to vary to fit it to the data  data (
Marray
object) – The data to fit  minimizer ((optional)
Minimizer
) – The minimizer to use to do the fit  random_subset (float (optional)) – Fit only a randomly selected fraction of the data points in data
Returns: result – an object containing the best fit parameters and information about the fit
Return type:  model (

make_subset_data
(data, random_subset=None, pixels=None, return_selection=False)¶
holopy.fitting.minimizer module¶
Interfaces to minimizers. Classes here provide a common interface to a variety of third party minimizers.

class
Minimizer
¶ Bases:
holopy.core.holopy_object.HoloPyObject
Common interface to all minimizers holopy supports

minimize
(parameters, cost_func)¶ Find the best solution to an optimization problem
Parameters:  parameters (list of
Parameter
objects) – Parameters to vary in the model  cost_func (function) – A function taking parameters as arguments that returns the residual for the minimization problem
 parameters (list of

pars_from_minimizer
(parameters, values)¶


class
Nmpfit
(quiet=False, ftol=1e10, xtol=1e10, gtol=1e10, damp=0, maxiter=100)¶ Bases:
holopy.fitting.minimizer.Minimizer
LevenbergMarquardt minimizer, from Numpy/Python translation of Craig Markwardt’s mpfit.pro.
Parameters:  quiet (Boolean) – If True, suppress output on minimizer convergence.
 ftol (float) – Convergence criterion for minimizer: converges if actual and predicted relative reductions in chi squared <= ftol
 xtol (float) – Convergence criterion for minimizer: converges if relative error between two LevenbergMarquardt iterations is <= xtol
 gtol (float) – Convergence criterion for minimizer: converges if absolute value of cosine of angle between vector of cost function evaluated at current solution for minimized parameters and any column of the Jacobian is <= gtol
 damp (float) – If nonzero, residuals larger than damp will be replaced by tanh. See nmpfit documentation.
 maxiter (int) – Maximum number of LevenbergMarquardt iterations to be performed.
Notes
See nmpfit documentation for further details. Not all functionalities of nmpfit are implemented here: in particular, we do not allow analytical derivatives of the residual function, which is impractical and/or impossible to calculate for holograms. If you want to weight the residuals, you need to supply a custom residual function.

minimize
(parameters, cost_func, debug=False)¶ Find the best solution to an optimization problem
Parameters:  parameters (list of
Parameter
objects) – Parameters to vary in the model  cost_func (function) – A function taking parameters as arguments that returns the residual for the minimization problem
 parameters (list of
holopy.fitting.model module¶
Classes for defining models of scattering for fitting

class
BaseModel
(scatterer, medium_index=None, illum_wavelen=None, illum_polarization=None, theory='auto')¶ Bases:
holopy.core.holopy_object.HoloPyObject

get_par
(name, pars, schema=None, default=None)¶

get_pars
(names, pars, schema=None)¶

par
(name, schema=None, default=None)¶

parameters
¶


class
Model
(scatterer, calc_func, medium_index=None, illum_wavelen=None, illum_polarization=None, theory='auto', alpha=None, use_random_fraction=None, constraints=[])¶ Bases:
holopy.fitting.model.BaseModel
Representation of a model to fit to data
Parameters:  alpha (float or Parameter) – Extra scaling parameter, hopefully this will be removed by improvements in our theory soon.
 constraints (function or list of functions) – One or a list of constraint functions. A constraint function should take a scaterer as an argument and return False if you wish to disallow that scatterer (usually because it is unphysical for some reason)

get_alpha
(pars=None)¶

guess
¶

guess_dict
¶

residual
(pars, data)¶

class
ParameterizedObject
(obj)¶ Bases:
holopy.fitting.model.Parametrization
Specify parameters for a fit by including them in an object
Parameters are named automatically from their position in the object
Parameters: obj ( scatterer
) – Object containing parameters specifying any values vary in the fit. It can also include numbers for any fixed values
guess
¶

make_from
(parameters)¶


class
Parametrization
(make_scatterer, parameters)¶ Bases:
holopy.core.holopy_object.HoloPyObject
Description of free parameters and how to make a scatterer from them
Parameters:  make_scatterer (function) – A function which should take the Parametrization parameters by name as keyword arguments and return a scatterer
 parameters (list) – The list of parameters for this Parametrization

guess
¶

make_from
(parameters)¶

limit_overlaps
(fraction=0.1)¶ Generator for constraint prohibiting overlaps beyond a certain tolerance
Parameters: fraction (float) – Fraction of the sphere diameter that the spheres should be allowed to overlap by Returns: constraint – A function which tests scatterers to see if the exceed the specified tolerance Return type: function (scatterer > bool)

tied_name
(name1, name2)¶
holopy.fitting.parameter module¶
Classes for describing free parameters in fitting models

class
ComplexParameter
(real, imag, name=None)¶ Bases:
holopy.fitting.parameter.Parameter
A complex free parameter
ComplexParameters have a real and imaginary part which can (potentially) vary separately.
Parameters:  imag (real,) – The real and imaginary parts of this parameter. Assign floats to fix that portion or parameters to allow it to vary. The parameters must be purely real. You should omit name’s for the parameters; ComplexParameter will name them
 name (string) – Short descriptive name of the ComplexParameter. Do not provide this if using a ParameterizedScatterer, a name will be assigned based its position within the scatterer.

guess
¶

class
Parameter
(guess=None, limit=None, name=None, **kwargs)¶ Bases:
holopy.core.holopy_object.HoloPyObject
Describe a free parameter in a fitting model
Parameters:  guess ((optional) float) – Your inital guess for this parameter
 limit ((optional) float or (float, float)) – Describe how the minimizer should allow a parameter to vary. A single value here fixes the parameter to that value, a pair limits the parameter to vary between (min, max)
 name ((optional) string) – A short descriptive name of the parameter. Do not provide this if using
the parameter in a
ParameterizedScatterer
, it will assign a name based on the Parameter’s position within the scatterer  **kwargs (varies) – Additional information made available to the minimizer. This can be used for things like step sizes.

fixed
¶

scale
(physical)¶ Scales parameters to approximately unity
Parameters: physical (np.array(dtype=float)) – Returns: scaled Return type: np.array(dtype=float)

unscale
(scaled)¶ Inverts scale’s transformation
Parameters: scaled (np.array(dtype=float)) – Returns: physical Return type: np.array(dtype=float)
holopy.inference package¶
Module contents¶
holopy.inference.noise_model module¶

class
AlphaModel
(scatterer, noise_sd, alpha, medium_index=None, illum_wavelen=None, illum_polarization=None, theory='auto')¶

class
NoiseModel
(scatterer, noise_sd, medium_index=None, illum_wavelen=None, illum_polarization=None, theory='auto')¶ Bases:
holopy.fitting.model.BaseModel
Model probabilites of observing data
Compute probabilities that observed data could be explained by a set of scatterer and observation parameters.

lnlike
(par_vals, data)¶

lnposterior
(par_vals, data)¶

lnprior
(par_vals)¶

holopy.inference.prior module¶

class
BoundedGaussian
(mu, sd, lower_bound=inf, upper_bound=inf, name=None)¶ Bases:
holopy.inference.prior.Gaussian

lnprob
(p)¶

sample
(size=None)¶


class
Gaussian
(mu, sd, name=None)¶ Bases:
holopy.inference.prior.Prior

guess
¶

lnprob
(p)¶

prob
(p)¶

sample
(size=None)¶


class
Prior
(guess=None, limit=None, name=None, **kwargs)¶

class
Uniform
(lower_bound, upper_bound, name=None)¶ Bases:
holopy.inference.prior.Prior

guess
¶

interval
¶

lnprob
(p)¶

sample
(size=None)¶


make_center_priors
(im, z_range_extents=5, xy_uncertainty_pixels=1, z_range_units=None)¶ Make sensible default priors for the center of a sphere in a hologram
Parameters:  im (xarray) – The image you wish to make priors for
 z_range_extents (float (optional)) – What range to extend a uniform prior for z over, measured in multiples of the total extent of the image. The default is 5 times the extent of the image, a large range, but since tempering is quite good at refining this, it is safer to just choose a large range to be sure to include the correct value.
 xy_uncertainty_pixels (float (optional)) – The number of pixels of uncertainty to assume for the centerfinder. The default is 1 pixel, and this is probably correct for most images.
 z_range_units (float) – Specify the range of the z prior in your data units. If this is provided, z_range_extents is ignored.

updated
(prior, v, extra_uncertainty=0)¶ Update a prior from a posterior
Parameters:  v (UncertainValue) – The new value, usually from an mcmc result
 extra_uncertainty (float) – provide a floor for uncertainty (sd) of the new parameter
holopy.inference.result module¶
Results of sampling

class
SamplingResult
(dataset, model, strategy)¶ Bases:
holopy.core.holopy_object.HoloPyObject

MAP
¶

data
¶

lnprobs
¶

mean
¶

median
¶

samples
¶

sigma_intervals
(sigmas=[2, 1, 1, 2])¶

values
(sigma_interval=1)¶


class
TemperedSamplingResult
(end_result, stage_results, strategy)¶ Bases:
holopy.inference.result.SamplingResult

dataset
¶

model
¶


class
UncertainValue
(value, plus, minus=None, n_sigma=1, kind='MAP')¶ Bases:
holopy.core.holopy_object.HoloPyObject
Represent an uncertain value
Parameters:

autocorr_from_sentinal
(d)¶

autocorr_to_sentinal
(d)¶

get_stage_names
(inf)¶
holopy.inference.sample module¶
Sample posterior probabilities given model and data

class
EmceeStrategy
(nwalkers=100, pixels=2000, threads='auto', cleanup_threads=True, seed=None)¶ Bases:
holopy.core.holopy_object.HoloPyObject

make_guess
(parameters)¶

sample
(model, data, nsamples, walker_initial_pos=None)¶


class
TemperedStrategy
(next_initial_dist=<function sample_one_sigma_gaussian>, nwalkers=100, min_pixels=50, max_pixels=1000, threads='auto', stages=3, stage_len=30, seed=None)¶ Bases:
holopy.inference.sample.EmceeStrategy

sample
(model, data, nsamples)¶


autothreads
(threads='auto', quiet=False)¶

emcee_lnprobs_DataArray
(sampler)¶

emcee_samples_DataArray
(sampler, parameters)¶

get_acor
(sampler)¶

sample_emcee
(model, data, nwalkers, nsamples, walker_initial_pos, threads='auto', cleanup_threads=True, seed=None)¶

sample_emcee_autocorr
(model, data, nwalkers, independent_samples, walker_initial_pos, estimated_autocorr, threads='auto')¶

sample_one_sigma_gaussian
(result)¶

tempered_sample
(model, data, nwalkers=100, min_pixels=50, max_pixels=2000, samples=600, next_initial_dist=<function sample_one_sigma_gaussian>, stages=3, stage_len=30, seed=None, threads='auto')¶
holopy.propagation package¶
Module contents¶
holopy.propagation.convolution_propagation module¶
Code to propagate objects/waves using scattering models.

propagate
(data, d, medium_index=None, illum_wavelen=None, cfsp=0, gradient_filter=False)¶ Propagates a hologram along the optical axis
Parameters:  data (
Image
orVectorGrid
) – Hologram to propagate  d (float or list of floats) – Distance to propagate, in meters, or desired schema. A list tells to propagate to several distances and return the volume
 cfsp (integer (optional)) – cascaded freespace propagation factor. If this is an integer > 0, the transfer function G will be calculated at d/csf and the value returned will be G**csf. This helps avoid artifacts related to the limited window of the transfer function
 gradient_filter (float) – For each distance, compute a second propagation a distance gradient_filter away and subtract. This enhances contrast of rapidly varying features
Returns: data – The hologram progagated to a distance d from its current location.
Return type: Image
orVolume
 data (

trans_func
(schema, d, med_wavelen, cfsp=0, gradient_filter=0)¶ Calculates the optical transfer function to use in reconstruction
This routine uses the analytical form of the transfer function found in in Kreis [1]. It can optionally do cascaded freespace propagation for greater accuracy [2], although the code will run slightly more slowly.
Parameters:  shape ((int, int)) – maximum dimensions of the transfer function
 spacing ((float, float)) – the spacing between points is the grid to calculate
 wavelen (float) – the wavelength in the medium you are propagating through
 d (float or list of floats) – reconstruction distance. If list or array, this function will return an array of transfer functions, one for each distance
 cfsp (integer (optional)) – cascaded freespace propagation factor. If this is an integer > 0, the transfer function G will be calculated at d/csf and the value returned will be G**csf.
 gradient_filter (float (optional)) – Subtract a second transfer function a distance gradient_filter from each z
Returns: trans_func – The calculated transfer function. This will be at most as large as shape, but may be smaller if the frequencies outside that are zero
Return type: np.ndarray
References
[1] Kreis, Handbook of Holographic Interferometry (Wiley, 2005), equation 3.79 (page 116) [2] Kreis, Optical Engineering 41(8):1829, section 5
holopy.propagation.point_source_propagate module¶

interpolate2D
(data, i, j, fill=None)¶ Interpolates values from a 2D array (data) given noninteger indecies i and j. If [i,j] is outside of the shape of data, fill is returned. If fill=None, the value of the closest edge pixel to (i,j) is used.

ps_propagate
(data, d, L, beam_c, out_schema=None)¶ Propagates light back through a hologram that was taken using a diverging reference beam. Only propagation through media with refractive index 1 is supported. This is a wrapper function for ps_propagate_plane() This function can handle a single reconstruction plane or a volume.
Based on the algorithm described in Manfred H. Jericho and H. Jurgen Kreuzer, “Point Source Digital InLine Holographic Microscopy,” Chapter 1 of Coherent Light Microscopy, Springer, 2010. http://link.springer.com/chapter/10.1007%2F9783642158131_1
data is a holopy Xarray. It is the hologram to reconstruct. Must be square. The pixel spacing must also be square. d = distance from pinhole to reconstructed image, in meters (this is z in Jericho and Kreuzer). Can be a scalar or a 1D list or array. L = distance from screen to pinhole, in meters beam_c = [x,y] coodinates of beam center, in pixels out_schema = size of output image and pixel spacing, default is the schema of data.
returns an image(volume) corresponding to the reconstruction at plane(s) d.

ps_propagate_plane
(data, d, L, beam_c, out_schema=None, old_Ip=False)¶ Propagates light back through a hologram that was taken using a diverging reference beam. Propataion can be to one plane only. Only propagation through media with refractive index 1 is supported.
Based on the algorithm described in Manfred H. Jericho and H. Jurgen Kreuzer, “Point Source Digital InLine Holographic Microscopy,” Chapter 1 of Coherent Light Microscopy, Springer, 2010. http://link.springer.com/chapter/10.1007%2F9783642158131_1
data is a holopy Xarray. It is the hologram to reconstruct. Must be square. The pixel spacing must also be square. d = distance from pinhole to reconstructed image, in meters (this is z in Jericho and Kreuzer). Must be a scalar. L = distance from screen to pinhole, in meters beam_c = [x,y] coodinates of beam center, in pixels out_schema = size of output image and pixel spacing, default is the schema of data. if Ip == True, returns Ip to be used on calculations in the stack if Ip == False compute reconstructed image as normal if Ip is an image, use this to speed up calculations
returns an image(volume) corresponding to the reconstruction at plane(s) d.
holopy.scattering package¶
Subpackages¶
holopy.scattering.scatterer package¶
Modules for defining different types of scatterers, including scattering primitives such as Spheres, and more complex objects such as Clusters.
Defines cylinder scatterers.

class
Bisphere
(n=None, h=None, d=None, center=None, rotation=(0, 0, 0))¶ Bases:
holopy.scattering.scatterer.scatterer.CenteredScatterer
Scattering object representing bisphere scatterers
Parameters:  n (complex) – Index of refraction
 h (distance between centers) –
 d (diameter) –
 center (3tuple, list or numpy array) – specifies coordinates of center of the scatterer
 rotation (3tuple, list or numpy.array) – specifies the Euler angles (alpha, beta, gamma) in radians defined in adda manual section 8.1
Defines capsule scatterers.

class
Capsule
(n=None, h=None, d=None, center=None, rotation=(0, 0, 0))¶ Bases:
holopy.scattering.scatterer.scatterer.CenteredScatterer
A cylinder with semispherical caps.
A particle with no rotation has its long axis pointing along +z, specify other orientations by euler angle rotations from that reference.
Parameters:  n (complex) – Index of refraction
 h (height of cylinder) –
 d (diameter) –
 center (3tuple, list or numpy array) – specifies coordinates of center of the scatterer
 rotation (3tuple, list or numpy.array) – specifies the Euler angles (alpha, beta, gamma) in radians

indicators
¶
Defines Scatterers, a scatterer that consists of other scatterers, including scattering primitives (e.g. Sphere) or other Scatterers scatterers (e.g. two trimers).

class
Scatterers
(scatterers=None)¶ Bases:
holopy.scattering.scatterer.scatterer.Scatterer
Contains optical and geometrical properties of a a composite scatterer. A Scatterers can consist of multiple scattering primitives (e.g. Sphere) or other Scatterers scatterers.

scatterers
¶ list – List of scatterers that make up this object
Notes
Stores information about components in a tree. This is the most generic container for a collection of scatterers.

add
(scatterer)¶

classmethod
from_parameters
(parameters)¶

get_component_list
()¶

in_domain
(points)¶

index_at
(point)¶

parameters
¶

rotated
(ang1, ang2=None, ang3=None)¶

Do Constructive Solid Geometry (CSG) with scatterers. Currently only useful with the DDA th

class
CsgScatterer
(s1, s2)¶ Bases:
holopy.scattering.scatterer.scatterer.Scatterer

bounds
¶

rotated
(alpha, beta, gamma)¶


class
Difference
(s1, s2)¶ Bases:
holopy.scattering.scatterer.csg.CsgScatterer

bounds
¶

in_domain
(points)¶


class
Intersection
(s1, s2)¶ Bases:
holopy.scattering.scatterer.csg.CsgScatterer

in_domain
(points)¶


class
Union
(s1, s2)¶ Bases:
holopy.scattering.scatterer.csg.CsgScatterer

in_domain
(points)¶

Defines cylinder scatterers.

class
Cylinder
(n=None, h=None, d=None, center=None, rotation=(0, 0, 0))¶ Bases:
holopy.scattering.scatterer.scatterer.CenteredScatterer
Scattering object representing cylinder scatterers
Parameters:  n (complex) – Index of refraction
 h (height of cylinder) –
 d (diameter) –
 center (3tuple, list or numpy array) – specifies coordinates of center of the scatterer
 rotation (3tuple, list or numpy.array) – specifies the Euler angles (alpha, beta, gamma) in radians defined in adda manual section 8.1
Defines ellipsoidal scatterers.

class
Ellipsoid
(n=None, r=None, center=None, rotation=(0, 0, 0))¶ Bases:
holopy.scattering.scatterer.scatterer.CenteredScatterer
Scattering object representing ellipsoidal scatterers
Parameters:  n (complex) – Index of refraction
 r (float or (float, float, float)) – x, y, z semiaxes of the ellipsoid
 center (3tuple, list or numpy array) – specifies coordinates of center of the scatterer
 rotation (3tuple, list or numpy.array) – specifies the Euler angles (alpha, beta, gamma) in radians defined in adda manual section 8.1

indicators
¶ NOTE – Ellipsoid indicators does not currently apply rotations

all_numbers
(x)¶

isnumber
(x)¶
Defines two types of Janus (two faced) Spheres as scattering primitives.

class
JanusSphere_Tapered
(n=None, r=None, rotation=(0, 0), center=None)¶ Bases:
holopy.scattering.scatterer.scatterer.CenteredScatterer

indicators
¶


class
JanusSphere_Uniform
(n=None, r=None, rotation=(0, 0, 0), center=None)¶ Bases:
holopy.scattering.scatterer.scatterer.CenteredScatterer

indicators
¶

The abstract base class for all scattering objects

class
CenteredScatterer
(center=None)¶ Bases:
holopy.scattering.scatterer.scatterer.Scatterer

classmethod
from_parameters
(parameters)¶ Create a Scatterer from a dictionary of parameters
Parameters: parameters (dict or list) – Parameters for a scatterer. This should be of the form returned by Scatterer.parameters. Returns: scatterer – A scatterer with the given parameter values Return type: Scatterer class

classmethod

class
Indicators
(functions, bound=None)¶ Bases:
holopy.core.holopy_object.HoloPyObject
Class holding functions describing a scatterer
One or more functions (one per domain) that take Nx3 arrays of points and return a boolean array of membership in each domain. More than one indicator is allowed to return true for a given point, in that case the point is considered a member of the first domain with a true value.

class
Scatterer
(indicators, n, center)¶ Bases:
holopy.core.holopy_object.HoloPyObject
Base class for scatterers

bounds
¶

contains
(points)¶

guess
()¶

in_domain
(points)¶ Tell which domain of a scatterer points are in
Parameters: points (np.ndarray (Nx3)) – Point or list of points to evaluate Returns: domain – The domain of each point. Domain 0 means not in the particle Return type: np.ndarray (N)

index_at
(points, background=0)¶

num_domains
¶

translated
(coord1, coord2=None, coord3=None)¶ Make a copy of this scatterer translated to a new location
Parameters: y, z (x,) – Value of the translation along each axis Returns: translated – A copy of this scatterer translated to a new location Return type: Scatterer

voxelate
(spacing, medium_index=0)¶ Represent a scatterer by discretizing into voxels
Parameters: Returns: voxelation – An array with refractive index at every pixel
Return type: np.ndarray

voxelate_domains
(spacing)¶

x
¶

y
¶

z
¶


bound_union
(d1, d2)¶

find_bounds
(indicator)¶ Finds the bounds needed to contain an indicator function
Notes
Will probably determine incorrect bounds for functions which are not convex
Defines Sphere, a scattering primitive

class
LayeredSphere
(n=None, t=None, center=None)¶ Bases:
holopy.scattering.scatterer.sphere.Sphere
Alternative description of a sphere where you specify layer thicknesses instead of radii

n
¶ list of complex – Index of each each layer

t
¶ list of float – Thickness of each layer

center
¶ length 3 listlike – specifies coordinates of center of sphere

r
¶


class
Sphere
(n=None, r=0.5, center=None)¶ Bases:
holopy.scattering.scatterer.scatterer.CenteredScatterer
Contains optical and geometrical properties of a sphere, a scattering primitive.
This can be a multiple layered sphere by making r and n lists.

n
¶ complex or list of complex – index of refraction of each layer of the sphere

r
¶ float or list of float – radius of the sphere or outer radius of each sphere.

center
¶ length 3 listlike – specifies coordinates of center of sphere

guess
()¶

indicators
¶

like_me
(**overrides)¶

num_domains
¶

rotated
(alpha, beta, gamma)¶

Defines Spheres, a Scatterers scatterer consisting of Spheres

class
Spheres
(scatterers, warn=True)¶ Bases:
holopy.scattering.scatterer.composite.Scatterers
Contains optical and geometrical properties of a cluster of spheres.

spheres
¶ list of Spheres – Spheres which will make up the cluster
Notes

add
(scatterer)¶

center
¶

centers
¶

guess
()¶

largest_overlap
()¶

n
¶

n_imag
¶

n_real
¶

overlaps
¶

r
¶

x
¶

y
¶

z
¶

Defines spheroidal scatterers.

class
Spheroid
(n=None, r=None, rotation=(0, 0, 0), center=None)¶ Bases:
holopy.scattering.scatterer.scatterer.CenteredScatterer
Scattering object representing spheroidal scatterers

n
¶ complex – Index of refraction

r
¶ (float, float) – length of xy and z semiaxes of the spheroid

rotation
¶ 3tuple, list or numpy array – specifies the Euler angles (alpha, beta, gamma) in radians

center
¶ 3tuple, list or numpy array – specifies coordinates of center of the scatterer

indicators
¶

holopy.scattering.theory package¶
Fortran extension module for calculating cluster holograms using tmatrix scattering theory.
Compute special functions needed for the computation of scattering coefficients in the LorenzMie scattering solution and related problems such as layered spheres.
These functions are not to be used for calculations at each field point. Rather, they should be used once for the calculation of scattering coefficients, which then get passed to faster Fortran code for field calculations.
Papers referenced herein:
D. W. Mackowski, R. A. Altenkirch, and M. P. Menguc, “Internal absorption cross sections in a stratified sphere,” Applied Optics 29, 15511559, (1990).
Yang, “Improved recursive algorithm for light scattering by a multilayered sphere,” Applied Optics 42, 17101720, (1993).

Qratio
(z1, z2, nstop, dns1=None, dns2=None, eps1=0.001, eps2=1e16)¶ Calculate ratio of RiccatiBessel functions defined in [Yang2003] eq. 23 by up recursion.
Notes
Logarithmic derivatives calculated automatically if not specified. Lentz continued fraction algorithm used to start downward recursion for logarithmic derivatives.

R_psi
(z1, z2, nmax, eps1=0.001, eps2=1e16)¶ Calculate ratio of RiccatiBessel function psi: psi(z1)/psi(z2).
Notes
See [Mackowski1990] eqns. 6566. Uses Lentz continued fraction algorithm for logarithmic derivatives.

log_der_1
(z, nmx, nstop)¶ Computes logarithmic derivative of RiccatiBessel function psi_n(z) by downward recursion as in BHMIE.
Parameters:  z (complex argument) –
 nmx (order from which downward recursion begins.) –
 nstop (integer, maximum order) –
Notes
psi_n(z) is related to the spherical Bessel function j_n(z). Consider implementing Lentz’s continued fraction method.

log_der_13
(z, nstop, eps1=0.001, eps2=1e16)¶ Calculate logarithmic derivatives of RiccatiBessel functions psi and xi for complex arguments. RiccatiBessel conventions follow Bohren & Huffman.
See Mackowski et al., Applied Optics 29, 1555 (1990).
Parameters:  z (complex number) –
 nstop (maximum order of computation) –
 eps1 (underflow criterion for Lentz continued fraction for Dn1) –
 eps2 (convergence criterion for Lentz continued fraction for Dn1) –
MieScatLib.py
Library of code to do Mie scattering calculations.

asymmetry_parameter
(al, bl)¶ Calculate asymmetry parameter of scattered field.
Parameters: bn (an,) – coefficient arrays from Mie solution Returns: Return type: float Notes
See discussion on Bohren & Huffman p. 120. The output of this function omits the prefactor of 4/(x^2 Q_sca).

cross_sections
(al, bl)¶ Calculates scattering and extinction cross sections given arrays of Mie scattering coefficients an and bn.
Parameters: bn (an,) – coefficient arrays from Mie solution Returns: Scattering, extinction, and radar backscattering cross sections Return type: ndarray(3) Notes
See Bohren & Huffman eqns. 4.61 and 4.62. The output omits a scaling prefactor of 2 * pi / k^2.

internal_coeffs
(m, x, n_max, eps1=0.001, eps2=1e16)¶ Calculate internal Mie coefficients c_n and d_n given relative index, size parameter, and maximum order of expansion.
Parameters: docstring for scatcoeffs (See) – Returns: Internal coefficients c_n and d_n Return type: ndarray(2,n) complex Notes
Follow Bohren & Huffman’s convention. Note that van de Hulst and Kerker have different conventions (labeling of c_n and d_n and factors of m) for their internal coefficients.

nstop
(x)¶ Calculate maximum expansion order of LorenzMie solution.
Parameters: x (float) – Particle size parameter Returns: nstop Return type: int Notes
Criterion taken from [Wiscombe1980].

scatcoeffs
(m, x, nstop, eps1=0.001, eps2=1e16)¶ Calculate expansion coefficients for scattered field in LorenzMie solution.
Parameters:  m (complex) – Sphere relative refractive index (n_sphere / n_medium)
 x (float) – Sphere size parameter (k_med * a)
 nstop (int) – Maximum order of scattered field expansion
 eps1 (float, optional) – In Lentz continued fraction algorithm for logarithmic derivative D_n(z), value of continued fraction numerator or denominator triggering illconditioning workaround.
 eps2 (float, optional) – Convergence criterion for Lentz continued fraction algorithm
Returns: Scattering coefficients a_n and b_n
Return type: Notes
Uses formula for scattering coefficients based on logarithmic derivative D_n(z) of spherical Bessel function psi_n(z). See [Bohren1983] eq. 4.88.
Following BHMIE, calculates D_n for complex argument using downward recursion, and RiccatiBessel functions psi and xi for real argument using upward recursion.
Initializes downward recursion for D_n using Lentz continued fraction algorithm [Lentz1976].
multilayer_sphere_lib.py
Author: Jerome Fung (fung@physics.harvard.edu)
Functions to calculate the scattering from a spherically symmetric particle with an arbitrary number of layers with different refractive indices.
Key reference for multilayer algorithm: Yang, “Improved recursive algorithm for light scattering by a multilayered sphere,” Applied Optics 42, 17101720, (1993).

scatcoeffs_multi
(marray, xarray, eps1=0.001, eps2=1e16)¶ Calculate scattered field expansion coefficients (in the Mie formalism) for a particle with an arbitrary number of spherically symmetric layers.
Parameters:  marray (array_like, complex128) – array of layer indices, innermost first
 xarray (array_like, real) – array of layer size parameters (k * outer radius), innermost first
 eps1 (float, optional) – underflow criterion for Lentz continued fraction for Dn1
 eps2 (float, optional) – convergence criterion for Lentz continued fraction for Dn1
Returns: scat_coeffs – Scattering coefficients
Return type: ndarray (complex)
Theories to compute scattering from objects.
All theories have a common interface defined by
holopy.scattering.theory.scatteringtheory.ScatteringTheory
.
Compute holograms using the discrete dipole approximation (DDA). Currently uses ADDA (http://code.google.com/p/adda/) to do DDA calculations. .. moduleauthor:: Thomas G. Dimiduk <tdimiduk@physics.harvard.edu>

class
DDA
(n_cpu=1, max_dpl_size=None, use_indicators=False, keep_raw_calculations=False, addacmd=[])¶ Bases:
holopy.scattering.theory.scatteringtheory.ScatteringTheory
Computes scattering using the the Discrete Dipole Approximation (DDA). It can (in principle) calculate scattering from any arbitrary scatterer. The DDA uses a numerical method that represents arbitrary scatterers as an array of point dipoles and then selfconsistently solves Maxwell’s equations to determine the scattered field. In practice, this model can be extremely computationally intensive, particularly if the size of the scatterer is larger than the wavelength of light. This model requires an external scattering code: adda .. attribute:: n_cpu
int (optional) – Number of threads to use for the DDA calculation
max_dpl_size
¶ float (optional) – Force a maximum dipole size. This is useful for forcing extra dipoles if necessary to resolve features in an object. This may make dda calculations take much longer.

use_indicators
¶ bool – If true, a scatterer’s indicators method will be used instead of its builtin adda definition

keep_raw_calculations
¶ bool – If true, do not delete the temporary file we run ADDA in, instead print its path so you can inspect its raw results
Notes
Does not handle near fields. This introduces ~5% error at 10 microns. This can in principle handle any scatterer, but in practice it will need excessive memory or computation time for particularly large scatterers.

required_spacing
(medium_wavelen, medium_index, n)¶

Calculates holograms of spheres using Fortran implementation of Mie theory. Uses superposition to calculate scattering from multiple spheres. Uses full radial dependence of spherical Hankel functions for scattered field.

class
Mie
(compute_escat_radial=True, full_radial_dependence=True, eps1=0.01, eps2=1e16)¶ Bases:
holopy.scattering.theory.scatteringtheory.ScatteringTheory
Compute scattering using the LorenzMie solution.
This theory calculates exact scattering for single spheres and approximate results for groups of spheres. It does not account for multiple scattering, hence the approximation in the case of multiple spheres. Neglecting multiple scattering is a good approximation if the particles are sufficiently separated.
This model can also calculate the exact scattered field from a spherically symmetric particle with an arbitrary number of layers with differing refractive indices, using Yang’s recursive algorithm ([Yang2003]).
By default, calculates radial component of scattered electric fields, which is nonradiative.
Currently, in calculating the LorenzMie scattering coefficients, the maximum size parameter x = ka is limited to 1000.
Defines Multisphere theory class, which calculates scattering for multiple spheres using the (exact) superposition method implemented in modified version of Daniel Mackowski’s SCSMFO1B.FOR. Uses full radial dependence of spherical Hankel functions for the scattered field.

class
Multisphere
(niter=200, eps=1e06, meth=1, qeps1=1e05, qeps2=1e08, compute_escat_radial=False, suppress_fortran_output=True)¶ Bases:
holopy.scattering.theory.scatteringtheory.ScatteringTheory
Exact scattering from a cluster of spheres.
Calculate the scattered field of a collection of spheres through a numerical method that accounts for multiple scattering and nearfield effects (see [Fung2011], [Mackowski1996]). This approach is much more accurate than Mie superposition, but it is also more computationally intensive. The Multisphere code can handle any number of spheres; see notes below for details.

niter
¶ integer (optional) – maximum number of iterations to use in solving the interaction equations

meth
¶ integer (optional) – method to use to solve interaction equations. Set to 0 for biconjugate gradient; 1 for orderofscattering

eps
¶ float (optional) – relative error tolerance in solution for interaction equations

qeps1
¶ float (optional) – error tolerance used to determine at what order the singlesphere spherical harmonic expansion should be truncated

qeps2
¶ float (optional) – error tolerance used to determine at what order the cluster spherical harmonic expansion should be truncated
Notes
According to Mackowski’s manual for SCSMFO1B.FOR [1] and later papers [2], the biconjugate gradient is generally the most efficient method for solving the interaction equations, especially for dense arrays of identical spheres. Orderofscattering may converge better for nonidentical spheres.
Multisphere does not check for overlaps becaue overlapping spheres can be useful for getting fits to converge. The results to be sensible for small overlaps even though mathemtically speaking they are not xstrictly valid.
Currently, Multisphere does not calculate the radial component of scattered electric fields. This is a good approximation for large kr, since the radial component falls off as 1/kr^2.
 scfodim.for contains three parameters, all integers:
 nod: Maximum number of spheres
 nod: Maximum order of individual sphere expansions. Will depend on
 size of largest sphere in cluster.
 notd: Maximum order of clustercentered expansion. Will depend on
 overall size of cluster.
Changing these values will require recompiling Fortran extensions.
The maximum size parameter of each individual sphere in a cluster is currently limited to 1000, indepdently of the above scfodim.for parameters.
References
[1] Daniel W. Mackowski, SCSMFO.FOR: Calculation of the Scattering Properties for a Cluster of Spheres, ftp://ftp.eng.auburn.edu/pub/dmckwski/scatcodes/scsmfo.ps [2] D.W. Mackowski, M.I. Mishchenko, A multiple sphere Tmatrix Fortran code for use on parallel computer clusters, Journal of Quantitative Spectroscopy and Radiative Transfer, In Press, Corrected Proof, Available online 11 March 2011, ISSN 00224073, DOI: 10.1016/j.jqsrt.2011.02.019. 

normalize_polarization
(illum_polarization)¶
Base class for scattering theories. Implements pythonbased calc_intensity and calc_holo, based on subclass’s calc_field .. moduleauthor:: Jerome Fung <jerome.fung@post.harvard.edu> .. moduleauthor:: Vinothan N. Manoharan <vnm@seas.harvard.edu> .. moduleauthor:: Thomas G. Dimiduk <tdimiduk@physics.harvard.edu>

class
ScatteringTheory
¶ Bases:
holopy.core.holopy_object.HoloPyObject
Defines common interface for all scattering theories. .. rubric:: Notes
A subclasses that do the work of computing scattering should do it by implementing _raw_fields and/or _raw_scat_matrs and (optionally) _raw_cross_sections. _raw_cross_sections is needed only for calc_cross_sections. Either of _raw_fields or _raw_scat_matrs will give you calc_holo, calc_field, and calc_intensity. Obviously calc_scat_matrix will only work if you implement _raw_cross_sections. So the simplest thing is to just implement _raw_scat_matrs. You only need to do _raw_fields there is a way to compute it more efficently and you care about that speed, or if it is easier and you don’t care about matrices.

stack_spherical
(a)¶

wavevec
(a)¶
Compute holograms using Mishchenko’s Tmatrix method for axisymmetric scatterers. Currently uses

class
Tmatrix
(delete=True)¶ Bases:
holopy.scattering.theory.scatteringtheory.ScatteringTheory
Computes scattering using the axisymmetric Tmatrix solution by Mishchenko with extended precision.
It can calculate scattering from axisymmetric scatterers such as cylinders and spheroids. Calculations for particles that are very large or have high aspect ratios may not converge.

delete
¶ bool (optional) – If true (default), delete the temporary directory where we store the input and output file for the fortran executable
Notes
Does not handle near fields. This introduces ~5% error at 10 microns.

Module contents¶
Scattering calculations
The scattering package provides objects and methods to define scatterer geometries, and theories to compute scattering from specified geometries. Scattering depends on holopy.core, and certain scattering theories may require external scattering codes.
The HoloPy scattering module is used to:
holopy.scattering.calculations module¶
Base class for scattering theories. Implements pythonbased calc_intensity and calc_holo, based on subclass’s calc_field

calc_cross_sections
(scatterer, medium_index=None, illum_wavelen=None, illum_polarization=None, theory='auto')¶ Calculate scattering, absorption, and extinction cross sections, and asymmetry parameter <cos heta>.
Parameters:  scatterer (
scatterer
object) – (possibly composite) scatterer for which to compute scattering  medium_index (float or complex) – Refractive index of the medium in which the scatter is imbedded
 illum_wavelen (float or ndarray(float)) – Wavelength of illumination light. If illum_wavelen is an array result will add a dimension and have all wavelengths
 theory (
theory
object (optional)) – Scattering theory object to use for the calculation. This is optional if there is a clear choice of theory for your scatterer. If there is not a clear choice, calc_intensity will error out and ask you to specify a theory
Returns: cross_sections – Dimensional scattering, absorption, and extinction cross sections, and <cos theta>
Return type: array (4)
 scatterer (

calc_field
(schema, scatterer, medium_index=None, illum_wavelen=None, illum_polarization=None, theory='auto')¶ Calculate hologram formed by interference between scattered fields and a reference wave
Parameters:  scatterer (
scatterer
object) – (possibly composite) scatterer for which to compute scattering  medium_index (float or complex) – Refractive index of the medium in which the scatter is imbedded
 illum_wavelen (float or ndarray(float)) – Wavelength of illumination light. If illum_wavelen is an array result will add a dimension and have all wavelengths
 theory (
theory
object (optional)) – Scattering theory object to use for the calculation. This is optional if there is a clear choice of theory for your scatterer. If there is not a clear choice, calc_intensity will error out and ask you to specify a theory
Returns: e_field – Calculated hologram from the given distribution of spheres
Return type: Vector
object scatterer (

calc_holo
(schema, scatterer, medium_index=None, illum_wavelen=None, illum_polarization=None, theory='auto', scaling=1.0)¶ Calculate hologram formed by interference between scattered fields and a reference wave
Parameters:  scatterer (
scatterer
object) – (possibly composite) scatterer for which to compute scattering  medium_index (float or complex) – Refractive index of the medium in which the scatter is imbedded
 illum_wavelen (float or ndarray(float)) – Wavelength of illumination light. If illum_wavelen is an array result will add a dimension and have all wavelengths
 theory (
theory
object (optional)) – Scattering theory object to use for the calculation. This is optional if there is a clear choice of theory for your scatterer. If there is not a clear choice, calc_intensity will error out and ask you to specify a theory  scaling (scaling value (alpha) for amplitude of reference wave) –
Returns: holo – Calculated hologram from the given distribution of spheres
Return type: Image
object scatterer (

calc_intensity
(schema, scatterer, medium_index=None, illum_wavelen=None, illum_polarization=None, theory='auto')¶ Calculate intensity at a location or set of locations
Parameters:  scatterer (
scatterer
object) – (possibly composite) scatterer for which to compute scattering  medium_index (float or complex) – Refractive index of the medium in which the scatter is imbedded
 illum_wavelen (float or ndarray(float)) – Wavelength of illumination light. If illum_wavelen is an array result will add a dimension and have all wavelengths
 theory (
theory
object (optional)) – Scattering theory object to use for the calculation. This is optional if there is a clear choice of theory for your scatterer. If there is not a clear choice, calc_intensity will error out and ask you to specify a theory
Returns: inten – scattered intensity
Return type: Image
 scatterer (

calc_scat_matrix
(schema, scatterer, medium_index=None, illum_wavelen=None, theory='auto')¶ Compute farfield scattering matrices for scatterer
Parameters:  scatterer (
holopy.scattering.scatterer
object) – (possibly composite) scatterer for which to compute scattering  medium_index (float or complex) – Refractive index of the medium in which the scatter is imbedded
 illum_wavelen (float or ndarray(float)) – Wavelength of illumination light. If illum_wavelen is an array result will add a dimension and have all wavelengths
 theory (
theory
object (optional)) – Scattering theory object to use for the calculation. This is optional if there is a clear choice of theory for your scatterer. If there is not a clear choice, calc_intensity will error out and ask you to specify a theory
Returns: scat_matr – Scattering matrices at specified positions
Return type: Marray
 scatterer (

check_schema
(schema, pol=True)¶

determine_theory
(scatterer)¶

finalize
(schema, result)¶

interpret_theory
(scatterer, theory='auto')¶

prep_schema
(schema, medium_index, illum_wavelen, illum_polarization)¶

scattered_field_to_hologram
(scat, ref, normals)¶ Calculate a hologram from an Efield
Parameters:
holopy.scattering.geometry module¶
Routines for common calculations and transformations of groups of spheres.
This code is in need of significant refactoring and simplification, refactoring which may break code that depends on it.

angles
(cluster, degrees=True)¶ calculate the angles between one particle and every pair of other particles
Parameters:  cluster (
holopy.scattering.scatterer.Scatterer
) – A sphere cluster to determine the interparticle distances of.  degrees (bool) – Whether to return angles in degrees (True) or in radians (False).
Notes
Angle abc is the acute angle formed by edges conecting points ab and bc. If a, b, and c are locations of particles (vertices), the returned 3D array has nonzero values for angles abc, zeros for angles aba, and NAN’s for “angles” aab.
 cluster (

distances
(cluster, gaponly=False)¶ calculate the distances between each sphere in a cluster and each of the others
Parameters:  cluster (
holopy.scattering.scatterer.Scatterer
) – A sphere cluster to determine the interparticle distances of.  gaponly (bool) – Whether to calculate the distances between particle centers or between particle surfaces (gap distances).
Notes
The returned array of distances includes redundant information. The identical distances between sphere 1 and sphere 2 and between sphere 2 and sphere 1 are both in the returned array. Calculating and returning the full array makes it easy for the user to access all the interparticle distances starting from any sphere of interest.
 cluster (

make_cubecluster
(index, radius, gap, xcom=0, ycom=0, zcom=0)¶ Returns a sphere cluster of eight particles forming a cube centered on a given center of mass.
Parameters:  index – Index of refraction of particles.
 radius – Radius if particles.
 gap – Space to add between the particles.
 xcom – Center of mass xcoordinate
 ycom – Center of mass ycoordinate
 zcom – Center of mass zcoordinate

make_octacluster
(index, radius, gap, xcom=0, ycom=0, zcom=0)¶ Returns a sphere cluster of six particles forming an octahedron centered on a given center of mass.
Parameters:  index – Index of refraction of particles.
 radius – Radius if particles.
 gap – Space to add between the particles.
 xcom – Center of mass xcoordinate
 ycom – Center of mass ycoordinate
 zcom – Center of mass zcoordinate

make_polytetracluster
(index, radius, gap, xcom=0, ycom=0, zcom=0)¶ Returns a sphere cluster of six particles forming a polytetrahedron centered on a given center of mass of the middle tetrahedron.
Parameters:  index – Index of refraction of particles.
 radius – Radius if particles.
 gap – Space to add between the particles.
 xcom – Center of mass of the middle tetrahedron xcoordinate
 ycom – Center of mass of the middle tetrahedron xcoordinate
 zcom – Center of mass of the middle tetrahedron xcoordinate

make_sqcluster
(index, radius, gap, xcom=0, ycom=0, zcom=0)¶ Returns a sphere cluster of four particles forming a square centered on a given center of mass.
Parameters:  index – Index of refraction of particles.
 radius – Radius if particles.
 gap – Space to add between the particles.
 xcom – Center of mass xcoordinate
 ycom – Center of mass ycoordinate
 zcom – Center of mass zcoordinate

make_tetracluster
(index, radius, gap, xcom=0, ycom=0, zcom=0)¶ Returns a sphere cluster of four particles forming a tetrahedron centered on a given center of mass.
Parameters:  index – Index of refraction of particles.
 radius – Radius if particles.
 gap – Space to add between the particles.
 xcom – Center of mass xcoordinate
 ycom – Center of mass ycoordinate
 zcom – Center of mass zcoordinate

make_tribipyrcluster
(index, radius, gap, xcom=0, ycom=0, zcom=0)¶ Returns a sphere cluster of five particles forming a triagonal bipyramid centered on a given center of mass.
Parameters:  index – Index of refraction of particles.
 radius – Radius if particles.
 gap – Space to add between the particles.
 xcom – Center of mass xcoordinate
 ycom – Center of mass ycoordinate
 zcom – Center of mass zcoordinate

make_tricluster
(index, radius, gap, xcom=0, ycom=0, zcom=0)¶ Returns a sphere cluster of three particles forming an equilateral triangle centered on a given center of mass.
Parameters:  index – Index of refraction of particles.
 radius – Radius if particles.
 gap – Space to add between the particles.
 xcom – Center of mass xcoordinate
 ycom – Center of mass ycoordinate
 zcom – Center of mass zcoordinate
holopy.vis package¶
Module contents¶
Visualize HoloPy objects
Uses Matplotlib and Mayavi to visualize holopy objects.
Image
,Volume
, orSpheres
object => plot or rendering
This module does not import plotting packages until they are actually needed so that holopy does not have a hard dependency on mayavi or matplotlib. Because of this you may see a small lag on your first plot.
holopy.vis.show module¶
A general show method that can display most holopy and scatterpy objects in a sensible way.

exception
VisualizationNotImplemented
(o)¶ Bases:
Exception

show
(o, color=(0.5, 0.5, 0.5))¶ Visualize a scatterer, hologram, or reconstruction
Parameters: o ( Image
,Volume
orSpheres
) – Object to visualizeNotes
Loads plotting library the first time it is required (so that we don’t have to import all of matplotlib or mayavi just to load holopy)

test_disp
()¶
holopy.vis.vis2d module¶
New custom display functions for holograms and reconstructions.

class
plotter
(im, plane_axes, slice_axis, starting_index, t)¶ Bases:
object

click
(event)¶

draw
()¶

location
(x, y)¶

pixel
(x, y)¶


show2d
(im, plane_axes=None, slice_axis=None, starting_index=0, t=0, phase=False)¶ Display a hologram or reconstruction
Allows scrolling through multidimensional holograms or reconstructions. Defaults to showing magnitude of complex images
Parameters:

show_scatterer_slices
(scatterer, spacing)¶ Show slices of a scatterer voxelation
 scatterer : .Scatterer
 scatterer to visualize
 spacing : float or (float, float, float)
 voxel spacing for the visualization
References and credits¶
The following references describe applications of HoloPy and technical advances. If you use HoloPy, we ask that you cite the articles that are relevant to your application.
[Dimiduk2016]  Dimiduk, Thomas G., and Vinothan N. Manoharan. “Bayesian Approach to Analyzing Holograms of Colloidal Particles.” Optics Express 24, no. 21 (October 17, 2016): 24045–60. doi:10.1364/OE.24.024045. 
[Wang2016]  Wang, Anna, Rees F. Garmann, and Vinothan N. Manoharan. “Tracking E. Coli Runs and Tumbles with Scattering Solutions and Digital Holographic Microscopy.” Optics Express 24, no. 21 (October 17, 2016): 23719–25. doi:10.1364/OE.24.023719. 
[Dimiduk2014]  Dimiduk, Thomas G., Rebecca W. Perry, Jerome Fung, and Vinothan N. Manoharan. “RandomSubset Fitting of Digital Holograms for Fast ThreeDimensional Particle Tracking.” Applied Optics 53, no. 27 (September 20, 2014): G177–83. doi:10.1364/AO.53.00G177. 
[Wang2014]  Wang, Anna, Thomas G. Dimiduk, Jerome Fung, Sepideh Razavi, Ilona Kretzschmar, Kundan Chaudhary, and Vinothan N. Manoharan. “Using the Discrete Dipole Approximation and Holographic Microscopy to Measure Rotational Dynamics of NonSpherical Colloidal Particles.” Journal of Quantitative Spectroscopy and Radiative Transfer 146 (October 2014): 499–509. doi:10.1016/j.jqsrt.2013.12.019. 
[Fung2013]  Fung, Jerome, and Vinothan N. Manoharan. “Holographic Measurements of Anisotropic ThreeDimensional Diffusion of Colloidal Clusters.” Physical Review E 88, no. 2 (August 30, 2013): 020302. doi:10.1103/PhysRevE.88.020302. 
[Fung2012]  Fung, Jerome, Rebecca W. Perry, Thomas G. Dimiduk, and Vinothan N. Manoharan. “Imaging Multiple Colloidal Particles by Fitting Electromagnetic Scattering Solutions to Digital Holograms.” Journal of Quantitative Spectroscopy and Radiative Transfer 113, no. 18 (December 2012): 2482–89. doi:10.1016/j.jqsrt.2012.06.007. 
[Kaz2012]  Kaz, David M., Ryan McGorty, Madhav Mani, Michael P. Brenner, and Vinothan N. Manoharan. “Physical Ageing of the Contact Line on Colloidal Particles at Liquid Interfaces.” Nature Materials 11, no. 2 (February 2012): 138–42. doi:10.1038/nmat3190. 
[Perry2012]  Perry, Rebecca W., Guangnan Meng, Thomas G. Dimiduk, Jerome Fung, and Vinothan N. Manoharan. “RealSpace Studies of the Structure and Dynamics of SelfAssembled Colloidal Clusters.” Faraday Discussions 159, no. 1 (June 7, 2012): 211–34. doi:10.1039/C2FD20061A. 
[Fung2011]  Fung, Jerome, K. Eric Martin, Rebecca W. Perry, David M. Kaz, Ryan McGorty, and Vinothan N. Manoharan. “Measuring Translational, Rotational, and Vibrational Dynamics in Colloids with Digital Holographic Microscopy.” Optics Express 19, no. 9 (April 25, 2011): 8051–65. doi:10.1364/OE.19.008051. 
Ovryn and Izen and Lee and coworkers were the first to develop methods to fit scattering models to digital holograms:
[Ovryn2000]  Ovryn, Ben, and Steven H. Izen. “Imaging of Transparent Spheres through a Planar Interface Using a HighNumericalAperture Optical Microscope.” Journal of the Optical Society of America A 17, no. 7 (July 1, 2000): 1202–13. doi:10.1364/JOSAA.17.001202. 
[Lee2007]  Lee, SangHyuk, Yohai Roichman, GiRa Yi, ShinHyun Kim, SeungMan Yang, Alfons van Blaaderen, Peter van Oostrum, and David G. Grier. “Characterizing and Tracking Single Colloidal Particles with Video Holographic Microscopy.” Optics Express 15, no. 26 (December 24, 2007): 18275–82. doi:10.1364/OE.15.018275. 
The following papers describe different methods for calculating scattering and various algorithms that HoloPy uses in its calculations:
[Yurkin2011]  Yurkin, Maxim A., and Alfons G. Hoekstra. “The DiscreteDipoleApproximation Code ADDA: Capabilities and Known Limitations.” Journal of Quantitative Spectroscopy and Radiative Transfer 112, no. 13 (September 2011): 2234–47. doi:10.1016/j.jqsrt.2011.01.031. 
[Mackowski1990]  Mackowski, D. W., Altenkirch, R. A., and Menguc, M. P. “Internal absorption cross sections in a stratified sphere.” Applied Optics 29, no. 10 (1990). do:10.1364/AO.29.001551. 
[Mackowski1996]  Mackowski, Daniel W., and Michael I. Mishchenko. “Calculation of the T Matrix and the Scattering Matrix for Ensembles of Spheres.” Journal of the Optical Society of America A 13, no. 11 (November 1, 1996): 2266. doi:10.1364/JOSAA.13.002266. 
[Wiscombe1996]  Wiscombe, J. “Mie Scattering Calculations: Advances in Technique and Fast, VectorSpeed Computer Codes,” 1979. doi:10.5065/D6ZP4414. 
[Wiscombe1980]  Wiscombe, W. J. “Improved Mie scattering algorithms.” Applied Optics 19, no. 9 (May 1, 1980): 1505. doi:10.1364/AO.19.001505. 
[Yang2003]  Yang, Wen. “Improved Recursive Algorithm for Light Scattering by a Multilayered Sphere.” Applied Optics 42, no. 9 (March 20, 2003): 1710–20. doi:10.1364/AO.42.001710. 
[Lentz1976]  Lentz, William J. “Generating Bessel Functions in Mie Scattering Calculations Using Continued Fractions.” Applied Optics 15, no. 3 (March 1, 1976): 668–71. doi:10.1364/AO.15.000668. 
For scattering calculations and formalism, we draw heavily on the treatise of Bohren & Huffman. We generally follow their conventions except where noted.
[Bohren1983]  C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles, Wiley (1983). 
The package includes code from several sources. We thank Daniel Mackowski for allowing us to include his TMatrix code, which computes scattering from clusters of spheres: SCSMFO1B.
We also make use of a modified version of the Python version of mpfit, originally developed by Craig Markwardt. The modified version we use is drawn from the stsci_python package.
We thank A. Ross Barnett for permitting us to use his routine SBESJY.FOR, which computes spherical Bessel functions.
HoloPy is based upon work supported by the National Science Foundation under grant numbers CBET0747625, DMR0820484, DMR1306410, and DMR1420570.