Utility Functions

Overview

This module provides various utility functions for radiometry calculations. Functions are provided for a maximally flat spectral filter, a simple photon detector spectral response, effective value calculation, conversion of spectral domain variables between [um], [cm^-1] and [Hz], conversion of spectral density quantities between [um], [cm^-1] and [Hz] and spectral convolution.

See the __main__ function for examples of use.

This package was partly developed to provide additional material in support of students and readers of the book Electro-Optical System Analysis and Design: A Radiometry Perspective, Cornelius J. Willers, ISBN 9780819495693, SPIE Monograph Volume PM236, SPIE Press, 2013. http://spie.org/x648.html?product_id=2021423&origin_id=x646

Module classes

class VersionInformation(line='')[source]

Bases: object

Gathers and presents version information for Python and installed modules.

Produces output suitable for IPython/Jupyter notebooks in plaintext, HTML, JSON, and LaTeX formats.

Example

print(ryutils.VersionInformation(‘matplotlib,numpy’))

version_information(line='')[source]

Show information about versions of modules.

Collects Python runtime info and, optionally, the versions of any additional modules supplied as a comma-separated string.

Parameters:

line (str) – Comma-separated list of module names (optional).

Returns:

self, to allow chaining.

Return type:

VersionInformation

class Spectral(ID, value, wl=None, wn=None, desc=None)[source]

Bases: object

Generic spectral can be used for any spectral vector

vecalign(other)[source]

returns two spectral values properly interpolated and aligned to same base

it is not intended that the function will be called directly by the user

Args:
other (Spectral): the other Spectral to be used in addition
Returns:
wl, wn, s, o
Raises:
No exception is raised.
plot(filename=None, heading=None, ytitle='')[source]

Do a simple plot of spectral variable(s)

Parameters:

  • filename (|) – filename for png graphic

  • heading (|) – graph heading

  • ytitle (|) – graph y-axis title

Returns:

Nothing, writes png file to disk

Raises:

| No exception is raised. –

class Atmo(ID, distance=None, wl=None, wn=None, tran=None, atco=None, prad=None, desc=None)[source]

Bases: object

Atmospheric spectral such as transittance or attenuation coefficient

tauR(distance)[source]

Calculates the transmittance at distance

Distance is in m

Args:
distance (scalar or np.array (M,)): distance in m if transmittance, or None if att coeff
Returns:
transmittance (np.array (N,M) ): transmittance along N at distance along M
Raises:
No exception is raised.
pathR(distance)[source]

Calculates the path radiance at distance

Distance is in m

Args:
distance (scalar or np.array (M,)): distance in m if transmittance, or None if att coeff
Returns:
transmittance (np.array (N,M) ): transmittance along N at distance along M
Raises:
No exception is raised.
class Sensor(ID, fno, detarea, inttime, tauOpt=1, quantEff=1, pfrac=1, desc='')[source]

Bases: object

Sensor characteristics

tauOpt()[source]

Returns scaler or np.array for optics transmittance

Parameters:

None (|)

Returns:

str

Raises:

| No exception is raised. –

QE()[source]

Returns scaler or np.array for detector quantEff

Parameters:

None (|)

Returns:

str

Raises:

| No exception is raised. –

class Target(ID, tmprt, emis, refl=1, cosTarg=1, taumed=1, scale=1, desc='')[source]

Bases: Spectral

Target / Source characteristics

emis()[source]

Returns scaler or np.array for emissivity

Parameters:

None (|)

Returns:

str

Raises:

| No exception is raised. –

refl()[source]

Returns scaler or np.array for reflectance

Parameters:

None (|)

Returns:

str

Raises:

| No exception is raised. –

taumed()[source]

Returns scaler or np.array for atmospheric transmittance to illuminating source

Parameters:

None (|)

Returns:

str

Raises:

| No exception is raised. –

radiance(units='el')[source]

Returns radiance spectral for target

The type of spectral is one of the following:

type=’el [W/(m$^2$.$mu$m)] type=’ql’ [q/(s.m$^2$.$mu$m)] type=’en’ [W/(m$^2$.cm$^{-1}$)] type=’qn’ [q/(s.m$^2$.cm$^{-1}$)]

Args:
None
Returns:
str
Raises:
No exception is raised.

Module functions

sfilter(spectral, center, width, exponent=6, taupass=1.0, taustop=0.0, filtertype='bandpass')[source]

Calculate a symmetrical filter response of shape exp(-x^n)

Given a number of parameters, calculates maximally flat, symmetrical transmittance. The function parameters controls the width, pass-band and stop-band transmittance and sharpness of cutoff. This function is not meant to replace the use of properly measured filter responses, but rather serves as a starting point if no other information is available. This function does not calculate ripple in the pass-band or cut-off band.

Filter types supported include band pass, high (long) pass and low (short) pass filters. High pass filters have maximal transmittance for all spectral values higher than the central value. Low pass filters have maximal transmittance for all spectral values lower than the central value.

Parameters:

  • spectral (|) – spectral vector in [um] or [cm-1].

  • center (|) – central value for filter passband

  • width (|) – proportional to width of filter passband

  • exponent (|) – even integer, define the sharpness of cutoff.

  • gaussian (| If exponent=2 then)

  • square (| If exponent=infinity then)

  • taupass (|) – the transmittance in the pass band (assumed constant)

  • taustop (|) – peak transmittance in the stop band (assumed constant)

  • filtertype (|) – filter type, one of ‘bandpass’, ‘lowpass’ or ‘highpass’

Returns:

transmittances at “spectral” intervals.

Return type:

