codebase / stimlib

#ifndef STIM_SCALARFIELD_H
#define STIM_SCALARFIELD_H


#include "../math/rect.h"
#include "../math/complex.h"
#include "../math/fft.h"

#ifdef CUDA_FOUND
#include "../cuda/crop.cuh"
#endif

namespace stim{

	template<typename T>
	__global__ void cuda_abs(T* img, stim::complex<T>* field, size_t N){
		size_t i = blockIdx.x * blockDim.x + threadIdx.x;
		if(i >= N) return;

		img[i] = field[i].abs();
	}

	template<typename T>
	__global__ void cuda_real(T* img, stim::complex<T>* field, size_t N){
		size_t i = blockIdx.x * blockDim.x + threadIdx.x;
		if(i >= N) return;

		img[i] = field[i].real();
	}

	template<typename T>
	__global__ void cuda_imag(T* img, stim::complex<T>* field, size_t N){
		size_t i = blockIdx.x * blockDim.x + threadIdx.x;
		if(i >= N) return;

		img[i] = field[i].imag();
	}

	template<typename T>
	__global__ void cuda_intensity(T* img, stim::complex<T>* field, size_t N){
		size_t i = blockIdx.x * blockDim.x + threadIdx.x;
		if(i >= N) return;

		img[i] = pow(field[i].abs(), 2);
	}

	template<typename T>
	__global__ void cuda_sum_intensity(T* img, stim::complex<T>* field, size_t N){
		size_t i = blockIdx.x * blockDim.x + threadIdx.x;
		if(i >= N) return;

		img[i] += pow(field[i].abs(), 2);
	}

	/// Perform a k-space transform of a scalar field (FFT). The given field has a width of x and the calculated momentum space has a
	///		width of kx (in radians).
	/// @param K is a pointer to the output array of all plane waves in the field
	/// @param kx is the width of the frame in momentum space
	/// @param ky is the height of the frame in momentum space
	/// @param E is the field to be transformed
	/// @param x is the width of the field in the spatial domain
	/// @param y is the height of the field in the spatial domain
	/// @param nx is the number of pixels representing the field in the x (and kx) direction
	/// @param ny is the number of pixels representing the field in the y (and ky) direction
	template<typename T>
	void cpu_scalar_to_kspace(stim::complex<T>* K, T& kx, T& ky, stim::complex<T>* E, T x, T y, size_t nx, size_t ny){

		kx = stim::TAU * nx / x;			//calculate the width of the momentum space
		ky = stim::TAU * ny / y;

		stim::complex<T>* dev_FFT;
		HANDLE_ERROR( cudaMalloc(&dev_FFT, sizeof(stim::complex<T>) * nx * ny) );		//allocate space on the CUDA device for the output array

		stim::complex<T>* dev_E;
		HANDLE_ERROR( cudaMalloc(&dev_E, sizeof(stim::complex<T>) * nx * ny) );		//allocate space for the field
		HANDLE_ERROR( cudaMemcpy(dev_E, E, sizeof(stim::complex<T>) * nx * ny, cudaMemcpyHostToDevice) );	//copy the field to GPU memory

		cufftResult result;
		cufftHandle plan;
		result = cufftPlan2d(&plan, nx, ny, CUFFT_C2C);
		if(result != CUFFT_SUCCESS){
			std::cout<<"Error creating cuFFT plan."<<std::endl;
			exit(1);
		}

		result = cufftExecC2C(plan, (cufftComplex*)dev_E, (cufftComplex*)dev_FFT, CUFFT_FORWARD);
		if(result != CUFFT_SUCCESS){
			std::cout<<"Error using cuFFT to perform a forward Fourier transform of the field."<<std::endl;
			exit(1);
		}
		cufftDestroy(plan);

		stim::complex<T>* fft = (stim::complex<T>*) malloc(sizeof(stim::complex<T>) * nx * ny);
		HANDLE_ERROR( cudaMemcpy(fft, dev_FFT, sizeof(stim::complex<T>) * nx * ny, cudaMemcpyDeviceToHost) );

		stim::cpu_fftshift(K, fft, nx, ny);
		
		HANDLE_ERROR( cudaFree(dev_FFT) );			//free GPU memory
		HANDLE_ERROR( cudaFree(dev_E) );
		free(fft);									//free CPU memory
	}

