~brenthuisman/cpp_headers

0a10511c2ab7bbc9af1a7534ca71a5391678a4ad — Brent Huisman 1 year, 10 months ago master
init
6 files changed, 1909 insertions(+), 0 deletions(-)

A include/ct.h
A include/fluence.h
A include/image.h
A include/rt.h
A include/settings.h
A include/tools.h
A  => include/ct.h +243 -0
@@ 1,243 @@
#pragma once

#include <istream>
#include <assert.h>

#include "pystring.h"
using namespace pystring;
using std::string;
#include "tools.h"
using namespace vect;
#include "image.h"

#include "Phantom.h"
#include "settings.h"


/*
 * Sets up and returns a (GPUMCD) phantom.
 */


class CT;

class CT{
public:
	//members needed by gpumcd
	Phantom phantom;
	vector<string> materials;

	//ctors
	CT() = default;
	CT(const CT &) = default;
	~CT() = default;

	CT(DosiaSettings &, BeamMetaData &);

	//methods
	int num_vox(){ return image.nvox(); };
	Image generate_image(const vector<float> &);

private:
	//members
	BeamMetaData beamMetaData;
	DosiaSettings sett;
	Image image;

	//methods
	Phantom generate_phantom(const Image &, const string &, const string &);
	template <typename T>
	void set_hu2density(const string &, const vector<T> &, vector<float> &);
	//void set_hu2material(const string & = "hu2mat.ini"); //use with Schneider data from Gate
	void set_density2material(const string &, const vector<float> &, vector<float> &);
};


CT::CT(DosiaSettings &_sett, BeamMetaData &_beamMetaData) : beamMetaData(_beamMetaData), sett(_sett){

	string &rt_files = sett.rt_files;
	image = Image(os::path::join(rt_files, "ct.xdr"));
	if (sett.verbose > 1) cerr << "phantom file loaded: " << os::path::join(rt_files, "ct.xdr") << "\n";
	
	phantom = generate_phantom(image, os::path::join(sett.hounsfield_conversion_dir, "hu2dens.ini"), os::path::join(sett.hounsfield_conversion_dir, "dens2mat.ini"));

	if (sett.verbose > 1) {
		fprintf(stderr, "voxelSizes: %.2f,%.2f,%.2f\n", phantom.voxelSizes.x, phantom.voxelSizes.y, phantom.voxelSizes.z);
		fprintf(stderr, "numVoxels: %i,%i,%i\n", phantom.numVoxels.x, phantom.numVoxels.y, phantom.numVoxels.z);
		fprintf(stderr, "phantomCorner: %.2f,%.2f,%.2f\n", phantom.phantomCorner.x, phantom.phantomCorner.y, phantom.phantomCorner.z);
	}

	if (sett.dbgoutput){
		Image mediumIndex = generate_image(phantom.mediumIndexArray);
		Image massDensity = generate_image(phantom.massDensityArray);
		mediumIndex.write(os::path::join(rt_files, "mediumIndex.xdr"));
		massDensity.write(os::path::join(rt_files, "massDensityArray.xdr"));
	}
};


Phantom CT::generate_phantom(const Image &im, const string &hu2dens_fname, const string &dens2mat_fname){
	//think of this function as a constructor for Phantom structures
	assert(im.ndim() == 3);

	Phantom phantom;
	vector<float> ct_voxels; // gpumcd cares only about floats

	//setup coords,dimensions

	phantom.numVoxels.x = im.dim_size[0];
	phantom.numVoxels.y = im.dim_size[1];
	phantom.numVoxels.z = im.dim_size[2];

	phantom.phantomCorner.x = im.min_ext[0];
	phantom.phantomCorner.y = im.min_ext[1];
	phantom.phantomCorner.z = im.min_ext[2];

	phantom.voxelSizes.x = im.voxel_sizes[0];
	phantom.voxelSizes.y = im.voxel_sizes[1];
	phantom.voxelSizes.z = im.voxel_sizes[2];

	//now, move phantomCorner a half voxel

	phantom.phantomCorner.x -= phantom.voxelSizes.x / 2.;
	phantom.phantomCorner.y -= phantom.voxelSizes.y / 2.;
	phantom.phantomCorner.z -= phantom.voxelSizes.z / 2.;

	//convert image to HU units
	ct_voxels.resize(im.nvox());
	for (int i = 0; i < im.nvox(); i++){
		ct_voxels[i] = im.imdata[i] * beamMetaData.hu_slope + beamMetaData.hu_intercept;
	}
	//convert HU units to mass density and materials.
	set_hu2density(hu2dens_fname, ct_voxels, phantom.massDensityArray);
	set_density2material(dens2mat_fname, phantom.massDensityArray, phantom.mediumIndexArray);

	if (sett.in_aqua_vivo || sett.score_and_transport_in_water || sett.score_dose_to_water) {
		// set ref medium to water
		materials.push_back("Water"); //preserve existing materials, because GPUMCD still checks presence.
		sett.physicsSettings.referenceMedium = materials.size() - 1;
	}
	if (sett.in_aqua_vivo || sett.score_and_transport_in_water) {
		// set all voxels to water (we added it above at the end of 'materials')
		std::fill(phantom.mediumIndexArray.begin(), phantom.mediumIndexArray.end(), materials.size() - 1);
	}
	if (sett.in_aqua_vivo) {
		// set all voxels density 1. FIXME: currently sets all to 1, must set 1 only to patient mask
		
		//std::fill(phantom.massDensityArray.begin(), phantom.massDensityArray.end(), 1.f);
	}
	
	return phantom;
}


Image CT::generate_image(const vector<float> &new_voxels){
	assert(new_voxels.size() == image.nvox());
	Image ret(image);
	ret.imdata = new_voxels;
	//ret.imdata.insert(ret.imdata.begin(), new_voxels.begin(), new_voxels.end());
	return ret;
}


template <typename T>
void CT::set_hu2density(const string &fname, const vector<T> &ct_voxels, vector<float> &massDensityArray) {
	io::isfile(fname, 43);

	vector<float> density_hu;
	vector<float> density;

	std::ifstream is(fname);
	string str;
	while (getline(is, str))
	{
		split(str);
		density_hu.push_back( stoi(split(str)[0]) );
		density.push_back(stof(split(str)[1]));
	}

	massDensityArray.resize(num_vox());

	for (int i = 0; i < num_vox(); i++){
		massDensityArray[i] = interpolate(density_hu, density, ct_voxels[i], true);
		if (massDensityArray[i] < 0) massDensityArray[i] = 0.f;
	}
};


void CT::set_density2material(const string &fname, const vector<float> &massDensityArray, vector<float> &mediumIndexArray) {
	assert(massDensityArray.size()>0);//this ensures the densities are available.
	io::isfile(fname, 45);
	
	//use dens2mat with data from AvL clinic and Monaco defaults in appendix C of research manual.
	//NOTE: indices are continuous: a matindex of 1.4 will be a mix of 40% material 1 and 60% material 2 by weight.
	//gpumcd uses continuous material indices to mix materials. this vector helps.
	vector<int> continuous_material_index_axis;
	int i = 0;

	//load file into following. 'matireals' is a class member, because gpumcd later needs it
	vector<std::pair<int, int>> material_hu;
	vector<float> material_dens; //if used, indices must map to 'material' member

	std::ifstream is(fname);
	string str;
	while (getline(is, str))
	{
		split(str);
		material_dens.push_back(stof(split(str)[0]));
		materials.push_back(split(str)[1]);
		continuous_material_index_axis.push_back(i++);
	}

	mediumIndexArray.resize(num_vox());

	if (sett.continous_materials){
		for (int i = 0; i < num_vox(); i++){
			mediumIndexArray[i] = interpolate(material_dens, continuous_material_index_axis, massDensityArray[i], false);
			//if (phantom.mediumIndexArray[i] < 1) phantom.mediumIndexArray[i] = 0.f; //clip first materials (probably air)
		}
	}
	else { //integer materials
		for (int i = 0; i < num_vox(); i++){
			mediumIndexArray[i] = 0; //any densities lower than minimum found in file are assumed to be of material at first line.
			for (int j = 0; j < material_dens.size(); j++){
				if (material_dens[j]>massDensityArray[i]){
					break;
				}
				mediumIndexArray[i] = j;
			}
		}
	}
};



/*
void CT::set_hu2material(const string &fname){
	//not used, for compat with Gate
	io::isfile(fname, 44);

	vector<std::pair<int, int>> material_hu;
	vector<float> material_dens; //if used, indices must map to 'material' member


	std::ifstream is(fname);
	string str;
	while (getline(is, str))
	{
		split(str);
		material_hu.push_back({ stoi(split(str)[0]), stoi(split(str)[1]) });
		materials.push_back(split(str)[2]);
	}

	phantom.mediumIndexArray.resize(num_vox());

	for (int i = 0; i < num_vox(); i++){
		for (int j = 0; j < material_hu.size(); j++){
			if (material_hu[j].first < ct_voxels[i] && ct_voxels[i] < material_hu[j].second){
				phantom.mediumIndexArray[i] = j;
			}
		}
	}
};*/


A  => include/fluence.h +53 -0
@@ 1,53 @@
#pragma once

//#include <set>
#include <assert.h>

#include "CalculationInformation.h"

#include <string>
#include "pystring.h"
using namespace pystring;
using std::string;
#include "tools.h"
using namespace parse;
using namespace vect;
using std::vector;
#include "settings.h"

#include "image.h"

class FluenceImage;

class FluenceBeam;


