/* parse_func.cpp */

/* Checks the validity of the arguments input to the convolution routine and 
transfers the Matlab-style arguments to the C structures defined in defs.h. */

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#include <string.h>
#include "mex.h"
#include "defs.h"

void check_MATLAB_KERNELS(const mxArray *Kernels)
/* Ensures that the MATLAB_KERNELS structure has the following format:

MATLAB_KERNEL = 

               radii: [1xNradii float]
              angles: [1xNangles float]
            energies: [1xNenergies float]
             primary: [Nradii x Nangles x Nenergies float]
       first_scatter: [Nradii x Nangles x Nenergies float]
      second_scatter: [Nradii x Nangles x Nenergies float]
    multiple_scatter: [Nradii x Nangles x Nenergies float]
          brem_annih: [Nradii x Nangles x Nenergies float]
               total: [Nradii x Nangles x Nenergies float]
          mean_radii: [Nradii x Nangles x Nenergies float] (not used)
         mean_angles: [Nradii x Nangles x Nenergies float] (not used)
           helpfield: [any x any char] 
             fluence: [1xNenergies float]
                  mu: [1xNenergies float]
               mu_en: [1xNenergies float]
*/
{
	mxArray *tmp;
	int Nenergies, Nangles, Nradii, length;
	const int *dims;
	char err_msg[100];

	// check 'energies' field
	if ( (tmp=mxGetField(Kernels,0,"energies")) == NULL)
		mexErrMsgTxt("Missing 'energies' field in Kernels structure.");
	else
		Nenergies = MAXX(mxGetM(tmp),mxGetN(tmp));

	// check fields that are dependent on 'energies' field ('fluence','mu','mu_en')
	if ( (tmp=mxGetField(Kernels,0,"fluence")) == NULL)
		mexErrMsgTxt("Missing 'fluence' field in Kernels structure.");
	else
	{
		length = MAXX(mxGetM(tmp),mxGetN(tmp));
		if (length != Nenergies)
			mexErrMsgTxt("'energies' field and 'fluence' field in Kernels structure must have the same length");
	}

	if ( (tmp=mxGetField(Kernels,0,"mu")) == NULL)
		mexErrMsgTxt("Missing 'mu' field in Kernels structure.");
	else
	{
		length = MAXX(mxGetM(tmp),mxGetN(tmp));
		if (length != Nenergies)
			mexErrMsgTxt("'energies' field and 'mu' field in Kernels structure must have the same length");
	}

	if ( (tmp=mxGetField(Kernels,0,"mu_en")) == NULL)
		mexErrMsgTxt("Missing 'mu_en' field in Kernels structure.");
	else
	{
		length = MAXX(mxGetM(tmp),mxGetN(tmp));
		if (length != Nenergies)
			mexErrMsgTxt("'energies' field and 'mu_en' field in Kernels structure must have the same length");
	}

	// check 'angles' field
	if ( (tmp=mxGetField(Kernels,0,"angles")) == NULL)
		mexErrMsgTxt("Missing 'angles' field in Kernels structure.");
	else
		Nangles = MAXX(mxGetM(tmp),mxGetN(tmp));

	// check 'radii' field
	if ( (tmp=mxGetField(Kernels,0,"radii")) == NULL)
		mexErrMsgTxt("Missing 'radii' field in Kernels structure.");
	else
		Nradii = MAXX(mxGetM(tmp),mxGetN(tmp));

	// check kernels, which are dependent on 'radii', 'angles', and 'energies' fields

	// check primary kernel
	if ( (tmp=mxGetField(Kernels,0,"primary")) == NULL)
		mexErrMsgTxt("Missing 'primary' kernel in Kernels structure.");
	else
	{
		length = mxGetNumberOfDimensions(tmp);
		dims = mxGetDimensions(tmp);
		if ((length != 3) || (dims[0] != Nradii) || (dims[1] != Nangles)
			|| (dims[2] != Nenergies))
		{ 
			sprintf(err_msg,"'primary' kernel in Kernels structure must be a %d x %d x %d matrix", 
			    Nradii, Nangles, Nenergies);
			mexErrMsgTxt(err_msg);
		}
	}

	// check first scatter kernel
	if ( (tmp=mxGetField(Kernels,0,"first_scatter")) == NULL)
		mexErrMsgTxt("Missing 'first_scatter' kernel in Kernels structure.");
	else
	{
		length = mxGetNumberOfDimensions(tmp);
		dims = mxGetDimensions(tmp);
		if ((length != 3) || (dims[0] != Nradii) || (dims[1] != Nangles)
			|| (dims[2] != Nenergies))
		{
			sprintf(err_msg,"'first_scatter' kernel in Kernels structure must be a %d x %d x %d matrix", 
			     Nradii, Nangles, Nenergies);
			mexErrMsgTxt(err_msg);
		}
	}

    // check second scatter kernel
	if ( (tmp=mxGetField(Kernels,0,"second_scatter")) == NULL)
		mexErrMsgTxt("Missing 'second_scatter' kernel in Kernels structure.");
	else
	{
		length = mxGetNumberOfDimensions(tmp);
		dims = mxGetDimensions(tmp);
		if ((length != 3) || (dims[0] != Nradii) || (dims[1] != Nangles)
			|| (dims[2] != Nenergies))
		{
			sprintf(err_msg,"'second_scatter' kernel in Kernels structure must be a %d x %d x %d matrix", 
			     Nradii, Nangles, Nenergies);
			mexErrMsgTxt(err_msg);
		}
	}

	// check multiple scatter kernel
	if ( (tmp=mxGetField(Kernels,0,"multiple_scatter")) == NULL)
		mexErrMsgTxt("Missing 'multiple_scatter' kernel in Kernels structure.");
	else
	{
		length = mxGetNumberOfDimensions(tmp);
		dims = mxGetDimensions(tmp);
		if ((length != 3) || (dims[0] != Nradii) || (dims[1] != Nangles)
			|| (dims[2] != Nenergies))
		{
			sprintf(err_msg,"'multiple_scatter' kernel in Kernels structure must be a %d x %d x %d matrix", 
			     Nradii, Nangles, Nenergies);
			mexErrMsgTxt(err_msg);
		}
	}

	// check bremstrahlung & annihilation in flight kernel
	if ( (tmp=mxGetField(Kernels,0,"brem_annih")) == NULL)
		mexErrMsgTxt("Missing 'brem_annih' kernel in Kernels structure.");
	else
	{
		length = mxGetNumberOfDimensions(tmp);
		dims = mxGetDimensions(tmp);
		if ((length != 3) || (dims[0] != Nradii) || (dims[1] != Nangles)
			|| (dims[2] != Nenergies))
		{
			sprintf(err_msg,"'brem_annih' kernel in Kernels structure must be a %d x %d x %d matrix", 
			     Nradii, Nangles, Nenergies);
			mexErrMsgTxt(err_msg);
		}
	}

	// check total kernel
	if ( (tmp=mxGetField(Kernels,0,"total")) == NULL)
		mexErrMsgTxt("Missing 'total' kernel in Kernels structure.");
	else
	{
		length = mxGetNumberOfDimensions(tmp);
		dims = mxGetDimensions(tmp);
		if ((length != 3) || (dims[0] != Nradii) || (dims[1] != Nangles)
			|| (dims[2] != Nenergies))
		{
			sprintf(err_msg,"'total' kernel in Kernels structure must be a %d x %d x %d matrix", 
			     Nradii, Nangles, Nenergies);
			mexErrMsgTxt(err_msg);
		}
	} 
 
	// the mean_radii and mean_angles fields are not needed, so don't check for them
}