transmittance (np.array[N,] or [N,1])

Raises:

  • | No exception is raised. –

  • | If an invalid filter type is specified, return None. –

  • | If negative spectral is specified, return None. –

responsivity(wavelength, lwavepeak, cuton=1, cutoff=20, scaling=1.0)[source]

Calculate a photon detector wavelength spectral responsivity

Given a number of parameters, calculates a shape that is somewhat similar to a photon detector spectral response, on wavelength scale. The function parameters controls the cutoff wavelength and shape of the response. This function is not meant to replace the use of properly measured spectral responses, but rather serves as a starting point if no other information is available.

Parameters:

  • wavelength (|) – vector in [um].

  • lwavepeak (|) – approximate wavelength at peak response

  • cutoff (|) – cutoff strength beyond peak, 5 < cutoff < 50

  • cuton (|) – cuton sharpness below peak, 0.5 < cuton < 5

  • scaling (|) – scaling factor

Returns:

responsivity at wavelength intervals.

Return type:

responsivity (np.array[N,] or [N,1])

Raises:

| No exception is raised. –

effectiveValue(spectraldomain, spectralToProcess, spectralBaseline)[source]

Normalise a spectral quantity to a scalar, using a weighted mapping by another spectral quantity.

Effectivevalue = integral(spectralToProcess * spectralBaseline) / integral( spectralBaseline)

The data in spectralToProcess and spectralBaseline must both be sampled at the same domain values as specified in spectraldomain.

The integral is calculated with numpy/scipy trapezoid trapezoidal integration function.

Parameters:

  • inspectraldomain (|) – spectral domain in wavelength, frequency or wavenumber.

  • spectralToProcess (|) – spectral quantity to be normalised

  • spectralBaseline (|) – spectral serving as baseline for normalisation

Returns:

effective value | Returns None if there is a problem

Return type:

(float)

Raises:

| No exception is raised. –

convertSpectralDomain(inspectraldomain, type='')[source]

Convert spectral domains, i.e. between wavelength [um], wavenummber [cm^-1] and frequency [Hz]

In string variable type, the ‘from’ domain and ‘to’ domains are indicated each with a single letter: ‘f’ for temporal frequency, ‘l’ for wavelength and ‘n’ for wavenumber The ‘from’ domain is the first letter and the ‘to’ domain the second letter.

Note that the ‘to’ domain vector is a direct conversion of the ‘from’ domain to the ‘to’ domain (not interpolated or otherwise sampled.

Parameters:

  • inspectraldomain (|) – spectral domain in wavelength, frequency or wavenumber.

  • [um] (| wavelength vector in)

  • [Hz] (| frequency vector in)

  • [cm^-1] (| wavenumber vector in)

  • type (|) – specify from and to domains:

  • frequency (| 'nf' convert from wavenumber to per)

  • wavenumber (| 'fn' convert from frequency to per)

  • wavelength (| 'nl' convert from wavenumber to per)

  • wavenumber

  • wavelength

  • frequency

Returns:

outspectraldomain | Returns zero length array if type is illegal, i.e. not one of the expected values

Return type:

[N,1]

Raises:

| No exception is raised. –

convertSpectralDensity(inspectraldomain, inspectralquantity, type='', outspecdomainFix=False, outspecdomain=None)[source]

Convert spectral density quantities, i.e. between W/(m^2.um), W/(m^2.cm^-1) and W/(m^2.Hz).

In string variable type, the ‘from’ domain and ‘to’ domains are indicated each with a single letter: ‘f’ for temporal frequency, ‘w’ for wavelength and ‘’n’ for wavenumber The ‘from’ domain is the first letter and the ‘to’ domain the second letter.

The return values from this function are always positive, i.e. not mathematically correct, but positive in the sense of radiance density.

The spectral density quantity input is given as a two vectors: the domain value vector and the density quantity vector. The output of the function is also two vectors, i.e. the ‘to’ domain value vector and the ‘to’ spectral density. Note that the ‘to’ domain vector is a direct conversion of the ‘from’ domain to the ‘to’ domain (not interpolated or otherwise sampled).

Parameters:

  • inspectraldomain (|) – spectral domain in wavelength, frequency or wavenumber.

  • inspectralquantity (|) – spectral density in same domain as domain vector above.

  • [um] (| wavelength vector in)

  • [Hz] (| frequency vector in)

  • [cm^-1] (| wavenumber vector in)

  • type (|) – specify from and to domains:

  • density (| 'nf' convert from per wavenumber interval density to per frequency interval)

  • density

  • density

  • density

  • density

  • density

  • outspecdomainFix (|) – if true resphape the first return vec into (-1,1)

  • outspecdomain (|) – if not None interpolate output to this vector

Returns:

outspectraldomain and outspectralquantity | Returns zero length arrays is type is illegal, i.e. not one of the expected values

Return type:

([N,1],[N,1])

Raises:

| No exception is raised. –

convolve(inspectral, samplingresolution, inwinwidth, outwinwidth, windowtype=<function bartlett>)[source]

Convolve (non-circular) a spectral variable with a window function, given the input resolution and input and output window widths.

This function is normally used on wavenumber-domain spectral data. The spectral data is assumed sampled at samplingresolution wavenumber intervals. The inwinwidth and outwinwidth window function widths are full width half-max (FWHM) for the window functions for the inspectral and returned spectral variables, respectively. The Bartlett function is used as default, but the user can use a different function. The Bartlett function is a triangular function reaching zero at the ends. Window function width is correct for Bartlett and only approximate for other window functions.

Spectral convolution is best done in frequency domain ([cm-1] units) because the filter or emission line shapes have better symmetry in frequency domain than in wavelength domain.

The input spectral vector must be in spectral density units of cm-1.

Parameters:

  • inspectral (|) – spectral variable input vector (e.g., radiance or transmittance).

  • samplingresolution (|) – wavenumber interval between inspectral samples

  • inwinwidth (|) – FWHM window width used to obtain the input spectral vector (e.g., spectroradiometer window width)

  • outwinwidth (|) – FWHM window width of the output spectral vector after convolution

  • windowtype (|) – name of a numpy/scipy function for the window function

Returns:

input vector, filtered to new window width. | windowfn (np.array[N,]): The window function used.

Return type:

outspectral (np.array[N,])

Raises:

| No exception is raised. –

savitzkyGolay1D(y, window_size, order, deriv=0, rate=1)[source]

Smooth (and optionally differentiate) data with a Savitzky-Golay filter.

Source: https://scipy-cookbook.readthedocs.io/items/SavitzkyGolay.html

The Savitzky Golay filter is a particular type of low-pass filter, well adapted for data smoothing. For further information see: https://scipy-cookbook.readthedocs.io/items/SavitzkyGolay.html

The Savitzky-Golay filter removes high frequency noise from data. It has the advantage of preserving the original shape and features of the signal better than other types of filtering approaches, such as moving averages techniques.

The Savitzky-Golay is a type of low-pass filter, particularly suited for smoothing noisy data. The main idea behind this approach is to make for each point a least-square fit with a polynomial of high order over a odd-sized window centered at the point.

Examples

t = np.linspace(-4, 4, 500) y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape) ysg = savitzkyGolay1D(y, window_size=31, order=4) import matplotlib.pyplot as plt plt.plot(t, y, label=’Noisy signal’) plt.plot(t, np.exp(-t**2), ‘k’, lw=1.5, label=’Original signal’) plt.plot(t, ysg, ‘r’, label=’Filtered signal’) plt.legend() plt.show()