class FluenceImage : private Image {
	FluenceImage(const std::string &);// : Image(const string &);
	std::vector<BeamInformation> beaminfos;
};

FluenceImage::FluenceImage(const string &fn) : Image(fn) {
	assert(ndim() == 2);
	beaminfos = std::vector<BeamInformation>(nvox());

	//x-coord panel ligt tussen min_ext[0]+voxel_sizes[0]*x waar x<max_ext[0]
	/*float relativeWeight;
	Float3 isoCenter;
	std::pair<float, float> gantryAngle;
	std::pair<float, float> couchAngle;
	std::pair<float, float> collimatorAngle;
	std::pair<float, float> fieldMin;
	std::pair<float, float> fieldMax;*/

	for (int y = 0; y < dim_size[1]; y++){
		for (int x = 0; x < dim_size[0]; x++){
			// todo convert x,y afmetingen in isoc
			// plus and minus halfpixel
			beaminfos[x + y*x].relativeWeight = imdata[x + y*x];
			beaminfos[x + y*x].fieldMin = { min_ext[0] + (x - 0.5)*voxel_sizes[0], min_ext[1] + (y - 0.5)*voxel_sizes[1] };
			beaminfos[x + y*x].fieldMax = { min_ext[0] + (x + 0.5)*voxel_sizes[0], min_ext[1] + (y + 0.5)*voxel_sizes[1] };
		}
	}

};
\ No newline at end of file

A  => include/image.h +323 -0
@@ 1,323 @@
#pragma once

#include <iostream>
#include "pystring.h" //pystring
#include "tools.h" //vect,std::vector
#include <cstring> //std::memcpy
#include "gpumcd/Phantom.h"
using namespace vect;
using namespace pystring;

class Image{
public:
	vector<int> dim_size;
	vector<float> voxel_sizes;
	vector<float> min_ext;
	vector<float> max_ext;
	vector<float> imdata; // watch out! I convert all to floats! my world is simple!

	Image() = default;
	~Image() = default;
	Image(const std::string &);

	void write(const std::string &);
	Image copy_with_new_voxels(const vector<float> &);

	int ndim() const { return dim_size.size(); };
	int nvox() const { return mul(dim_size); };

private:
	void read_xdr(const std::string &);
	void read_mhd(const std::string &);

	void write_xdr(const std::string &);
	void write_mhd(const std::string &);
};


Image::Image(const std::string &fname){
	if (pystring::endswith(fname, ".xdr")){
		read_xdr(fname);
	}
	else if (pystring::endswith(fname, ".mhd")){
		read_mhd(fname);
	}
	else {
		printf("Unkown file extension encountered, aborting write out...");
	}
}


void Image::write(const std::string &fname){
	if (pystring::endswith(fname, ".xdr")){
		write_xdr(fname);
	}
	else if (pystring::endswith(fname, ".mhd")){
		write_mhd(fname);
	}
	else {
		printf("Unkown file extension encountered, aborting write out...");
	}
}


void Image::read_xdr(const std::string &xdrfile) {
	//Phantom phantom;

	auto xdr = fromfile<char>(xdrfile);

	std::string header;

	//calc offset of magic bytes
	char lasti = ' ';
	long imdata_offset = 0;
	for (char i : xdr){
		imdata_offset++;
		if (i == 0x0c && lasti == 0x0c){
			break;
		}
		lasti = i;
		header += i;
	}
	header.pop_back(); //one magic byte was added.

	int type = -1; //2 = >i2, 4 = >f4

	for (const auto &line : parse::parse_dump(header)){
		if (startswith(line.first, "ndim")) {
			int _ndim = stoi(line.second);
			assert(_ndim > 0 && _ndim < 4);
			dim_size.resize(_ndim);
			voxel_sizes.resize(_ndim);
			min_ext.resize(_ndim);
			max_ext.resize(_ndim);
			continue;
		}
		if (startswith(line.first, "field")) {
			assert(startswith(line.second, "uniform"));
			continue;
		}
		if (startswith(line.first, "data")) {
			if (startswith(line.second, "xdr_short")){
				type = 2;
			}
			else if (startswith(line.second, "xdr_real") || startswith(line.second, "xdr_float")){
				type = 4;
			}
			else{
				assert(type != -1); //blow up
			}
			continue;
		}
		if (startswith(line.first, "dim1")) {
			dim_size[0] = stoi(line.second);
			continue;
		}
		if (startswith(line.first, "dim2")) {
			dim_size[1] = stoi(line.second);
			continue;
		}
		if (startswith(line.first, "dim3")) {
			dim_size[2] = stoi(line.second);
			continue;
		}
		/*if (startswith(line.first, "min_ext")) {
			min_ext.push_back(stof(split(line.second)[0]));
			min_ext.push_back(stof(split(line.second)[1]));
			min_ext.push_back(stof(split(line.second)[2]));
			continue;
			}
			if (startswith(line.first, "max_ext")) {
			max_ext.push_back(stof(split(line.second)[0]));
			max_ext.push_back(stof(split(line.second)[1]));
			max_ext.push_back(stof(split(line.second)[2]));
			continue;
			}*/
	}

	//header loaded.

	long imdata_bytes = nvox()*type;
	long ext_offset = imdata_offset + imdata_bytes;
	long ext_bytes = ndim() * 2 * sizeof(float);

	//check that the size of the .xdr corresponds to the header+voxels*voxeltype+exts:
	assert(xdr.size() == imdata_offset + imdata_bytes + ext_bytes);

	vector<char> _imdata(xdr.begin() + imdata_offset, xdr.begin() + imdata_offset + imdata_bytes);

	if (type == 2){
		types::swap_endianness<short>(_imdata);
		vector<short> shorts = types::reinterpret<short>(_imdata);
		imdata.insert(imdata.begin(), shorts.begin(), shorts.end()); //this upcasts shorts
	}
	else if (type == 4){
		types::swap_endianness<float>(_imdata);
		imdata = types::reinterpret<float>(_imdata);
	}

	//now the extents in the final ndim*2*sizeof(float) bytes
	//write extents, looped pairwise over axis
	//xmin, xmax, ymin, ymax, zmin, zmax

	vector<char> _exts(xdr.begin() + ext_offset, xdr.begin() + ext_offset + ext_bytes);
	types::swap_endianness<float>(_exts);
	vector<float> exts = types::reinterpret<float>(_exts);

	for (size_t i = 0; i < ndim(); i++) {
		min_ext[i] = exts[2 * i];
		max_ext[i] = exts[2 * i + 1];
		//calc binsize
		voxel_sizes[i] = (max_ext[i] - min_ext[i]) / (dim_size[i] - 1);
	}
}


void Image::write_mhd(const std::string &fn){
	FILE* ffile = fopen(fn.c_str(), "wb");
	if (ffile == nullptr){
		throw std::pair<int, std::string>(70, "Problem writing file '" + fn + "'.");
	}
	// TODO:header
	fprintf(ffile, "ObjectType = Image\n");
	fprintf(ffile, "NDims=%i\n", ndim());
	fprintf(ffile, "BinaryData = True\n");
	fprintf(ffile, "BinaryDataByteOrderMSB = False\n");
	fprintf(ffile, "CompressedData = False\n");
	fprintf(ffile, "Offset = ");
	for (int i = 0; i < ndim(); i++) {
		fprintf(ffile, "%f ", min_ext[i]*10);
	}
	fprintf(ffile, "\n");
	fprintf(ffile, "ElementSpacing = ");
	for (int i = 0; i < ndim(); i++) {
		fprintf(ffile, "%f ", voxel_sizes[i]*10);
	}
	fprintf(ffile, "\n");
	fprintf(ffile, "DimSize = ");
	for (int i = 0; i < ndim(); i++) {
		fprintf(ffile, "%i ", dim_size[i]);
	}
	fprintf(ffile, "\n");
	fprintf(ffile, "ElementType = MET_FLOAT\n");
	std::string rawfile;
	std::string ext;
	os::path::splitext(rawfile,ext,fn);
	rawfile += ".raw";
	fprintf(ffile, "ElementDataFile = %s\n",rawfile.c_str());
	fclose(ffile);

	//write rawfile
	tofile<float>(imdata,rawfile);
}


void Image::write_xdr(const std::string &fn){
    FILE* ffile = fopen(fn.c_str(), "wb");
    if (ffile == nullptr){
        throw std::pair<int, std::string>(70, "Problem writing file '" + fn + "'.");
    }
    fprintf(ffile, "# AVS WRITER BY BRENT\n");
    fprintf(ffile, "ndim=%i\n",ndim());
	for (int i = 0; i < ndim(); i++) {
        fprintf(ffile, "dim%i=%i\n",i+1,dim_size[i]);
    }
    fprintf(ffile, "nspace=%i\n",ndim()); //dont know if used
    fprintf(ffile, "veclen=1\n"); //dont know if used
    fprintf(ffile, "data=xdr_real\n"); //dont know if used
    fprintf(ffile, "field=uniform\n"); //dont know if used
    fprintf(ffile, "%c%c",0x0c,0x0c); //magic bytes

    //convert imdata to bytes
	vector<char> imdata_swapped = types::reinterpret<char>(imdata);
	types::swap_endianness<float>(imdata_swapped);
    fwrite(imdata_swapped.data(), sizeof(char), imdata_swapped.size(), ffile);

	//extents
	//xmin, xmax, ymin, ymax, zmin, zmax
	vector<float> exts(ndim() * 2);
	for (size_t i = 0; i < ndim(); i++) {
		exts[2 * i]=min_ext[i];
		exts[2 * i + 1]=max_ext[i];
	}
	vector<char> exts_swapped = types::reinterpret<char>(exts);
	types::swap_endianness<float>(exts_swapped);
    fwrite(exts_swapped.data(), sizeof(char), exts_swapped.size(), ffile);
    fclose(ffile);
}