	template<typename T>
	void cpu_scalar_from_kspace(stim::complex<T>* E, T& x, T& y, stim::complex<T>* K, T kx, T ky, size_t nx, size_t ny){

		x = stim::TAU * nx / kx;			//calculate the width of the momentum space
		y = stim::TAU * ny / ky;
		
		stim::complex<T>* fft = (stim::complex<T>*) malloc(sizeof(stim::complex<T>) * nx * ny);
		stim::cpu_ifftshift(fft, K, nx, ny);
		//memcpy(fft, K, sizeof(stim::complex<T>) * nx * ny);

		stim::complex<T>* dev_FFT;
		HANDLE_ERROR( cudaMalloc(&dev_FFT, sizeof(stim::complex<T>) * nx * ny) );		//allocate space on the CUDA device for the output array
		HANDLE_ERROR( cudaMemcpy(dev_FFT, fft, sizeof(stim::complex<T>) * nx * ny, cudaMemcpyHostToDevice) );	//copy the field to GPU memory

		stim::complex<T>* dev_E;
		HANDLE_ERROR( cudaMalloc(&dev_E, sizeof(stim::complex<T>) * nx * ny) );		//allocate space for the field

		cufftResult result;
		cufftHandle plan;
		result = cufftPlan2d(&plan, nx, ny, CUFFT_C2C);
		if(result != CUFFT_SUCCESS){
			std::cout<<"Error creating cuFFT plan."<<std::endl;
			exit(1);
		}

		result = cufftExecC2C(plan, (cufftComplex*)dev_FFT, (cufftComplex*)dev_E, CUFFT_INVERSE);
		if(result != CUFFT_SUCCESS){
			std::cout<<"Error using cuFFT to perform a forward Fourier transform of the field."<<std::endl;
			exit(1);
		}
		cufftDestroy(plan);

		HANDLE_ERROR( cudaMemcpy(E, dev_E, sizeof(stim::complex<T>) * nx * ny, cudaMemcpyDeviceToHost) );

		HANDLE_ERROR( cudaFree(dev_FFT) );			//free GPU memory
		HANDLE_ERROR( cudaFree(dev_E) );
		free(fft);									//free CPU memory
		
	}

	

	/// Propagate a field slice along its orthogonal direction by a given distance z
	/// @param Enew is the resulting propogated field
	/// @param E is the field to be propogated
	/// @param sx is the size of the field in the lateral x direction
	/// @param sy is the size of the field in the lateral y direction
	/// @param z is the distance to be propagated
	/// @param k is the wavenumber 2*pi/lambda
	/// @param nx is the number of samples in the field along the lateral x direction
	/// @param ny is the number of samples in the field along the lateral y direction
	template<typename T>
	void cpu_scalar_propagate(stim::complex<T>* Enew, stim::complex<T>* E, T sx, T sy, T z, T k, size_t nx, size_t ny){
		
		stim::complex<T>* K = (stim::complex<T>*) malloc( sizeof(stim::complex<T>) * nx * ny );

		T Kx, Ky;											//width and height in k space
		cpu_scalar_to_kspace(K, Kx, Ky, E ,sx, sy, nx, ny);

		//T* mag = (T*) malloc( sizeof(T) * nx * ny );
		//stim::abs(mag, K, nx * ny);
		//stim::cpu2image<float>(mag, "kspace_pre_shift.bmp", nx, ny, stim::cmBrewer);
		
		size_t kxi, kyi;
		size_t i;
		T kx, kx_sq, ky, ky_sq, k_sq;
		T kz;
		stim::complex<T> shift;
		T min_kx = -Kx / 2;
		T dkx = Kx / (nx);

		T min_ky = -Ky / 2;
		T dky = Ky / (ny);

		for(kyi = 0; kyi < ny; kyi++){						//for each plane wave in the ky direction
			for(kxi = 0; kxi < nx; kxi++){					//for each plane wave in the ky direction
				i = kyi * nx + kxi;

				kx = min_kx + kxi * dkx;					//calculate the position of the current plane wave
				ky = min_ky + kyi * dky;

				kx_sq = kx * kx;
				ky_sq = ky * ky;
				k_sq = k*k;
				
				if(kx_sq + ky_sq < k_sq){
					kz = sqrt(k_sq - kx_sq - ky_sq);			//estimate kz using the Fresnel approximation				
					shift = -exp(stim::complex<T>(0, kz * z));
					K[i] *= shift;
					K[i] /= (nx*ny);							//normalize the DFT
				}
				else{
					K[i] = 0;
				}
			}
		}
		
		//stim::abs(mag, K, nx * ny);
		//stim::cpu2image<float>(mag, "kspace_post_shift.bmp", nx, ny, stim::cmBrewer);
		
		cpu_scalar_from_kspace(Enew, sx, sy, K, Kx, Ky, nx, ny);
		free(K);
	}