References

[1] A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of

Data by Simplified Least Squares Procedures. Analytical Chemistry, 1964, 36 (8), pp 1627-1639.

[2] Numerical Recipes 3rd Edition: The Art of Scientific Computing

W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery Cambridge University Press ISBN-13: 9780521880688

Parameters:

  • y (|) – array_like, shape (N,) the values of the time history of the signal.

  • window_size (|) – int the length of the window. Must be an odd integer number.

  • order (|) – int the order of the polynomial used in the filtering. Must be less then window_size - 1.

  • deriv (|) – int the order of the derivative to compute (default = 0 means only smoothing)

Returns:

ndarray, shape (N) the smoothed signal (or it’s n-th derivative).

Return type:

ys
Raises:
Exception raised for window size errors.

abshumidity(T, equationSelect=1)[source]

Atmopsheric absolute humidity [g/m3] for temperature in [K] between 248 K and 342 K.

This function provides two similar equations, but with different constants.

Parameters:

  • temperature (|) – in [K].

  • equationSelect (|) – select the equation to be used.

Returns:

abs humidity in [g/m3]

Return type:

absolute humidity (np.array[N,] or [N,1])

Raises:

| No exception is raised. –

TFromAbshumidity(AH, equationSelect=1)[source]

temperature in [K] between 248 K and 342 K, given atmopsheric absolute humidity [g/m3], assuming 100% RH

This function uses two similar equations, but with different constants.

Parameters:

  • AH (|) – absolute humidity in g/m3.

  • equationSelect (|) – select the equation to be used.

Returns:

in K

Return type:

temperature (float)

Raises:

| No exception is raised. –

rangeEquation(Intensity, Irradiance, rangeTab, tauTab, rangeGuess=1, n=2)[source]

Solve the range equation for arbitrary transmittance vs range.

This function solve for the range R in the range equation

E = \frac{I\tau_a(R)}{R^n}

where E is the threshold irradiance in [W/m2], and I is the intensity in [W/sr]. This range equation holds for the case where the target is smaller than the field of view.

The range R must be in [m], and \tau_a(R) is calculated from a lookup table of atmospheric transmittance vs. range. The transmittance lookup table can be calculated from the simple Bouguer law, or it can have any arbitrary shape, provided it decreases with increasing range. The user supplies the lookup table in the form of an array of range values and an associated array of transmittance values. The range values need not be on constant linear range increment.

The parameter n

  • n=2 (default value) the general case of a radiating source smaller than the field of view.

  • n=4 the special case of a laser range finder illuminating a target smaller than the field of view, viewed against the sky. In this case there is an R^2 attenuation from the laser to the source and another R^2 attenuation from the source to the receiver, hence R^4 overall.

If the range solution is doubtful (e.g. not a trustworthy solution) the returned value is made negative.

Parameters:

  • Intensity (|) – in [W/sr].

  • Irradiance (|) – in [W/m2].

  • rangeTab (|) – range vector for tauTab lookup in [m]

  • tauTab (|) – transmittance vector for lookup in [m]

  • rangeGuess (|) – starting value range estimate in [m] (optional)

  • n (|) – range power (2 or 4) (optional)

Returns:

Solution to the range equation in [m].

Value is negative if calculated range exceeds the top value in range table, or if calculated range is too near the lower resolution limit.

Return type:

range (float or np.array[N,] or [N,1])

Raises:

| No exception is raised. –

detectThresholdToNoiseTpFAR(pulseWidth, FAR)[source]

Solve for threshold to noise ratio, given pulse width and FAR, for matched filter.

