Back End Code

Using Directives

The following is the complete set of Using Directives I used for all back-end code for this site. Not all code examples use all Using Directives, but if you are trying to duplicate this code in your code editor, use all of these Using Directives at first until you become more familiar with each one's usages.

using System;
using System.Web;
using System.Numerics;
using Accord.Math;
using Accord.Math.Transforms;
using Accord.Math.Integration;
using System.Web.Script.Serialization;
using System.Drawing;
using System.Drawing.Imaging;
using System.Runtime.InteropServices;

arraySize

Nearly every function uses a static readonly variable called arraySize. It is used to define the dimensions of the Fourier Transform, and thus the size of the resulting images. I chose 512 because FFTs operate on powers of 2, and 512 is the best mix of meaningful resolution and computation speed.

private static readonly int arraySize = 512;

Accord.NET Framework

I use functions from the Accord.NET Framework extensively throughout my code, especially Accord.Math.Transforms.FourierTransform2

FFT Shift

// Description  :   This function performs a 2D FFT shift equivalent to MATLAB's fftshift function.
// Inputs       :   shiftme - the array to be shifted
// Outputs      :   Results stored in referenced array shiftme

public static void fftshift(ref Complex[][] shiftme)
{
    Complex[,] quad1 = new Complex[arraySize / 2, arraySize / 2];
    Complex[,] quad2 = new Complex[arraySize / 2, arraySize / 2];
    Complex[,] quad3 = new Complex[arraySize / 2, arraySize / 2];
    Complex[,] quad4 = new Complex[arraySize / 2, arraySize / 2];

    for (int i = 0; i < (arraySize / 2); i++)
    {
        for (int j = 0; j < (arraySize / 2); j++)
        {
            quad1[i, j] = shiftme[i][j];
            quad2[i, j] = shiftme[i + (arraySize / 2)][j];
            quad3[i, j] = shiftme[i][j + (arraySize / 2)];
            quad4[i, j] = shiftme[i + (arraySize / 2)][j + (arraySize / 2)];
            shiftme[i][j] = quad4[i, j];
            shiftme[i + (arraySize / 2)][j] = quad3[i, j];
            shiftme[i][j + (arraySize / 2)] = quad2[i, j];
            shiftme[i + (arraySize / 2)][j + (arraySize / 2)] = quad1[i, j];
        }
    }
}

Mesh Grid

// Description  :   This function generates a 2D meshgrid equivalent to MATLAB's meshgrid function.
// Inputs       :   f - array of numbers defining meshgrid space
//              :   X - matrix of f in X-direction
//              :   Y - matrix of f in Y-direction
// Outputs      :   Results stored in referenced arrays X and Y

public static void meshgrid(double[] f, ref double[,] X, ref double[,] Y)
{
    for (int i = 0; i < arraySize; i++)
    {
        for (int j = 0; j < arraySize; j++)
        {
            Y[i, j] = f[i];
            X[j, i] = f[i];
        }
    }
}

Circular Shift

// Description  :   This function performs a 2D circular shift on a 
//                  matrix equivalent to MATLAB's circshift function.
// Inputs       :   old_array - matrix to be circularly shifted
//              :   xshift - amount to circularly shift in x-direction
//              :   yshift - amount to circularly shift in y-direction
// Outputs      :   new_array - the circularly shifted matrix

public static double[,] circshift(ref double[,] old_array, int xshift, int yshift)
{
    double[,] temp_array = new double[arraySize, arraySize];
    double[,] new_array = new double[arraySize, arraySize];

    for (int i = 0; i < arraySize; i++)
    {
        if (xshift > 0)
        {
            for (int j = 0; j < arraySize; j++)
            {
                if (j >= xshift)
                {
                    temp_array[i, j] = old_array[i, j - xshift];
                }
                else
                {
                    temp_array[i, j] = old_array[i, j - xshift + arraySize];
                }
            }
        }
        else
        {
            for (int j = 0; j < arraySize; j++)
            {
                if (j <= (arraySize + xshift - 1))
                {
                    temp_array[i, j] = old_array[i, j - xshift];
                }
                else
                {
                    temp_array[i, j] = old_array[i, j - xshift - arraySize];
                }
            }
        }
    }

    for (int i = 0; i < arraySize; i++)
    {
        if (yshift > 0)
        {
            for (int j = 0; j < arraySize; j++)
            {
                if (i >= yshift)
                {
                    new_array[i, j] = temp_array[i - yshift, j];
                }
                else
                {
                    new_array[i, j] = temp_array[i - yshift + arraySize, j];
                }
            }
        }
        else
        {
            for (int j = 0; j < arraySize; j++)
            {
                if (i <= (arraySize + yshift - 1))
                {
                    new_array[i, j] = temp_array[i - yshift, j];
                }
                else
                {
                    new_array[i, j] = temp_array[i - yshift - arraySize, j];
                }
            }
        }
    }
    return new_array;
}

Matrix Multiplication

// Description  :   This function returns the element-wise product of two square matricies
// Inputs       :   a,b - the matricies to multiply
// Outputs      :   returnValue - the element-wise product of a and b

public static Complex[][] matrixMultiply(Complex[][] a, Complex[][] b)
{
    int size = a.GetLength(0);
    Complex[][] returnValue = new Complex[size][];
    for (int i = 0; i < size; i++)
    {
        returnValue[i] = new Complex[size];
        for (int j = 0; j < size; j++)
        {
            returnValue[i][j] = a[i][j] * b[i][j];
        }
    }
    return returnValue;
}

