10 Image Utilities

This notebook forms part of a series on computational optical radiometry. The notebooks can be downloaded from Github. These notebooks are constantly revised and updated, please revisit from time to time.

The date of this document and module versions used in this document are given at the end of the file.
Feedback is appreciated: neliswillers at gmail dot com.

Overview

The pyradi module is not intended to provide a full-function image processing functionality. If you require a comprehensive image processing capability consider scikit-image, PIL (not actively developed anymore), or Pillow (a PIL fork, actively developed). In fact, the functions described may use scikit-image internally.

Pyradi provides several productivity-enhancing image utilities and a few unique functions not available elsewhere. These functions are supplied to support the broader pyradi objectives.

Loading multi-frame two-dimensional arrays from a raw binary data files

The ryfiles.readRawFrames function reads any number of two-dimensional array frames from a binary file with known data type.

The file must consist of multiple frames, all with the same number of rows and columns. Frames of different data types can be read, according to the user specification. The user can specify which frames must be loaded (if not the whole file).

The function signature is:

`readRawFrames(fname, rows, cols, vartype, loadFrames=[])`

• fname (string) filename.
• rows (int) number of rows in each frame.
• cols (int) number of columns in each frame.
• vartype (np.dtype) numpy data type of data to be read: int8, int16, int32, int64, uint8, uint16, uint32, uint64, float16, float32, float64.
• loadFrames ([int]) list of frames to load, zero-based , empty list (default) loads all frames (optional).

The function returns:

• frames (int) number of frames in the returned data set, 0 if error occurred.
• rawShaped (np.ndarray) vartype numpy array of dimensions (frames,rows,cols), None if error occurred.

Saving a two-dimensional array as an image

The ryfiles.rawFrameToImageFile function saves the data in a two-dimensional array as an image.

The file type must be given, e.g. png or jpg. The image need not be scaled beforehand, it is done prior to writing out the image. Could be one of BMP, JPG, JPEG, PNG, PPM, TIFF, XBM, XPM) but the file types available depends on the imsave plugin in use.

The function signature is:

`rawFrameToImageFile(image, filename)`

• image (np.ndarray) two-dimensional array representing an image.
• filename (string) name of file to be written to, with extension.

Reading a binary file and saving as image

In the following example a data file is downloaded from the internet. The data file contains 8 frames, each of 100 rows and 100 colummns. The data type is a four-byte unsigned long. In the first case all the frames in the imge is loaded, in the second case only selected frames are loaded. The loadFrames list must provide a zero-based list of frames to be loaded.

Image sequence ordering

Image sequences are stored in three-dimensional arrays, in rows, columns and frames. Not all libraries share the same sequencing, some store frames along axis=0 and others store frames along axis=2.

The ryutils.framesFirst function reorders an image sequence with frames along axis=2 to an image sequence with frames along axis=0. The ryutils.framesFirst function uses np.transpose(imageSequence, (2,0,1)). The function signature is framesFirst(imageSequence) where imageSequence is the image sequence. A view or a copy of the reordered sequence is returned.

The ryutils.framesLast function reorders an image sequence with frames along axis=0 to an image sequence with frames along axis=2. The ryutils.framesLast function uses np.transpose(imageSequence, (1,2,0)). The function signature is framesLast(imageSequence) where imageSequence is the image sequence. A view or a copy of the reordered sequence is returned.

Image histogram equalisation

pyradi is not meant to be an image processing package, but a histogram equalisation capability is included as a support for plotting tasks. The ryplot.ProcessImage class contains the function ryplot.ProcessImage.compressEqualizeImage to do histogram equalisation. This function has a unique feature in that it returns a set of colourbar values and labels, which retains the original data set levels, so that the true original input values can be read off on the colourbar.

Compress an image (and then inversely expand the color bar values), prior to histogram equalisation to ensure that the two keep in step, we store the compression function names as pairs, and invoke the compression function as follows: linear, log. sqrt. Note that the image is histogram equalised in all cases.

There are three compression functions available: Linear, Natural Log, and Square Root. Different datasets equalise differently with the different compressors. Experiment and find the most appropriate for your data set.

compressSet = [
           [lambda x : x , lambda x : x, 'Linear'],
           [np.log,  np.exp, 'Natural Log'],
           [np.sqrt,  np.square, 'Square Root']
           ]

The function signature is compressEqualizeImage(image, selectCompressSet=2, numCbarlevels=20, cbarformat='.3f')

• image (np.ndarray) the image to be processed.
• selectCompressSet (int) compression selection [0,1,2] (optional).
• numCbarlevels (int) number of labels in the colourbar (optional).
• cbarformat (string) colourbar label format, e.g., '10.3f', '.5e' (optional).

The ryplot.ProcessImage.compressEqualizeImage function returns a tuple, where the first element is the compressed image and the second element is a zip of colourbar levels and associated colourbar labels.

