/* 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 #include #include #include #include #include #include #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; }
Click to toggle
does not end with </html> tag
does not end with </body> tag
The output has ended thus: Kernel((float *)nin->data, (int)nrrdElementNumber(nin),inverseMatrix, sAxis, out); return 0; }