void Image::read_mhd(const std::string &header) {
	std::string rawfile;
	int type = -1; //2 = >i2, 4 = >f4
	for (const auto &line : parse::load_dump(header)) {
		if (startswith(line.first, "NDims")) {
			int _ndim = stoi(line.second);
			assert(_ndim > 0 && _ndim < 4);
			dim_size.resize(_ndim);
			voxel_sizes.resize(_ndim);
			min_ext.resize(_ndim);
			max_ext.resize(_ndim);
			continue;
		}
		if (startswith(line.first, "BinaryData") && !startswith(line.first, "BinaryDataByteOrderMSB")) {
			assert(startswith(line.second, "True"));
			continue;
		}
		if (startswith(line.first, "BinaryDataByteOrderMSB")) {
			assert(startswith(line.second, "False"));
			continue;
		}
		if (startswith(line.first, "CompressedData")) {
			assert(startswith(line.second, "False"));
			continue;
		}
		if (startswith(line.first, "ElementSpacing")) {
			voxel_sizes = types::split<float>(line.second);
			for (auto &i : voxel_sizes){
				i /= 10;
			}
			continue;
		}
		if (startswith(line.first, "DimSize")) {
			dim_size = types::split<int>(line.second);
			continue;
		}
		if (startswith(line.first, "Offset")) {
			min_ext = types::split<float>(line.second);
			for (auto &i : min_ext){
				i /= 10;
			}
			continue;
		}
		if (startswith(line.first, "ElementDataFile")) {
			rawfile = line.second;
			continue;
		}
		if (startswith(line.first, "ElementType")) {
			if (startswith(line.second, "MET_SHORT")){
				type = 2;
			}
			else if (startswith(line.second, "MET_FLOAT")){
				type = 4;
			}
			else{
				assert(type != -1); //blow up
			}
			continue;
		}
	}

	if (type == 2){
		vector<short> shorts = fromfile<short>(rawfile);
		imdata.insert(imdata.begin(), shorts.begin(), shorts.end());
	}
	else if (type == 4){
		imdata = fromfile<float>(rawfile);
	}
	assert(imdata.size() == nvox());

	for (size_t i = 0; i < ndim(); i++) {
		max_ext[i] = min_ext[i] + voxel_sizes[i] * (dim_size[i] -1);
	}
}
\ No newline at end of file

A  => include/rt.h +812 -0
@@ 1,812 @@
#pragma once

//#include <set>
//#include <string>
//using std::string;
#include <assert.h>

#include "CalculationInformation.h"

#include "pystring.h"
using namespace pystring;
using std::string;
#include "tools.h"
using namespace parse;
using namespace vect;
#include "settings.h"

class RTBeam;
class Parser;
class dicomParser;
class pinnacleParser;

struct ControlPoint {
	// gpumcd computes one ModifierInformation+BeamInformation pair at a time, and a full description of a controlpoint consists of both
	ModifierInformation collimator; //ASMX == parallelJaw
	BeamInformation beamInfo;
};

enum class AcceleratorType { UNKNOWN, EMPTY, MLCi80, Agility, MRLinac };
enum class Energy { UNKNOWN, MV6, MV7, MV10 };
enum class Filter { UNKNOWN, FF, NoFF }; //YES is the normal condition.

class Accelerator {
public:
	AcceleratorType type;
	int leafs_per_bank;
	Energy energy;
	Filter filter = Filter::FF;//default is WITH flattening filter, is overridden to FFF if encountered

	Accelerator(AcceleratorType _type){
		if (_type == AcceleratorType::MLCi80){
			leafs_per_bank = 40;
		}
		else if (_type == AcceleratorType::Agility){
			leafs_per_bank = 80;
		}
		else if (_type == AcceleratorType::MRLinac){
			leafs_per_bank = 80;
		}
		type = _type;
	}
	Accelerator(){ type = AcceleratorType::EMPTY; };
};

enum class BeamType { UNKNOWN, IMRT, VMAT };

struct BeamMetaData {
	//plan variables
	BeamType beamtype;
	float weight;
	float mu_per_fraction;// = 1.f;
	int nr_fractions;// = 1;
	Accelerator accelerator;
	float prescriptiondose;
	string isocentername;
	Float3 isoc;
	Float3 bfield;

	//ct metadata
	float hu_slope;
	float hu_intercept;
	float outsidepatientairthreshold;
	int outsidepatientisctnumber;
	string patient_position;
	float couchremovalycoordinate;
	float couch_height;
	
	//params
	bool dose_per_fraction;
	float fieldMargin;
	bool pinnacle_vmat_interpolation;
	bool table_type;
};


class RTBeam{ //base
public:
	//members
	BeamMetaData metaData;
	vector<ControlPoint> controlPoints;

	//methods
	void printInfo();
	void printFirstLeaf();
	int num_cps(){ return controlPoints.size(); };

	//ctors
	RTBeam() = default;
	RTBeam(const RTBeam &) = default;
	~RTBeam() = default;

	RTBeam(DosiaSettings &);

private:
	DosiaSettings sett;
};


class Parser{
protected:
	//members
	BeamMetaData metaData;
	vector<ControlPoint> controlPoints;

	//params
	bool debug;
	
	//methods
	virtual void parseTrial(const string &){}; //default implementations do nothing
	virtual void parsePlan(const string &){};
	virtual void parseBeam(const string &){};
	virtual void parseScan(const string &){};
	virtual void parseDose(const string &){};
	virtual void setCPIs();
	virtual int num_cps(){ return controlPoints.size(); };

	//ctors, which cant be used by anyone
	Parser() = default;
	Parser(const Parser &) = default;
	~Parser() = default;
	Parser(BeamMetaData, bool = false); //metaData, debug (metaData moet pinnacleVMATmode, dose_per_fraction, fieldMargin geset hebben)

	friend class RTBeam; //only RTBeam can touch or construct this class.
};


class pinnacleParser : virtual private Parser{
public:

private:
	pinnacleParser(BeamMetaData, bool = false);
	void parseTrial(const string &);
	void parsePlan(const string &);
	void parseBeam(const string &);
	//void parseScan(const string &){};//not needed
	void parseDose(const string &);//prescription dose
	//void setCPIs();

	//ctors

	friend class RTBeam; //only RTBeam can touch or construct this class.
};


class dicomParser : virtual private Parser{
public:

private:
	dicomParser(BeamMetaData, bool = false);
	//void parseTrial(const string &){};//doesnt exist
	//void parsePlan(const string &); //TODO NOT IMPLEMENTED, SUSPICION: NO RELEVANT INFO HERE
	void parseBeam(const string &); //TODO
	void parseScan(const string &); //intercept, slope
	//void parseDose(const string &){};//prescription dose
	void setCPIs();

	/*
	'NumberOfControlPoints',   //BEAM_NRCONTROLPOINTS
		'BeamLimitingDeviceAngle', //BEAM_COLLIMATOR
		'GantryAngle',             //BEAM_GANTRY,
		'PatientSupportAngle',     //BEAM_COUCH
		'NominalBeamEnergy[0]'     //BEAM_ENERGY
		*/

	//ctors

	friend class RTBeam; //only RTBeam can touch or construct this class.
};


Parser::Parser(BeamMetaData _metaData, bool _debug) : debug(_debug), metaData(_metaData){
}


pinnacleParser::pinnacleParser(BeamMetaData _metaData, bool _debug) : Parser(_metaData, _debug){
}


dicomParser::dicomParser(BeamMetaData _metaData, bool _debug) : Parser(_metaData, _debug){
}


void pinnacleParser::parseTrial(const string &inFile) {
	auto dump = load_dump(inFile);
	assert(!dump.empty());

	for (auto &line : dump) {
		/*if (startswith(line.first, "requestedmonitorunitsperfraction[")) {	//better take this from dose.dump
			metaData.mu_per_fraction = stof(line.second);
			continue;
		}*/
		if (startswith(line.first, "numberoffractions[")) {
			metaData.nr_fractions = stoi(line.second);
			continue;
		}
		if (startswith(line.first, "negativemupenalty")) {
			metaData.hu_intercept = -stoi(line.second); //NEGATIVE!!!
			metaData.hu_slope = 1.f; //There is no slope in pinnacle.
			continue;
		}
		if (startswith(line.first, "outsidepatientairthreshold")) {
			metaData.outsidepatientairthreshold = stof(line.second);
			continue;
		}
		if (startswith(line.first, "outsidepatientisctnumber")) {
			metaData.outsidepatientisctnumber = stoi(line.second); // this is pre-offset
			continue;
		}
		if (startswith(line.first, "patient_position")) {
			metaData.patient_position = line.second;
			continue;
		}
		if (startswith(line.first, "couchremovalycoordinate")) {
			metaData.couchremovalycoordinate = stof(line.second);
			continue;
		}
	}

	if (debug) fprintf(stderr,"RTPLan: dicom intercept, slope: %.2f,%.2f\n", metaData.hu_intercept, metaData.hu_slope);
	if (debug) fprintf(stderr,"RTPLan: mu_per_fraction: %.2f.\nnr_fractions: %i\n", metaData.mu_per_fraction, metaData.nr_fractions);
}

