deprecated / bimsim

//AnyOption for command-line processing
//#include "anyoption.h"

#include "rts/optics/material.h"

#include "nearfield.h"
#include "microscope.h"
#include "rts/graphics/colormap.h"
#include "fileout.h"
//extern nearfieldStruct* NF;
extern microscopeStruct* SCOPE;
extern fileoutStruct gFileOut;

//default values
#include "defaults.h"

#include <string>
#include <sstream>
#include <fstream>
#include <limits>
using namespace std;

#include <boost/program_options.hpp>
namespace po = boost::program_options;

extern bool verbose;



static void lNearfield(po::variables_map vm)
{
	//test to see if we are simulating a plane wave
	bool planeWave = DEFAULT_PLANEWAVE;
	if(vm.count("plane-wave"))
		planeWave = !planeWave;
	SCOPE->nf.planeWave = planeWave;

	//get the incident field amplitude
	SCOPE->nf.A = vm["amplitude"].as<ptype>();

	//get the condenser parameters
    SCOPE->nf.condenser[0] = DEFAULT_CONDENSER_MIN;
    SCOPE->nf.condenser[1] = DEFAULT_CONDENSER_MAX;

    if(vm.count("condenser"))
    {
        vector<ptype> cparams = vm["condenser"].as< vector<ptype> >();

        if(cparams.size() == 1)
            SCOPE->nf.condenser[1] = cparams[0];
        else
        {
            SCOPE->nf.condenser[0] = cparams[0];
            SCOPE->nf.condenser[1] = cparams[1];
        }
    }


	//get the focal rtsPoint position
    SCOPE->nf.focus[0] = DEFAULT_FOCUS_X;
    SCOPE->nf.focus[1] = DEFAULT_FOCUS_Y;
    SCOPE->nf.focus[2] = DEFAULT_FOCUS_Z;
    if(vm.count("focus"))
    {
        vector<ptype> fpos = vm["focus"].as< vector<ptype> >();
        if(fpos.size() != 3)
        {
            cout<<"BIMSIM Error - the incident focal point is incorrectly specified; it must have three components."<<endl;
            exit(1);
        }
        SCOPE->nf.focus[0] = fpos[0];
        SCOPE->nf.focus[1] = fpos[1];
        SCOPE->nf.focus[2] = fpos[2];
    }

	//get the incident light direction (k-vector)
	bsVector spherical(1, 0, 0);

    //if a k-vector is specified
    if(vm.count("k"))
    {
        vector<ptype> kvec = vm["k"].as< vector<ptype> >();
        if(kvec.size() != 2)
        {
            cout<<"BIMSIM Error - k-vector is not specified correctly: it must contain two elements"<<endl;
            exit(1);
        }
        spherical[1] = kvec[0];
        spherical[2] = kvec[1];
    }
	SCOPE->nf.k = spherical.sph2cart();


    //incident field order
    SCOPE->nf.m = vm["field-order"].as<int>();

    //number of Monte-Carlo samples
    SCOPE->nf.nWaves = vm["samples"].as<int>();

	//random number seed for Monte-Carlo samples
	if(vm.count("seed"))
		srand(vm["seed"].as<unsigned int>());



}


static void loadOutputParams(po::variables_map vm)
{
    //append simulation results to previous binary files
    gFileOut.append = DEFAULT_APPEND;
    if(vm.count("append"))
        gFileOut.append = true;

	//image parameters
	//component of the field to be saved
	std::string fieldStr;
    fieldStr = vm["output-type"].as<string>();

    if(fieldStr == "magnitude")
        gFileOut.field = fileoutStruct::fieldMag;
    else if(fieldStr == "intensity")
        gFileOut.field = fileoutStruct::fieldIntensity;
    else if(fieldStr == "polarization")
        gFileOut.field = fileoutStruct::fieldPolar;
    else if(fieldStr == "imaginary")
        gFileOut.field = fileoutStruct::fieldImag;
    else if(fieldStr == "real")
        gFileOut.field = fileoutStruct::fieldReal;
    else if(fieldStr == "angular-spectrum")
        gFileOut.field = fileoutStruct::fieldAngularSpectrum;


	//image file names
	gFileOut.intFile = vm["intensity"].as<string>();
	gFileOut.absFile = vm["absorbance"].as<string>();
	gFileOut.transFile = vm["transmittance"].as<string>();
	gFileOut.nearFile = vm["near-field"].as<string>();
	gFileOut.farFile = vm["far-field"].as<string>();

	//colormap
	std::string cmapStr;
    cmapStr = vm["colormap"].as<string>();
    if(cmapStr == "brewer")
        gFileOut.colormap = rts::cmBrewer;
    else if(cmapStr == "gray")
        gFileOut.colormap = rts::cmGrayscale;
    else
        cout<<"color-map value not recognized (using default): "<<cmapStr<<endl;
}

