SimulateSpectrum.cpp 19.5 KB
``````#include <math.h>
#include <complex>
#include <iostream>
#include <fstream>
#include "globals.h"
//#include "cufft.h"
using namespace std;

#define pi 3.14159

typedef complex<float> scComplex;

int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);
int bessjyv(double v,double x,double &vm,double *jv,double *yv,
double *djv,double *dyv);

complex<double> Jl_neg(complex<double> x)
{
//this function computes the bessel function of the first kind Jl(x) for l = -0.5
return ( sqrt(2.0/pi) * cos(x) )/sqrt(x);
}

double Jl_neg(double x)
{
//this function computes the bessel function of the first kind Jl(x) for l = -0.5
return ( sqrt(2.0/pi) * cos(x) )/sqrt(x);
}

double Yl_neg(double x)
{
//this function computes the bessel function of the second kind Yl(x) for l = -0.5;
return ( sqrt(2.0/pi) * sin(x) )/sqrt(x);
}

void computeB(complex<float>* B, float radius, complex<double> refIndex, float lambda, int Nl)
{
double k = (2*pi)/lambda;
int b = 2;

//allocate space for the real bessel functions
double* jv = (double*)malloc(sizeof(double)*(Nl+b));
double* yv = (double*)malloc(sizeof(double)*(Nl+b));
double* jvp = (double*)malloc(sizeof(double)*(Nl+b));
double* yvp = (double*)malloc(sizeof(double)*(Nl+b));

//allocate space for the complex bessel functions
complex<double>* cjv = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
complex<double>* cyv = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
complex<double>* cjvp = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
complex<double>* cyvp = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));

double kr = k*radius;
complex<double> knr = k*refIndex*(double)radius;
complex<double> n = refIndex;

//compute the bessel functions for k*r
double vm;// = Nl - 1;
bessjyv((Nl)+0.5, kr, vm, jv, yv, jvp, yvp);
//cout<<"Nl: "<<Nl<<"  vm: "<<vm<<endl;
//printf("Nl: %f, vm: %f\n", (float)Nl, (float)vm);

//compute the bessel functions for k*n*r
cbessjyva((Nl)+0.5, knr, vm, cjv, cyv, cjvp, cyvp);

//scale factor for spherical bessel functions
double scale_kr = sqrt(pi/(2.0*kr));
complex<double> scale_knr = sqrt(pi/(2.0*knr));

complex<double> numer, denom;
double j_kr;
double y_kr;
complex<double> j_knr;
complex<double> j_d_knr;
double j_d_kr;
complex<double> h_kr;
complex<double> h_d_kr;
complex<double> h_neg;
complex<double> h_pos;

//cout<<"B coefficients:"<<endl;
for(int l=0; l<Nl; l++)
{
//compute the spherical bessel functions
j_kr = jv[l] * scale_kr;
y_kr = yv[l] * scale_kr;
j_knr = cjv[l] * scale_knr;

//compute the Hankel function
h_kr = complex<double>(j_kr, y_kr);

//compute the derivatives
if(l == 0)
{
//spherical bessel functions for l=0
j_d_kr = scale_kr * (Jl_neg(kr) - (jv[l] + kr*jv[l+1])/kr )/2.0;
j_d_knr = scale_knr * ( Jl_neg(knr) - (cjv[l] + knr*cjv[l+1])/knr )/2.0;
h_neg = complex<double>(scale_kr*Jl_neg(kr), scale_kr*Yl_neg(kr));
h_pos = complex<double>(scale_kr*jv[l+1], scale_kr*yv[l+1]);
h_d_kr = (h_neg - (h_kr + kr*h_pos)/kr)/2.0;
}
else
{
//spherical bessel functions
j_d_kr = scale_kr * (jv[l-1] - (jv[l] + kr*jv[l+1])/kr )/2.0;
j_d_knr = scale_knr * ( cjv[l-1] - (cjv[l] + knr*cjv[l+1])/knr )/2.0;
h_neg = complex<double>(scale_kr*jv[l-1], scale_kr*yv[l-1]);
h_pos = complex<double>(scale_kr*jv[l+1], scale_kr*yv[l+1]);
h_d_kr = (h_neg - (h_kr + kr*h_pos)/kr)/2.0;
}

numer = j_kr*j_d_knr*n - j_knr*j_d_kr;
denom = j_knr*h_d_kr - h_kr*j_d_knr*n;
complex<double> temp =  numer/denom;

B[l] = scComplex(temp.real(), temp.imag());
//cout<<B[l]<<endl;
}

