// build a weighting function based on distance to the original location
@Override
	protected void execute(CalculationMonitor monitor){
		// import the image data into 1D arrays : TO DO
		ImageDataFloat	layersImg = new ImageDataFloat(layersImage.getImageData());
		ImageDataFloat	intensImg = new ImageDataFloat(intensityImage.getImageData());
		int nx = layersImg.getRows();
		int ny = layersImg.getCols();
		int nz = layersImg.getSlices();
		int nlayers = layersImg.getComponents()-1;
		int nxyz = nx*ny*nz;
		float rx = layersImg.getHeader().getDimResolutions()[0];
		float ry = layersImg.getHeader().getDimResolutions()[1];
		float rz = layersImg.getHeader().getDimResolutions()[2];
		float[][] layers = new float[nlayers+1][nxyz];
		float[][][][] buffer4 = layersImg.toArray4d();
		for (int x=0;x<nx;x++) for (int y=0;y<ny;y++) for (int z=0;z<nz;z++) for (int l=0;l<=nlayers;l++) {
			int xyz = x+nx*y+nx*ny*z;
			layers[l][xyz] = buffer4[x][y][z][l];
		}
		buffer4 = null;
		layersImg = null;
		float[] intensity = new float[nxyz];
		float[][][] buffer3 = intensImg.toArray3d();
		for (int x=0;x<nx;x++) for (int y=0;y<ny;y++) for (int z=0;z<nz;z++) {
			int xyz = x+nx*y+nx*ny*z;
			intensity[xyz] = buffer3[x][y][z];
		}
		buffer3 = null;
		intensImg = null;
		// create a mask for all the regions outside of the area where layer 1 is > 0 and layer 2 is < 0
		boolean[] ctxmask = new boolean[nxyz];
		if (maskImage.getImageData()!=null) {
			ImageDataUByte	maskImg = new ImageDataUByte(maskImage.getImageData());
			byte[][][] bufferbyte = maskImg.toArray3d();
			for (int x=0;x<nx;x++) for (int y=0;y<ny;y++) for (int z=0;z<nz;z++) {
				int xyz = x+nx*y+nx*ny*z;
				ctxmask[xyz] = (layers[0][xyz]>=0.0 && layers[nlayers][xyz]<=0.0 && bufferbyte[x][y][z]>0);
			}
			bufferbyte = null;
			maskImg = null;
		} else {
			for (int xyz=0;xyz<nxyz;xyz++) {
				ctxmask[xyz] = (layers[0][xyz]>=0.0 && layers[nlayers][xyz]<=0.0);
			}
		}
		// main algorithm
		// 1. define partial voume for each layer, each voxel
		float[][] pvol = new float[nlayers+1][nxyz];
		for (int x=0; x<nx; x++) for (int y=0; y<ny; y++) for (int z = 0; z<nz; z++) {
			int xyz = x + nx*y + nx*ny*z;
			if (ctxmask[xyz]) {
				for (int l=0;l<=nlayers;l++) {
					pvol[l][xyz] = partialVolumeFromSurface(x, y, z, layers[l], nx, ny, nz);
				}
			}
		}
		// 2. build and invert the GLM for each profile / voxel?
		float delta = paramExtent.getValue().floatValue()/Numerics.min(rx,ry,rz);
		float stdev = paramStdev.getValue().floatValue()/Numerics.min(rx,ry,rz);
		CorticalProfile profile = new CorticalProfile(nlayers, nx, ny, nz, rx, ry, rz);
		float maskval = 1e13f;
		float[][][][] mapping = new float[nx][ny][nz][nlayers];
		float[][][] residual = new float[nx][ny][nz];
		BitSet sampled = new BitSet(nx*ny*nz);
		float[] profiledist = new float[nx*ny*nz];
		for (int x=0; x<nx; x++) for (int y=0; y<ny; y++) for (int z = 0; z<nz; z++) {
			int xyz = x + nx*y + nx*ny*z;
			if (ctxmask[xyz]) {
				findFastMarchingProfileNeighborhood(sampled, profiledist, x,y,z, delta, layers, profile, ctxmask, nx, ny, nz, nlayers);
				int nsample = sampled.size();
				if (nsample>=nlayers) {
					double[][] glm = new double[nlayers][nsample];
					double[][] data = new double[nsample][1];
					int idx = 0;
					for (int n=0;n<nsample;n++) {
						// get the next non-zero value
						idx = sampled.nextSetBit(idx);
						// build a weighting function based on distance to the original location
						double weight = FastMath.exp(-0.5*Numerics.square(profiledist[idx]/stdev));
						for (int l=0;l<nlayers;l++) {
							glm[l][n] = weight*(pvol[l+1][idx]-pvol[l][idx]);
						}
						data[n][0] = weight*intensity[idx];
					}
					// invert the linear model
					Matrix mtx = new Matrix(glm);
					Matrix smp = new Matrix(data);
					Matrix val = mtx.solve(smp);
					for (int l=0;l<nlayers;l++) {
						mapping[x][y][z][l] = (float)val.get(l,0);
					}
					Matrix res = mtx.times(val).minus(smp);
					residual[x][y][z] = (float)res.normInf();
				}
			}
		}
		// output
		String imgname = intensityImage.getImageData().getName();
		ImageDataFloat mapData = new ImageDataFloat(mapping);		
		mapData.setHeader(layersImage.getImageData().getHeader());
		mapData.setName(imgname+"_glmprofiles");
		mappedImage.setValue(mapData);
		mapData = null;
		mapping = null;
		ImageDataFloat resData = new ImageDataFloat(residual);		
		resData.setHeader(layersImage.getImageData().getHeader());
		resData.setName(imgname+"_glmresidual");
		residualImage.setValue(resData);
		resData = null;
		residual = null;
	}
