optics / imie

  #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 "qwtSpectrumDisplay.h"

  #include "globals.h"

  #include "rtsGUIConsole.h"

  #include "PerformanceData.h"

  #include <complex>

  #include <sstream>

  //#include <direct.h>

  

  PerformanceData PD;

  

  qtSpectrumDisplay* gpSpectrumDisplay;

  //qwtSpectrumDisplay* SpectrumDisplay;

  

  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 nMag = 1.0;

  double kMax = 1.0;

  

  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 = -1;

  

  //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();

  	//SpectrumDisplay->replot();

  }

  

  void LoadMaterial(string fileName, string materialName)

  {

  	//open the file

  	ifstream inFile(fileName.c_str());

  	if(!inFile)

  	{

  		cout<<"Error loading material: "<<fileName<<endl;

  		return;

  	}

  

  	Material newMaterial;

  	newMaterial.validN = false;

  	newMaterial.validK = false;

  	//read the header

  	string units;

  	getline(inFile, units, '\t');

  	char c = 0;

  	while(c != '\n')

  	{

  		inFile.get(c);

  		if(c == 'n') newMaterial.validN = true;

  		if(c == 'k') newMaterial.validK = true;

  	}

  	//read the entire refractive index (both real and imaginary)

  	float nu;

  	float n;

  	float k;

  	while(inFile>>nu)

  	{

  		n = 1.0;

  		k = 0.0;

  		if(newMaterial.validN)

  			inFile>>n;

  		if(newMaterial.validK)

  			inFile>>k;

  		//ignore the rest of the line

  		inFile.ignore();

  

  		newMaterial.nu.push_back(nu);

  		newMaterial.eta.push_back(complex<double>(n, k));

  	}

  	//iterate through the material to compute the necessary bounds for display

  	float minN = 9999;

  	float maxN = -9999;

  	float sumN = 0;

  	float maxK = 0;

  	for(int i = 0; i<newMaterial.nu.size(); i++)

  	{

  		n = newMaterial.eta[i].real();

  		k = newMaterial.eta[i].imag();

  

  		if(n < minN) minN = n;

  		if(n > maxN) maxN = n;

  		if(k > maxK) maxK = k;

  		sumN += n;

  	}

  	float meanN = sumN / newMaterial.nu.size();

  	nMag = max(maxN - meanN, meanN - minN);

  	kMax = maxK;

  	baseIR = meanN;

  

  	//set the name of the material object

  	newMaterial.name = materialName;

  	//add it to the material list

  	MaterialList.push_back(newMaterial);

  	currentMaterial = MaterialList.size() - 1;

  }

  

  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();

  	if(currentMaterial == -1) return;

  

  	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("etaToluene.txt", "Toluene");

  	LoadMaterial("kPMMA.txt", "PMMA");

  	LoadMaterial("kPolyethylene.txt", "Polyethylene");

  	LoadMaterial("kPTFE.txt", "Teflon");

  

  

  	//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);

  

  	//SpectrumDisplay = new qwtSpectrumDisplay();

  

  	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();

  }