Complex Magnitude

// Description  :   This function returns the complex magnitudes of a
//                  complex number array, raised to the specified power
// Inputs       :   a - a 2D array containing containing complex numbers
//              :   power - element-wise raise the complex magnitude aray
//                  by this number
// Outputs      :   returnValue - a matrix of complex magnitudes

public static double[,] complexMag(Complex[][] a, double power)
{
    double[,] returnValue = new double[arraySize, arraySize];
    for (int i = 0; i < arraySize; i++)
    {
        for (int j = 0; j < arraySize; j++)
        {
            returnValue[i, j] = Math.Pow(a[i][j].Magnitude, power);
        }
    }
    return returnValue;
}

Double To Complex Conversion

// Description  :   This function converts a 2D matrix of type double[,] to a 2D array of type Complex[][]
// Inputs       :   a = a matrix of doubles
// Outputs      :   returnValue = a 2D array of complex numbers where the real part of each
//                  element is the corresponding element of a 

public static Complex[][] double2Complex(double[,] a)
{
    int size = a.GetLength(0);
    Complex[][] returnValue = new Complex[size][];
    for (int i = 0; i < size; i++)
    {
        returnValue[i] = new Complex[size];
        for (int j = 0; j < size; j++)
        {
            returnValue[i][j] = new Complex(a[i, j], 0);
        }
    }
    return returnValue;
}

Real Part of Complex

// Description  :   This function extracts the real part of a 2D Complex array.
// Inputs       :   a - a 2D array of complex numbers
// Outputs      :   returnValue - a matrix of doubles where each element is the real
//                  part of the corresponding element of a

public static double[,] realPart(Complex[][] a)
{
    int size = a.GetLength(0);
    double[,] returnValue = new double[size, size];
    for (int i = 0; i < size; i++)
    {
        for (int j = 0; j < size; j++)
        {
            returnValue[i, j] = a[i][j].Real;
        }
    }
    return returnValue;
}

Circular Pupil

// Description  :   This function generates a circular pupil.
// Inputs       :   L1 - the pupil field side length
//              :   diameter - the diameter of the pupil
// Outputs      :   Uin - the circular pupil

public static Complex[][] circ(double L1, double diameter)
{
    double dx1 = L1 / arraySize;
    double[] x = new double[arraySize];
    for (int i = 0; i < arraySize; i++)
    {
        x[i] = -L1 / 2.0 + i * dx1;
    }
    double[,] x1 = new double[arraySize, arraySize];
    double[,] y1 = new double[arraySize, arraySize];
    meshgrid(x, ref x1, ref y1);
    Complex[][] Uin = new Complex[arraySize][];
    for (int i = 0; i < arraySize; i++)
    {
        Uin[i] = new Complex[arraySize];
        for (int j = 0; j < arraySize; j++)
        {
            if (Math.Sqrt(Math.Pow(x1[i,j],2) + Math.Pow(y1[i,j],2)) < (diameter/2))
            {
                Uin[i][j] = new Complex(1.0, 0.0);
            }
            else
            {
                Uin[i][j] = new Complex(0.0, 0.0);
            }
        }
    }
    return Uin;
}

Rectangular Pupil

// Description  :   This function generates a rectangular pupil.
// Inputs       :   L1 - the pupil field side length
//              :   height - the height of the pupil
//              :   width - the width of the pupil
// Outputs      :   Uin - the rectangular pupil

public static Complex[][] rect(double L1, double height, double width)
{
    double dx1 = L1 / arraySize;
    double[] x = new double[arraySize];
    for (int i = 0; i < arraySize; i++)
    {
        x[i] = -L1 / 2.0 + i * dx1;
    }
    double[,] x1 = new double[arraySize, arraySize];
    double[,] y1 = new double[arraySize, arraySize];
    meshgrid(x, ref x1, ref y1);
    Complex[][] Uin = new Complex[arraySize][];
    for (int i = 0; i < arraySize; i++)
    {
        Uin[i] = new Complex[arraySize];
        for (int j = 0; j < arraySize; j++)
        {
            if (x1[i,j] > -width/2 && x1[i,j] < width/2 && y1[i,j] > -height/2 && y1[i,j] < height/2)
            {
                Uin[i][j] = new Complex(1.0, 0.0);
            }
            else
            {
                Uin[i][j] = new Complex(0.0, 0.0);
            }
        }
    }
    return Uin;
}

Obscured Circular Pupil

// Description  :   This function generates a circularly obscured pupil.
// Inputs       :   L1 - the pupil field side length
//              :   d - outer diameter of the pupil
//              :   d2 - obscuration ratio of pupil
// Outputs      :   Uin - the circularly obscured pupil