void check_GEOMETRY(const mxArray *Geometry)
/* Ensure that the GEOMETRY structure has the following format:
  Geometry = 

         start: [3 element float]
    voxel_size: [3 element float]
          data: [M x N x Q float]
          beam: [1xB cell array] 
		  
 Geometry.beam = 

    ip: [3 element float]
    jp: [3 element float]
    kp: [3 element float]
     y_vec: [3 element float]
        xp: [float scalar]
        yp: [float scalar]
    del_xp: [float scalar]
    del_yp: [float scalar]
	   SAD: [float scalar]
	
*/
{
	mxArray *tmp1, *tmp2, *cell;
	int length, B, b;
	// beam's eye view axis unit vectors, must be orthonormal and form a right-handed
	// coordinate system.
	float *ip, *jp, *kp;
	const int *dims;
	char errstr[200];

	// check that Geometry.start is a 3 element vector
	if ( (tmp1 = mxGetField(Geometry,0,"start")) == NULL)
		mexErrMsgTxt("Missing 'start' field in Geometry structure.");
	else
		if (mxGetNumberOfDimensions(tmp1) != 2)
			mexErrMsgTxt("Geometry.start must have 2 dimensions but be a vector.");
		else
			if (mxGetM(tmp1)*mxGetN(tmp1) != 3)
				mexErrMsgTxt("Geometry.start must be a 3 element vector.");

	// check that Geometry.voxel_size is a 3 element vector
	if ( (tmp1 = mxGetField(Geometry,0,"voxel_size")) == NULL)
		mexErrMsgTxt("Missing 'voxel_size' field in Geometry structure.");
	else
		if (mxGetNumberOfDimensions(tmp1) != 2)
			mexErrMsgTxt("Geometry.voxel_size must have 2 dimensions but be a vector.");
		else
			if (mxGetM(tmp1)*mxGetN(tmp1) != 3)
				mexErrMsgTxt("Geometry.voxel_size must be a 3 element vector.");

	// check that all parts of 'beam' field are present
	if ( (tmp1 = mxGetField(Geometry,0,"beam")) == NULL)
		mexErrMsgTxt("Missing 'beam' field in Geometry structure.");
	else
	{
		if (!mxIsCell(tmp1) || (mxGetNumberOfDimensions(tmp1) != 2))
			mexErrMsgTxt("'beam' field in Geometry structure must be a 2-D cell array.");
		else
		{
			dims = mxGetDimensions(tmp1);  // dimensions of cell array
			B = dims[0]*dims[1];           // total number of cells
			
			for (b=0;b<B;b++) // check all of the beam fields
			{
				if((cell = mxGetCell(tmp1,b)) == NULL || !mxIsStruct(cell))
				{
					sprintf(errstr,"Geometry.beam{%d} either doesn't exist or isn't a structure.",b+1);
					mexErrMsgTxt(errstr);
				}
				
				// start checking the beam fields
				if((tmp2 = mxGetField(cell,0,"ip")) == NULL)
				{
					sprintf(errstr,"Geometry.beam{%d} must contain a ip field.",b+1);
					mexErrMsgTxt(errstr);
				}
				else if (mxGetNumberOfDimensions(tmp2) != 2)
				{
					sprintf(errstr,"Geometry.beam{%d}.ip must have 2 dimensions but be a vector.",b+1);
					mexErrMsgTxt(errstr);
				}
				else if (mxGetM(tmp2)*mxGetN(tmp2) != 3)
				{
					sprintf(errstr,"Geometry.beam{%d}.ip must be a 3 element vector.",b+1);
					mexErrMsgTxt(errstr);
				}
				else
					ip = (float *)mxGetData(tmp2);

				if((tmp2 = mxGetField(cell,0,"jp")) == NULL)
				{
					sprintf(errstr,"Geometry.beam{%d} must contain a jp field.",b+1);
					mexErrMsgTxt(errstr);
				}
				else if (mxGetNumberOfDimensions(tmp2) != 2)
				{
					sprintf(errstr,"Geometry.beam{%d}.jp must have 2 dimensions but be a vector.",b+1);
					mexErrMsgTxt(errstr);
				}
				else if (mxGetM(tmp2)*mxGetN(tmp2) != 3)
				{
					sprintf(errstr,"Geometry.beam{%d}.jp must be a 3 element vector.",b+1);
					mexErrMsgTxt(errstr);
				}
				else
					jp = (float *)mxGetData(tmp2);

				if((tmp2 = mxGetField(cell,0,"kp")) == NULL)
				{
					sprintf(errstr,"Geometry.beam{%d} must contain a ip field.",b+1);
					mexErrMsgTxt(errstr);
				}
				else if (mxGetNumberOfDimensions(tmp2) != 2)
				{
					sprintf(errstr,"Geometry.beam{%d}.kp must have 2 dimensions but be a vector.",b+1);
					mexErrMsgTxt(errstr);
				}
				else if (mxGetM(tmp2)*mxGetN(tmp2) != 3)
				{
					sprintf(errstr,"Geometry.beam{%d}.kp must be a 3 element vector.",b+1);
					mexErrMsgTxt(errstr);
				}
				else
					kp = (float *)mxGetData(tmp2);
				
				// ensure that ip,jp,kp are orthonormal and form a right-handed coordinate system
				if (   (fabs(ip[0]*ip[0] + ip[1]*ip[1] + ip[2]*ip[2] - 1.0) > 1e-6)
					|| (fabs(jp[0]*jp[0] + jp[1]*jp[1] + jp[2]*jp[2] - 1.0) > 1e-6)
					|| (fabs(kp[0]*kp[0] + kp[1]*kp[1] + kp[2]*kp[2] - 1.0) > 1e-6)
					|| (fabs(ip[0]*jp[0] + ip[1]*jp[1] + ip[2]*jp[2]) > 1e-6)
					|| (fabs(ip[0]*kp[0] + ip[1]*kp[1] + ip[2]*kp[2]) > 1e-6)
					|| (fabs(jp[0]*kp[0] + jp[1]*kp[1] + jp[2]*kp[2]) > 1e-6))
				{
					sprintf(errstr,"Beam's eye view unit vectors ip, jp, and kp for beam %d must be orthonormal.",b+1);
					mexErrMsgTxt(errstr);
				}

				// check that cross(ip,jp) = kp
				if (   (fabs(ip[1]*jp[2] - ip[2]*jp[1] - kp[0]) > 1e-6)
					|| (fabs(ip[2]*jp[0] - ip[0]*jp[2] - kp[1]) > 1e-6)
					|| (fabs(ip[0]*jp[1] - ip[1]*jp[0] - kp[2]) > 1e-6))
				{
					sprintf(errstr,"Must have cross(ip,jp) = kp for beam %d in beam's eye view.",b+1);
					mexErrMsgTxt(errstr);
				}

				if((tmp2 = mxGetField(cell,0,"y_vec")) == NULL)
				{
					sprintf(errstr,"Geometry.beam{%d} must contain a y_vec field.",b+1);
					mexErrMsgTxt(errstr);
				}
				else
					if (mxGetNumberOfDimensions(tmp2) != 2)
					{
						sprintf(errstr,"Geometry.beam{%d}.y_vec must have 2 dimensions but be a vector.",b+1);
						mexErrMsgTxt(errstr);
					}
					else
						if (mxGetM(tmp2)*mxGetN(tmp2) != 3)
						{
							sprintf(errstr,"Geometry.beam{%d}.y_vec must be a 3 element vector.",b+1);
							mexErrMsgTxt(errstr);
						}

				if((tmp2 = mxGetField(cell,0,"xp")) == NULL)
				{
					sprintf(errstr,"Geometry.beam{%d} must contain a xp field.",b+1);
					mexErrMsgTxt(errstr);
				}
				else
					if (mxGetNumberOfDimensions(tmp2) != 2)
					{
						sprintf(errstr,"Geometry.beam{%d}.xp must have 2 dimensions but be a scalar.",b+1);
						mexErrMsgTxt(errstr);
					}
					else
						if (mxGetM(tmp2)*mxGetN(tmp2) != 1)
						{
							sprintf(errstr,"Geometry.beam{%d}.xp must be a scalar.",b+1);
							mexErrMsgTxt(errstr);
						}

				if((tmp2 = mxGetField(cell,0,"yp")) == NULL)
				{
					sprintf(errstr,"Geometry.beam{%d} must contain a yp field.",b+1);
					mexErrMsgTxt(errstr);
				}
				else
					if (mxGetNumberOfDimensions(tmp2) != 2)
					{
						sprintf(errstr,"Geometry.beam{%d}.yp must have 2 dimensions but be a scalar.",b+1);
						mexErrMsgTxt(errstr);
					}
					else
						if (mxGetM(tmp2)*mxGetN(tmp2) != 1)
						{
							sprintf(errstr,"Geometry.beam{%d}.yp must be a scalar.",b+1);
							mexErrMsgTxt(errstr);
						}

				if((tmp2 = mxGetField(cell,0,"del_xp")) == NULL)
				{
					sprintf(errstr,"Geometry.beam{%d} must contain a del_xp field.",b+1);
					mexErrMsgTxt(errstr);
				}
				else
					if (mxGetNumberOfDimensions(tmp2) != 2)
					{
						sprintf(errstr,"Geometry.beam{%d}.del_xp must have 2 dimensions but be a scalar.",b+1);
						mexErrMsgTxt(errstr);
					}
					else
						if (mxGetM(tmp2)*mxGetN(tmp2) != 1)
						{
							sprintf(errstr,"Geometry.beam{%d}.del_xp must be a scalar.",b+1);
							mexErrMsgTxt(errstr);
						}

				if((tmp2 = mxGetField(cell,0,"del_yp")) == NULL)
				{
					sprintf(errstr,"Geometry.beam{%d} must contain a del_yp field.",b+1);
					mexErrMsgTxt(errstr);
				}
				else
					if (mxGetNumberOfDimensions(tmp2) != 2)
					{
						sprintf(errstr,"Geometry.beam{%d}.del_yp must have 2 dimensions but be a scalar.",b+1);
						mexErrMsgTxt(errstr);
					}
					else
						if (mxGetM(tmp2)*mxGetN(tmp2) != 1)
						{
							sprintf(errstr,"Geometry.beam{%d}.del_yp must be a scalar.",b+1);
							mexErrMsgTxt(errstr);
						}

				if((tmp2 = mxGetField(cell,0,"SAD")) == NULL)
				{
					sprintf(errstr,"Geometry.beam{%d} must contain an SAD field.",b+1);
					mexErrMsgTxt(errstr);
				}
				else
					if (mxGetNumberOfDimensions(tmp2) != 2)
					{
						sprintf(errstr,"Geometry.beam{%d}.SAD must have 2 dimensions but be a scalar.",b+1);
						mexErrMsgTxt(errstr);
					}
					else
						if (mxGetM(tmp2)*mxGetN(tmp2) != 1)
						{
							sprintf(errstr,"Geometry.beam{%d}.SAD must be a scalar.",b+1);
							mexErrMsgTxt(errstr);
						}
			}
		}
	}

	// check that data matrix is present and 3-D 
    //(3-D part maybe not necessary for calculation itself, I'm looking into this)
	if ( (tmp1 = mxGetField(Geometry,0,"data")) == NULL)
		mexErrMsgTxt("Missing 'data' field in Geometry structure.");
	else
	{
		length = mxGetNumberOfDimensions(tmp1);
		if (length != 3)
			mexErrMsgTxt("Geometry.data must be a 3-D matrix.");
	}
}