free(jv);
free(yv);
free(jvp);
free(yvp);
free(cjv);
free(cyv);
free(cjvp);
free(cyvp);
}

void Legendre(float* P, float x, int Nl)
{
//computes the legendre polynomials from orders 0 to Nl-1
P[0] = 1;
if(Nl == 1) return;
P[1] = x;
for(int l = 2; l < Nl; l++)
{
P[l] = ((2*l - 1)*x*P[l-1] - (l - 1)*P[l-2])/l;
}

}

complex<float> integrateUi(float cAngleI, float cAngleO, float oAngleI, float oAngleO, float M = 2*pi)
{
/*This function integrates the incident field of magnitude M in the far zone
in order to evaluate the field at the central pixel of a detector.
cNAi = condenser inner angle
cNAo = condenser outer angle
oNAi = objective inner angle
oNAo = objective outer angle
M = field magnitude*/

float alphaIn = max(cAngleI, oAngleI);
float alphaOut = min(cAngleO,oAngleO);

complex<float> Ui;
if(alphaIn > alphaOut)
Ui = complex<float>(0.0, 0.0);
else
Ui = complex<float>(M * 2 * pi * (cos(alphaIn) - cos(alphaOut)), 0.0f);

return Ui;

}

void computeCondenserAlpha(float* alpha, int Nl, float cAngleI, float cAngleO)
{
/*This function computes the condenser integral in order to build the field of incident light
alpha = list of Nl floating point values representing the condenser alpha as a function of l
Nl = number of orders in the incident field
cAngleI, cAngleO = inner and outer condenser angles (inner and outer NA)*/

//compute the Legendre polynomials for the condenser aperature
float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1));
Legendre(PcNAo, cos(cAngleO), Nl+1);
float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1));
Legendre(PcNAi, cos(cAngleI), Nl+1);

for(int l=0; l<Nl; l++)
{
//integration term
if(l == 0)
alpha[l] = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
else
alpha[l] = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];

alpha[l] *= 2 * pi;
}

}

complex<float> integrateUs(float r, float lambda, complex<float> eta,
float cAngleI, float cAngleO, float oAngleI, float oAngleO, float M = 2*pi)
{
/*This function integrates the incident field of magnitude M in the far zone
in order to evaluate the field at the central pixel of a detector.
r = sphere radius
lambda = wavelength
eta = index of refraction
cNAi = condenser inner NA
cNAo = condenser outer NA
oNAi = objective inner NA
oNAo = objective outer NA
M = field magnitude*/

//compute the required number of orders
float k = 2*pi/lambda;
int Nl = ceil( k + 4 * exp(log(k*r)/3) + 3 );

//compute the material coefficients B
complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);
//compute the Legendre polynomials for the condenser and objective aperatures
float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1));
Legendre(PcNAo, cos(cAngleO), Nl+1);
float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1));
Legendre(PcNAi, cos(cAngleI), Nl+1);

float* PoNAo = (float*)malloc(sizeof(float)*(Nl+1));
Legendre(PoNAo, cos(oAngleO), Nl+1);
float* PoNAi = (float*)malloc(sizeof(float)*(Nl+1));
Legendre(PoNAi, cos(oAngleI), Nl+1);

//store the index of refraction;
complex<double> IR(eta.real(), eta.imag());

//compute the scattering coefficients
computeB(B, r, IR, lambda, Nl);

//aperature terms for the condenser (alpha) and objective (beta)
float alpha;
float beta;
float c;
complex<float> Us(0.0, 0.0);

for(int l=0; l<Nl; l++)
{
//integration term
if(l == 0)
{
alpha = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
beta = -PoNAo[l+1] + PoNAo[0] + PoNAi[l+1] - PoNAi[0];
}
else
{
alpha = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];
beta = -PoNAo[l+1] + PoNAo[l-1] + PoNAi[l+1] - PoNAi[l-1];
}
c = (2*pi)/(2.0 * l + 1.0);
Us += c * alpha * beta * B[l] * M;

}
free(PcNAo);
free(PcNAi);
free(PoNAo);
free(PoNAi);
free(B);

return Us;

}