public static Complex[][] obscuredCirc(double L1, double d, double d2)
{
    double dx1 = L1 / arraySize;
    double[] x = new double[arraySize];
    for (int i = 0; i < arraySize; i++)
    {
        x[i] = -L1 / 2.0 + i * dx1;
    }
    double[,] x1 = new double[arraySize, arraySize];
    double[,] y1 = new double[arraySize, arraySize];
    meshgrid(x, ref x1, ref y1);
    Complex[][] Uin = new Complex[arraySize][];
    for (int i = 0; i < arraySize; i++)
    {
        Uin[i] = new Complex[arraySize];
        for (int j = 0; j < arraySize; j++)
        {
            if ((Math.Sqrt(Math.Pow(x1[i, j], 2) + Math.Pow(y1[i, j], 2)) > (d / 2.0)) || (Math.Sqrt(Math.Pow(x1[i, j], 2) + Math.Pow(y1[i, j], 2)) < (d2 * d / 2.0)))
            {
                Uin[i][j] = new Complex(0.0, 0.0);
            }
            else
            {
                Uin[i][j] = new Complex(1.0, 0.0);
            }
        }
    }
    return Uin;
}

Fresnel Transfer Function Method

The Fresnel diffraction pattern using the Transfer Function method is computed by $$U_2(x,y) = \left | \Im ^{-1}\left ( \Im \left ( U_1(x,y) \right )H(x,y) \right ) \right |^2$$ where $H$ is the transfer function $$H(x,y) = e^{ikz}e^{-i \pi \lambda z (x^2 + y^2)}$$ The Transfer Function method is implemented using the following function:

// Description  :   This function generates a Fresnel Diffraction pattern
//                  using the Fresnel Transfer function method
// Inputs       :   u1 - the pupil
//              :   L1 - pupil plane side length
//              :   lambda - wavelength
//              :   z - pupil-to-observation plane distance
// Outputs      :   u2 - observation plane diffraction pattern

public static Complex[][] propTF(Complex[][] u1, double L1, double lambda, double z)
{
    double dx = L1 / Convert.ToDouble(arraySize);
    double[] fx = new double[arraySize];
    double[,] FX = new double[arraySize, arraySize];
    double[,] FY = new double[arraySize, arraySize];
    Complex[][] H = new Complex[arraySize][];
    for (int i = 0; i < arraySize; i++)
    {
        fx[i] = -1.0 / (2.0 * dx) + i / L1;
    }
    meshgrid(fx, ref FX, ref FY);

    for (int i = 0; i < arraySize; i++)
    {
        H[i] = new Complex[arraySize];
        for (int j = 0; j < arraySize; j++)
        {
            H[i][j] = Complex.Exp(-Complex.ImaginaryOne * Math.PI * lambda * z * (FX[i, j] * FX[i, j] + FY[i, j] * FY[i, j]));
        }
    }

    fftshift(ref H);
    fftshift(ref u1);
    FourierTransform2.FFT2(u1, FourierTransform.Direction.Forward);
    Complex[][] u2 = matrixMultiply(H, u1);
    FourierTransform2.FFT2(u2, FourierTransform.Direction.Backward);
    fftshift(ref u2);
    return u2;
}

Fresnel Impulse Response Method

The Fresnel diffraction pattern using the Impulse Response method is computed by $$U_2(x,y) = \left | \Im ^{-1}\left ( \Im \left ( U_1\left (x,y\right ) \right )\Im \left ( h\left (x,y\right ) \right ) \right ) \right |^2$$ where $h$ is the impulse response $$h(x,y) = \frac{e^{ikz}}{i \lambda z}e^{\frac{ik}{2z}(x^2+y^2)}$$ The Impulse Response method is implemented using the following function:

// Description  :   This function generates a Fresnel Diffraction pattern
//                  using the Fresnel Impulse Response method
// Inputs       :   u1 - the pupil
//              :   L1 - pupil plane side length
//              :   lambda - wavelength
//              :   z - pupil-to-observation plane distance
// Outputs      :   u2 - observation plane diffraction pattern

public static Complex[][] propIR(Complex[][] u1, double L, double lambda, double z)
{
    double dx = L / Convert.ToDouble(arraySize);
    double k = 2 * Math.PI / lambda;
    double[] x = new double[arraySize];
    double[,] X = new double[arraySize, arraySize];
    double[,] Y = new double[arraySize, arraySize];
    Complex[][] H = new Complex[arraySize][];

    for (int i = 0; i < arraySize; i++)
    {
        x[i] = -L / 2.0 + i * dx;
    }

    meshgrid(x, ref X, ref Y);

    for (int i = 0; i < arraySize; i++)
    {
        H[i] = new Complex[arraySize];
        for (int j = 0; j < arraySize; j++)
        {
            H[i][j] = 1 / (Complex.ImaginaryOne * lambda * z) * Complex.Exp(Complex.ImaginaryOne * k / (2 * z) * (X[i, j] * X[i, j] + Y[i, j] * Y[i, j]));
        }
    }

    fftshift(ref H);
    FourierTransform2.FFT2(H, FourierTransform.Direction.Forward);

    for (int i = 0; i < arraySize; i++)
    {
        for (int j = 0; j < arraySize; j++)
        {
            H[i][j] = H[i][j] * dx * dx;
        }
    }

    fftshift(ref u1);
    FourierTransform2.FFT2(u1, FourierTransform.Direction.Forward);
    Complex[][] u2 = matrixMultiply(H, u1);
    FourierTransform2.FFT2(u2, FourierTransform.Direction.Backward);
    fftshift(ref u2);
    return u2;
}

Fraunhofer FFT Method

The Fraunhofer diffraction pattern using the FFT method is computed by $$U_2(x_2,y_2) = \left | \Im\left ( U_1\left ( x_1,y_1\right )\right ) h\left ( x_2,y_2\right ) \right |^2$$ where $h$ is the impulse response $$h(x,y) = \frac{e^{ikz}}{i \lambda z}e^{\frac{ik}{2z}(x^2+y^2)}$$ The Fraunhofer FFT method is implemented using the following function:

// Description  :   This function generates a Fraunhofer Diffraction pattern
//                  using the Fraunhofer FFT method
// Inputs       :   Uin - the pupil
//              :   L2 - Observation plane side length
//              :   lambda - wavelength
//              :   z - pupil-to-observation plane distance
// Outputs      :   matrixMultiply(Uin, c) - observation plane diffraction pattern

public static Complex[][] propFF(Complex[][] Uin, double L2, double lambda, double z)
{
    double k = 2 * Math.PI / lambda;
    double L1 = lambda * z * arraySize / L2;
    double dx2 = lambda * z / L1;
    double[] x2 = new double[arraySize];
    Complex[][] c = new Complex[arraySize][];

    for (int i = 0; i < arraySize; i++)
    {
        x2[i] = -L2 / 2.0 + i * dx2;
    }
    
    double[,] X2 = new double[arraySize, arraySize];
    double[,] Y2 = new double[arraySize, arraySize];
    meshgrid(x2, ref X2, ref Y2);

    for (int i = 0; i < arraySize; i++)
    {
        c[i] = new Complex[arraySize];
        for (int j = 0; j < arraySize; j++)
        {
            c[i][j] = Complex.One / Complex.ImaginaryOne * lambda * z * Complex.Exp(Complex.ImaginaryOne * k / (2 * z) * (X2[i, j] * X2[i, j] + Y2[i, j] * Y2[i, j]));
        }
    }

    fftshift(ref Uin);
    FourierTransform2.FFT2(Uin, FourierTransform.Direction.Forward);
    fftshift(ref Uin);
    return matrixMultiply(Uin, c);
}

Fraunhofer Exact Analytic Solution

The Fraunhofer diffraction pattern is computed using known exact analytic solutions. The solution used depends on the shape of the pupil.

This is the function that implements the exact analytic solution for a circular pupil $$I(\rho ) \propto \left (\frac{2J_1(\pi \rho D/\lambda z)}{\pi \rho D/\lambda z} \right )^{2}$$

// Description  :   Computes exact analytic solution for the
//                  Fraunhofer Diffraction Pattern of a circle
// Inputs       :   L - pupil plane side length
//              :   lambda - wavelength
//              :   z - pupil-to-observation plane distance
//              :   d - pupil diameter
// Outputs      :   u2 - the observation plane diffraction pattern

private Complex[][] exact(double L, double lambda, double z, double d)
{
    double dx = L / Convert.ToDouble(fftsize);
    double[] x = new double[fftsize];
    double[,] X = new double[fftsize, fftsize];
    double[,] Y = new double[fftsize, fftsize];
    Complex[][] u2 = new Complex[fftsize][];
    for (int i = 0; i < fftsize; i++)
    {
        x[i] = -L / 2.0 + i * dx;
    }

    OTWFunctions.meshgrid(x, ref X, ref Y);
    for (int i = 0; i < fftsize; i++)
    {
        u2[i] = new Complex[fftsize];
        for (int j = 0; j < fftsize; j++)
        {
            u2[i][j] = new Complex(Math.Pow(2 * Bessel.J(Math.PI * d * Math.Sqrt(Math.Pow(X[i, j], 2) + Math.Pow(Y[i, j], 2)) / (z * lambda)) / (Math.PI * d * Math.Sqrt(Math.Pow(X[i, j], 2) + Math.Pow(Y[i, j], 2)) / (z * lambda)), 2), 0.0);
        }
    }
    u2[fftsize/2][fftsize/2] = new Complex(1.0, 0.0);
    return u2;
}

This is the function that implements the exact analytic solution for a rectangular pupil $$I(x,y) \propto sinc^{2} \left ( \frac{\pi Wx}{\lambda z} \right ) sinc^{2} \left ( \frac{\pi Hy}{\lambda z} \right )$$

// Description  :   Computes exact analytic solution for the
//                  Fraunhofer Diffraction Pattern of a rectangular pupil
// Inputs       :   L - pupil plane side length
//              :   lambda - wavelength
//              :   z - pupil-to-observation plane distance
//              :   h - pupil height
//              :   w - pupil width
// Outputs      :   u2 - the observation plane diffraction pattern

private Complex[][] exact(Complex[][] u1, double L, double lambda, double z, double h, double w)
{
    double dx = L / Convert.ToDouble(fftsize);
    double[] x = new double[fftsize];
    double[,] X = new double[fftsize, fftsize];
    double[,] Y = new double[fftsize, fftsize];
    Complex[][] u2 = new Complex[fftsize][];

    for (int i = 0; i < fftsize; i++)
    {
        x[i] = -L / 2.0 + i * dx;
    }

    meshgrid(x, ref X, ref Y);

    for (int i = 0; i < fftsize; i++)
    {
        u2[i] = new Complex[fftsize];
        for (int j = 0; j < fftsize; j++)
        {
            u2[i][j] = new Complex(Math.Pow(sinc(Math.PI * w * X[i, j] / (lambda * z)), 2) * Math.Pow(sinc(Math.PI * h * Y[i, j] / (lambda * z)), 2), 0.0);
        }
    }
    return u2;
}