void pinnacleParser::parseDose(const string &inFile) {
	auto dump = load_dump(inFile);
	assert(!dump.empty());

	for (auto &line : dump) {
		if (startswith(line.first, "prescriptiondose")) {
			metaData.prescriptiondose = stof(line.second);
			continue;
		}
		if (startswith(line.first, "requestedmonitorunitsperfraction")) {
			metaData.mu_per_fraction = stof(line.second);
			continue;
		}
	}
}

void dicomParser::parseScan(const string &inFile) {
	auto dump = load_dump(inFile);
	assert(!dump.empty());

	for (auto &line : dump) {
		if (startswith(line.first, "RescaleIntercept")) {
			metaData.hu_intercept = stof(line.second);
			continue;
		}
		if (startswith(line.first, "RescaleSlope")) {
			metaData.hu_slope = stof(line.second);
			continue;
		}
		if (startswith(line.first, "ImagePositionPatient")) {
			/*
			auto pos = split_vec<float>(line.second, "\\");
			fprintf(stderr, "RTPLan: dicom ImagePositionPatient : %.2f,%.2f,%.2f\n", x,y,z);
			*/
			continue;
		}

		
	}
	if (debug) fprintf(stderr, "RTPLan: dicom intercept, slope: %.2f,%.2f\n", metaData.hu_intercept, metaData.hu_slope);
}


void pinnacleParser::parsePlan(const string &inFile) {
	auto dump = load_dump(inFile);
	assert(!dump.empty());

	assert(!metaData.isocentername.empty());

	for (auto &line : dump) {
		if (startswith(lower(line.first), lower(metaData.isocentername) + "_x")) { //normalize case for isoc name
			metaData.isoc.x = stof(line.second);
			continue;
		}
		if (startswith(lower(line.first), lower(metaData.isocentername) + "_y")) { //minus!
			metaData.isoc.y = -stof(line.second);
			continue;
		}
		if (startswith(lower(line.first), lower(metaData.isocentername) + "_z")) { //minus!
			metaData.isoc.z = -stof(line.second);
			continue;
		}
	}
	if (debug) fprintf(stderr,"RTPLan: isoc %.2f %.2f %.2f\n", metaData.isoc.x, metaData.isoc.y, metaData.isoc.z);
}

void dicomParser::parseBeam(const string &inFile) {
	auto dump = load_dump(inFile);
	assert(!dump.empty());

	vector<int> ReferencedBeamNumber;
	vector<float> BeamDose;
	vector<float> BeamMeterset;

	vector<string> BeamLimitingDeviceSequence(3); //should usually have 2 or 3.

	for (auto &line : dump) {
		if (startswith(line.first, "FractionGroupSequence")) { //here we look for the beam weight, which is not yet known here
			vector<string> line_breakdown = split(line.first, ".");
			if (line_breakdown.size() <= 1) continue; //one or less: nothing to split. ignore
			
			if (line_breakdown[1] == "NumberOfFractionsPlanned"){
				metaData.nr_fractions = stoi(line.second);
			}

			int blds = get_index(line_breakdown[1]); //centerpiece has beam index
			if (blds < 0) continue; // can be no number, then skip

			if (line_breakdown[2] == "BeamDose"){
				BeamDose.push_back(stof(line.second));
			}
			if (line_breakdown[2] == "BeamMeterset"){
				BeamMeterset.push_back(stof(line.second));
			}
			continue;
		}
		if (line.first == "TreatmentMachineName") {
			if (startswith(line.second, "MLC160")){
				metaData.accelerator = Accelerator(AcceleratorType::Agility);
				continue;
			}
			else if (startswith(line.second, "M160")){
				metaData.accelerator = Accelerator(AcceleratorType::Agility);
				continue;
			}
			else if (startswith(line.second, "MLC80")){
				metaData.accelerator = Accelerator(AcceleratorType::MLCi80);
				continue;
			}
			else if (startswith(line.second, "M80")){
				metaData.accelerator = Accelerator(AcceleratorType::MLCi80);
				continue;
			}
			else {
				metaData.accelerator = Accelerator(AcceleratorType::UNKNOWN);
				continue;
			}
		}
		if (startswith(line.first, "BeamLimitingDeviceSequence")) {
			vector<string> line_breakdown = split(line.first, ".");
			if (line_breakdown.size() <= 1) continue; //one or less: nothing to split. ignore

			int blds = get_index(line_breakdown[0]);
			if (blds < 0) continue; // can be no number, then skip
			if (blds > 2) throw std::pair<int, string>(24, "Unexpected number of BeamLimitingDevices.");

			if (line_breakdown[1] == "RTBeamLimitingDeviceType"){
				BeamLimitingDeviceSequence[blds] = line.second;
			}
			if (line_breakdown[1] == "NumberOfLeafJawPairs"){
				if (BeamLimitingDeviceSequence[blds] == "MLCX"){
					if (metaData.accelerator.leafs_per_bank != stoi(line.second)) throw std::pair<int, string>(25, "Number of Leafs does not correspond to specified Accelerator.");
				}
			}
		}
		if (line.first == "BeamNumber") {
			assert(len(BeamDose) == len(BeamMeterset) == len(ReferencedBeamNumber));
			assert(len(ReferencedBeamNumber) > 0);
			//metaData.weight = BeamDose[stoi(line.first)]; // TODO: not sure which
			metaData.weight = BeamMeterset[index(ReferencedBeamNumber, stoi(line.first))];
			continue;
		}
		if (line.first == "NumberOfControlPoints") {
			controlPoints.resize(stoi(line.second));
			for (int i = 0; i < stoi(line.second); i++){
				controlPoints[i].collimator.mlc.leftLeaves.resize(metaData.accelerator.leafs_per_bank); // leaf_per_bank should be set before
				controlPoints[i].collimator.mlc.rightLeaves.resize(metaData.accelerator.leafs_per_bank);
			}
			continue;
		}
		if (startswith(line.first, "BeamType")) {
			if (startswith(line.second, "DYNAMIC")){
				metaData.beamtype = BeamType::VMAT;
			}
			//following copied from pinnacle::parseBeam
			else if (startswith(line.second, "Dynamic Arc")){
				metaData.beamtype = BeamType::VMAT;
			}
			else if (startswith(line.second, "Step & Shoot MLC")){
				metaData.beamtype = BeamType::IMRT;
			}
			else if (startswith(line.second, "Static")){
				metaData.beamtype = BeamType::IMRT;
			}
			else {
				metaData.beamtype = BeamType::UNKNOWN;
			}
			continue;
		}
		if (startswith(line.first, "PatientSetupSequence")) {
			vector<string> line_breakdown = split(line.first, ".");
			if (line_breakdown[1] == "PatientPosition"){
				metaData.patient_position = line.second;
			}
			continue;
		}
		if (startswith(line.first, "ControlPointSequence")) {
			vector<string> line_breakdown = split(line.first, ".");
			if (line_breakdown.size() <= 1) continue; //one or less: nothing to split. ignore

			//get CPI
			int cpi = get_index(line_breakdown[0]);
			if (cpi < 0) continue; // can be no number, then skip
			if (cpi > num_cps() - 1) throw std::pair<int, string>(26, "CPI out of range.");

			if (line_breakdown[1] == "ControlPointIndex") {
				assert(cpi == stoi(line.second));
				continue;
			}
			if (line_breakdown[1] == "NominalBeamEnergy") {
				if (startswith(line.second, "6")){
					metaData.accelerator.energy = Energy::MV6;
				}
				else if (startswith(line.second, "10")){
					metaData.accelerator.energy = Energy::MV10;
				}
				else if (startswith(line.second, "7")){ //cant happen I guess
					metaData.accelerator.energy = Energy::MV7;
				}
				else {
					metaData.accelerator.energy = Energy::UNKNOWN;
				}
				continue;
			}
			if (line_breakdown[1] == "NumberOfCompensators") { //TODO
				/*if (endswith(line.second, "FFF")){
				metaData.accelerator.filter = Filter::NoFF;
				//TODO: imrtfilter="Compensator" == FFF?
				}*/
				continue;
			}
			if (line_breakdown[1] == "GantryAngle") {
				controlPoints[cpi].beamInfo.gantryAngle = { stof(line.second), stof(line.second) };
				continue;
			}
			if (line_breakdown[1] == "PatientSupportAngle") {
				controlPoints[cpi].beamInfo.couchAngle = { stof(line.second), stof(line.second) };
				continue;
			}
			if (line_breakdown[1] == "BeamLimitingDeviceAngle") {
				controlPoints[cpi].beamInfo.collimatorAngle = { stof(line.second), stof(line.second) };
				continue;
			}
			if (line_breakdown[1] == "IsocenterPosition") {
				auto pos = types::split<float>(line.second, "\\");
				controlPoints[cpi].beamInfo.isoCenter.x = pos[0];
				controlPoints[cpi].beamInfo.isoCenter.y = pos[1];
				controlPoints[cpi].beamInfo.isoCenter.z = pos[2];
				continue;
			}
			if (line_breakdown[1] == "CumulativeMetersetWeight") {
				if (cpi == 0) {
					controlPoints[cpi].beamInfo.relativeWeight = stof(line.second);
				}
				else {
					controlPoints[cpi].beamInfo.relativeWeight = stof(line.second) - controlPoints[cpi - 1].beamInfo.relativeWeight;
				}
				continue;
			}
			if (startswith(line_breakdown[1], "BeamLimitingDevicePositionSequence")) {
				int bldi = get_index(line_breakdown[1]);
				if (bldi < 0) continue; // can be no number, then skip

				if (line_breakdown[2] == "LeafJawPositions") {
					auto pos = types::split<float>(line.second,"\\");
					if (BeamLimitingDeviceSequence[bldi] == "ASMX"){
						controlPoints[get_index(line.first)].collimator.parallelJaw.j1 = { pos[0], pos[0] }; //J1 always most negative coord according to doc.
						controlPoints[get_index(line.first)].beamInfo.fieldMin.first = pos[0] - metaData.fieldMargin;

						controlPoints[get_index(line.first)].collimator.parallelJaw.j2 = { pos[1], pos[1] };
						controlPoints[get_index(line.first)].beamInfo.fieldMax.first = pos[1] + metaData.fieldMargin;
					}
					if (BeamLimitingDeviceSequence[bldi] == "ASMY"){
						controlPoints[get_index(line.first)].collimator.perpendicularJaw.j1 = { pos[0], pos[0] };
						controlPoints[get_index(line.first)].beamInfo.fieldMin.second = pos[0] - metaData.fieldMargin;

						controlPoints[get_index(line.first)].collimator.perpendicularJaw.j2 = { pos[1], pos[1] };
						controlPoints[get_index(line.first)].beamInfo.fieldMax.second = pos[1] + metaData.fieldMargin;
					}
					if (BeamLimitingDeviceSequence[bldi] == "MLCX"){
						//pos heeft right bank eerst
						for (int p = 0; p < metaData.accelerator.leafs_per_bank; p++){
							controlPoints[cpi].collimator.mlc.rightLeaves[p] = { pos[p], pos[p] };
							controlPoints[cpi].collimator.mlc.leftLeaves[p] = { pos[p + metaData.accelerator.leafs_per_bank], pos[p + metaData.accelerator.leafs_per_bank] };
						}
						//NOOT: for square fields dicom does not necesarily give MLCX positions. Only jaw positions may be encountered!!
					}
				}
				continue;
			}
			continue;
		}
	}
}