In the following example a data set is histogram-equalised and compared to the original input dataset. On the second figure, note the non-linear scale on the colourbar; ryplot.ProcessImage.compressEqualizeImage keeps track of the original levels and even though the image is compressed, the colourbar still shows the original input values.

Three dimensional noise analysis

The pyradi.ry3dnoise module provides a set of functions to aid in the calculation of 3D noise parameters from noise images. The functions are based on the work done by John D'Agostino and Curtis Webb. For details see "3-D Analysis Framwork and Measurement Methodology for Imaging System Nioise" p110-121 in "Infrared Imaging Systems: Design, Analysis, Modelling, and Testing II", Holst, G. C., ed., Volume 1488, SPIE (1991), DOI 10.1117/12.45794.

The three-dimensional noise analysis requires a sequence of images. The functions all receive an numpy.ndarray[frames][rows][cols] as first argument, where the array is a three-dimensional sequence of images. The functions in this module are:

• getNT(imgSeq) Average for all pixels as a function of time/frames. Represents noise which consists of fluctuations in the temporal direction affecting the mean of each frame. Returns the noise (double), frame-to-frame intensity variation.
• getNVH(imgSeq) Average over all frames, for each pixel. Represents non-uniformity spatial noise that does not change from frame-to-frame. Returns the noise (double), fixed spatial noise.
• getNTV(imgSeq) Average for each row and frame over all columns. Represents variations in row averages that change from frame-to-frame. Returnsthe noise (double), row temporal noise.
• getNTH(imgSeq) Average for each column and frame over all rows. Represents variations in column averages that change from frame-to-frame. Returnsthe noise (double), column temporal noise.
• def getNV(imgSeq) Average for each column over all frames and rows. Represents variations in row averages that are fixed in time. Returns the noise (double), fixed row noise.
• getNH(imgSeq) Average for each row over all frames and cols. Represents variations in column averages that are fixed in time. Returns the noise (double), fixed column noise.
• getNTVH(imgSeq)Noise for each row, frame and column. Represents random noise in the detector and electronics. Returns the noise (double), temporal pixel noise.
• getTotal(imgSeq) Total system noise. Returns the noise (double), total system noise.

The following example reads a file with an image sequence, and proceeds to calculate the three-dimensional noise parameters.

Image warping: cartesian to polar

ryplot.ProcessImage.reprojectImageIntoPolar function wraps or reprojects an image sequence from cartesian to polar coordinates, relative to some origin.

The origin of the new coordinate system defaults to the center of the image, unless the user supplies a new origin. The data format can be data.shape = (rows, cols, frames) or data.shape = (frames, rows, cols), the format of which is indicated by the framesFirst parameter.

The function signature is: reprojectImageIntoPolar(self, data, origin=None, framesFirst=True)

• data (np.array) 3-D array to which transformation must be applied.
• origin ( (x-orig, y-orig) ) data-coordinates of where origin should be placed
• framesFirst (bool) True if data.shape is (frames, rows, cols), False if data.shape is (rows, cols, frames)

the function returns

• output (float np.array) transformed images/array data in the same sequence as input sequence.
• r_i (np.array[N,]) radial values for returned image.
• theta_i (np.array[M,]) angular values for returned image.

This code in this function is based on code by Joe Kington

The following example converts Siemens Star and a bullseye images from cartesian to polar coordinates. In both cases are the centre of the object, near or on the centre of the image. You can experiment with this code by moving the origin to other locations.

MTF from a distorted wavefront

The optical transfer function is the Fourier transform of the point spread function, and the point spread function is the square absolute of the inverse Fourier transformed pupil function. The optical transfer function can also be calculated directly from the pupil function. From the convolution theorem it can be seen that the optical transfer function is the autocorrelation of the pupil function https://en.wikipedia.org/wiki/Optical_transfer_function.

The pupil function comprises a masking shape (the binary shape of the pupil) and a transmittance and spatial phase delay inside the mask. A perfect aperture has unity transmittance and zero phase delay in the mask. Some pupils have irregular pupil functions/shapes and hence the diffraction MTF has to be calculated numerically using images (masks) of the pupil function.

From the OSA Handbook of Optics, Vol II, p 32.4:
For an incoherent optical system, the OTF is proportional to the two-dimensional autocorrelation of the exit pupil. This calculation can account for any phase factors across the pupil, such as those arising from aberrations or defocus. A change of variables is required for the identification of an autocorrelation (a function of position in the pupil) as a transfer function (a function of image-plane spatial frequency). The change of variables is

where $x$ is the autocorrelation shift distance in the pupil, $\lambda$ is the wavelength, and $d_i$ is the distance from the exit pupil to the image. A system with an exit pupil of full width $D$ has an image-space cutoff frequency (at infinite conjugates) of

In this analysis we assume that

1. the sensor is operating at infinite conjugates.
2. the mask falls in the entrance pupil shape.

The MTF is calculated as follows:

1. Read in the pupil function mask and create an image of the mask.
2. Calculate the two-dimensional autocorrelation function of the binary image (using the SciPy two-dimensional correlation function signal.correlate2d).
3. Scale the magnitude and $(x,y)$ dimensions according to the dimensions of the physical pupil.

Pupil shape

For this analysis the pupil has a circular shape, the size of the mirror. If the mirror has a central obscuration the pupil must also model the obscuration.

Pupil transmittance

For the purpose of this analysis the transmittance is assumed spatially constant at unity.

Pupil phase function

The mirror deviation from the ideal shape introduces a phase error in the optical wave front. The phase error, expressed in radians is given by

where $\Delta$ is the figure error in m, the factor 2 converts the mirror figure error into wavefront displacement, $\lambda$ the wavelength in m converts the wavefront displacement into cycles, and the $2\pi$ converts the displacement in cycles into displacement in radians.

The following graphs plot the wavefront error in m and in radian for the spectral bands defined in specdef, plus the reference case with no wavefront error. In the following example the maximum wavefront error is of the order of 0.25 $\mu$m. The maximum phase error is in the visible band and is approximately 3 rad, which is about half a cycle, which corresponds to the maximum wavefront error.

The following graph shows the two dimensional MTF for the entries in specdef.

The following graph shows the two dimensional ratios of MTF: true MTF divided by the diffraction limit MTF. In other words, the graphs show the degradation in MTF resulting from the wavefront error. In this example, the MTF is not significantly affected in the IR band, but quite poor in the visual band.

The next graph shows the radial MTF for a large number of rotational angles (amongst which are tangential and sagital). The IR MTF does not show significant degradation and variation around rotation, whereas the visible MTF does show significant rotational variation.

The final set of graphs show the mean MTF across all rotational angles (left-side graphs) and the ratio of mean MTF (across all rotational angles) to diffraction MTF (right-side graphs). The graphs on the right side shows that in the visible band the mirror's figure error degrades the MTF by about 40% at most spatial frequencies, whereas the IR MTF is not affected significantly.

Threshold image by interactive widgets

https://ipywidgets.readthedocs.io/en/latest/user_install.html

Introduction

IPython 2.0 introduced interactive widgets that allows the creation of GUI widgets in the notebook. The widgets are relatively simple to implement in notebook code. The interactive widgets API support different layers of detail, providing increasing control at lower levels.

This notebook explores a very simple scenario where a slider is used to specify a threshold used to segment an image.

An excellent introduction to the IPython display system and interactive widgets was given at a recent PyData conference by Granger and Frederic. The introduction provides IPython notebooks as well as videos of the presentations at the conference.

Videos are here:
https://www.youtube.com/watch?v=VaV10VNZCLA
https://www.youtube.com/watch?v=vE_CJTen15M
https://www.youtube.com/watch?v=o7Tb7YhJZR0

IPython notebooks are also available. Clone or download the complete repository at
https://github.com/ipython/ipython-in-depth
and then open the file ipython-in-depth/examples/Interactive Widgets/Index.ipynb

Read at least the first two tutorials Using Interact and Widget Basics. There is much more to read, but these two should get you going.

Segmenting an image

The following example loads an image file, sets up a slider and then display the segmented image and its statistical properties.

Play with the two images to observe the effect of the different gray level distributions.

The following function does the image segmentation and will be called from widgets.interactive. The function is given the image and the threshold value to be used in the segmentation. The image is segmented and displayed inside this function - required because the image must be displayed every time the segmentation threshold is changed. The function returns the segmented image and threshold value for the image (these returns can be used subsequently elsewhere). The function displays the histogram of the segmented image, excluding the first bin, which contains the discarded elements, as well as the cumulative probability of the original input image. All this processing takes a while, so the software is sluggish for large images.

The widgets.interactive function is given three types of parameters:

• The first parameter is the name of the function to be called when the widget's states changes.
• Any input data to be passed to the function (any number of parameters). These inputs must be flagged as inputs by using widgets.fixed() function. If so flagged the parameters will not be interpreted as widgets, but as input data.
• The definition of the widgets to be drawn and interrogated. All parameters not flagged as widgets.fixed() are interpreted as widget definitions, either in shorthand form or by explicit definition (see the Granger/Frederic notebooks).

In this case the function is threshold_image, and only one data item is passed img, and there a float slider widget ranging from 0 to 255, with step increment of 1.0.

The interactive function is used and assigned to a variable. This is necessary because later we want to use the variable to read the return values from the segmentation function. The display function must be called on the widget variable to display the widget. Once the widget is displayed, the user can change the slide value, which in turn calls the threshold function, which does the segmentation and displays the image.

The segmentation function returns a segmented image and the value of the segmentation threshold. These values can be accessed after the slider has been moved, as the result attribute of the widget. The inputs to the segmentation function is available in the kwargs attribute.

Python and module versions, and dates