/* mip_opencl.c
 *
 * COPYRIGHT (c) 2010 The Diderot Project (http://diderot.cs.uchicago.edu)
 * All rights reserved.
 *
 * An OpenCL mip implementation.
 *
 * To View the Image 
 * =========================
 * unu reshape -i mip.txt -s 200 200 | unu quantize -b 8 -o new.png
 */
#include <OpenCL/OpenCl.h> 
#include <assert.h> 
#include <stdio.h> 
#include <stdlib.h> 
#include <sys/sysctl.h>
#include <sys/stat.h>
#include <teem/nrrd.h> 

#define SIZE 200

//Loads the Kernel from a file 
char * loadKernel (const char * filename) 
{  
	struct stat statbuf;
	FILE *fh; 
	char *source; 
	
	fh = fopen(filename, "r");
	if (fh == 0)
		return 0; 
	
	stat(filename, &statbuf);
	source = (char *) malloc(statbuf.st_size + 1);
	fread(source, statbuf.st_size, 1, fh);
	source[statbuf.st_size] = '\0'; 
	
	return source; 
} 
void saveResults (float * matrix, int size) 
{ 
	int i; 
	float max = -INFINITY; 
	FILE * out_file; 
	out_file = fopen("mip.txt", "w");
    if (out_file == NULL) {
        fprintf(stderr,"Can not open output file\n");
        exit (8);
    }
    
    for(i = 0; i < size; i++) 
    {
    	if(matrix[i] == -INFINITY || matrix[i] < 0)
    	   fprintf(out_file,"%.4f\n",0.0f); 
    	else
    	   fprintf(out_file,"%.4f\n",matrix[i]); 
    	   
    	 if(matrix[i] > max) 
    	 	max = matrix[i]; 
   
   }
   fclose(out_file);


} 
float det3x3(float a, float b, float c, float d, float e, float f, float g, float h, float i) 
{
   return ( (a)*(e)*(i) 
   		  + (d)*(h)*(c) 
   		  + (g)*(b)*(f) 
   		  - (g)*(e)*(c) 
   		  - (d)*(b)*(i) 
   		  - (a)*(h)*(f)); 
}
float det4x4(cl_float16 m) 
{ 
	return	(m[ 0] * det3x3(m[ 5], m[ 6], m[ 7], 
						   m[ 9], m[10], m[11], 
						   m[13], m[14], m[15]) 
						   
	- m[ 1] * det3x3(m[ 4], m[ 6], m[ 7], 
						   m[ 8], m[10], m[11], 
						   m[12], m[14], m[15]) 
	+ m[ 2] * det3x3(m[ 4], m[ 5], m[ 7], 
						   m[ 8], m[ 9], m[11], 
						   m[12], m[13], m[15]) 
						   
	- m[ 3] * det3x3(m[ 4], m[ 5], m[ 6], 
						   m[ 8], m[ 9], m[10], 
						   m[12], m[13], m[14])); 



} 
void invMatrix(cl_float16 m, cl_float16 i) 
{  
	float det = det4x4(m); 

	
	i[0] = det3x3(m[5],m[ 6],m[ 7], 
					   m[ 9],m[10],m[11], 
					   m[13],m[14],m[15])/det; 
	
	i[ 1] = -det3x3(m[ 1],m[ 2],m[ 3], 
					   m[ 9],m[10],m[11], 
					   m[13],m[14],m[15])/det; 
					   
	i[ 2] =  det3x3(m[ 1],m[ 2],m[ 3], 
					   m[ 5],m[ 6],m[ 7], 
					   m[13],m[14],m[15])/det; 
					   
	i[ 3] = -det3x3(m[ 1],m[ 2],m[ 3], 
					   m[ 5],m[ 6],m[ 7], 
					   m[ 9],m[10],m[11])/det; 
					   
	i[ 4] = -det3x3(m[ 4],m[ 6],m[ 7], 
					   m[ 8],m[10],m[11], 
					   m[12],m[14],m[15])/det; 
					   
	i[ 5] =  det3x3(m[ 0],m[ 2],m[ 3], 
					   m[ 8],m[10],m[11], 
					   m[12],m[14],m[15])/det; 
					   
	i[ 6] = -det3x3(m[ 0],m[ 2],m[ 3], 
					   m[ 4],m[ 6],m[ 7], 
					   m[12],m[14],m[15])/det; 
					   
	i[ 7] =  det3x3(m[ 0],m[ 2],m[ 3], 
					   m[ 4],m[ 6],m[ 7], 
					   m[ 8],m[10],m[11])/det; 
					   
	i[ 8] =  det3x3(m[ 4],m[ 5],m[ 7], 
					   m[ 8],m[ 9],m[11], 
					   m[12],m[13],m[15])/det; 
					   
	i[ 9] = -det3x3(m[ 0],m[ 1],m[ 3], 
					   m[ 8],m[ 9],m[11], 
					   m[12],m[13],m[15])/det; 
					   
	i[10] =  det3x3(m[ 0],m[ 1],m[ 3], 
					   m[ 4],m[ 5],m[ 7], 
					   m[12],m[13],m[15])/det; 
					   
	i[11] = -det3x3(m[ 0],m[ 1],m[ 3], 
					   m[ 4],m[ 5],m[ 7], 
					   m[ 8],m[ 9],m[11])/det; 
					   
	i[12] = -det3x3(m[ 4],m[ 5],m[ 6], 
					   m[ 8],m[ 9],m[10], 
					   m[12],m[13],m[14])/det; 
					   
	i[13] =  det3x3(m[ 0],m[ 1],m[ 2], 
					   m[ 8],m[ 9],m[10],
					   m[12],m[13],m[14])/det; 
					   
	i[14] = -det3x3(m[ 0],m[ 1],m[ 2], 
					   m[ 4],m[ 5],m[ 6], 
					   m[12],m[13],m[14])/det; 
					   
	i[15] =  det3x3(m[ 0],m[ 1],m[ 2], 
					   m[ 4],m[ 5],m[ 6], 
					   m[ 8],m[ 9],m[10])/det; 
} 
void printMatrix(float * matrix, int rowSize) 
{
	int index = 0, end = 1, arraySize = rowSize * rowSize; 
	
	for(index = 0; index < arraySize; index++) 
	{ 
		if(end == rowSize) 
		{ 
			printf(" %.2f\n",matrix[index]);
			end = 1;
		} 
		else
		{
			printf(" %.2f ",matrix[index]); 
			end++; 
		} 
	}
	printf("\n"); 
}
void loadTransformMatrix(Nrrd * nin, cl_float16 transformMatrix) 
{ 
	int i,j, size =  nin->spaceDim; 
	NrrdAxisInfo axisInfo; 

	//Image axis Scaling and Rotation 
	for(i = 0; i < size; i++) 
	{ 
		axisInfo = nin->axis[i];
		for(j = 0; j < size; j++) 
		{
	     	transformMatrix[ (size+ 1) * j + i] = axisInfo.spaceDirection[j]; 
	    }
	  
		//Image Location 
	    transformMatrix[ (i * (size + 1)) + size] = nin->spaceOrigin[i];  
	    
	    //Bottom row of the Transform Matrix 
	    transformMatrix[((size + 1) * (size)) + i ] = 0; 
	}
	transformMatrix[((size + 1) * (size)) + size ] = 1;

} 
Nrrd * loadNrrdFile(char * filename) 
{
	/* create a nrrd; at this point this is just an empty container */
  	Nrrd * nin; 
  	
  	nin = nrrdNew();
	char *err;

  	/* read in the nrrd from file */
  	if (nrrdLoad(nin, filename, NULL)) {
    	err = biffGetDone(NRRD);
    	fprintf(stderr, "Mip: trouble reading \"%s\":\n%s", filename, err);
    	free(err);
   	    return NULL;
   }

  /* say something about the array 
  printf("Mip: \"%s\" is a %d-dimensional nrrd of type %d (%s)\n", 
         filename, nin->dim, nin->type,
         airEnumStr(nrrdType, nin->type));
  printf("Mip: the array contains %d elements, each %d bytes in size\n",
         (int)nrrdElementNumber(nin), (int)nrrdElementSize(nin));*/ 
         
  return nin; 
  
} 
int exe_MIP_Kernel(float * img, int imageDataSize, float inverseMatrix[16], int sAxis[3], float * out)  
{ 

	cl_program program; 
	cl_kernel kernel; 

	cl_command_queue queue; 
	cl_context context; 
	
	cl_device_id cpu = NULL, device = NULL; 
	
	cl_int err = 0; 

	
	cl_mem imageData_mem, out_mem,sAxis_mem; 
	
	/** Setup Device **/ 
	err = clGetDeviceIDs(NULL,CL_DEVICE_TYPE_CPU,1,&cpu,NULL); 
	assert(err==CL_SUCCESS); 
	
	err = clGetDeviceIDs(NULL,CL_DEVICE_TYPE_GPU,1,&device,NULL); 
	//if(err != CL_SUCCESS)  
	device = cpu; 
	
	assert(device); 
	
	/* Setup Context and Command Queue */ 
	context = clCreateContext(0,1,&device,NULL,NULL,&err); 
	assert(err == CL_SUCCESS); 
	
	queue = clCreateCommandQueue(context,device,0,NULL); 
	
	/** Load the Kernel and Program **/ 
	const char * filename = "mip.cl"; 
	char * kernel_source = loadKernel(filename); 
	
	assert(kernel_source != 0); 
	
	program = clCreateProgramWithSource(context,1,(const char **)&kernel_source,NULL,&err); 
	
	assert(err == CL_SUCCESS); 
			
	err = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);


   /** Retrieve information about the program build to check for any possible errors **/ 
   char * build_log;
   size_t log_size;
   
   // First call to know the proper size
   clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
   build_log = (char *) malloc(log_size+1);
   // Second call to get the log
   clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, log_size, build_log, NULL);
   build_log[log_size] = '\0';
   printf("\nBuild Log:\n%s\n",build_log); 
   free(build_log); 
   
    assert(err == CL_SUCCESS); 
   
	kernel = clCreateKernel(program,"mip",&err); 
		
	assert(err == CL_SUCCESS); 
	
	/** Memory Allocation for the Matrices **/    
	
	imageData_mem = clCreateBuffer(context,CL_MEM_READ_ONLY,sizeof(float) * imageDataSize,NULL,NULL); 
	err |= clEnqueueWriteBuffer(queue,imageData_mem,CL_TRUE,0,sizeof(float) * imageDataSize, 
								(void *)img ,0,NULL,NULL); 	
	
	sAxis_mem = clCreateBuffer(context,CL_MEM_READ_ONLY,sizeof(int) * 3,NULL,NULL); 
	err |= clEnqueueWriteBuffer(queue,sAxis_mem,CL_TRUE,0,sizeof(int) * 3, 
								(void *)sAxis ,0,NULL,NULL); 	
	
	assert(err == CL_SUCCESS); 
	
	out_mem = clCreateBuffer(context,CL_MEM_READ_WRITE, sizeof(float) * (SIZE *SIZE),NULL,NULL); 
	
	clFinish(queue); 
	
	size_t global_work_size[2], local_work_size[2]; 
	
	global_work_size[0] = 256; 
	global_work_size[1] = 256; 
	
	local_work_size[0] = 1; 
	local_work_size[1] = 1; 
	
	cl_int2 workDim = {SIZE,SIZE}; 
	
	int index = 0; 
	
	err  =clSetKernelArg(kernel,index++,sizeof(cl_mem), &imageData_mem); 
	err |=clSetKernelArg(kernel,index++,sizeof(cl_mem), &out_mem); 
	err |=clSetKernelArg(kernel,index++,sizeof(cl_float16), inverseMatrix); 
	err |=clSetKernelArg(kernel,index++,sizeof(cl_int2), &workDim); 
	err |=clSetKernelArg(kernel,index++,sizeof(cl_mem), &sAxis_mem);

	assert(err == CL_SUCCESS); 
		
	
	err = clEnqueueNDRangeKernel(queue,kernel,2,NULL,global_work_size,
								local_work_size,0,NULL,NULL); 
	
	assert(err == CL_SUCCESS); 
		
	clFinish(queue); 
	
	err = clEnqueueReadBuffer(queue,out_mem,CL_TRUE,0, sizeof(float) * (SIZE *SIZE),out,0,NULL,NULL); 
	 
	saveResults(out,SIZE * SIZE); 
	  
	clReleaseKernel(kernel); 
	clReleaseProgram(program); 
	clReleaseCommandQueue(queue); 
	clReleaseContext(context); 
	clReleaseMemObject(imageData_mem); 
	clReleaseMemObject(out_mem); 
	
	return CL_SUCCESS; 	
}
int main (int argc, char ** argv) 
{ 
	char * dataFile = "txs.nrrd"; 
	float transformMatrix[16]; 
	float inverseMatrix[16]; 

	float * out; 
	out = (float *) malloc(sizeof(float) * (SIZE * SIZE)); 
	Nrrd * nin = loadNrrdFile(dataFile); 

	int sAxis[] = {nin->axis[0].size, nin->axis[1].size,nin->axis[2].size};  
	
	loadTransformMatrix(nin,transformMatrix); 
	invMatrix(transformMatrix,inverseMatrix); 
	

	exe_MIP_Kernel((float *)nin->data, (int)nrrdElementNumber(nin),inverseMatrix, sAxis, out);  
	
	return 0; 
}

<script type="text/javascript">//<![CDATA[
		function toggle_ffErrors() {
			var errorsblock = document.getElementById("ffErrorsBlock");
			var errorsgroup = document.getElementById("ffErrors");
			if (errorsblock.style.display == "none") {
				errorsblock.style.display = "block";
				errorsgroup.style.right = "10px";
			} else {
				errorsblock.style.display = "none";
				errorsgroup.style.right = "300px";
			}
		}
//]]></script>
<div id="ffErrors">
<a href="javascript:toggle_ffErrors();">Click to toggle</a>
<div id="ffErrorsBlock">
	<div class="error">does not end with &lt;/html&gt; tag</div>
	<div class="error">does not end with &lt;/body&gt; tag</div>
	<div class="info">The output has ended thus: Kernel((float *)nin-&gt;data, (int)nrrdElementNumber(nin),inverseMatrix, sAxis, out);  
	
	return 0; 
}</div>
</div></div>
</body></html>