void pinnacleParser::parseBeam(const string &inFile) {
	auto dump = load_dump(inFile);
	assert(!dump.empty());

	int skip_n_lines = 0;
	double total_weight = 0.;

	for (auto &line : dump) {
		if (skip_n_lines > 0){
			//it was set, so decrement and skip.
			skip_n_lines--;
			if (debug) fprintf(stderr,"RTPLan: skipping line: %s\n", line.second);
			continue;
		}
		if (line.first == "isocentername"){
			metaData.isocentername = line.second;
			continue;
		}
		if (line.first == "machinenameandversion") {
			if (startswith(line.second, "MLC160")){
				metaData.accelerator = Accelerator(AcceleratorType::Agility);
				continue;
			}
			else if (startswith(line.second, "M160")){
				metaData.accelerator = Accelerator(AcceleratorType::Agility); //FFF?
				continue;
			}
			else if (startswith(line.second, "MLC80")){
				metaData.accelerator = Accelerator(AcceleratorType::MLCi80);
				continue;
			}
			else if (startswith(line.second, "M80")){
				metaData.accelerator = Accelerator(AcceleratorType::MLCi80);
				continue;
			}
			else {
				metaData.accelerator = Accelerator(AcceleratorType::UNKNOWN);
				continue;
			}
		}
		if (line.first == "machineenergyname") {
			if (startswith(line.second, "6")){
				metaData.accelerator.energy = Energy::MV6;
			}
			else if (startswith(line.second, "10")){
				metaData.accelerator.energy = Energy::MV10;
			}
			else if (startswith(line.second, "7")){ //cant happen I guess
				metaData.accelerator.energy = Energy::MV7;
			}
			else {
				metaData.accelerator.energy = Energy::UNKNOWN;
				metaData.accelerator.filter = Filter::UNKNOWN;
			}
			if (endswith(line.second, "FFF")){
				metaData.accelerator.filter = Filter::NoFF;
				//TODO: imrtfilter="Compensator" == FFF?
			}
			continue;
		}
		if (line.first == "numberofcontrolpoints") {
			controlPoints.resize(stoi(line.second));
			for (int i = 0; i < stoi(line.second); i++){
				controlPoints[i].collimator.mlc.leftLeaves.resize(metaData.accelerator.leafs_per_bank); // leaf_per_bank should be set before
				controlPoints[i].collimator.mlc.rightLeaves.resize(metaData.accelerator.leafs_per_bank);
			}
			continue;
		}
		if (startswith(line.first, "gantry[")) {
			controlPoints[get_index(line.first)].beamInfo.gantryAngle = { stof(line.second), stof(line.second) };
			continue;
		}
		if (startswith(line.first, "couch[")) {
			controlPoints[get_index(line.first)].beamInfo.couchAngle = { 360-stof(line.second), 360-stof(line.second) };
			// we doen 360 - waarde, want Pinnacle...
			continue;
		}
		if (startswith(line.first, "collimator[")) {
			controlPoints[get_index(line.first)].beamInfo.collimatorAngle = { stof(line.second), stof(line.second) };
			continue;
		}
		if (startswith(line.first, "setbeamtype")) {
			if (startswith(line.second, "Dynamic Arc")){
				metaData.beamtype = BeamType::VMAT;
			}
			else if (startswith(line.second, "Step & Shoot MLC")){
				metaData.beamtype = BeamType::IMRT;
			}
			else if (startswith(line.second, "Static")){
				metaData.beamtype = BeamType::IMRT;
			}
			else {
				metaData.beamtype = BeamType::UNKNOWN;
			}
			continue;
		}
		if (startswith(line.first, "numberofpoints[")) {
			//use as check.
			if (stoi(line.second) == metaData.accelerator.leafs_per_bank){
				continue; //all ok
			}
			else {
				fprintf(stderr,"RTPLan: Invalid CPI detected: incorrect nr leafs: %s. Skipping CPI...\n", line.second);
				skip_n_lines = stoi(line.second) * 2; //skip this nr*2 (pairs) of lines
				continue;
			}
		}
		if (startswith(line.first, "points_element[")) {
			//loop backwards!
			int cpi = get_index(line.first);
			int leafindex = get_index(line.first, 2);
			float pos = stof(line.second);
			//fprintf(stderr,"points element %i %i %f\n", cpi, leafindex, pos);
			if (leafindex % 2 == 0){ //even, dus leftbank
				//std::cout << "left[" << cpi << "][" << leafindex / 2 << "] " << -pos;
				controlPoints[cpi].collimator.mlc.leftLeaves[metaData.accelerator.leafs_per_bank - 1 - (leafindex / 2)] = { -pos, -pos }; //minus want vanaf het midden!
			}
			else{
				//std::cout << "right[" << cpi << "][" << (leafindex - 1) / 2 << "] " << pos;
				controlPoints[cpi].collimator.mlc.rightLeaves[metaData.accelerator.leafs_per_bank - 1 - ((leafindex - 1) / 2)] = { pos, pos };
			}
			continue;
		}
		if (startswith(line.first, "leftjawposition[")) {
			controlPoints[get_index(line.first)].collimator.parallelJaw.j1 = { -stof(line.second), -stof(line.second) }; //J1 always most negative coord according to doc.
			controlPoints[get_index(line.first)].beamInfo.fieldMin.first = -stof(line.second) - metaData.fieldMargin; // pinnacle gives distance to center of colli, so add minus
			continue;
		}
		if (startswith(line.first, "rightjawposition[")) {
			controlPoints[get_index(line.first)].collimator.parallelJaw.j2 = { stof(line.second), stof(line.second) };
			controlPoints[get_index(line.first)].beamInfo.fieldMax.first = stof(line.second) + metaData.fieldMargin;
			continue;
		}
		if (startswith(line.first, "topjawposition[")) { //TRF reader: y axis swapped ? increases downward
			controlPoints[get_index(line.first)].collimator.perpendicularJaw.j2 = { stof(line.second), stof(line.second) }; //J1 always most negative coord according to doc.
			controlPoints[get_index(line.first)].beamInfo.fieldMax.second = stof(line.second) + metaData.fieldMargin;
			continue;
		}
		if (startswith(line.first, "bottomjawposition[")) {
			controlPoints[get_index(line.first)].collimator.perpendicularJaw.j1 = { -stof(line.second), -stof(line.second) };
			controlPoints[get_index(line.first)].beamInfo.fieldMin.second = -stof(line.second) - metaData.fieldMargin;
			continue;
		}
		if (startswith(line.first, "weight[")) { //weight sums to one per beam.
			controlPoints[get_index(line.first)].beamInfo.relativeWeight = stof(line.second); //relweight is relative to computation, so total should be 1, and per beam it is.
			total_weight += stof(line.second);
			continue;
		}
		if (startswith(line.first, "weight") && !startswith(line.first, "weight[")) { // beam weight, one per beam. summed over beams should be one.
			metaData.weight = stof(line.second);
			continue;
		}
	}

	if (debug) fprintf(stderr,"RTPLan: One beam loaded with a beam weight of %f, and cps weights sum to %f.\n", metaData.weight, total_weight);
}


