main.cpp

#include <fstream>
#include <iostream>
using namespace std;
#include "interactivemie.h"
#include <QtGui/QApplication>
#include <QGraphicsScene>
#include <QGraphicsView>
#include <QGraphicsPixmapItem>
#include <QLayout>
#include "qtSpectrumDisplay.h"
#include "globals.h"
#include "rtsGUIConsole.h"
#include "PerformanceData.h"
#include <complex>
//#include <direct.h>

PerformanceData PD;

qtSpectrumDisplay* gpSpectrumDisplay;

QGraphicsScene* distortionScene = NULL;
QGraphicsView* distortionWindow = NULL;
QGraphicsPixmapItem* pixmapItem = NULL;

vector<vector<SpecPair> > RefSpectrum;
vector<SpecPair> SimSpectrum;
vector<SpecPair> EtaK;
vector<SpecPair> EtaN;

//source spectrum
bool useSourceSpectrum = false;
vector<SpecPair> SourceSpectrum;

//resample the source based on the material samples (EtaK and EtaN)
vector<SpecPair> SourceResampled;

int currentSpec = 0;

double nuMin = 800;
double nuMax = 4000;
double dNu = 2;

double aMin = 0;
double aMax = 1;

double scaleI0 = 1.0;
double refSlope = 0.0;

bool dispRefSpec = true;
bool dispSimSpec = true;
bool dispSimK = true;
bool dispMatK = true;
bool dispSimN = true;
bool dispMatN = true;
double dispScaleK = 1.0;
double dispScaleN = 1.0;
SpecType dispSimType = AbsorbanceSpecType;
bool dispNormalize = false;
double dispNormFactor = 1.0;


//material parameters
double radius = 4.0f;
double baseIR = 1.49f;
double cA = 1.0;
//vector<SpecPair> KMaterial;
//vector<SpecPair> NMaterial;
bool applyMaterial = true;
vector<Material> MaterialList;
int currentMaterial = 0;

//optical parameters
double cNAi = 0.0;
double cNAo = 0.6;
double oNAi = 0.0;
double oNAo = 0.6;
OpticsType opticsMode = TransmissionOpticsType;
bool pointDetector = false;
int objectiveSamples = 200;

//fitting parameters
double minMSE = 0.00001;
int maxFitIter = 20;

void TempSimSpectrum()
{
	SpecPair temp;
	for(int i=800; i<4000; i++)
	{
		temp.nu = i;
		temp.A = sin((double)i/200);
		SimSpectrum.push_back(temp);
	}
}

void UpdateDisplay(){
	gpSpectrumDisplay->updateGL();
}

void LoadMaterial(string fileNameK, string fileNameN, string materialName)
{
	Material newMaterial;
	newMaterial.name = materialName;

	vector<SpecPair> KMaterial = LoadSpectrum(fileNameK.c_str());
	vector<SpecPair> NMaterial = LoadSpectrum(fileNameN.c_str());

	//make sure that the sizes are the same
	if(KMaterial.size() != NMaterial.size()){
		cout<<"Error, material properties don't match."<<endl;
		exit(1);
	}

	complex<double> eta;
	//int j;
	for(unsigned int i=0; i<KMaterial.size(); i++){
		newMaterial.nu.push_back(KMaterial[i].nu);
		eta = complex<double>(NMaterial[i].A, KMaterial[i].A);
		newMaterial.eta.push_back(eta);
	}
	MaterialList.push_back(newMaterial);
}

void LoadMaterial(string fileNameK, string materialName){

	//load the material absorbance
	vector<SpecPair> KMaterial = LoadSpectrum(fileNameK.c_str());
	vector<SpecPair> NMaterial;
	//KMaterial = LoadSpectrum("eta_TolueneK.txt");

	//compute the real IR using Kramers Kronig
	//copy the absorbance values into a linear array
	double* k = (double*)malloc(sizeof(double) * KMaterial.size());
	double* n = (double*)malloc(sizeof(double) * KMaterial.size());
	for(unsigned int i=0; i<KMaterial.size(); i++)
		k[i] = KMaterial[i].A;

	//use Kramers Kronig to determine the real part of the index of refraction
	cudaKramersKronig(n, k, KMaterial.size(), KMaterial[0].nu, KMaterial.back().nu, baseIR);
	SpecPair temp;
	for(unsigned int i=0; i<KMaterial.size(); i++)
	{
		temp.nu =  KMaterial[i].nu;
		temp.A = n[i];
		NMaterial.push_back(temp);
	}

	//create the material
	Material newMaterial;
	newMaterial.name = materialName;
	complex<double> eta;
	for(unsigned int i=0; i<KMaterial.size(); i++){
		newMaterial.nu.push_back(KMaterial[i].nu);
		eta = complex<double>(NMaterial[i].A, KMaterial[i].A);
		newMaterial.eta.push_back(eta);
	}

	MaterialList.push_back(newMaterial);
}