Using the theory of matched filter design, calculate the threshold to noise ratio, to achieve a required false alarm rate.

References:

“Electro-optics handbook,” Tech. Rep. EOH-11, RCA, 1974. RCA Technical Series Publication.

    1. Hippenstiel, Detection Theory: Applications and Digital Signal Processing, CRC Press, 2002

Parameters:

  • pulseWidth (|) – the signal pulse width in [s].

  • FAR (|) – the false alarm rate in [alarms/s]

Returns:

threshold to noise ratio

Return type:

range (float)

Raises:

| No exception is raised. –

detectSignalToNoiseThresholdToNoisePd(ThresholdToNoise, pD)[source]

Solve for the signal to noise ratio, given the threshold to noise ratio and probability of detection.

Using the theory of matched filter design, calculate the signal to noise ratio, to achieve a required probability of detection.

References:

“Electro-optics handbook,” Tech. Rep. EOH-11, RCA, 1974. RCA Technical Series Publication.

    1. Hippenstiel, Detection Theory: Applications and Digital Signal Processing, CRC Press, 2002

Parameters:

  • ThresholdToNoise (|) – the threshold to noise ratio [-]

  • pD (|) – the probability of detection [-]

Returns:

signal to noise ratio

Return type:

range (float)

Raises:

| No exception is raised. –

detectThresholdToNoiseSignalToNoisepD(SignalToNoise, pD)[source]

Solve for the threshold to noise ratio, given the signal to noise ratio and probability of detection.

References:

“Electro-optics handbook,” Tech. Rep. EOH-11, RCA, 1974. RCA Technical Series Publication.

    1. Hippenstiel, Detection Theory: Applications and Digital Signal Pro-cessing, CRC Press, 2002

Parameters:

  • SignalToNoise (|) – the signal to noise ratio [-]

  • pD (|) – the probability of detection [-]

Returns:

signal to noise ratio

Return type:

range (float)

Raises:

| No exception is raised. –

detectProbabilityThresholdToNoiseSignalToNoise(ThresholdToNoise, SignalToNoise)[source]
Solve for the probability of detection, given the signal to noise ratio and

threshold to noise ratio

References:

“Electro-optics handbook,” Tech. Rep. EOH-11, RCA, 1974. RCA Technical Series Publication.

    1. Hippenstiel, Detection Theory: Applications and Digital Signal Pro-cessing, CRC Press, 2002

Parameters:

  • ThresholdToNoise (|) – the threshold to noise ratio [-]

  • SignalToNoise (|) – the signal to noise ratio [-]

Returns:

probability of detection

Return type:

range (float)

Raises:

| No exception is raised. –

detectFARThresholdToNoisepulseWidth(ThresholdToNoise, pulseWidth)[source]

Solve for the FAR, given the threshold to noise ratio and pulse width, for matched filter.

References:

“Electro-optics handbook,” Tech. Rep. EOH-11, RCA, 1974. RCA Technical Series Publication.

    1. Hippenstiel, Detection Theory: Applications and Digital Signal Processing, CRC Press, 2002

Parameters:

  • ThresholdToNoise (|) – the threshold to noise ratio.

  • pulseWidth (|) – the signal pulse width in [s].

Returns:

the false alarm rate in [alarms/s]

Return type:

FAR (float)

Raises:

| No exception is raised. –

buildLogSpace(Vmin, Vmax, nDec, patn=False)[source]

Calculate a log space given low, high and number samples per decade

If patn is True, the upper limit is adjusted to obtain a repeat numeric pattern in each dcade.

Parameters:

  • Vmin (|)

  • Vmax (|)

  • nDec (|)

  • patn (|)

Returns:

vector with equal spacing in log

Raises:

| No exception is raised. –

update_progress(progress, bar_length=20)[source]

Simple text-based progress bar for Jupyter notebooks.

Note that clear_output, and hence this function wipes the entire cell output, including previous output and widgets.

Usage:

import pyradi.ryutils as ryutils
import time
print('before')
#Replace this with a real computation
number_of_elements = 100
for i in range(number_of_elements):
    time.sleep(0.1)
    # progress must be a float between 0 and 1
    ryutils.update_progress((i+1) / number_of_elements,bar_length=40)
print('after')

source:

https://mikulskibartosz.name/how-to-display-a-progress-bar-in-jupyter-notebook-47bd4c2944bf https://ipython.org/ipython-doc/3/api/generated/IPython.display.html#IPython.display.clear_output Wait to clear the output until new output is available to replace it.

solidAngleSquare(width, breadth, height, stype, numsamples)[source]

Calculate the solid angle of a rectagular plate from a point on the normal at its centre

The solid angle of a rectangular flat surface, with dimensions $W$ and $D$, as seen from a reference point centered above the surface, is determined by the integral of the projected area of a small elemental area $costheta,dd,dw$ across the full size of the surface: $$ omega_{rm s}=int_Wint_Dfrac{dw,dd,cos^{n-2}theta}{R^2} $$ $$ omega_{rm s}=int_Wint_Dfrac{dw,dd,cos^ntheta}{H^2} $$ $$ omega_{rm s}=int_Wint_Dfrac{dw,dd,}{H^2}left(frac{H}{R}right)^n $$ $$omega_{rm s}=int_Wint_Dfrac{dw,dd,}{H^2}left(frac{H}{sqrt{w^2+d^2+H^2}}right)^n, $$ where $H$ is the reference point height above the surface, and $n=3$ for the geometrical solid angle and $n=4$ for the projected solid angle. The integral is performed along the $W$ and $D$ dimensions with increments of $dw$ and $dd$. The slant range between the reference point and the elemental area $ddtimes dw$ is $R=H/costheta$.