This is the function that implements the exact analytic solution for an obscured circular pupil $$I(\rho) \propto \frac{1}{\left ( 1-\varepsilon ^2 \right )^2}\left(\left ( \frac{2J_1(\pi D \rho / \lambda z)}{\pi D \rho / \lambda z} \right ) - \varepsilon^2 \left (\frac{2J_1(\varepsilon \pi D \rho / \lambda z)}{\varepsilon \pi D \rho / \lambda z} \right )\right )^2$$

// Description  :   Computes exact analytic solution for the
//                  Fraunhofer Diffraction Pattern of an obscured circular pupil
// Inputs       :   L - pupil plane side length
//              :   lambda - wavelength
//              :   z - pupil-to-observation plane distance
//              :   d - pupil outer diameter
//              :   d2 - pupil obscuration ratio
// Outputs      :   u2 - the observation plane diffraction pattern

private Complex[][] exact(Complex[][] u1, double L, double lambda, double z, double d, double d2)
{
    double dx = L / Convert.ToDouble(fftsize);
    double[] x = new double[fftsize];
    double[,] X = new double[fftsize, fftsize];
    double[,] Y = new double[fftsize, fftsize];
    Complex[][] u2 = new Complex[fftsize][];
    for (int i = 0; i < fftsize; i++)
    {
        x[i] = -L / 2.0 + i * dx;
    }

    meshgrid(x, ref X, ref Y);

    for (int i = 0; i < fftsize; i++)
    {
        u2[i] = new Complex[fftsize];
        for (int j = 0; j < fftsize; j++)
        {
            u2[i][j] = new Complex(Math.Pow((2 * Bessel.J(Math.PI * d * Math.Sqrt(Math.Pow(X[i, j], 2) + Math.Pow(Y[i, j], 2)) / (z * lambda)) / (Math.PI * d * Math.Sqrt(Math.Pow(X[i, j], 2) + Math.Pow(Y[i, j], 2)) / (z * lambda))) - Math.Pow(d2, 2) * (2 * Bessel.J(Math.PI * d2 * d * Math.Sqrt(Math.Pow(X[i, j], 2) + Math.Pow(Y[i, j], 2)) / (z * lambda)) / (Math.PI * d2 * d * Math.Sqrt(Math.Pow(X[i, j], 2) + Math.Pow(Y[i, j], 2)) / (z * lambda))), 2) / Math.Pow(1 - Math.Pow(d2, 2), 2), 0.0);
        }
    }
    u2[fftsize / 2][fftsize / 2] = new Complex(1.0, 0.0);
    return u2;
}

Read In Original Test Image

When performing a diffraction limited imaging simulation, the first thing we need to do is read in the original test image. This is accomplished using the following function call:

public static double[,] originalImage(string path)
{
    return complexMag(readInTestImg(path),1);
}

The innermost function call, which actually does most of the work, is:

// Description  :   Load graphic file from disk and store in a 2D numeric array
// Inputs       :   path - the location of the graphic file on the disk
// Outputs      :   grayValues - a 2D array containing the graphic data, converted
//                  to grayscale, normalized to 1, and assigned to the real part
//                  of each element

private static Complex[][] readInTestImg(string path)
{
    Bitmap bmp = new Bitmap(path);
    Rectangle rect = new Rectangle(0, 0, bmp.Width, bmp.Height);
    BitmapData bmpData = bmp.LockBits(rect, ImageLockMode.ReadWrite, bmp.PixelFormat);
    IntPtr ptr = bmpData.Scan0;
    int bytes = Math.Abs(bmpData.Stride) * bmp.Height;
    byte[] rgbValues = new byte[bytes];
    Marshal.Copy(ptr, rgbValues, 0, bytes);
    bmp.UnlockBits(bmpData);
    Complex[][] grayValues = new Complex[arraySize][];
    for (int i = 0; i < arraySize; i++)
    {
        grayValues[i] = new Complex[arraySize];
        for (int j = 0; j < arraySize; j++)
        {
            grayValues[i][j] = new Complex(Math.Sqrt(Convert.ToDouble(rgbValues[i * arraySize + j]) / 255), 0);
        }
    }
    return grayValues;
}

Coherent Imaging Simulation

A diffraction-limited coherent camera system with the specified aperture produces an image based on the equation $$U_2(x_2,y_2) = |\Im ^{-1}\left (H(x_2,y_2) \Im \left (I_g(x_2,y_2) \right ) \right)|^2$$ where $I_g(x_2,y_2)$ is the ideal image field and $H(x_2,y_2)$ is the transfer function $$H(x_2,y_2)=P(-\lambda z x_2,-\lambda z y_2)$$ where $P$ is the pupil function. This algorithm is implemented using the following function call:

// Description  :   Computs a diffraction-limited image generated by a simulated
//                  coherent camera system
// Inputs       :   lambda - wavelength
//              :   z - pupil-to-observation plane distance
//              :   L2 - observation plane side length
//              :   path - location of observation plane graphic on disk
//              :   H - the pupil
// Outputs      :   complexMag(Gi,1) - diffraction-limited image

public static double[,] coherentNoSpeckle(double lambda, double z, double L2, string path, Complex[][] H)
{
    Complex[][] grayValues = readInTestImg(path);
    fftshift(ref grayValues);
    FourierTransform2.FFT2(grayValues, FourierTransform.Direction.Forward);
    Complex[][] Gi = matrixMultiply(grayValues, H);
    FourierTransform2.FFT2(Gi, FourierTransform.Direction.Backward);
    fftshift(ref Gi);
    return complexMag(Gi,1);
}