void LoadSource(string fileNameSource)
{
	SourceSpectrum = LoadSpectrum(fileNameSource);  
}

void ResampleSource()
{
	//clear the current resampled spectrum
	SourceResampled.clear();
	
	//get the number of source and material samples
	int nMatSamples = EtaK.size();
	int nSourceSamples = SourceSpectrum.size();
	
	float nu, I;
	for(int i=0; i<nMatSamples; i++)
	{
		SpecPair newSample;
		newSample.nu = EtaK[i].nu;
		
		//iterate through the SourceSpectrum to find the bounding wavelengths
		if(newSample.nu < SourceSpectrum[0].nu)
			newSample.A = 0.0;
		else if(newSample.nu > SourceSpectrum[nSourceSamples-1].nu)
			newSample.A = 0.0;
		else
		{
			int k=1;
			while(SourceSpectrum[k].nu < newSample.nu)
				k++;
			
			//interpolate
			float a = (newSample.nu - SourceSpectrum[k-1].nu)/(SourceSpectrum[k].nu - SourceSpectrum[k-1].nu);
			newSample.A = a * SourceSpectrum[k].A + (1.0 - a) * SourceSpectrum[k-1].A;
		}
		
		//insert the new spectral point into the resampled spectrum
		SourceResampled.push_back(newSample);
		//cout<<newSample.nu<<"   "<<newSample.A<<endl;
	}
  
}

void FitDisplay(){
	double minA = 99999.0;
	double maxA = -99999.0;
	double k, n;

	if(dispSimSpec)
		for(unsigned int i=0; i<SimSpectrum.size(); i++)
		{
			if(SimSpectrum[i].A < minA)
				minA = SimSpectrum[i].A;
			if(SimSpectrum[i].A > maxA)
				maxA = SimSpectrum[i].A;
		}

	if(dispRefSpec && RefSpectrum.size() > 0)
		for(unsigned int i=0; i<RefSpectrum[currentSpec].size(); i++)
		{
			if(RefSpectrum[currentSpec][i].A < minA)
				minA = RefSpectrum[currentSpec][i].A;
			if(RefSpectrum[currentSpec][i].A > maxA)
				maxA = RefSpectrum[currentSpec][i].A;
		}
	if(dispMatK)
		for(unsigned int i=0; i<EtaK.size(); i++)
		{
			k = MaterialList[currentMaterial].eta[i].imag() * dispScaleK;
			if(k < minA)
				minA = k;
			if(k > maxA)
				maxA = k;
		}
	if(dispSimK)
		for(unsigned int i=0; i<EtaK.size(); i++)
		{
			k = EtaK[i].A * dispScaleK;
			if(k < minA)
				minA = k;
			if(EtaK[i].A > maxA)
				maxA = k;
		}
	if(dispMatN)
		for(unsigned int i=0; i<EtaN.size(); i++)
		{
			n = (MaterialList[currentMaterial].eta[i].real() - baseIR) * dispScaleN;
			if(n < minA)
				minA = n;
			if(n > maxA)
				maxA = n;
		}
	if(dispSimN)
		for(unsigned int i=0; i<EtaN.size(); i++)
		{
			n = (EtaN[i].A - baseIR) * dispScaleN;
			if(n < minA)
				minA = n;
			if(n > maxA)
				maxA = n;
		}

	aMin = minA;
	aMax = maxA;
	UpdateDisplay();
}