Parameters:

  • width (|) – size along one edge of rectangle

  • breadth (|) – size along the second edge of rectangle

  • height (|) – distance along normal to the rect to reference point

  • stype (|) – type of solid angle can be one of (‘g’ or ‘p’) for (‘geometric’,’projected’)

  • numsamples (|) – number of samples along edges

Returns:

solid angle (float) or None if incorrect type

Raises:

| No exception is raised. –

getFHWM(wl, tau, normaliseMax=False)[source]

Determine the full-width half-maximum (FWHM) of a spectral curve.

Finds the two wavelengths where the curve crosses the 50% level and returns the distance between them, as well as the crossing wavelengths themselves.

Parameters:

  • wl (np.ndarray) – Spectral domain vector (e.g. wavelength in µm).

  • tau (np.ndarray) – Spectral quantity vector (same length as wl).

  • normaliseMax (bool) – If True, normalise tau to its maximum before computing the half-max crossing. Default False.

Returns:

(fwhm, lower_wl, upper_wl) where fwhm is

the full-width half-maximum and lower_wl / upper_wl are the two 50%-crossing wavelengths.

Return type:

tuple[float, float, float]

Raises:

No exception is raised. –

upMu(uprightMu=True, textcomp=False)[source]

Returns a LaTeX micron symbol, either an upright version or the normal symbol.

The upright symbol requires that the siunitx LaTeX package be installed on the computer running the code. This function also changes the Matplotlib rcParams file.

Parameters:

  • uprightMu (|) – signals upright (True) or regular (False) symbol (optional).

  • textcomp (|) – if True use the textcomp package, else use siunitx package (optional).

Returns:

LaTeX code for the micro symbol.

Return type:

range (string)

Raises:

| No exception is raised. –

cart2polar(x, y)[source]

Converts from cartesian to polar coordinates, given (x,y) to (r,theta).

Parameters:

  • x (|) – x values in array format.

  • y (|) – y values in array format.

Returns:

radial component for given (x,y). | theta (float np.array): angular component for given (x,y).

Return type:

r (float np.array)

Raises:

| No exception is raised. –

original code by Joe Kington https://stackoverflow.com/questions/3798333/image-information-along-a-polar-coordinate-system

polar2cart(r, theta)[source]

Converts from polar to cartesian coordinates, given (r,theta) to (x,y).

Parameters:

  • r (|) – radial values in array format.

  • theta (|) – angular values in array format.

Returns:

x component for given (r, theta). | y (float np.array): y component for given (r, theta).

Return type:

x (float np.array)

Raises:

| No exception is raised. –

original code by Joe Kington https://stackoverflow.com/questions/3798333/image-information-along-a-polar-coordinate-system

polar2cartesian(xycoords, inputshape, origin)[source]

Converting polar (r,theta) array indices to Cartesian (x,y) array indices.

This function is called from scipy.ndimage.geometric_transform which calls this function as follows: polar2cartesian: A callable object that accepts a tuple of length equal to the output array rank, and returns the corresponding input coordinates as a tuple of length equal to the input array rank.

theta goes from 0 to 2pi. the x,y coords maps from -r to +r.

For an example application, see the function warpPolarImageToCartesianImage below

http://stackoverflow.com/questions/2164570/reprojecting-polar-to-cartesian-grid note that we changed the code from the original

Parameters:

  • xycoords (|) – x,y values for which r,theta coords must be found

  • inputshape (|) – shape of the polar input array

  • origin (|) – x and y indices of where the origin should be in the output array

Returns:

indices into the the r,theta array corresponding to xycoords

Return type:

r,theta_index (tuple)

Raises:

| No exception is raised. –

cartesian2polar(rtcoords, inputshape, origin)[source]

Converting Cartesian (x,y) to polar (r,theta) array indices.

This function is called from scipy.ndimage.geometric_transform which calls this function as follows: cartesian2polar: A callable object that accepts a tuple of length equal to the output array rank, and returns the corresponding input coordinates as a tuple of length equal to the input array rank.

the x,y coords maps from -r to +r. theta goes from 0 to 2pi.

Parameters:

  • rtcoords (|) – r, theta values for which xy coords must be found

  • inputshape (|) – shape of the cartesian input array

  • origin (|) – x and y indices of where the origin should be in the output array

Returns:

indices into the the r,theta array corresponding to xycoords

Return type:

r,theta_index (tuple)

Raises:

| No exception is raised. –

index_coords(data, origin=None, framesFirst=True)[source]

Creates (x,y) zero-based coordinate arrrays for a numpy array indices, relative to some origin.

This function calculates two meshgrid arrays containing the coordinates of the input array. 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.

Parameters:

  • data (|) – array for which coordinates must be calculated.

  • origin (|) – data-coordinates of where origin should be

  • framesFirst (|) – True if data.shape is (frames, rows, cols), False if data.shape is (rows, cols, frames)

Returns:

x coordinates in array format. | y (float np.array): y coordinates in array format.

Return type:

x (float np.array)

Raises:

| No exception is raised. –

original code by Joe Kington https://stackoverflow.com/questions/3798333/image-information-along-a-polar-coordinate-system

framesFirst(imageSequence)[source]

Image sequence with frames along axis=2 (last index), reordered such that frames are along axis=0 (first index).

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. This function reorders an image sequence with frames along axis=2 to an image sequence with frames along axis=0. The function uses np.transpose(imageSequence, (2,0,1))

Parameters:

imageSequence (|) – image sequence in three-dimensional array, frames along axis=2