void Parser::setCPIs() {
	if (metaData.accelerator.type == AcceleratorType::EMPTY || metaData.accelerator.type == AcceleratorType::UNKNOWN) throw std::pair<int, string>(27, "Undefined accelerator encountered.");
	if (metaData.accelerator.energy == Energy::UNKNOWN) throw std::pair<int, string>(21, "Unknown energy encountered.");
	if (metaData.accelerator.filter == Filter::UNKNOWN) throw std::pair<int, string>(23, "Unknown filter encountered.");
	if (metaData.beamtype == BeamType::UNKNOWN) throw std::pair<int, string>(22, "Unknown beamtype encountered.");

	metaData.prescriptiondose *= (metaData.dose_per_fraction ? 1 : metaData.nr_fractions);

	double total_weight = 0.f;
	// zet orientaties, isoc, correct weights
	for (auto &cp : controlPoints){
		cp.beamInfo.isoCenter.x = metaData.isoc.x;
		cp.beamInfo.isoCenter.y = metaData.isoc.y;
		cp.beamInfo.isoCenter.z = metaData.isoc.z;
		cp.collimator.mlc.orientation = ModifierOrientation::IECY;
		cp.collimator.perpendicularJaw.orientation = ModifierOrientation::IECX;

		if (metaData.accelerator.type == AcceleratorType::MLCi80){
			cp.collimator.parallelJaw.orientation = ModifierOrientation::IECY;
		}
		else {
			cp.collimator.parallelJaw.orientation = ModifierOrientation::NOT_PRESENT;
		}
		total_weight += cp.beamInfo.relativeWeight;
		//(beam)weight given in perc, to factor (/100), and then from Gy to cGy (*100)
		cp.beamInfo.relativeWeight *= metaData.weight * metaData.mu_per_fraction * (metaData.dose_per_fraction ? 1 : metaData.nr_fractions);
		//								^-beam weight		^- mu/frac						^- times nb frac or not.
	}
	assert((total_weight >= 0.9999) && (total_weight <= 1.0001)); //cp weigths within beam must sum to 1. only machine error is tolerated.

	//VMAT to dynamic arcs
	if (metaData.beamtype == BeamType::VMAT && !metaData.pinnacle_vmat_interpolation){
		for (int i = 0; i < (num_cps() - 1); i++){ //skip last CPI
			// van N CPIs naar N-1 segmenten.
			// kijk ACHTERUIT: skippen dus N=0. Want, laatste CP heeft weight nul
			// Let op: segment fieldMin,fieldMax moeten omhullende van naastgelegen CPIs worden

			controlPoints[i].beamInfo.collimatorAngle.second = controlPoints[i + 1].beamInfo.collimatorAngle.first; //first is already correct
			controlPoints[i].beamInfo.couchAngle.second = controlPoints[i + 1].beamInfo.couchAngle.first;
			controlPoints[i].beamInfo.gantryAngle.second = controlPoints[i + 1].beamInfo.gantryAngle.first;

			//fieldMin/Max is static!!!, take most outward of the two CPIs.
			controlPoints[i].beamInfo.fieldMin.first = std::min(controlPoints[i].beamInfo.fieldMin.first, controlPoints[i + 1].beamInfo.fieldMin.first);
			controlPoints[i].beamInfo.fieldMin.second = std::min(controlPoints[i].beamInfo.fieldMin.second, controlPoints[i + 1].beamInfo.fieldMin.second);
			controlPoints[i].beamInfo.fieldMax.first = std::max(controlPoints[i].beamInfo.fieldMax.first, controlPoints[i + 1].beamInfo.fieldMax.first);
			controlPoints[i].beamInfo.fieldMax.second = std::max(controlPoints[i].beamInfo.fieldMax.second, controlPoints[i + 1].beamInfo.fieldMax.second);

			//j1 == negative, so minimum
			controlPoints[i].collimator.perpendicularJaw.j1.second = controlPoints[i + 1].collimator.perpendicularJaw.j1.first;
			controlPoints[i].collimator.perpendicularJaw.j2.second = controlPoints[i + 1].collimator.perpendicularJaw.j2.first;
			controlPoints[i].collimator.parallelJaw.j1.second = controlPoints[i + 1].collimator.parallelJaw.j1.first;
			controlPoints[i].collimator.parallelJaw.j2.second = controlPoints[i + 1].collimator.parallelJaw.j2.first;

			for (int j = 0; j < metaData.accelerator.leafs_per_bank; j++){
				controlPoints[i].collimator.mlc.leftLeaves[j].second = controlPoints[i + 1].collimator.mlc.leftLeaves[j].first;
				controlPoints[i].collimator.mlc.rightLeaves[j].second = controlPoints[i].collimator.mlc.rightLeaves[j].first;
			}
		}
		controlPoints.pop_back(); //we moved every CP forward, so last one is now superfluous
	}
	//VMAT pinnacle interpol. no dynamic thingies, but only reweighting
	else if (metaData.beamtype == BeamType::VMAT && metaData.pinnacle_vmat_interpolation){

		for (int i = 0; i < (num_cps() - 1); i++){ // loop over CPIs except last
			// we take average weight of i and i + 1, because of the pinnacle shifted weights wrt to dicom and tabel 5 in NLCOMvStralingsdosimetrie Rapport 26
			// i=1 has weight zero to we get half weight of next CP anyway.
			controlPoints[i].beamInfo.relativeWeight = controlPoints[i].beamInfo.relativeWeight / 2. + controlPoints[i + 1].beamInfo.relativeWeight / 2.;
		}

		controlPoints.pop_back(); //weights have shifted forward, laatste mag weg
	}

}


void dicomParser::setCPIs() {
	//first, convert CUMULATIVE mu to mu.
	// kijk vooruit: skippen dus N=N. Want, eerste CP heeft weight 0
	// skip last, because lookahead on next line
	for (int i = 0; i < num_cps() - 1; i++){
		controlPoints[i].beamInfo.relativeWeight = controlPoints[i + 1].beamInfo.relativeWeight - controlPoints[i].beamInfo.relativeWeight; //i+1 - i
	}
	controlPoints[num_cps() - 1].beamInfo.relativeWeight = 0; //last one can be zeroed out.
	Parser::setCPIs(); //then regular procedure
}


//RTBeam::RTBeam(const string &rt_files, float _fieldMargin, bool _debug, bool _pinnacleVMATmode) {
RTBeam::RTBeam(DosiaSettings &_sett) : sett(_sett){
	string &rt_files = sett.rt_files;
	io::isfile(rt_files + "/dbtype.dump", 20);
	auto dbtype = load_dump(sett.rt_files + "/dbtype.dump");
	for (auto &line : dbtype) {
		if (line.first == "pinnacle"){
			auto parser = pinnacleParser(metaData, (sett.verbose > 1)?true:false);
			parser.parseTrial(rt_files + "/trialname.dump");	//mu, nr_fracties, HU-intercept. NIET beschikbaar bij dicom
			parser.parseBeam(rt_files + "/beam.dump");
			parser.parsePlan(rt_files + "/plan.dump");	//isoc, requires info from beam!
			parser.parseScan(rt_files + "/scan.dump");	//not needed
			parser.parseDose(rt_files + "/dose.dump");

			parser.setCPIs();
			metaData = parser.metaData;
			controlPoints = parser.controlPoints;
		}
		else if (line.first == "dicom"){
			throw std::pair<int, string>(28, "Dicom not fully implemented.");
			auto parser = dicomParser(metaData, (sett.verbose > 1) ? true : false);
			parser.parseTrial(rt_files + "/trialname.dump");	//NIET beschikbaar bij dicom
			parser.parsePlan(rt_files + "/plan.dump");	//no required info?
			parser.parseBeam(rt_files + "/beam.dump");
			parser.parseScan(rt_files + "/scan.dump");	//dicom slope/intercept, ImagePositionPatient
			parser.parseDose(rt_files + "/dose.dump"); //not needed

			parser.setCPIs();
			metaData = parser.metaData;
			controlPoints = parser.controlPoints;
		}
	}

	if (sett.verbose > 2) {
		printInfo();
		printFirstLeaf();
	}
};


void RTBeam::printInfo(){
	fprintf(stderr,"RTPLan: Number of fractions is %i, %.2f MU per fraction and a prescription dose of %.2f.\n", metaData.nr_fractions, metaData.mu_per_fraction, metaData.prescriptiondose);
	fprintf(stderr,"RTPLan: This beam of weight %.2f has %i controlpoints. Per CPI:\n", metaData.weight, num_cps());
	for (auto &cp : controlPoints){
		fprintf(stderr,"fieldMin/Max                 \t %.2f \t %.2f \t %.2f \t %.2f \n", cp.beamInfo.fieldMin.first, cp.beamInfo.fieldMin.second, cp.beamInfo.fieldMax.first, cp.beamInfo.fieldMax.second);
		fprintf(stderr,"parallelJaw.j1/j2/orient     \t %.2f \t %.2f \t %i \n", cp.collimator.parallelJaw.j1.first, cp.collimator.parallelJaw.j2.first, cp.collimator.parallelJaw.orientation);
		fprintf(stderr,"perpendicularJaw.j1/j2/orient\t %.2f \t %.2f \t %i \n", cp.collimator.perpendicularJaw.j1.first, cp.collimator.perpendicularJaw.j2.first, cp.collimator.perpendicularJaw.orientation);
		fprintf(stderr,"gantryangle %.2f\n", cp.beamInfo.gantryAngle.first);
		fprintf(stderr,"relativeWeight (in MU) %.2f\n", cp.beamInfo.relativeWeight);
	}
}


void RTBeam::printFirstLeaf(){
	fprintf(stderr,"RTPLan: This beam has %i controlpoints. First CPI:\n", num_cps());
	for (int i = 0; i < controlPoints[0].collimator.mlc.leftLeaves.size(); i++){
		fprintf(stderr,"%i th leafbank (l, r): %.2f, %.2f\n", i, controlPoints[0].collimator.mlc.leftLeaves[i].first, controlPoints[0].collimator.mlc.rightLeaves[i].first);
	}
}