// get the next non-zero value
idx = sampled.nextSetBit(idx);
// build a weighting function based on distance to the original location
double weight = FastMath.exp(-0.5*Numerics.square(profiledist[idx]/stdev));
for (int l=0;l<nlayers;l++) {
if (ctxmask[xyz]) {
findFastMarchingProfileNeighborhood(sampled, profiledist, x,y,z, delta, layers, profile, ctxmask, nx, ny, nz, nlayers);
int nsample = sampled.size();
if (nsample>=nlayers) {
double[][] glm = new double[nlayers][nsample];
double[][] data = new double[nsample][1];
int idx = 0;
for (int n=0;n<nsample;n++) {
// get the next non-zero value
idx = sampled.nextSetBit(idx);
// build a weighting function based on distance to the original location
double weight = FastMath.exp(-0.5*Numerics.square(profiledist[idx]/stdev));
for (int l=0;l<nlayers;l++) {
glm[l][n] = weight*(pvol[l+1][idx]-pvol[l][idx]);
}
data[n][0] = weight*intensity[idx];
}
// invert the linear model
Matrix mtx = new Matrix(glm);
Matrix smp = new Matrix(data);
Matrix val = mtx.solve(smp);
float delta = paramExtent.getValue().floatValue()/Numerics.min(rx,ry,rz);
float stdev = paramStdev.getValue().floatValue()/Numerics.min(rx,ry,rz);
CorticalProfile profile = new CorticalProfile(nlayers, nx, ny, nz, rx, ry, rz);
float maskval = 1e13f;
float[][][][] mapping = new float[nx][ny][nz][nlayers];
float[][][] residual = new float[nx][ny][nz];
BitSet sampled = new BitSet(nx*ny*nz);
float[] profiledist = new float[nx*ny*nz];
for (int x=0; x<nx; x++) for (int y=0; y<ny; y++) for (int z = 0; z<nz; z++) {
int xyz = x + nx*y + nx*ny*z;
if (ctxmask[xyz]) {
findFastMarchingProfileNeighborhood(sampled, profiledist, x,y,z, delta, layers, profile, ctxmask, nx, ny, nz, nlayers);
int nsample = sampled.size();
if (nsample>=nlayers) {
double[][] glm = new double[nlayers][nsample];
double[][] data = new double[nsample][1];
int idx = 0;
for (int n=0;n<nsample;n++) {
// get the next non-zero value
idx = sampled.nextSetBit(idx);
// build a weighting function based on distance to the original location
double weight = FastMath.exp(-0.5*Numerics.square(profiledist[idx]/stdev));
for (int l=0;l<nlayers;l++) {
glm[l][n] = weight*(pvol[l+1][idx]-pvol[l][idx]);
}
data[n][0] = weight*intensity[idx];
}
// invert the linear model
Matrix mtx = new Matrix(glm);
Matrix smp = new Matrix(data);
Matrix val = mtx.solve(smp);
for (int l=0;l<nlayers;l++) {
mapping[x][y][z][l] = (float)val.get(l,0);
}
Matrix res = mtx.times(val).minus(smp);
residual[x][y][z] = (float)res.normInf();
}
}
}
// output
String imgname = intensityImage.getImageData().getName();