Incoherent Imaging Simulation

A diffraction-limited incoherent camera system with the specified aperture produces an image based on the equation $$U_2(x_2,y_2) = \Im ^{-1}\left (\Im(|H(x_2,y_2)|^2) \Im \left (I_g(x_2,y_2) \right ) \right)$$ where $I_g(x_2,y_2)$ is the ideal image field and $H(x_2,y_2)$ is the transfer function $$H(x_2,y_2)=P(-\lambda z x_2,-\lambda z y_2)$$ normalized by the area of the transfer function, where $P$ is the pupil function. This algorithm is implemented using the following function call:

// Description  :   Computs a diffraction-limited image generated by a simulated
//                  coherent camera system
// Inputs       :   lambda - wavelength
//              :   z - pupil-to-observation plane distance
//              :   L2 - observation plane side length
//              :   path - location of observation plane graphic on disk
//              :   H - the pupil
// Outputs      :   complexMag(Gi,1) - diffraction-limited image

public static double[,] incoherent(double lambda, double z, double L2, string path, Complex[][] H)
{
    FourierTransform2.FFT2(H, FourierTransform.Direction.Forward);
    Complex[][] OTF = double2Complex(complexMag(H, 2));
    FourierTransform2.FFT2(OTF, FourierTransform.Direction.Backward);
    Complex normalization = OTF[0][0];
    for (int i = 0; i < arraySize; i++)
    {
        for (int j = 0; j < arraySize; j++)
        {
            OTF[i][j] = OTF[i][j] / normalization;
        }
    }

    Complex[][] grayValues = readInTestImg(path);
    fftshift(ref grayValues);
    FourierTransform2.FFT2(grayValues, FourierTransform.Direction.Forward);
    Complex[][] Gi = matrixMultiply(grayValues, OTF);
    FourierTransform2.FFT2(Gi, FourierTransform.Direction.Backward);
    fftshift(ref Gi);
    return realPart(Gi);
}

Seidel Wavefront

This function computes the Seidel power series expansion $$W_{ijk} = H^{i}\rho ^{j}\cos ^{k}\theta$$

// Description  :   Computes the Seidel power series expansion
// Inputs       :   beta - angle of observation plane PSF location vector
//              :   u0r - length of observation plane PSF location vector
//              :   X, Y - observation plane location
//              :   w*** - Seidel coefficients
// Outputs      :   Seidel power series expansion total value

public static double seidel_5(double beta, double u0r, double X, double Y, double w020, double w111, double w040, double w131, double w222, double w220, double w311, double w060, double w240, double w420, double w151, double w331, double w511, double w242, double w422, double w333)
{
    double Xr = X * Math.Cos(beta) + Y * Math.Sin(beta);
    double Yr = -X * Math.Sin(beta) + Y * Math.Cos(beta);
    double rho2 = Math.Pow(Xr, 2) + Math.Pow(Yr, 2);
    return 
        (w020 * rho2) + 
        (w111 * u0r * Xr) + 
        (w040 * Math.Pow(rho2, 2)) + 
        (w131 * u0r * rho2 * Xr) + 
        (w222 * Math.Pow(u0r, 2) * Math.Pow(Xr, 2)) + 
        (w220 * Math.Pow(u0r, 2) * rho2) + 
        (w311 * Math.Pow(u0r, 3) * Xr) + 
        (w060 * Math.Pow(rho2, 3)) + 
        (w240 * Math.Pow(u0r, 2) * Math.Pow(rho2, 2)) + 
        (w420 * Math.Pow(u0r, 4) * rho2) + 
        (w151 * u0r * Math.Pow(rho2, 2) * Xr) + 
        (w331 * Math.Pow(u0r, 3) * rho2 * Xr) + 
        (w511 * Math.Pow(u0r, 5) * Xr) + 
        (w242 * Math.Pow(u0r, 2) * rho2 * Math.Pow(Xr, 2)) + 
        (w422 * Math.Pow(u0r, 4) * Math.Pow(Xr, 2)) + 
        (w333 * Math.Pow(u0r, 3) * Math.Pow(Xr, 3));
}

Aberrated Pupil Function

This function computes the aberrated pupil function $$H(f_U, f_V) = P(-\lambda zf_U, -\lambda z f_V)$$ where $$P(x_p, y_p) = circ\left ( \frac{\sqrt{x_p^2 + y_p^2}}{r_{xp}} \right )e^{-ikW(x_p,y_p)}$$

// Description  :   Computes the aberration pupil function in the Seidel domain
// Inputs       :   l - observation plane side length
//              :   lambda - wavelength
//              :   z - pupil-to-observation plane distance
//              :   d - pupil diameter
//              :   u0, v0 - normalized observation plane PSF centroid location
//              :   w*** - Seidel coefficients
// Outputs      :   H - aberrated pupil function