	template<typename T>
	void gpu_scalar_propagate(stim::complex<T>* Enew, stim::complex<T>* E, T sx, T sy, T z, T k, size_t nx, size_t ny){
		
		size_t field_bytes = sizeof(stim::complex<T>) * nx * ny;
		stim::complex<T>* host_E = (stim::complex<T>*) malloc( field_bytes);
		HANDLE_ERROR( cudaMemcpy(host_E, E, field_bytes, cudaMemcpyDeviceToHost) );

		stim::complex<T>* host_Enew = (stim::complex<T>*) malloc(field_bytes);

		cpu_scalar_propagate(host_Enew, host_E, sx, sy, z, k, nx, ny);

		HANDLE_ERROR( cudaMemcpy(Enew, host_Enew, field_bytes, cudaMemcpyHostToDevice) );
		free(host_E);
		free(host_Enew);
	}

	/// Apply a lowpass filter to a field slice
	/// @param Enew is the resulting propogated field
	/// @param E is the field to be propogated
	/// @param sx is the size of the field in the lateral x direction
	/// @param sy is the size of the field in the lateral y direction
	/// @param highest is the highest spatial frequency that can pass through the filter
	/// @param nx is the number of samples in the field along the lateral x direction
	/// @param ny is the number of samples in the field along the lateral y direction
	template<typename T>
	void cpu_scalar_lowpass(stim::complex<T>* Enew, stim::complex<T>* E, T sx, T sy, T highest, size_t nx, size_t ny){
		
		stim::complex<T>* K = (stim::complex<T>*) malloc( sizeof(stim::complex<T>) * nx * ny );

		T Kx, Ky;											//width and height in k space
		cpu_scalar_to_kspace(K, Kx, Ky, E ,sx, sy, nx, ny);

		//T* mag = (T*) malloc( sizeof(T) * nx * ny );
		//stim::abs(mag, K, nx * ny);
		//stim::cpu2image<float>(mag, "kspace_pre_lowpass.bmp", nx, ny, stim::cmBrewer);
		
		size_t kxi, kyi;
		size_t i;
		T kx, kx_sq, ky, ky_sq, k_sq;
		T kz;
		stim::complex<T> shift;
		T min_kx = -Kx / 2;
		T dkx = Kx / (nx);

		T min_ky = -Ky / 2;
		T dky = Ky / (ny);

		T highest_sq = highest * highest;

		for(kyi = 0; kyi < ny; kyi++){						//for each plane wave in the ky direction
			for(kxi = 0; kxi < nx; kxi++){					//for each plane wave in the ky direction
				i = kyi * nx + kxi;

				kx = min_kx + kxi * dkx;					//calculate the position of the current plane wave
				ky = min_ky + kyi * dky;

				kx_sq = kx * kx;
				ky_sq = ky * ky;
				
				if(kx_sq + ky_sq > highest_sq){
					K[i] = 0;
				}
				else
					K[i] /= nx * ny;						//normalize the DFT
			}
		}
		
		//stim::abs(mag, K, nx * ny);
		//stim::cpu2image<float>(mag, "kspace_post_lowpass.bmp", nx, ny, stim::cmBrewer);
		
		cpu_scalar_from_kspace(Enew, sx, sy, K, Kx, Ky, nx, ny);
		free(K);
	}

	enum locationType {CPUmem, GPUmem};

	/// Class represents a scalar optical field.

	/// In general, this class is designed to operate between the CPU and GPU. So, make sure all functions have an option to create the output on either.
	///		The field is stored *either* on the GPU or host memory, but not both. This enforces that there can't be different copies of the same field.
	///		This class is designed to be included in all of the other scalar optics classes, allowing them to render output data so make sure to keep it general and compatible.

template<typename T>
class scalarfield : public rect<T>{

protected:
	stim::complex<T>* E;
	size_t R[2];
	locationType loc;
	using rect<T>::X;
	using rect<T>::Y;

	T* p[3];											//scalar position for each point in E