void pointSpectrum()
{
PD.StartTimer(SIMULATE_SPECTRUM);
//clear the previous spectrum
SimSpectrum.clear();

float dNu = 2.0f;
float lambda;

//compute the angles based on NA
float cAngleI = asin(cNAi);
float cAngleO = asin(cNAo);
float oAngleI = asin(oNAi);
float oAngleO = asin(oNAo);

//implement a reflection-mode system if necessary
if(opticsMode == ReflectionOpticsType){

//set the condenser to match the objective
cAngleI = oAngleI;
cAngleO = oAngleO;

//invert the objective
oAngleO = pi - cAngleI;
oAngleI = pi - cAngleO;
}

//integrate the incident field at the detector position
complex<float> Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
float I0 = Ui.real() * Ui.real() + Ui.imag() * Ui.imag();
I0 *= scaleI0;

float I;
SpecPair temp;
float nu;
complex<float> eta;
complex<float> Us, U;

float vecLen = 0.0;
for(int i=0; i<EtaK.size(); i++)
{
nu = EtaK[i].nu;
lambda = 10000.0f/nu;
if(applyMaterial)
eta = complex<float>(EtaN[i].A, EtaK[i].A);
else
eta = complex<float>(baseIR, 0.0);

//integrate the scattered field at the detector position
Us = integrateUs(radius, lambda, eta, cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
U = Us + Ui;
float I = U.real() * U.real() + U.imag() * U.imag();

temp.nu = nu;

//set the spectrum value based on the current display type
if(dispSimType == AbsorbanceSpecType)
temp.A = -log10(I/I0);
else
temp.A = I;

if(dispNormalize)
vecLen += temp.A * temp.A;

SimSpectrum.push_back(temp);
}
vecLen = sqrt(vecLen);

if(dispNormalize)
for(int i=0; i<SimSpectrum.size(); i++)
SimSpectrum[i].A = (SimSpectrum[i].A / vecLen) * dispNormFactor;

PD.EndTimer(SIMULATE_SPECTRUM);
}

/*
complex<float> sampleUs(complex<float>* B, float* Alpha, int Nl, float r,
float cAngleI, float cAngleO, float theta, float M = 2*pi)
{
/*This function takes a point sample of the scattered field in the far zone
in order to evaluate the field at the central pixel of a detector.
r = sphere radius
lambda = wavelength
eta = index of refraction
cNAi = condenser inner NA
cNAo = condenser outer NA
theta = angle of the sample
M = field magnitude*/

/*
//compute the material coefficients B
//complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);

//compute the Legendre polynomials for theta (at the objective)
float* Ptheta = (float*)malloc(sizeof(float)*(Nl+1));
Legendre(Ptheta, cos(theta), Nl+1);

//complex<double> IR(eta.real(), eta.imag());

//aperature terms for the condenser (alpha) and objective (beta)
float beta;
float c;
complex<float> Us(0.0, 0.0);

for(int l=0; l<Nl; l++)
{

complex<float> numer(0.0, -(l*pi)/2.0);
Us += B[l] * exp(numer) * Ptheta[l] * Alpha[l] * pow(complex<float>(0.0, 1.0), l);

}
//printf("Ptheta: %f\n", Ptheta[Nl-1]);

return Us;

}*/

/*
void detectorSpectrum(int numSamples)
{
//integrate across the objective aperature and calculate the resulting intensity on a detector
PD.StartTimer(SIMULATE_SPECTRUM);
//clear the previous spectrum
SimSpectrum.clear();

float dNu = 2.0f;
float lambda;

//compute the angles based on NA
float cAngleI = asin(cNAi);
float cAngleO = asin(cNAo);
float oAngleI = asin(oNAi);
float oAngleO = asin(oNAo);

//implement a reflection-mode system if necessary
if(opticsMode == ReflectionOpticsType){

//set the condenser to match the objective
cAngleI = oAngleI;
cAngleO = oAngleO;

//invert the objective
oAngleO = pi - cAngleI;
oAngleI = pi - cAngleO;
}

//compute Nl (maximum order of the spectrum)
//****************************************************************************
float maxNu = EtaK.back().nu;
float maxLambda = 10000.0f/maxNu;
float k = 2*pi/maxLambda;
int Nl = ceil( k + 4 * exp(log(k*radius)/3) + 3 );
int nLambda = EtaK.size();

//compute alpha (condenser integral)
//****************************************************************************
//compute the Legendre polynomials for the condenser aperature
float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1));
Legendre(PcNAo, cos(cAngleO), Nl+1);
float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1));
Legendre(PcNAi, cos(cAngleI), Nl+1);

//allocate space for the alpha array
float* alpha = (float*)malloc(sizeof(float)*(Nl + 1));
computeCondenserAlpha(alpha, Nl, cAngleI, cAngleO);

for(int l=0; l<Nl; l++)
{
//integration term
if(l == 0)
alpha[l] = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
else
alpha[l] = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];

alpha[l] *= 2 * pi;
}

//compute the information based on wavelength and sphere

//evaluate the incident field intensity
float I0 = 0.0;
float theta;
float dTheta = (oAngleO - oAngleI)/numSamples;
complex<float> Ui;

Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
I0 = Ui.real()*2*pi;

float I;
SpecPair temp;
float nu;
complex<float> eta;
complex<float> Us, U;

float vecLen = 0.0;
for(int i=0; i<EtaK.size(); i++)
{
//compute information based on wavelength and material
nu = EtaK[i].nu;
lambda = 10000.0f/nu;
if(applyMaterial)
eta = complex<float>(EtaN[i].A, EtaK[i].A);
else
eta = complex<float>(baseIR, 0.0);

//allocate memory for the scattering coefficients
complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);

complex<double> IR(eta.real(), eta.imag());
computeB(B, radius, IR, lambda, Nl);

//integrate the scattered field at the detector position
I = 0.0;

for(int iTheta = 0; iTheta < numSamples; iTheta++)
{
theta = oAngleI + iTheta * dTheta;
Us = sampleUs(B, alpha, Nl, radius, cAngleI, cAngleO, theta, 2*pi);

//calculate the intensity and add
if(theta >= cAngleI && theta <= cAngleO)
U = Us + 2*(float)pi;
else
U = Us;

I += (U.real() * U.real() + U.imag() * U.imag()) * sin(theta) * 2 * pi * dTheta;

}

temp.nu = nu;

if(i == 0)
printf("I: %f        I0: %f\n", I, I0);

//set the spectrum value based on the current display type
if(dispSimType == AbsorbanceSpecType)
temp.A = -log10(I/I0);
else
temp.A = I;

if(dispNormalize)
vecLen += temp.A * temp.A;
//temp.A = Us.real();
SimSpectrum.push_back(temp);

free(B);
}
vecLen = sqrt(vecLen);

if(dispNormalize)
for(int i=0; i<SimSpectrum.size(); i++)
SimSpectrum[i].A = (SimSpectrum[i].A / vecLen) * dispNormFactor;

free(alpha);

PD.EndTimer(SIMULATE_SPECTRUM);
}*/

void updateSpectrum(float* I, float I0, int n)
{
SimSpectrum.clear();
SpecPair temp;

//update the displayed spectrum based on the computed intensity I
for(int i=0; i<n; i++)
{
temp.nu = EtaK[i].nu;

//set the spectrum value based on the current display type
if(dispSimType == AbsorbanceSpecType)
temp.A = -log10(I[i]/I0);
else
temp.A = I[i];

SimSpectrum.push_back(temp);
}
}

void computeCassegrainAngles(float& cAngleI, float& cAngleO, float& oAngleI, float& oAngleO)
{
//compute the angles based on NA
cAngleI = asin(cNAi);
cAngleO = asin(cNAo);
oAngleI = asin(oNAi);
oAngleO = asin(oNAo);

//implement a reflection-mode system if necessary
if(opticsMode == ReflectionOpticsType){

//set the condenser to match the objective
cAngleI = oAngleI;
cAngleO = oAngleO;

//invert the objective
oAngleO = pi - cAngleI;
oAngleI = pi - cAngleO;
}

}

int computeNl()
{
float maxNu = EtaK.back().nu;
float maxLambda = 10000.0f/maxNu;
float k = 2*pi/maxLambda;
int Nl = ceil( k + 4 * exp(log(k*radius)/3) + 3 );

return Nl;
}

void computeBArray(complex<float>* B, int Nl, int nLambda)
{
float nu;
complex<float> eta;
float* Lambda = (float*)malloc(sizeof(float) * nLambda);

//for each wavenumber nu
for(int i=0; i<EtaK.size(); i++)
{
//compute information based on wavelength and material
nu = EtaK[i].nu;
Lambda[i] = 10000.0f/nu;
if(applyMaterial)
eta = complex<float>(EtaN[i].A, EtaK[i].A);
else
eta = complex<float>(baseIR, 0.0);

//allocate memory for the scattering coefficients
//complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);

complex<double> IR(eta.real(), eta.imag());
computeB(&B[i * Nl], radius, IR, Lambda[i], Nl);
}
}

void computeOpticalParameters(float& cAngleI, float& cAngleO, float& oAngleI, float& oAngleO, float& I0, float* alpha, int Nl)
{
computeCassegrainAngles(cAngleI, cAngleO, oAngleI, oAngleO);

//evaluate the incident field intensity
I0 = 0.0;
float theta;
complex<float> Ui;

Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
I0 = Ui.real()*2*pi;

//compute alpha (condenser integral)
computeCondenserAlpha(alpha, Nl, cAngleI, cAngleO);
}

void gpuDetectorSpectrum(int numSamples)
{
//integrate across the objective aperature and calculate the resulting intensity on a detector
PD.StartTimer(SIMULATE_SPECTRUM);
//clear the previous spectrum
SimSpectrum.clear();

//compute Nl (maximum order of the spectrum)
int Nl = computeNl();

float* alpha = (float*)malloc(sizeof(float)*(Nl + 1));
float cAngleI, cAngleO, oAngleI, oAngleO, I0;
computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl);

//allocate space for a list of wavelengths
int nLambda = EtaK.size();

//allocate space for the 2D array (Nl x nu) of scattering coefficients
complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>) * Nl * nLambda);
computeBArray(B, Nl, nLambda);