A  => include/settings.h +117 -0
@@ 1,117 @@
#pragma once

#include "INIreader.h"
#include "gpumcd/Settings.h"

class DosiaSettings{
public:
	//from cmd args
	string rt_files;

	//from ini
	string gpumcd_material_data_dir;
	string gpumcd_machine_dir;
	string hounsfield_conversion_dir;

	int verbose;
	bool dbgoutput;
	
	float field_margin;
	bool dose_per_fraction;
	bool continous_materials;
	bool pinnacle_vmat_interpolation;
	bool monte_carlo_high_precision;
	bool score_dose_to_water;
	bool score_and_transport_in_water;
	bool in_aqua_vivo;
	
	bool gamma_comparison;
	bool gamma_global_dose;
	float gamma_isodose_region;
	float gamma_dd;
	float gamma_dta;

	string MRLinac_MV7;
	string Agility_MV6_FF;
	string Agility_MV6_NoFF;
	string Agility_MV10_FF;
	string Agility_MV10_NoFF;

	PhysicsSettings physicsSettings;
	PlanSettings planSettings;

	DosiaSettings() = default;
	DosiaSettings(const string &, const string &);
	DosiaSettings(const DosiaSettings &) = default;
	~DosiaSettings() = default;
};

DosiaSettings::DosiaSettings(const string &ini_file, const string & _rt_files) : rt_files(_rt_files){
	INIReader ini(ini_file);
	if (ini.ParseError() < 0) {
		throw std::pair<int, string>(1, "Can't load " + ini_file + "\n");
	}
	
	gpumcd_material_data_dir = ini.Get("directories", "gpumcd_material_data_dir", ".");
	gpumcd_machine_dir = ini.Get("directories", "gpumcd_machine_dir", ".");
	hounsfield_conversion_dir = ini.Get("directories", "hounsfield_conversion_dir", ".");

	MRLinac_MV7 = ini.Get("gpumcd_machines", "MRLinac_MV7", ".");
	Agility_MV6_FF = ini.Get("gpumcd_machines", "Agility_MV6_FF", ".");
	Agility_MV6_NoFF = ini.Get("gpumcd_machines", "Agility_MV6_NoFF", ".");
	Agility_MV10_FF = ini.Get("gpumcd_machines", "Agility_MV10_FF", ".");
	Agility_MV10_NoFF = ini.Get("gpumcd_machines", "Agility_MV10_NoFF", ".");

	verbose = ini.GetInteger("debug", "verbose", 0);
	dbgoutput = ini.GetBoolean("debug", "output", false);

	field_margin = ini.GetReal("dose", "field_margin", 5.f);
	dose_per_fraction = ini.GetBoolean("dose", "dose_per_fraction", true);
	continous_materials = ini.GetBoolean("dose", "continous_materials", true);
	pinnacle_vmat_interpolation = ini.GetBoolean("dose", "pinnacle_vmat_interpolation", false);
	monte_carlo_high_precision = ini.GetBoolean("dose", "monte_carlo_high_precision", false);
	score_dose_to_water = ini.GetBoolean("dose", "score_dose_to_water", false);
	score_and_transport_in_water = ini.GetBoolean("dose", "score_and_transport_in_water", false);
	in_aqua_vivo = ini.GetBoolean("dose", "in_aqua_vivo", false);

	gamma_comparison = ini.GetBoolean("gamma", "comparison", false);
	gamma_global_dose = ini.GetBoolean("gamma", "global_dose", true);
	gamma_isodose_region = ini.GetReal("gamma", "isodose_region", 10);
	gamma_dd = ini.GetReal("gamma", "dd", 3);
	gamma_dta = ini.GetReal("gamma", "dta", 2);
	
	// Load settings gpumcd physics, defaults from gpumcd/Settings.h
	physicsSettings.photonTransportCutoff = ini.GetReal("gpumcd_physicssettings", "photonTransportCutoff", 0.01f);
	physicsSettings.electronTransportCutoff = ini.GetReal("gpumcd_physicssettings", "electronTransportCutoff", 0.189f);
	physicsSettings.inputMaxStepLength = ini.GetReal("gpumcd_physicssettings", "inputMaxStepLength", 0.75f);
	physicsSettings.referenceMedium = ini.GetInteger("gpumcd_physicssettings", "referenceMedium", -1);
	physicsSettings.useElectronInAirSpeedup = ini.GetBoolean("gpumcd_physicssettings", "useElectronInAirSpeedup", true);
	physicsSettings.electronInAirSpeedupDensityThreshold = ini.GetReal("gpumcd_physicssettings", "electronInAirSpeedupDensityThreshold", 0.002f);

	// Load settings gpumcd plan, NOT ALL defaults from gpumcd/Settings.h, parts may be overridden by the RTBeam
	planSettings.goalSfom = ini.GetReal("gpumcd_plansettings", "goalSfom", 54.32f); //54.32f is a flag value, will be used to identify UNSET
	planSettings.statThreshold = ini.GetReal("gpumcd_plansettings", "statThreshold", 0.5f);
	planSettings.maxNumParticles = static_cast<uint64_t>(ini.GetReal("gpumcd_plansettings", "maxNumParticles", 1e10)); //GetReal so that we can cope with scientific notation
	planSettings.densityThresholdSfom = ini.GetReal("gpumcd_plansettings", "densityThresholdSfom", 0.2f);
	planSettings.densityThresholdOutput = ini.GetReal("gpumcd_plansettings", "densityThresholdOutput", 0.f);
	planSettings.useApproximateStatistics = ini.GetBoolean("gpumcd_plansettings", "useApproximateStatistics", true); // is AUTOMATICALLY set to FALSE for SEGMENT calcs.

	// Done. Print settings for log?
	if (verbose > 0){
		cerr << "Verbose level set at " << verbose << ".\n";

		cerr << "field_margin = " << field_margin << ".\n";
		cerr << "dose_per_fraction = " << dose_per_fraction << ".\n";
		cerr << "pinnacle_vmat_interpolation = " << pinnacle_vmat_interpolation << ".\n";
		cerr << "monte_carlo_high_precision = " << monte_carlo_high_precision << ".\n";

		if (gamma_comparison) cerr << "Gamma comparison enabled.\n";
		if (dbgoutput) cerr << "Debug outputs will be written to disk.\n";
		//if (in_aqua_vivo) cerr << "Forcing all densities inside patient threshold to 1.0g/cm3 (as EpidTrial.py in Pinnacle).\n";
		if (in_aqua_vivo) cerr << "in_aqua_vivo currently not correctly implemented. Will be removed. Dosia dump should fix this.\n";
		if (in_aqua_vivo) in_aqua_vivo = false;

		if (score_and_transport_in_water) { cerr << "Computing dose and transport in water instead of medium.\n"; }
		else if (score_dose_to_water) { cerr << "Computing dose to water instead of medium.\n"; };
	}
};
\ No newline at end of file

A  => include/tools.h +361 -0
@@ 1,361 @@
#pragma once

#include <iostream>
#include <vector>
#include <algorithm>    // std::min_element, std::max_element, std::sort, std::reverse, std::transform
#include <functional> //std::plus, std::minus
#include <numeric>	// std::iota
#include <random>
#include <fstream>
#include <assert.h>
#include <sstream> //stringstream
#include <cstring> //std::memcpy
#include "pystring.h"
#pragma warning(disable : 4996) // disable fopen warning vs

namespace vect {
	using std::vector;

	template <typename T, typename U, typename V>
	double interpolate(const vector<T> &xData, const vector<U> &yData, const V &x, const bool extrapolate) {
		//Can call with any combo of input types, but result is always double.
		//extrapolate sets gradient outside xData to zero or continue from nearest bin.
		size_t size = xData.size();

		size_t i = 0;// find left end of interval for interpolation
		if (x >= xData[size - 2]) {// special case: beyond right end
			i = size - 2;
		}
		else {
			while (x > xData[i + 1]) i++;
		}
		double xL = xData[i], yL = yData[i], xR = xData[i + 1], yR = yData[i + 1];
		if (!extrapolate) {	// if beyond ends of array and not extrapolating, set gradient to zero
			if (x < xL) yR = yL;
			if (x > xR) yL = yR;
		}
		double gradient = (yR - yL) / (xR - xL);

		return yL + gradient * (x - xL); // linear interpolation
	};