	/// Convert the field to a k-space representation (do an FFT)
	void to_kspace(T& kx, T& ky){
		cpu_scalar_to_kspace(E, kx, ky, E, X.len(), Y.len(), R[0], R[1]);
	}

	void from_kspace(){
		kx = stim::TAU * R[0] / X.len();			//calculate the width of the momentum space
		ky = stim::TAU * R[1] / Y.len();
		T x, y;
		cpu_scalar_from_kspace(E, x, y, E, kx, ky, R[0], R[1]);
	}

public:

	/// Returns the number of values in the field
	CUDA_CALLABLE size_t size(){
		return R[0] * R[1];
	}

	CUDA_CALLABLE size_t grid_bytes(){
		return sizeof(stim::complex<T>) * R[0] * R[1];
	}

	scalarfield(size_t X, size_t Y, T size = 1, T z_pos = 0) : rect<T>::rect(size, z_pos){
		R[0] = X;											//set the field resolution
		R[1] = Y;

		E = (stim::complex<T>*) malloc(grid_bytes());		//allocate in CPU memory
		memset(E, 0, grid_bytes());
		loc = CPUmem;

		p[0] = p[1] = p[2] = NULL;							//set the position vector to NULL

	}

	~scalarfield(){
		if(loc == CPUmem) free(E);
		else cudaFree(E);
	}	

	/// Calculates the distance between points on the grid
	T spacing(){
		T du = rect<T>::X.len() / R[0];
		T dv = rect<T>::Y.len() / R[1];
		return min<T>(du, dv);
	}

	/// Copy the field array to the GPU, if it isn't already there
	void to_gpu(){
		if(loc == GPUmem) return;
		else{
			stim::complex<T>* dev_E;
			HANDLE_ERROR( cudaMalloc(&dev_E, grid_bytes()) );								//allocate GPU memory
			HANDLE_ERROR( cudaMemcpy(dev_E, E, grid_bytes(), cudaMemcpyHostToDevice) );	//copy the field to the GPU
			free(E);																	//free the CPU memory
			E = dev_E;																	//swap pointers

			if(p[0]){
				size_t meshgrid_bytes = size() * sizeof(T);								//calculate the number of bytes in each meshgrid
				T* dev_X;																//allocate variables to store the device meshgrid
				T* dev_Y;
				T* dev_Z;
				HANDLE_ERROR( cudaMalloc(&dev_X, meshgrid_bytes) );						//allocate space for the meshgrid on the device
				HANDLE_ERROR( cudaMalloc(&dev_Y, meshgrid_bytes) );
				HANDLE_ERROR( cudaMalloc(&dev_Z, meshgrid_bytes) );

				HANDLE_ERROR( cudaMemcpy(dev_X, p[0], meshgrid_bytes, cudaMemcpyHostToDevice) );	//copy from the host to the device
				HANDLE_ERROR( cudaMemcpy(dev_Y, p[1], meshgrid_bytes, cudaMemcpyHostToDevice) );
				HANDLE_ERROR( cudaMemcpy(dev_Z, p[2], meshgrid_bytes, cudaMemcpyHostToDevice) );

				free(p[0]);																//free device memory
				free(p[1]);
				free(p[2]);

				p[0] = dev_X;															//swap in the new pointers to device memory
				p[1] = dev_Y;
				p[2] = dev_Z;
			}
			loc = GPUmem;																//set the location flag
		}

	}

	/// Copy the field array to the CPU, if it isn't already there
	void to_cpu(){
		if(loc == CPUmem) return;
		else{
			stim::complex<T>* host_E = (stim::complex<T>*) malloc(grid_bytes());			//allocate space in main memory
			HANDLE_ERROR( cudaMemcpy(host_E, E, grid_bytes(), cudaMemcpyDeviceToHost) );	//copy from GPU to CPU
			HANDLE_ERROR( cudaFree(E) );												//free device memory
			E = host_E;																	//swap pointers

			//copy a meshgrid has been created
			if(p[0]){
				size_t meshgrid_bytes = size() * sizeof(T);								//move X to the CPU
				T* host_X = (T*) malloc( meshgrid_bytes );
				T* host_Y = (T*) malloc( meshgrid_bytes );
				T* host_Z = (T*) malloc( meshgrid_bytes );

				HANDLE_ERROR( cudaMemcpy(host_X, p[0], meshgrid_bytes, cudaMemcpyDeviceToHost) );
				HANDLE_ERROR( cudaMemcpy(host_Y, p[1], meshgrid_bytes, cudaMemcpyDeviceToHost) );
				HANDLE_ERROR( cudaMemcpy(host_Z, p[2], meshgrid_bytes, cudaMemcpyDeviceToHost) );

				HANDLE_ERROR( cudaFree(p[0]) );
				HANDLE_ERROR( cudaFree(p[1]) );
				HANDLE_ERROR( cudaFree(p[2]) );

				p[0] = host_X;
				p[1] = host_Y;
				p[2] = host_Z;
			}
			loc = CPUmem;
		}
	}