//allocate temporary space for the spectrum
float* I = (float*)malloc(sizeof(float) * EtaK.size());

//compute the spectrum on the GPU
PD.StartTimer(SIMULATE_GPU);
cudaComputeSpectrum(I, (float*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, numSamples);
PD.EndTimer(SIMULATE_GPU);

PD.EndTimer(SIMULATE_SPECTRUM);

}

void SimulateSpectrum()
{
if(pointDetector)
pointSpectrum();
else
gpuDetectorSpectrum(objectiveSamples);
//detectorSpectrum(objectiveSamples);
}

float absorbanceDistortion(){

//compute the mean of the spectrum
float sumSim = 0.0;
for(int i=0; i<SimSpectrum.size(); i++)
{
sumSim += SimSpectrum[i].A;
}
float meanSim = sumSim/SimSpectrum.size();

//compute the distortion (MSE from the mean)
float sumSE = 0.0;
for(int i=0; i<SimSpectrum.size(); i++)
{
sumSE += pow(SimSpectrum[i].A - meanSim, 2);
}
float MSE = sumSE / SimSpectrum.size();

return MSE;
}

float intensityDistortion(){

//compute the magnitude of the spectrum
float sumSim = 0.0;
for(int i=0; i<SimSpectrum.size(); i++)
{
sumSim += SimSpectrum[i].A * SimSpectrum[i].A;
}
float magSim = sqrt(sumSim);

//compute the distortion (MSE from the mean)
float sumSE = 0.0;
for(int i=0; i<SimSpectrum.size(); i++)
{
sumSE += (SimSpectrum[i].A/magSim) * (1.0/SimSpectrum.size());
}
float MSE = sumSE;

return MSE;
}