Returns:

reordered three-dimensional array (view or copy)

Return type:

((3-D np.array)

Raises:

| No exception is raised. –

framesLast(imageSequence)[source]

Image sequence with frames along axis=0 (first index), reordered such that frames are along axis=2 (last index).

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. This function reorders an image sequence with frames along axis=0 to an image sequence with frames along axis=2. The function uses np.transpose(imageSequence, (1,2,0))

Parameters:

imageSequence (|) – image sequence in three-dimensional array, frames along axis=0

Returns:

reordered three-dimensional array (view or copy)

Return type:

((3-D np.array)

Raises:

| No exception is raised. –

rect(x, y, sx=1, sy=1)[source]

Generation of a rectangular aperture.

Parameters:

  • x (|) – x-grid, metres

  • y (|) – x-grid, metres

  • sx (|) – full size along x.

  • sy (|) – full size along y.

Returns:

Nothing.

Raises:

| No exception is raised. –

Author: CJ Willers

Original source: http://arxiv.org/pdf/1412.4031.pdf

circ(x, y, d=1)[source]

Generation of a circular aperture.

Parameters:

  • x (|) – x-grid, metres

  • y (|) – y-grid, metres

  • d (|) – diameter in metres.

  • comment (|) – the symbol used to comment out lines, default value is None.

  • delimiter (|) – delimiter used to separate columns, default is whitespace.

Returns:

z-grid, 1’s inside radius, meters/pixels.

Return type:

z (np.array[N,M])

Raises:

| No exception is raised. –

Author: Prof. Jason Schmidt, revised/ported by CJ Willers

Original source: http://arxiv.org/pdf/1412.4031.pdf

poissonarray(inp, seedval=None, tpoint=1000)[source]

This routine calculates a Poisson random variable for an array of input values with potentially very high event counts.

At high mean values the Poisson distribution calculation overflows. For mean values exceeding 1000, the Poisson distribution may be approximated by a Gaussian distribution.

The function accepts a two-dimensional array and calculate a separate random value for each element in the array, using the element value as the mean value. A typical use case is when calculating shot noise for image data.

From http://en.wikipedia.org/wiki/Poisson_distribution#Related_distributions For sufficiently large values of \lambda, (say \lambda>1000), the normal distribution with mean \lambda and variance \lambda (standard deviation \sqrt{\lambda}) is an excellent approximation to the Poisson distribution. If \lambda is greater than about 10, then the normal distribution is a good approximation if an appropriate continuity correction is performed, i.e., P(X \le x), where (lower-case) x is a non-negative integer, is replaced by P(X\le\,x+0.5).

F_\mathrm{Poisson}(x;\lambda)\approx\,F_\mathrm{normal}(x;\mu=\lambda,\sigma^2=\lambda)

This function returns values of zero when the input is zero.

Parameters:

  • inp (|) – array with mean value

  • seedval (|) – seed for random number generator, None means use system time.

  • tpoint (|) – Threshold when to switch over between Poisson and Normal distributions

Returns:

Poisson random variable for given mean value

Return type:

outp (np.array[N,M])

Raises:

| No exception is raised. –

Author: CJ Willers

draw_siemens_star(outfile, n, dpi)[source]

Siemens star chart generator

by Libor Wagner, http://cmp.felk.cvut.cz/~wagnelib/utils/star.html

Parameters:

  • outfile (|) – output image filename (monochrome only)

  • n (|) – number of spokes in the output image.

  • dpi (|) – dpi in output image, determines output image size.

Returns:

Nothing, creates a monochrome siemens star image

Raises:

| No exception is raised. –

Author: Libor Wagner, adapted by CJ Willers

gen_siemens_star(origin, radius, n)[source]

Generate a Siemens star patch collection for use with matplotlib.

Creates a PatchCollection of alternating black wedges arranged in a circle, which form the spokes of a Siemens star test pattern.

Parameters:

  • origin (tuple[float, float]) – (x, y) centre of the star in data coordinates.

  • radius (float) – Outer radius of the star in data units.

  • n (int) – Number of spoke pairs (the star will have 2*n sectors).

Returns:

Collection of black wedge patches

ready to be added to an Axes with ax.add_collection().

Return type:

matplotlib.collections.PatchCollection

Raises:

No exception is raised. –

drawCheckerboard(rows, cols, numPixInBlock, imageMode, colour1, colour2, imageReturnType='image', datatype=<class 'numpy.uint8'>)[source]

Draw checkerboard with 8-bit pixels

From http://stackoverflow.com/questions/2169478/how-to-make-a-checkerboard-in-numpy

Parameters:

  • rows (|) – number or rows in checkerboard

  • cols (|) – number of columns in checkerboard

  • numPixInBlock (|) – number of pixels to be used in one block of the checkerboard

  • imageMode (|) – PIL image mode [e.g. L (8-bit pixels, black and white),

  • RGB (|)

  • colour1 (|) – colour 1 specified according to the imageMode

  • colour2 (|) – colour 2 specified according to the imageMode

  • imageReturnType (|) – ‘image’ for PIL image, ‘nparray’ for numpy array

  • datatype (|) – numpy data type for the returned np.array

makemotionsequence(imgfilename, mtnfilename, postfix, intTime, frmTim, outrows, outcols, imgRatio, pixsize, numsamples, fnPlotInput=None)[source]

Builds a video from a still image and a displacement motion file.

The objective with this function is to create a video sequence from a still image, as if the camera moved minutely during the sensor integration time.

A static image is moved according to the (x,y) displacement motion in an input file. The input file must be at least ten times plus a bit larger than the required output file. The image input file is sampled with appropriate displacement for each point in the displacement file and pixel vlaues are accumulated in the output image. All of this temporal displacement and accumulation takes place in the context of a frame integration time and frame frequency.

The key requirements for accuracy in this method is an input image with much higher resolution than the output image, plus a temporal displacement file with much higher temporal sampling than the sensor integration time.

The function creates a sequence of images that can be used to create a video. Images are numbered in sequence, using the same base name as the input image. The sequence is generated in the current working directory.

The function currently processes only monochrome images (M,N) arrays.

The motion data file must be a compressed numpy npz or text file, with three columns: First column must be time, then movement along rows, then movement along columns. The units and scale of the motion columns must be the same units and scale as the pixel size in the output image.

imgRatio x imgRatio number of pixels in the input (hires) image are summed together and stored in one output image pixel. In other words if imgRatio is ten, each pixel in the output image will be the sum of 100 pixels in the imput image. During one integration time period the hires input image will be sampled at slightly different offsets (according to the motion file) and accumulated in an intermediate internal hires file. This intermediate internal file is collapsed as described above.

The function creates a series-numbered sequence if images that can be used to construct a video. One easy means to create the video is to use VirtualDub, available at www.virtualdub.org/index. In VirtualDub open the first image file in the numbered sequence, VirtualDub will then recognise the complete sequence as a video. Once loaded in VirtualDub, save the video as avi.

Parameters:

  • imgfilename (|) – static image filename (monochrome only)

  • mtnfilename (|) – motion data filename.

  • postfix (|) – add this string to the end of the output filename.

  • intTime (|) – sensor integration time.

  • frmTim (|) – sensor frame time.

  • outrows (|) – number of rows in the output image.

  • outcols (|) – number of columns in the output image.

  • imgRatio (|) – hires image pixel count block size of one output image pixel

  • pixsize (|) – pixel size in same units as motion file.

  • numsamples (|) – number of motion input samples to be processed (-1 for all).

  • fnPlotInput (|) – output plot filename (None for no plot).

Returns:

True if successful, message otherwise, creates numbered images in current working directory

Raises:

| No exception is raised. –

Author: CJ Willers

warpPolarImageToCartesianImage(mesh, order=0)[source]

Convert an image in (r,theta) format to (x,y) format

The 0th and 1st axes correspond to r and theta, respectively. theta goes from 0 to 2pi, and r’s units are just its indices. output image is twice the size of r length in both x and y

http://stackoverflow.com/questions/2164570/reprojecting-polar-to-cartesian-grid

Example code: size = 100 dset = np.random.random((size,size)) mesh_cart = warpPolarImageToCartesianImage(dset) p = ryplot.Plotter(1,1,2); p.showImage(1, dset); p.showImage(2, mesh_cart);

Parameters:

mesh (|) – array in r,theta coordinates

Returns:

array in x,y coordinates

Return type:

mesh_cart (np.array)

Raises:

| No exception is raised. –

warpCartesianImageToPolarImage(mesh, order=0)[source]

Convert an image in (x,y) format to (r,theta) format

The 0th and 1st axes correspond to x and y, respectively. size along x and y must be equal and the conversion origin is assumed in the center of the image. theta goes from 0 to 2pi, and r’s units are just its indices. output image is the same length as x and y in one axes and sqrt(2) bigger in the other axis.

Parameters:

mesh (|) – array in x,y coordinates

Returns:

array in r,theta coordinates

Return type:

mesh_Polar (np.array)

Raises:

| No exception is raised. –

extractGraph(filename, xmin, xmax, ymin, ymax, outfile=None, doPlot=False, xaxisLog=False, yaxisLog=False, step=None, value=None)[source]

Scan an image containing graph lines and produce (x,y,value) data.

This function processes an image, calculate the location of pixels on a graph line, and then scale the (r,c) or (x,y) values of pixels with non-zero values. The

Get a bitmap of the graph (scan or screen capture). Take care to make the graph x and y axes horizontal/vertical. The current version of the software does not work with rotated images. Bitmap edit the graph. Clean the graph to the maximum extent possible, by removing all the clutter, such that only the line to be scanned is visible. Crop only the central block that contains the graph box, by deleting the x and y axes notation and other clutter. The size of the cropped image must cover the range in x and y values you want to cover in the scan. The graph image/box must be cut out such that the x and y axes min and max correspond exactly with the edges of the bitmap. You must end up with nothing in the image except the line you want to digitize.

The current version only handles single lines on the graph, but it does handle vertical and horizontal lines.

The function can also write out a value associated with the (x,y) coordinates of the graph, as the third column. Normally these would have all the same value if the line represents an iso value.

The x,y axes can be lin/lin, lin/log, log/lin or log/log, set the flags.

Parameters:

  • filename (|) – name of the image file

  • xmin (|) – the value corresponding to the left side (column=0)

  • xmax (|) – the value corresponding to the right side (column=max)

  • ymin (|) – the value corresponding to the bottom side (row=bottom)

  • ymax (|) – the value corresponding to the top side (row=top)

  • outfile (|) – write the sampled points to this output file

  • doPlot (|) – plot the digitised graph for visual validation

  • xaxisLog (|) – x-axis is in log10 scale (min max are log values)

  • yaxisLog (|) – y-axis is in log10 scale (min max are log values)

  • step (|) – if not None only ouput every step values

  • value (|) – if not None, write this value as the value column

Returns:

a numpy array with columns (xval, yval, value) | side effect: a file may be written | side effect: a graph may be displayed

Return type:

outA

Raises:

| No exception is raised. –

Author: neliswillers@gmail.com

luminousEfficiency(vlamtype='photopic', wavelen=None, eqnapprox=False)[source]

Returns the photopic luminous efficiency function on wavelength intervals

Type must be one of:

photopic: CIE Photopic V(lambda) modified by Judd (1951) and Vos (1978) [also known as CIE VM(lambda)] scotopic: CIE (1951) Scotopic V’(lambda) CIE2008v2: 2 degree CIE “physiologically-relevant” luminous efficiency Stockman & Sharpe CIE2008v10: 10 degree CIE “physiologically-relevant” luminous efficiency Stockman & Sharpe

For the equation approximations (only photoic and scotopic), if wavelength is not given a vector is created 0.3-0.8 um.

For the table data, if wavelength is not given a vector is read from the table.

CIE Photopic V(l) modified by Judd (1951) and Vos (1978) [also known as CIE VM(l)] from http://www.cvrl.org/index.htm

Parameters:

  • vlamtype (|) – type of curve required

  • wavelen (|) – wavelength in um

  • eqnapprox (|) – if False read tables, if True use equation

Returns:

luminous efficiency | wavelen (np.array[]): wavelength in um

Return type:

luminousEfficiency (np.array[])

Raises:

| No exception is raised. –

Author: CJ Willers

calcMTFwavefrontError(sample, wfdisplmnt, xg, yg, specdef, samplingStride=1, clear='Clear')[source]

Given a mirror figure error, calculate MTF degradation from ideal

An aperture has an MTF determined by its shape. A clear aperture has zero phase delay and the MTF is determined only by the aperture shape. Any phase delay/error in the wavefront in the aperture will result in a

lower MTF than the clear aperture diffraction MTF.

This function calculates the MTF degradation attributable to a wavefront error, relative to the ideal aperture MTF.

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

xi = {x}/{lambda d_i}

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

xi_{cutoff} ={D}/{lambda f}

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.

The the array containing the wavefront displacement in the pupil must have np.nan values outside the pupil. The np.nan values are ignored and not included in the calculation. Obscurations can be modelled by placing np.nan in the obscuration.

The specdef dictionary has a string key to identify (name) the band, with a single float contents which is the wavelength associated with this band.

Parameters:

  • sample (|) – an identifier string to be used in the plots

  • wfdisplmnt (|) – wavefront displacement in m

  • xg (|) – x values from meshgrid, for wfdisplmnt

  • yg (|) – y values from meshgrid, for wfdisplmnt

  • specdef (|) – dictionary defining spectral wavelengths

  • samplingStride (|) – sampling stride to limit size and processing

  • clear (|) – defines the dict key for clear aperture reference

Returns:

dictionaries below have entries for all keys in specdef.
wfdev (dict): subsampled wavefront error in m
phase (dict): subsampled wavefront error in rad
pupfn (dict): subsampled complex pupil function
MTF2D (dict): 2D MTF in (x,y) format
MTFpol (dict): 2D MTF in (r,theta) format
specdef (): specdef dictionary as passed plus clear entry
MTFmean (dict): mean MTF across all rotation angles
rho (nd.array[M,]): spatial frequency scale in cy/mrad
fcrit (float): cutoff or critical spatial frequency cy/mrad
clear (string): key used to signify the clear aperture case.

Raises:

| No exception is raised. –

differcommonfiles(dir1, dir2, patterns='*')[source]

Find if common files in two dirs are different Directories are not recursed, only common files are binary compared

Parameters:

  • filename (|) – first directory name

  • filename – first directory name

Returns:

list of common files | different (list) : list of flags of different flags

Return type:

file (list)

Raises:

| No exception is raised. –

blurryextract(img, inputOrigin, outputSize, blocksize, sigma=0.0, filtermode='reflect')[source]

Slice region from 2d array, blur it with gaussian kernel, and coalesce to lower resolution

The image is blurred with a gaussian kernel of size sigma, using filtermode. Then the slice is calculated using the inputOrigin, required output size and the block size. The sliced image is then lowered in resolution by averaging blocks of values in the input image into the output image.

If the input image size/bounds are to be exceeded in the calculation the function returns None.

https://numpy.org/doc/stable/reference/generated/numpy.einsum.html
http://ajcr.net/Basic-guide-to-einsum/
https://obilaniu6266h16.wordpress.com/2016/02/04/einstein-summation-in-numpy/
https://stackoverflow.com/questions/26089893/understanding-numpys-einsum
https://stackoverflow.com/questions/26064594/how-does-numpy-einsum-work

http://stackoverflow.com/questions/16856788/slice-2d-array-into-smaller-2d-arrays
http://stackoverflow.com/questions/21220942/sums-of-subarrays
https://en.wikipedia.org/wiki/Einstein_notation
https://numpy.org/doc/stable/reference/generated/numpy.einsum.html
https://stackoverflow.com/questions/28652976/passing-array-range-as-argument-to-a-function
https://stackoverflow.com/questions/13706258/passing-python-slice-syntax-around-to-functions

https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html
Parameters:

  • img (|) – 2D input array

  • inputOrigin (|) – start coordinates in the input image

  • outputSize (|) – size of the output array

  • blocksize (|) – input image block size to be

  • image (| averaged into single output)

  • sigma (|) – gaussian kernel rms size in pixel counts

  • sigma – gaussian kernel filter mode

Returns:

smaller image or None if out of bounds

Return type:

imgo (nd.array of floats)

Raises:

| No exception is raised. –

intify_tuple(tup)[source]

Make tuple entries int type