	bool gpu(){
		if(loc == GPUmem) return true;
		else return false;
	}

	/// Propagate the field along its orthogonal direction by a distance d
	void propagate(T d, T k){
		if(loc == CPUmem){
			cpu_scalar_propagate(E, E, X.len(), Y.len(), d, k, R[0], R[1]);
		}
		else{
			gpu_scalar_propagate(E, E, X.len(), Y.len(), d, k, R[0], R[1]);
		}
	}

	/// Propagate the field along its orthogonal direction by a distance d
	void lowpass(T highest){
		cpu_scalar_lowpass(E, E, X.len(), Y.len(), highest, R[0], R[1]);
	}

	/// Crop an image based on a given padding parameter (crop out the center)
	void crop(size_t padding, stim::scalarfield<T>& cropped){
		size_t Cx = R[0] / (2 * padding + 1);										//calculate the size of the cropped image based on the padding value
		size_t Cy = R[1] / (2 * padding + 1);

		if(cropped.R[0] != Cx || cropped.R[1] != Cy){
			std::cout<<"Error: cropped field resolution ("<<cropped.R[0]<<" x "<<cropped.R[1]<<") does not match the required resolution ("<<Cx<<" x "<<Cy<<")."<<std::endl;
			exit(1);
		}

		if(loc == CPUmem){
			cropped.to_cpu();										//make sure that the cropped image is on the CPU
			size_t x, y;
			size_t sx, sy, si, di;
			for(y = 0; y < Cy; y++){
				sy = y + Cy * padding;								//calculate the y-index into the source image
				for(x = 0; x < Cx; x++){
					sx = x + Cx * padding;							//calculate the x-index into the source image
					si = sy * R[0] + sx;							//calculate the 1D index into the source image
					di = y * Cx + x;
					cropped.E[di] = E[si];
				}
			}
		}
		else{
			cropped.to_gpu();										//make sure that the cropped image is also on the GPU
			gpu_crop2d<stim::complex<T>>(cropped.E, E, R[0], R[1], Cx * padding, Cy * padding, Cx, Cy);
		}
	}

	std::string str(){
		std::stringstream ss;
		ss<<rect<T>::str()<<std::endl;
		ss<<"[ "<<R[0]<<" x "<<R[1]<<" ]"<<std::endl;
		ss<<"location: ";
		if(loc == CPUmem) ss<<"CPU";
		else ss<<"GPU";

		ss<<endl;
		return ss.str();
	}

	stim::complex<T>* ptr(){
		return E;
	}

	T* x(){ return p[0]; }
	T* y(){ return p[1]; }
	T* z(){ return p[2]; }

	/// Evaluate the cartesian coordinates of each point in the field. The resulting arrays are allocated in the same memory where the field is stored.
	void meshgrid(T* X, T* Y, T* Z, locationType location){
		//size_t array_size = sizeof(T) * R[0] * R[1];
		if(location == CPUmem){

			T du = (T)1.0 / (R[0] - 1);					//calculate the spacing between points in the grid
			T dv = (T)1.0 / (R[1] - 1);

			size_t ui, vi, i;
			stim::vec3<T> p;
			for(vi = 0; vi < R[1]; vi++){
				i = vi * R[0];
				for(ui = 0; ui < R[0]; ui++){
					p = rect<T>::p(ui * du, vi * dv);
					X[i] = p[0];
					Y[i] = p[1];
					Z[i] = p[2];
					i++;					
				}
			}
			//stim::cpu2image(X, "X.bmp", R[0], R[1], stim::cmBrewer);
			//stim::cpu2image(Y, "Y.bmp", R[0], R[1], stim::cmBrewer);
			//stim::cpu2image(Z, "Z.bmp", R[0], R[1], stim::cmBrewer);
		}
		else{
			std::cout<<"GPU allocation of a meshgrid isn't supported yet. You'll have to write kernels to do the calculation.";
			exit(1);
		}
	}