void load_kernels(const mxArray *Kernels, MONO_KERNELS *mono_kernels)
/* Turns Matlab structures into C structures. */
{
	int i, Nenergies, Nangles, Nradii;
	float *dataPtr;
	mxArray *tmp;

	/* find number of energies to consider in the spectrum */
	tmp = mxGetField(Kernels,0,"energies");
	Nenergies = MAXX(mxGetM(tmp),mxGetN(tmp));
	mono_kernels->energy = (float *)mxGetData(tmp);
	mono_kernels->nkernels = Nenergies;  // one kernel for each energy

	tmp = mxGetField(Kernels,0,"fluence");
	mono_kernels->fluence = (float *)mxGetData(tmp);

	tmp = mxGetField(Kernels,0,"mu");
	mono_kernels->mu = (float *)mxGetData(tmp);

	tmp = mxGetField(Kernels,0,"mu_en");
	mono_kernels->mu_en = (float *)mxGetData(tmp);

	/*  Load kernels in one at a time */
	for (i=0; i<Nenergies; i++)
	{
		tmp = mxGetField(Kernels,0,"radii");
		Nradii = MAXX(mxGetM(tmp),mxGetN(tmp));
		mono_kernels->kernel[i].nradii = Nradii;
		mono_kernels->kernel[i].radial_boundary = (float *)mxGetData(tmp); 

		tmp = mxGetField(Kernels,0,"angles");
		Nangles = MAXX(mxGetM(tmp),mxGetN(tmp));
		mono_kernels->kernel[i].ntheta = Nangles;
		mono_kernels->kernel[i].angular_boundary = (float *)mxGetData(tmp); 

		/* Assign data pointers. Data are in Nradii x Nangles x Nenergies matrices */
		tmp = mxGetField(Kernels,0,"primary");
		dataPtr = (float *)mxGetData(tmp);
		mono_kernels->kernel[i].matrix[0] = &dataPtr[Nradii*Nangles*i];

		tmp = mxGetField(Kernels,0,"first_scatter");
		dataPtr = (float *)mxGetData(tmp);
		mono_kernels->kernel[i].matrix[1] = &dataPtr[Nradii*Nangles*i];

		tmp = mxGetField(Kernels,0,"second_scatter");
		dataPtr = (float *)mxGetData(tmp);
		mono_kernels->kernel[i].matrix[2] = &dataPtr[Nradii*Nangles*i];

		tmp = mxGetField(Kernels,0,"multiple_scatter");
		dataPtr = (float *)mxGetData(tmp);
		mono_kernels->kernel[i].matrix[3] = &dataPtr[Nradii*Nangles*i];

		tmp = mxGetField(Kernels,0,"brem_annih");
		dataPtr = (float *)mxGetData(tmp);
		mono_kernels->kernel[i].matrix[4] = &dataPtr[Nradii*Nangles*i];

		tmp = mxGetField(Kernels,0,"total");
		dataPtr = (float *)mxGetData(tmp);
		mono_kernels->kernel[i].total_matrix = &dataPtr[Nradii*Nangles*i];
	}
}