void lFlags(po::variables_map vm, po::options_description desc)
{
    //display help and exit
	if(vm.count("help"))
	{
		cout<<desc<<endl;
		exit(1);
	}

    //flag for verbose output
	if(vm.count("verbose"))
        verbose = true;

    if(vm.count("recursive"))
    {
        SCOPE->nf.lut_us = false;
        SCOPE->nf.lut_uf = false;
    }
    else if(vm.count("recursive-us"))
    {
        SCOPE->nf.lut_us = false;
    }
    else if(vm.count("lut-uf"))
    {
        SCOPE->nf.lut_uf = true;
    }
}

void lWavelength(po::variables_map vm)
{
    //load the wavelength
	if(vm.count("nu"))
	{
		//wavelength is given in wavenumber - transform and flag
		SCOPE->nf.lambda = 10000/vm["nu"].as<ptype>();
		gFileOut.wavenumber = true;
	}
	//otherwise we are using lambda = wavelength
	else
	{
		SCOPE->nf.lambda = vm["lambda"].as<ptype>();
		gFileOut.wavenumber = false;
	}
}

static void lSpheres(string sphereList)
{
    /*This function loads a list of sphere given in the string sphereList
        The format is:
            x y z a m
        where
            x, y, z = sphere position (required)
            a = sphere radius (required)
            m = material ID (optional) */

    std::stringstream ss(sphereList);

    while(!ss.eof())
    {
        //create a new sphere
        sphere newS;

        //get the sphere data
        ss>>newS.p[0];
        ss>>newS.p[1];
        ss>>newS.p[2];
        ss>>newS.a;

        if(ss.peek() != '\n')
            ss>>newS.iMaterial;

        //add the new sphere to the sphere vector
        SCOPE->nf.sVector.push_back(newS);

        //ignore the rest of the line
        ss.ignore(std::numeric_limits<std::streamsize>::max(), '\n');

        //check out the next element (this should set the EOF error flag)
        ss.peek();
    }
}