	/// Create a local meshgrid
	void meshgrid(){
		if(p[0]) return;								//if the p[0] value is not NULL, a meshgrid has already been created
		if(loc == CPUmem){
			p[0] = (T*) malloc( size() * sizeof(T) );
			p[1] = (T*) malloc( size() * sizeof(T) );
			p[2] = (T*) malloc( size() * sizeof(T) );
		}
		else{
			std::cout<<"GPUmem meshgrid isn't implemented yet."<<std::endl;
			exit(1);
		}
		meshgrid(p[0], p[1], p[2], loc);
	}

	//clear the field, setting all values to zero
	void clear(){
		if(loc == GPUmem)
			HANDLE_ERROR(cudaMemset(E, 0, grid_bytes()));
		else
			memset(E, 0, grid_bytes());
	}

	void image(std::string filename, stim::complexComponentType type = complexMag, stim::colormapType cmap = stim::cmBrewer){

		if(loc == GPUmem){
			T* image;
			HANDLE_ERROR( cudaMalloc(&image, sizeof(T) * size()) );
			int threads = stim::maxThreadsPerBlock();												//get the maximum number of threads per block for the CUDA device
			dim3 blocks( R[0] * R[1] / threads + 1 );												//create a 1D array of blocks

			switch(type){
			case complexMag:
				cuda_abs<T><<< blocks, threads >>>(image, E, size());
				break;
			case complexReal:
				cuda_real<T><<< blocks, threads >>>(image, E, size());
				break;
			case complexImaginary:
				cuda_imag<T><<< blocks, threads >>>(image, E, size());
				break;
			case complexIntensity:
				cuda_intensity<T><<< blocks, threads >>>(image, E, size());
				break;
			}
			stim::gpu2image<T>(image, filename, R[0], R[1], stim::cmBrewer);
			HANDLE_ERROR( cudaFree(image) );
		}
		else{
			T* image = (T*) malloc( sizeof(T) * size() );				//allocate space for the real image

			switch(type){												//get the specified component from the complex value
			case complexMag:
				stim::abs(image, E, size());
				break;
			case complexReal:
				stim::real(image, E, size());
				break;
			case complexImaginary:
				stim::imag(image, E, size());
				break;
			case complexIntensity:
				stim::intensity(image, E, size());
				break;
			}
			stim::cpu2image(image, filename, R[0], R[1], cmap);			//save the resulting image
			free(image);												//free the real image
		}
	}

	void image(T* img, stim::complexComponentType type = complexMag){
		if(loc == GPUmem) to_cpu();									//if the field is in the GPU, move it to the CPU

		switch(type){												//get the specified component from the complex value
		case complexMag:
			stim::abs(img, E, size());
			break;
		case complexReal:
			stim::real(img, E, size());
			break;
		case complexImaginary:
			stim::imag(img, E, size());
			break;
		case complexIntensity:
			stim::intensity(img, E, size());
			break;
		}
		//stim::cpu2image(image, filename, R[0], R[1], cmap);			//save the resulting image
		//free(image);												//free the real image
	}

	//adds the field intensity to the output array (useful for calculating detector response to incoherent fields)
	void intensity(T* out){
		if(loc == GPUmem){
			//T* image;
			//HANDLE_ERROR( cudaMalloc(&image, sizeof(T) * size()) );
			int threads = stim::maxThreadsPerBlock();												//get the maximum number of threads per block for the CUDA device
			dim3 blocks( R[0] * R[1] / threads + 1 );												//create a 1D array of blocks
			cuda_sum_intensity<T><<< blocks, threads >>>(out, E, size());
		}
		else{
			T* image = (T*) malloc( sizeof(T) * size() );				//allocate space for the real image
			stim::intensity(image, E, size());							//calculate the intensity
		
			size_t N = size();											//calculate the number of elements in the field
			for(size_t n = 0; n < N; n++)								//for each point in the field
				out[n] += image[n];										//add the field intensity to the output image

			free(image);												//free the temporary intensity image
		}
	}

};				//end class scalarfield
}

//stream insertion operator
template<typename T>
std::ostream& operator<<(std::ostream& os, stim::scalarfield<T>& rhs){
	os<<rhs.str();
	return os;
}


#endif