	template <typename T>
	vector<size_t> sort_indexes(const vector<T> &v) {
		// src: https://stackoverflow.com/questions/1577475/c-sorting-and-keeping-track-of-indexes#12399290
		// as per std::sort convention, low to high. use std::reverse

		// initialize original index locations
		vector<size_t> idx(v.size());
		std::iota(idx.begin(), idx.end(), 0);

		// sort indexes based on comparing values in v
		std::sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) {return v[i1] < v[i2]; });

		return idx;
	}

	vector<double> uniform_random_vector(int length = 100) {
		std::random_device rng;
		std::uniform_real_distribution<double> distribution(0.0, 1.0);
		vector<double> result(length);
		for (auto& item : result) item = distribution(rng);
		return result;
	}

	template <typename T>
	T sum(vector<T> &v){
		T sum = 0;
		for (const auto &i : v) {
			sum += i;
		}
		return sum;
	}

	template <typename T>
	void reverse(vector<T> &v){ //REFERENCE, otherwise its copied!!! it reverses IN PLACE.
		std::reverse(v.begin(), v.end());
	}

	//element wise add two vectors
	template <typename T>
	vector<T> add(const vector<T> &v1, const vector<T> &v2){
		assert(v1.size() == v2.size());
		vector<T> v3(v1.size());
		std::transform(v1.begin(), v1.end(), v2.begin(), v3.begin(), std::plus<T>());
		return v3;
	}

	//element wise subtract two vectors
	template <typename T>
	vector<T> sub(const vector<T> &v1, const vector<T> &v2){
		assert(v1.size() == v2.size());
		vector<T> v3(v1.size());
		std::transform(v1.begin(), v1.end(), v2.begin(), v3.begin(), std::minus<T>());
		return v3;
	}

	//element wise multiply two vectors
	template <typename T>
	vector<T> mul(const vector<T> &v1, const vector<T> &v2){
		assert(v1.size() == v2.size());
		vector<T> v3(v1.size());
		std::transform(v1.begin(), v1.end(), v2.begin(), v3.begin(), std::multiplies<T>());
		return v3;
	}

	//multiply elements of vector
	template <typename T>
	T mul(const vector<T> &v1){
		T retval = 1;
		for (const auto &n : v1) {
			retval *= n;
		}
		return retval;
	}

	//element wise divide two vectors
	template <typename T>
	vector<T> div(const vector<T> &v1, const vector<T> &v2){
		assert(v1.size() == v2.size());
		vector<T> v3(v1.size());
		std::transform(v1.begin(), v1.end(), v2.begin(), v3.begin(), std::divides<T>());
		return v3;
	}

	//write vector to binary file
	template <typename T>
	void tofile(const vector<T> &v, const std::string &fn){
		FILE* ffile = fopen(fn.c_str(), "wb");
		if (ffile == nullptr){
			throw std::pair<int, std::string>(70, "Problem writing file '" + fn + "'.");
			//throw "Problem writing file '" + fn + "'.";
		}
		fwrite(reinterpret_cast<const char*>(&v[0]), sizeof(T), v.size(), ffile);
		fclose(ffile);
	}

	//template specialization for vector of std::strings
	template <>
	void tofile<std::string>(const vector<std::string> &v, const std::string &fn){
		FILE* ffile = fopen(fn.c_str(), "w");
		if (ffile == nullptr){
			throw std::pair<int, std::string>(71, "Problem writing file '" + fn + "'.");
			//throw "Problem writing file '" + fn + "'.";
		}
		for (const auto &line : v) {
			fprintf(ffile, "%s\n",line.c_str());
		}
		fclose(ffile);
	}

	//read vector from binary file, unknown size
	template <typename T>
	vector<T> fromfile(const std::string &fn){
		FILE* ffile = fopen(fn.c_str(), "rb");
		if (ffile == nullptr){
			throw std::pair<int, std::string>(72, "Problem reading file '" + fn + "'.");
		}

		//get filesize
		fseek(ffile, 0, SEEK_END);
		long fsize = ftell(ffile) / sizeof(T);
		rewind(ffile);

		vector<T> v(fsize);

		fread(&v[0], sizeof(T), v.size(), ffile);
		fclose(ffile);
		return v;
	}

	//fromfile specialization for textfile
	template <>
	vector<std::string> fromfile(const std::string &fn) {
		std::ifstream f(fn);// .c_str());
		if (!f.good()) throw std::pair<int, std::string>(73, "Problem reading file '" + fn + "'.");
		vector<std::string> ret;
		std::string str;
		while (std::getline(f, str)) {
			ret.push_back(str);
		}
		return ret;
	}

	//len()
	template <typename T>
	size_t len(const vector<T> &source) {
		return source.size();
	}

	//index()
	template <typename T>
	size_t index(const vector<T> &source, const T &item) {
		size_t index = std::distance(source.begin(), std::find(std::begin(source), std::end(source), item));
		if (index == len(source)) throw std::pair<int, std::string>(74, "Element not found");
		return index;
	}

	//printers voor (nested) vectors van any type
	void print(const vector<std::string> &v, std::ostream& dest = std::cout){
		for (const auto &i : v) {
			dest << i << std::endl;
		}
	}
	template <typename T>
	void print(const vector<vector<T>> &v, std::ostream& dest = std::cout){
		for (const auto &row : v) {
			for (const auto &col : row) {
				dest << col << "\t";
			}
			dest << std::endl;
		}
	}
	template <typename T>
	void print(const vector<T> &v, std::ostream& dest = std::cout){
		for (const auto &i : v) {
			dest << i << "\t";
		}
		dest << std::endl;
	}
	template <typename T>
	void print(const T &t, std::ostream& dest = std::cout){
		dest << t << std::endl;
	}
}

namespace io {
	bool isfile(const std::string& filename) {
		std::ifstream f(filename);// .c_str());
		return f.good();
	}

	bool isfile(const std::string& filename, int errcode) {
		if (!isfile(filename)) throw std::pair<int, std::string>(errcode, "Failed to load " + filename + ".");
		return true; //if no throw, then OK!
	}
}

namespace types {
	template <typename U,typename T>
	std::vector<U> reinterpret(const std::vector<T> &in){
		//change vector type but keep an exact copy of buffer
		std::vector<U> out;
		if (in.empty()) return out; //nothing to do
		size_t typesize = sizeof(T);
		size_t nbytes = in.size()*typesize;
		size_t out_nelems = nbytes / sizeof(U);
		//printf("%i %i %i %i\n",typesize,sizeof(U),nbytes,out_nelems);
		out.resize(out_nelems);
		std::memcpy(out.data(), in.data(), nbytes);
		return out;
	}

	template <typename T>
	void swap_endianness(std::vector<char> &in){ //in place
		if (in.empty()) return; //nothing to do
		size_t typesize = sizeof(T);
		if (typesize == 1) return; //no endianness with onebyte types

		std::vector<char> out(in);
		for (size_t i = 0; i < in.size(); i += typesize) {
			for (size_t j = 0; j < typesize; j++){
				in[i + j] = out[i + typesize - 1 - j];
			}
		}
		return;
	}

	// lexical_cast, default is fallback
	template<class T>
	T lexical_cast(const std::string &str)
	{
		static std::istringstream very_long_name_ss; /* reusing has severe (positive) impact on performance */
		T value;
		very_long_name_ss.str(str);
		very_long_name_ss >> value;
		very_long_name_ss.clear();
		return value;
	}

	// lexical_cast, trivial conversion
	template<> std::string lexical_cast(const std::string &str) { return str; }

	// lexical_cast, conversions that exist in stl
	template<> float lexical_cast(const std::string &str) { return std::stof(str); }
	template<> double lexical_cast(const std::string &str) { return std::stod(str); }
	template<> long double lexical_cast(const std::string &str) { return std::stold(str); }
	template<> int lexical_cast(const std::string &str) { return std::stoi(str); }
	template<> long lexical_cast(const std::string &str) { return std::stol(str); }
	template<> long long lexical_cast(const std::string &str) { return std::stoll(str); }
	template<> unsigned long lexical_cast(const std::string &str) { return std::stoul(str); }
	template<> unsigned long long lexical_cast(const std::string &str) { return std::stoull(str); }

	// cast, conversions that need to be truncated
	template<> short lexical_cast(const std::string &str) { return static_cast<short>(lexical_cast<long>(str)); }
	template<> unsigned short lexical_cast(const std::string &str) { return static_cast<unsigned short>(lexical_cast<unsigned long>(str)); }
	template<> unsigned int lexical_cast(const std::string &str) { return static_cast<unsigned int>(lexical_cast<unsigned long>(str)); }

	// cast while splitting std::string
	template <typename T>
	std::vector<T> split(const std::string &source, const std::string &delim = ""){
		std::vector<std::string> s;
		if (delim.size() == 0){
			s = pystring::split(source);
		}
		else {
			s = pystring::split(source, delim);
		}
		std::vector<T> ret;
		for (const auto &i : s){
			ret.push_back(lexical_cast<T>(i));
		}
		return ret;
	}
}

namespace parse {

	int get_index(const std::string &str, const int &occurance){
		if (str.find("[", 1) != std::string::npos) { // there are brackets!
			int index_first_bracket = 1; //skip first char
			for (int i = 0; i < occurance; i++){
				index_first_bracket = str.find("[", index_first_bracket) + 1;
			}
			int index_second_bracket = str.find("]", index_first_bracket);
			return stoi(str.substr(index_first_bracket, index_second_bracket - index_first_bracket));
		}
		else {
			return -1;
		}

	}
	int get_index(const std::string &str){ return get_index(str, 1); }

	std::vector<std::pair<std::string, std::string>> parse_dump(const std::vector<std::string> &source) {
		std::vector<std::pair<std::string, std::string>> ret;
		for (const auto &item : source){
			if (pystring::split(item, "=").size() == 2) { //cant have pairs otherwise
				std::string key(pystring::strip(pystring::split(item, "=")[0]));
				if (pystring::startswith(key, "\"")) continue;
				std::string value(pystring::strip(pystring::split(item, "=")[1]));
				if (pystring::startswith(value, "\"") && value.length() > 2) {
					value = value.substr(1, value.size() - 2); //strips quotation marks
				}
				ret.push_back({ key, value });
			}
			else { // no two parts
				ret.push_back({ item, "" });
			}
		}
		return ret;
	}

	std::vector<std::pair<std::string, std::string>> parse_dump(const std::string &source) {
		return parse_dump(pystring::split(source, "\n"));
	}

	std::vector<std::pair<std::string, std::string>> load_dump(const std::string &dumpfile) {
		return parse_dump(vect::fromfile<std::string>(dumpfile));
	}

}