void lSpheres(po::variables_map vm)
{
    //if a sphere is specified at the command line
    if(vm.count("spheres"))
    {
        //convert the sphere to a string
        vector<ptype> sdesc = vm["spheres"].as< vector<ptype> >();

        //compute the number of spheres specified
        unsigned int nS;
        if(sdesc.size() <= 5)
            nS = 1;
        else
        {
            //if the number of parameters is divisible by 4, compute the number of spheres
            if(sdesc.size() % 5 == 0)
                nS = sdesc.size() / 5;
            else
            {
                cout<<"BIMSIM Error: Invalid number of sphere parameters."<<endl;
                exit(1);
            }
        }

        stringstream ss;

        //for each sphere
        for(unsigned int s=0; s<nS; s++)
        {
            //compute the number of sphere parameters
            unsigned int nP;
            if(nS == 1) nP = sdesc.size();
            else nP = 5;

            //store each parameter as a string
            for(unsigned int i=0; i<nP; i++)
            {
                ss<<sdesc[s*5 + i]<<" ";
            }
            ss<<endl;
        }



        //convert the string to a sphere list
        lSpheres(ss.str());
    }

    //if a files are specified
    if(vm.count("sphere-file"))
    {

        vector<string> filenames = vm["sphere-file"].as< vector<string> >();
        //load each file
        for(int iS=0; iS<filenames.size(); iS++)
        {
            //load the file into a string
            std::ifstream ifs(filenames[iS].c_str());

            if(!ifs)
            {
                cout<<"Error loading sphere file."<<endl;
                exit(1);
            }

            std::string instr((std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>());

            //load the list of spheres from a string
            lSpheres(instr);
        }
    }

    //make sure the appropriate materials are loaded
    unsigned int nS = SCOPE->nf.sVector.size();

    //for each sphere
    for(unsigned int s = 0; s<nS; s++)
    {
        //make sure the corresponding material exists
        if(SCOPE->nf.sVector[s].iMaterial + 1 > SCOPE->nf.mVector.size())
        {
            //otherwise output an error
            cout<<"BIMSIM Error - A material is not loaded for sphere "<<s+1<<"."<<endl;
            exit(1);
        }
    }
}

static void lMaterials(po::variables_map vm)
{
	//if materials are specified at the command line
	if(vm.count("materials"))
	{
		vector<ptype> matVec = vm["materials"].as< vector<ptype> >();
		if(matVec.size() == 1)
		{
			rts::material<ptype> newM(SCOPE->nf.lambda, matVec[0], 0);
			SCOPE->nf.mVector.push_back(newM);
		}
		else if(matVec.size() %2 != 0)
		{
			cout<<"BIMSim Error: materials must be specified in n, k pairs"<<endl;
			exit(1);
		}
		else
		{
			for(int i=0; i<matVec.size(); i+=2)
			{
				rts::material<ptype> newM(SCOPE->nf.lambda, matVec[i], matVec[i+1]);
				SCOPE->nf.mVector.push_back(newM);
			}
		}
	}

	//if file names are specified, load the materials
	if(vm.count("material-file"))
	{
        vector<string> filenames = vm["material-file"].as< vector<string> >();
        for(int i=0; i<filenames.size(); i++)
        {
            //load the file into a string
            std::ifstream ifs(filenames[i].c_str());

            std::string instr((std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>());

            //load the list of spheres from a string
            rts::material<ptype> newM;
            newM.fromStr(instr, "");
            SCOPE->nf.mVector.push_back(newM);
        }
	}

}

static void lOptics(po::variables_map vm)
{
    SCOPE->objective[0] = DEFAULT_OBJECTIVE_MIN;
    SCOPE->objective[1] = DEFAULT_OBJECTIVE_MAX;
    if(vm.count("objective"))
    {
        vector<ptype> oparams = vm["objective"].as< vector<ptype> >();

        if(oparams.size() == 1)
            SCOPE->objective[1] = oparams[0];
        else
        {
            SCOPE->objective[0] = oparams[0];
            SCOPE->objective[1] = oparams[1];
        }
    }
}

static void lImagePlane(po::variables_map vm)
{
	bsPoint pMin(DEFAULT_PLANE_MIN_X, DEFAULT_PLANE_MIN_Y, DEFAULT_PLANE_MIN_Z);
	bsPoint pMax(DEFAULT_PLANE_MAX_X, DEFAULT_PLANE_MAX_Y, DEFAULT_PLANE_MAX_Z);
	bsVector normal(DEFAULT_PLANE_NORM_X, DEFAULT_PLANE_NORM_Y, DEFAULT_PLANE_NORM_Z);

	//set the default values for the slice position and orientation
	if(vm.count("plane-lower-left") && vm.count("plane-upper-right") && vm.count("plane-normal"))
	{
		vector<ptype> ll = vm["plane-lower-left"].as< vector<ptype> >();
		if(ll.size() != 3)
		{
			cout<<"BIMSIM Error - The lower-left corner of the image plane is incorrectly specified."<<endl;
			exit(1);
		}

		vector<ptype> ur = vm["plane-lower-left"].as< vector<ptype> >();
		if(ur.size() != 3)
		{
			cout<<"BIMSIM Error - The upper-right corner of the image plane is incorrectly specified."<<endl;
			exit(1);
		}

		vector<ptype> norm = vm["plane-lower-left"].as< vector<ptype> >();
		if(norm.size() != 3)
		{
			cout<<"BIMSIM Error - The normal of the image plane is incorrectly specified."<<endl;
			exit(1);
		}

		pMin = bsPoint(ll[0], ll[1], ll[2]);
		pMax = bsPoint(ur[0], ur[1], ur[2]);
		normal = bsVector(norm[0], norm[1], norm[2]);
	}
	else if(vm.count("xy"))
	{
		//default plane size in microns
		ptype s = DEFAULT_PLANE_SIZE;
		ptype pos = DEFAULT_PLANE_POSITION;

		vector<ptype> xy = vm["xy"].as< vector<ptype> >();
		if(xy.size() >= 1)
			s = xy[0];
		if(xy.size() >= 2)
			pos = xy[1];

		//calculate the plane corners and normal based on the size and position
		pMin = bsPoint(-s/2, -s/2, pos);
		pMax = bsPoint(s/2, s/2, pos);
		normal = bsVector(0, 0, 1);
	}
	else if(vm.count("xz"))
	{
		//default plane size in microns
		ptype size = DEFAULT_PLANE_SIZE;
		ptype pos = DEFAULT_PLANE_POSITION;

		vector<ptype> xz = vm["xz"].as< vector<ptype> >();
		if(xz.size() >= 1)
			size = xz[0];
		if(xz.size() >= 2)
			pos = xz[1];

		//calculate the plane corners and normal based on the size and position
		pMin = bsPoint(-size/2, pos, -size/2);
		pMax = bsPoint(size/2, pos, size/2);
		normal = bsVector(0, -1, 0);
	}
	else if(vm.count("yz"))
	{
		//default plane size in microns
		ptype size = DEFAULT_PLANE_SIZE;
		ptype pos = DEFAULT_PLANE_POSITION;

		vector<ptype> yz = vm["yz"].as< vector<ptype> >();
		if(yz.size() >= 1)
			size = yz[0];
		if(yz.size() >= 2)
			pos = yz[1];

		//calculate the plane corners and normal based on the size and position
		pMin = bsPoint(pos, -size/2, -size/2);
		pMax = bsPoint(pos, size/2, size/2);
		normal = bsVector(1, 0, 0);
	}
	SCOPE->setPos(pMin, pMax, normal);

	//resolution
	SCOPE->setRes(vm["resolution"].as<unsigned int>(),
				  vm["resolution"].as<unsigned int>(),
				  vm["padding"].as<unsigned int>(),
				  vm["supersample"].as<unsigned int>());





	SCOPE->setNearfield();
}

static void OutputOptions()
{
	cout<<SCOPE->toStr();

	cout<<"# of source points: "<<SCOPE->focalPoints.size()<<endl;

}

vector<ptype> test;
static void SetOptions(po::options_description &desc)
{
	desc.add_options()
		("help", "prints this help")
		("verbose", "verbose output\n")

		("intensity", po::value<string>()->default_value(DEFAULT_INTENSITY_FILE), "output measured intensity (filename)")
		("absorbance", po::value<string>()->default_value(DEFAULT_ABSORBANCE_FILE), "output measured absorbance (filename)")
		("transmittance", po::value<string>()->default_value(DEFAULT_TRANSMITTANCE_FILE), "output measured transmittance (filename)")
		("far-field", po::value<string>()->default_value(DEFAULT_FAR_FILE), "output far-field at detector (filename)")
		("near-field", po::value<string>()->default_value(DEFAULT_NEAR_FILE), "output field at focal plane (filename)")
		("extended-source", po::value<string>()->default_value(DEFAULT_EXTENDED_SOURCE), "image of source at focus (filename)\n")

		("spheres", po::value< vector<ptype> >()->multitoken(), "sphere position: x y z a m")
		("sphere-file", po::value< vector<string> >()->multitoken(), "sphere file:\n [x y z radius material]")
		("materials", po::value< vector<ptype> >()->multitoken(), "refractive indices as n, k pairs:\n ex. -m n0 k0 n1 k1 n2 k2")
		("material-file", po::value< vector<string> >()->multitoken(), "material file:\n [lambda n k]\n")

		("lambda", po::value<ptype>()->default_value(DEFAULT_LAMBDA), "incident wavelength")
		("nu", po::value<ptype>(), "incident frequency (in cm^-1)\n(if specified, lambda is ignored)")
		("k", po::value< vector<ptype> >()->multitoken(), "k-vector direction: -k theta phi\n theta = [0 2*pi], phi = [0 pi]")
		("amplitude", po::value<ptype>()->default_value(DEFAULT_AMPLITUDE), "incident field amplitude")
		("condenser", po::value< vector<ptype> >()->multitoken(), "condenser numerical aperature\nA pair of values can be used to specify an inner obscuration: -c NAin NAout")
		("objective", po::value< vector<ptype> >()->multitoken(), "objective numerical aperature\nA pair of values can be used to specify an inner obscuration: -c NAin NAout")
		("focus", po::value< vector<ptype> >()->multitoken(), "focal position for the incident point source\n (default = --focus 0 0 0)")
		("plane-wave", "simulates an incident plane wave\n")

		("resolution", po::value<unsigned int>()->default_value(DEFAULT_SLICE_RES), "resolution of the detector")
		("plane-lower-left", po::value< vector<ptype> >()->multitoken(), "lower-left position of the image plane")
		("plane-upper-right", po::value< vector<ptype> >()->multitoken(), "upper-right position of the image plane")
		("plane-normal", po::value< vector<ptype> >()->multitoken(), "normal for the image plane")
		("xy", po::value< vector<ptype> >()->multitoken(), "specify an x-y image plane\n (standard microscope)")
		("xz", po::value< vector<ptype> >()->multitoken(), "specify a x-z image plane\n (cross-section of the focal volume)")
		("yz", po::value< vector<ptype> >()->multitoken(), "specify a y-z image plane\n (cross-section of the focal volume)\n")

		("samples", po::value<int>()->default_value(DEFAULT_SAMPLES), "Monte-Carlo samples used to compute Us")
		("padding", po::value<unsigned int>()->default_value(DEFAULT_PADDING), "FFT padding for the objective bandpass")
		("supersample", po::value<unsigned int>()->default_value(DEFAULT_SUPERSAMPLE), "super-sampling rate for the detector field")
		("field-order", po::value<int>()->default_value(DEFAULT_FIELD_ORDER), "order of the incident field")
		("seed", po::value<unsigned int>(), "seed for the Monte-Carlo random number generator")
		("recursive", "evaluate all Bessel functions recursively\n")
		("recursive-us", "evaluate scattered-field Bessel functions recursively\n")
		("lut-uf", "evaluate the focused-field using a look-up table\n")

		("output-type", po::value<string>()->default_value(DEFAULT_FIELD_TYPE), "output field value:\n magnitude, polarization, real, imaginary, angular-spectrum")
		("colormap", po::value<string>()->default_value(DEFAULT_COLORMAP), "colormap: gray, brewer")
		("append", "append result to an existing file\n (binary files only)")
		;
}

static void LoadParameters(int argc, char *argv[])
{
	//create an option description
	po::options_description desc("BimSim arguments");

	//fill it with options
	SetOptions(desc);

    po::variables_map vm;
	po::store(po::parse_command_line(argc, argv, desc, po::command_line_style::unix_style ^ po::command_line_style::allow_short), vm);
	po::notify(vm);


    //load flags (help, verbose output)
    lFlags(vm, desc);

    //load the wavelength
    lWavelength(vm);

    //load materials
	//loadMaterials(vm);
	lMaterials(vm);

    //load the sphere data
    lSpheres(vm);

    //load the optics
    lOptics(vm);

	//load the position and orientation of the image plane
	lImagePlane(vm);

	//load spheres
	//loadSpheres(vm);



	lNearfield(vm);

	loadOutputParams(vm);

	//loadMicroscopeParams(vm);

	//loadSliceParams(vm);

    //if an extended source will be used
    if(vm["extended-source"].as<string>() != "")
    {
        //load the point sources
        string filename = vm["extended-source"].as<string>();
        SCOPE->LoadExtendedSource(filename);

    }





}