public static Complex[][] HSeidel(double l, double lambda, double z, double d, double u0, double v0, double w020, double w111,
        double w040, double w220, double w131, double w311, double w222, double w060, double w240, double w420, double w151,
        double w331, double w511, double w242, double w422, double w333)
{
    double du = l / Convert.ToDouble(arraySize);
    double k = 2.0 * Math.PI / lambda;
    double wxp = d / 2.0;
    double lz = lambda * z;

    double[] fu = new double[arraySize];
    for (int i = 0; i < arraySize; i++)
    {
        fu[i] = (-1.0 / (2.0 * du)) + (i / l);
    }

    double[,] Fu = new double[arraySize, arraySize];
    double[,] Fv = new double[arraySize, arraySize];

    meshgrid(fu, ref Fu, ref Fv);

    double beta = Math.Atan2(v0, u0);
    double u0r = Math.Sqrt(Math.Pow(u0, 2) + Math.Pow(v0, 2));
    double[,] W = new double[arraySize, arraySize];
    Complex[][] H = new Complex[arraySize][];

    for (int i = 0; i < arraySize; i++)
    {
        H[i] = new Complex[arraySize];
        for (int j = 0; j < arraySize; j++)
        {
            W[i, j] = seidel_5(beta, u0r, -lz * Fu[i, j] / wxp, -lz * Fv[i, j] / wxp, w020, w111, w040, w131, w222, w220, w311, w060, w240, w420, w151, w331, w511, w242, w422, w333);
            if ((Math.Sqrt(Math.Pow(Fu[i, j], 2) + Math.Pow(Fv[i, j], 2)) * lz / wxp) < 1.0)
            {
                H[i][j] = Math.Sqrt(Math.Pow(Fu[i, j], 2) + Math.Pow(Fv[i, j], 2)) * lz / wxp * Complex.Exp(-Complex.ImaginaryOne * k * W[i, j]);
            }
            else
            {
                H[i][j] = Complex.Zero;
            }
        }
    }
    return H;
}

9-point PSF Simulation

// Description  :   This function generates the 9 aberrated pupil functions in the 9-point PSF simulation
// Inputs       :   l - observation plane side length
//              :   lambda - wavelength
//              :   d - pupil diameter
//              :   z - pupil-to-observation plane distance
//              :   w*** - Seidel coefficients
// Outputs      :   I - observation plane intensity

public static double[,] image_simulation_49pt(double l, double lambda, double d, double z, double w020, double w111, double w040, double w131, double w222, double w220, double w311, double w060, double w240, double w420, double w151, double w331, double w511, double w242, double w422, double w333)
{
    double du = l / Convert.ToDouble(arraySize);
    double k = 2.0 * Math.PI / lambda;
    double wxp = d / 2.0;
    double lz = lambda * z;

    double[] fu = new double[arraySize];
    for (int i = 0; i < arraySize; i++)
    {
        fu[i] = (-1.0 / (2.0 * du)) + (i / l);
    }

    double[,] Fu = new double[arraySize, arraySize];
    double[,] Fv = new double[arraySize, arraySize];

    meshgrid(fu, ref Fv, ref Fu);
    
    double[,] I = new double[arraySize, arraySize];
    for (double u0 = -0.7; u0 <= 0.7; u0 = u0 + 0.7)
    {
        for (double v0 = -0.7; v0 <= 0.7; v0 = v0 + 0.7)
        {
            double beta = Math.Atan2(v0, u0);
            double u0r = Math.Sqrt(Math.Pow(u0, 2) + Math.Pow(v0, 2));
            double[,] W = new double[arraySize, arraySize];
            Complex[][] H = new Complex[arraySize][];

            for (int i = 0; i < arraySize; i++)
            {
                H[i] = new Complex[arraySize];
                for (int j = 0; j < arraySize; j++)
                {
                    W[i, j] = seidel_5(beta, u0r, -lz * Fu[i, j] / wxp, -lz * Fv[i, j] / wxp, -w020, -w111, -w040, -w131, -w222, -w220, -w311, -w060, -w240, -w420, -w151, -w331, -w511, -w242, -w422, -w333);
                    if ((Math.Sqrt(Math.Pow(Fu[i, j], 2) + Math.Pow(Fv[i, j], 2)) * lz / wxp) < 1.0)
                    {
                        H[i][j] = Math.Sqrt(Math.Pow(Fu[i, j], 2) + Math.Pow(Fv[i, j], 2)) * lz / wxp * Complex.Exp(-Complex.ImaginaryOne * k * W[i, j]);
                    }
                    else
                    {
                        H[i][j] = Complex.Zero;
                    }
                }
            }

            double[,] h2 = PSF(H);
            double[,] h2_new = circshift(ref h2, Convert.ToInt32(-v0 * arraySize / 2), Convert.ToInt32(-u0 * arraySize / 2));
            for (int i = 0; i < arraySize; i++)
            {
                for (int j = 0; j < arraySize; j++)
                {
                    I[i, j] = I[i, j] + h2_new[i, j];
                }
            }
        }
    }
    return I;
}

Zernike Wavefront

This function computes the Zernike series expansion $$\begin{eqnarray} W(\rho, \theta) = \overline{\Delta W} + \sum_{n=1}^{\infty}\left [ A_n \sum_{s=0}^{n}(-1)^s \frac{(2n -s)! \rho^{2(n-s)}}{s!(n-s)!(n-s)!} + \\ \sum_{m=1}^{n} \sum_{s=0}^{n-m}(-1)^s \frac{(2n -m -s)! \rho^{2(n-m-s)}}{s!(n-s)!(n-m-s)!} \rho^m (B_{nm} \cos m\theta' + C_{nm} \sin m\theta')\right ] \end{eqnarray}$$ where $\rho = \sqrt{x_p^2+y_p^2}$, $\rho \cos \theta = x_p$, $\overline{\Delta W}$ is the mean wavefront ODP and $A_n$, $B_{nm}$, and $C_{nm}$ are the Zernike coefficients.