void load_Geometry(const mxArray *Geometry, FLOAT_GRID *density)
// Loads the non-beam information of the Geometry Matab structure into a C structure.
// Note that the CT matrix is specified in units of density, not Hounsfield units.
// All Hounsfield unit conversions must be done before the convolution program is called.

/* Ensure that the GEOMETRY structure has the following format:
  Geometry = 

         start: [3 element float]
    voxel_size: [3 element float]
          data: [M x N x Q float]
          beam: [1x1 struct] 
		 
*/
{
	const int *dimensions;
	int i;
	float *dataPtr;
	mxArray *tmp1, *tmp2;

	/* Geometry is a Matlab structure with fields that correspond to the 
	C structures defined for this program.  */

	tmp1 = mxGetField(Geometry,0,"start");

	dataPtr = (float *)mxGetData(tmp1);
	density->start.x = dataPtr[0];
	density->start.y = dataPtr[1];
	density->start.z = dataPtr[2];

	/* read in voxel size */
	tmp1 = mxGetField(Geometry,0,"voxel_size");

	dataPtr = (float *)mxGetData(tmp1);
	density->inc.x = dataPtr[0];
	density->inc.y = dataPtr[1];
	density->inc.z = dataPtr[2];

	/* get pointer to CT data matrix */
	/* Data stored in an (x,y,z) matrix */
	tmp1 = mxGetField(Geometry,0,"data");
	density->matrix = (float *)mxGetData(tmp1);
	dimensions = mxGetDimensions(tmp1);
	
	density->x_count = dimensions[0];
	density->y_count = dimensions[1];
	density->z_count = dimensions[2];
}