void ChangeAbsorbance(){

	//compute the real part of the index of refraction

	//copy the absorbance values into a linear array
	int nSamples = MaterialList[currentMaterial].eta.size();
	double startNu = MaterialList[currentMaterial].nu.front();
	double endNu = MaterialList[currentMaterial].nu.back();
	double* k = (double*)malloc(sizeof(double) * nSamples);
	double* n = (double*)malloc(sizeof(double) * nSamples);
	for(int i=0; i<nSamples; i++)
		k[i] = MaterialList[currentMaterial].eta[i].imag() * cA;

	//NMaterial.clear();
	EtaK.clear();
	EtaN.clear();
	//use Kramers Kronig to determine the real part of the index of refraction
	cudaKramersKronig(n, k, nSamples, startNu, endNu, baseIR);

	//copy the real part of the index of refraction into the vector
	SpecPair temp;

	//load the imaginary IR from the absorbance data
	double nu;
	for(int i=0; i<nSamples; i++){
		nu = MaterialList[currentMaterial].nu[i];
		if(nu >= nuMin && nu <= nuMax){
			temp.nu = nu;
			temp.A = k[i];
			EtaK.push_back(temp);
			//temp.A = NMaterial[i].A;
			temp.A = n[i];
			EtaN.push_back(temp);
		}
	}

	free(k);
	free(n);
}

void SetMaterial()
{
	EtaK.clear();
	EtaN.clear();

	int nSamples = MaterialList[currentMaterial].eta.size();
	double nu;
	SpecPair temp;

	//initialize the current nuMin and nuMax values
	nuMin = MaterialList[currentMaterial].nu[0];
	nuMax = nuMin;
	for(int i=0; i<nSamples; i++){
		nu = MaterialList[currentMaterial].nu[i];
		//if(nu >= nuMin && nu <= nuMax){

		//update the min and max values for display
		if(nu < nuMin) nuMin = nu;
		if(nu > nuMax) nuMax = nu;

			temp.nu = nu;
			temp.A = MaterialList[currentMaterial].eta[i].imag();
			EtaK.push_back(temp);
			temp.A = MaterialList[currentMaterial].eta[i].real();
			EtaN.push_back(temp);
	}
	cA = 1.0;
	
	//resample the source spectrum
	if(SourceSpectrum.size() != 0)
		ResampleSource();

}


int main(int argc, char *argv[])
{


	//load the default project file (any previous optical settings)
	LoadState();

	//load the default materials
	LoadMaterial("eta_TolueneK.txt", "eta_TolueneN.txt", "Toluene");
	LoadMaterial("kPMMA.txt", "PMMA");
	LoadMaterial("eta_polystyreneK.txt", "Polystyrene");
	//LoadMaterial("../../../../data/materials/rtsSU8_k.txt", "../../../../data/materials/rtsSU8_n.txt", "SU8");
	SetMaterial();
	
	//load a mid-infrared source
	LoadSource("source_midIR.txt");
	ResampleSource();

	//compute the analytical solution for the Mie scattered spectrum
	SimulateSpectrum();

	QApplication a(argc, argv);
	
	InteractiveMie w;
	w.show() ;
	w.move(0, 0);
	QSize uiFrame = w.frameSize() + QSize(10, 10);
	cout<<"Frame: "<<uiFrame.width()<<endl;
	QSize uiNoFrame = w.size();
	cout<<"No Frame: "<<uiNoFrame.width()<<endl;
	int frameHeight = uiFrame.height() - uiNoFrame.height();

	//activate a console for output
	RedirectIOToConsole(0, uiFrame.height(), uiFrame.width(), 400);
	printf("Frame height: %d\n", frameHeight);

	//set the size and position of the spectrum window
	int visWinSize = uiFrame.height()/2 - frameHeight;
	//create the far field window
	gpSpectrumDisplay = new qtSpectrumDisplay();
	gpSpectrumDisplay->resize(visWinSize*2, visWinSize);
	gpSpectrumDisplay->move(uiFrame.width(), 0);
	gpSpectrumDisplay->show();

	//distortion dialog box
	distortionDialog = new qtDistortionDialog();
	distortionDialog->move(0, 0);

	//display the distortion map
	distortionScene = new QGraphicsScene();
	distortionWindow = new QGraphicsView(distortionScene);
	distortionWindow->move(uiFrame.width(), visWinSize);

	//refresh the UI
	w.refreshUI();

	return a.exec();
}