// Description  :   Computes the Zernike series expansion
// Inputs       :   x, y - observation plane locations
//              :   z** - Zernike coefficients
// Outputs      :   Zernike series expansion total value

public static double zernike(double x, double y, double z1_1, double z11, double z2_2, double z20, double z22, double z3_3, double z3_1, double z31, double z33, double z4_4, double z4_2, double z40, double z42, double z44)
{
    return 
        (z1_1 * 2.0 * x) + 
        (z11 * 2.0 * y) + 
        (z2_2 * Math.Sqrt(6.0) * 2.0 * x * y) + 
        (z20 * Math.Sqrt(3.0) * (2 * Math.Pow(x, 2) + 2 * Math.Pow(y, 2) + 1)) + 
        (z22 * Math.Sqrt(6.0) * (-Math.Pow(x, 2) + Math.Pow(y, 2))) + 
        (z3_3 * Math.Sqrt(6.0) * (-Math.Pow(x, 3) + 3 * x * Math.Pow(y, 2))) + 
        (z3_1 * Math.Sqrt(6.0) * (3 * Math.Pow(x, 3) + 3 * x * Math.Pow(y, 2) - 2 * x)) + 
        (z31 * Math.Sqrt(6.0) * (3 * Math.Pow(x, 2) * y + 3 * Math.Pow(y, 3) - 2 * y)) + 
        (z33 * Math.Sqrt(6.0) * (-3 * Math.Pow(x, 2) * y + Math.Pow(y, 3))) + 
        (z4_4 * Math.Sqrt(10.0) * (-4 * Math.Pow(x, 3) * y + 4 * x * Math.Pow(y, 3))) + 
        (z4_2 * Math.Sqrt(10) * (8 * Math.Pow(x, 3) * y + 8 * x * Math.Pow(y, 3) + 6 * x * y)) + 
        (z40 * Math.Sqrt(5.0) * (6 * Math.Pow(x, 4) + 12 * Math.Pow(x, 2) * Math.Pow(y, 2) + 6 * Math.Pow(y, 4) - 6 * Math.Pow(x, 2) - 6 * Math.Pow(y, 2) + 1)) + 
        (z42 * Math.Sqrt(10.0) * (-4 * Math.Pow(x, 4) + 4 * Math.Pow(y, 4) + 3 * Math.Pow(x, 2) - 3 * Math.Pow(y, 2))) + 
        (z44 * Math.Sqrt(10.0) * (Math.Pow(x, 4) - 6 * Math.Pow(x, 2) * Math.Pow(y, 2) + Math.Pow(y, 4)));
}

Aberrated Pupil Function

This function computes the aberrated pupil function $$H(f_U, f_V) = P(-\lambda zf_U, -\lambda z f_V)$$ where $$P(x_p, y_p) = circ\left ( \frac{\sqrt{x_p^2 + y_p^2}}{r_{xp}} \right )e^{-ikW(x_p,y_p)}$$

// Description  :   Computes the aberrated pupil function in the Zernike domain
// Inputs       :   l - observation plane side length
//              :   lambda - wavelength
//              :   z - pupil-to-observation plane distance
//              :   d - pupil diameter
//              :   z** - Zernike coefficients
// Outputs      :   H - aberrated pupil function

public static Complex[][] HZernike(double l, double lambda, double z, double d, double z1_1, double z11, double z2_2, double z20,
        double z22, double z3_3, double z3_1, double z31, double z33, double z4_4, double z4_2, double z40, double z42, double z44)
{
    double du = l / Convert.ToDouble(arraySize);

    double[] u = new double[arraySize];
    for (int i = 0; i < arraySize; i++)
    {
        u[i] = (-l / 2.0) + (i * du);
    }

    double k = 2.0 * Math.PI / lambda;
    double wxp = d / 2.0;
    double lz = lambda * z;

    double[] fu = new double[arraySize];
    for (int i = 0; i < arraySize; i++)
    {
        fu[i] = (-1.0 / (2.0 * du)) + (i / l);
    }

    double[,] Fu = new double[arraySize, arraySize];
    double[,] Fv = new double[arraySize, arraySize];

    meshgrid(fu, ref Fu, ref Fv);

    double[,] W = new double[arraySize, arraySize];
    Complex[][] H = new Complex[arraySize][];

    for (int i = 0; i < arraySize; i++)
    {
        H[i] = new Complex[arraySize];
        for (int j = 0; j < arraySize; j++)
        {
            W[i, j] = zernike(-lz * Fu[i, j] / wxp, -lz * Fv[i, j] / wxp, z1_1, z11, z2_2, z20, z22, z3_3, z3_1, z31, z33, z4_4, z4_2, z40, z42, z44);
            if ((Math.Sqrt(Math.Pow(Fu[i, j], 2) + Math.Pow(Fv[i, j], 2)) * lz / wxp) < 1.0)
            {
                H[i][j] = Math.Sqrt(Math.Pow(Fu[i, j], 2) + Math.Pow(Fv[i, j], 2)) * lz / wxp * Complex.Exp(-Complex.ImaginaryOne * k * W[i, j]);
            }
            else
            {
                H[i][j] = Complex.Zero;
            }
        }
    }

    return H;
}