void load_beam(const mxArray *Geometry, int beam_ind, BEAM *beam)
// Loads Geometry.beam{beam_ind} into a BEAM structure.  To avoid errors
// ensure that check_Geometry was run before this function.
/* Ensure that the 'beam' structure has the proper format:

  Geometry.beam = 

    ip: [3 element float]
    jp: [3 element float]
    kp: [3 element float]
     y_vec: [3 element float]
        xp: [float scalar]
        yp: [float scalar]
    del_xp: [float scalar]
    del_yp: [float scalar] */
{

	const int *dimensions;
	int i;
	float *dataPtr;
	mxArray *cell, *struc;
	
	/* read in beam information */
	cell = mxGetCell(mxGetField(Geometry,0,"beam"),beam_ind);

	struc = mxGetField(cell,0,"ip");
	dataPtr = (float *)mxGetData(struc);

	for (i=0;i<=2;i++) beam->ip[i] = dataPtr[i];

	struc = mxGetField(cell,0,"jp");
	dataPtr = (float *)mxGetData(struc);
	for (i=0;i<=2;i++) beam->jp[i] = dataPtr[i];

	struc = mxGetField(cell,0,"kp");
	dataPtr = (float *)mxGetData(struc);
	for (i=0;i<=2;i++) beam->kp[i] = dataPtr[i];

	struc = mxGetField(cell,0,"y_vec");
	dataPtr = (float *)mxGetData(struc);
	for (i=0;i<=2;i++) beam->y_vec[i] = dataPtr[i];

	struc = mxGetField(cell,0,"xp");
	beam->xp = (float)mxGetScalar(struc);

	struc = mxGetField(cell,0,"yp");
	beam->yp = (float)mxGetScalar(struc);

	struc = mxGetField(cell,0,"del_xp");
	beam->del_xp = (float)mxGetScalar(struc);

	struc = mxGetField(cell,0,"del_yp");
	beam->del_yp = (float)mxGetScalar(struc);

	struc = mxGetField(cell,0,"SAD");
	beam->SAD = (float)mxGetScalar(struc);
}