void MinimizeDistortion(){
ofstream outFile("distortion.txt");

//set the parameters for the distortion simulation
float step = 0.001;

oNAi = 0.2;
oNAo = 0.5;

//compute the optical parameters
//compute Nl (maximum order of the spectrum)
int Nl = computeNl();

float* alpha = (float*)malloc(sizeof(float)*(Nl + 1));
float cAngleI, cAngleO, oAngleI, oAngleO, I0;

//allocate space for a list of wavelengths
int nLambda = EtaK.size();

//allocate temporary space for the spectrum
float* I = (float*)malloc(sizeof(float) * EtaK.size());

//calculate the material parameters
//allocate space for the 2D array (Nl x nu) of scattering coefficients
complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>) * Nl * nLambda);
computeBArray(B, Nl, nLambda);

float D;
float e = 0.001;
for(float i=0.0; i<=oNAo-step; i+=step)
{

for(float o=oNAi+step; o<=1.0; o+=step)
{

//set the current optical parameters
cNAi = i;
cNAo = o;

//compute the optical parameters
computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl);

//simulate the spectrum
cudaComputeSpectrum(I, (float*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, objectiveSamples);

if(dispSimType == AbsorbanceSpecType)
{
if(i + e >= o || i + e >= oNAo || oNAi + e >= o || oNAi + e >= oNAo)
D = 0.0;
else
D = absorbanceDistortion();
}
else
{
if(i >= o || oNAi >= oNAo)
D=0;
else
D = intensityDistortion();
}
outFile<<D<<"       ";
}
outFile<<endl;
cout<<i<<endl;
}
outFile.close();
}``````