Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /trunk/test/MIP/mip_opencl.c
ViewVC logotype

Annotation of /trunk/test/MIP/mip_opencl.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 435 - (view) (download) (as text)

1 : lamonts 302 /* mip_opencl.c
2 : lamonts 203 *
3 : jhr 435 * COPYRIGHT (c) 2010 The Diderot Project (http://diderot-language.cs.uchicago.edu)
4 : lamonts 302 * All rights reserved.
5 :     *
6 :     * An OpenCL mip implementation.
7 :     *
8 : lamonts 203 * To View the Image
9 :     * =========================
10 : lamonts 303 * unu reshape -i mip.txt -s 200 200 | unu quantize -b 8 -o new.png
11 : lamonts 203 */
12 : lamonts 177 #include <OpenCL/OpenCl.h>
13 :     #include <assert.h>
14 :     #include <stdio.h>
15 :     #include <stdlib.h>
16 :     #include <sys/sysctl.h>
17 :     #include <sys/stat.h>
18 :     #include <teem/nrrd.h>
19 :    
20 : jhr 310 #define SIZE 640
21 : lamonts 177
22 :     //Loads the Kernel from a file
23 :     char * loadKernel (const char * filename)
24 :     {
25 :     struct stat statbuf;
26 :     FILE *fh;
27 :     char *source;
28 :    
29 :     fh = fopen(filename, "r");
30 :     if (fh == 0)
31 :     return 0;
32 :    
33 :     stat(filename, &statbuf);
34 :     source = (char *) malloc(statbuf.st_size + 1);
35 :     fread(source, statbuf.st_size, 1, fh);
36 :     source[statbuf.st_size] = '\0';
37 :    
38 :     return source;
39 :     }
40 : lamonts 203 void saveResults (float * matrix, int size)
41 :     {
42 :     int i;
43 :     float max = -INFINITY;
44 :     FILE * out_file;
45 :     out_file = fopen("mip.txt", "w");
46 :     if (out_file == NULL) {
47 :     fprintf(stderr,"Can not open output file\n");
48 :     exit (8);
49 :     }
50 :    
51 :     for(i = 0; i < size; i++)
52 :     {
53 :     if(matrix[i] == -INFINITY || matrix[i] < 0)
54 : lamonts 302 fprintf(out_file,"%.4f\n",0.0f);
55 : lamonts 203 else
56 : lamonts 302 fprintf(out_file,"%.4f\n",matrix[i]);
57 : lamonts 203
58 :     if(matrix[i] > max)
59 :     max = matrix[i];
60 :    
61 :     }
62 :     fclose(out_file);
63 :    
64 :    
65 :     }
66 : lamonts 177 float det3x3(float a, float b, float c, float d, float e, float f, float g, float h, float i)
67 :     {
68 :     return ( (a)*(e)*(i)
69 :     + (d)*(h)*(c)
70 :     + (g)*(b)*(f)
71 :     - (g)*(e)*(c)
72 :     - (d)*(b)*(i)
73 :     - (a)*(h)*(f));
74 :     }
75 :     float det4x4(cl_float16 m)
76 :     {
77 :     return (m[ 0] * det3x3(m[ 5], m[ 6], m[ 7],
78 :     m[ 9], m[10], m[11],
79 :     m[13], m[14], m[15])
80 :    
81 :     - m[ 1] * det3x3(m[ 4], m[ 6], m[ 7],
82 :     m[ 8], m[10], m[11],
83 :     m[12], m[14], m[15])
84 :     + m[ 2] * det3x3(m[ 4], m[ 5], m[ 7],
85 :     m[ 8], m[ 9], m[11],
86 :     m[12], m[13], m[15])
87 :    
88 :     - m[ 3] * det3x3(m[ 4], m[ 5], m[ 6],
89 :     m[ 8], m[ 9], m[10],
90 :     m[12], m[13], m[14]));
91 :    
92 :    
93 :    
94 :     }
95 :     void invMatrix(cl_float16 m, cl_float16 i)
96 :     {
97 :     float det = det4x4(m);
98 :    
99 :    
100 :     i[0] = det3x3(m[5],m[ 6],m[ 7],
101 :     m[ 9],m[10],m[11],
102 :     m[13],m[14],m[15])/det;
103 :    
104 :     i[ 1] = -det3x3(m[ 1],m[ 2],m[ 3],
105 :     m[ 9],m[10],m[11],
106 :     m[13],m[14],m[15])/det;
107 :    
108 :     i[ 2] = det3x3(m[ 1],m[ 2],m[ 3],
109 :     m[ 5],m[ 6],m[ 7],
110 :     m[13],m[14],m[15])/det;
111 :    
112 :     i[ 3] = -det3x3(m[ 1],m[ 2],m[ 3],
113 :     m[ 5],m[ 6],m[ 7],
114 :     m[ 9],m[10],m[11])/det;
115 :    
116 :     i[ 4] = -det3x3(m[ 4],m[ 6],m[ 7],
117 :     m[ 8],m[10],m[11],
118 :     m[12],m[14],m[15])/det;
119 :    
120 :     i[ 5] = det3x3(m[ 0],m[ 2],m[ 3],
121 :     m[ 8],m[10],m[11],
122 :     m[12],m[14],m[15])/det;
123 :    
124 :     i[ 6] = -det3x3(m[ 0],m[ 2],m[ 3],
125 :     m[ 4],m[ 6],m[ 7],
126 :     m[12],m[14],m[15])/det;
127 :    
128 :     i[ 7] = det3x3(m[ 0],m[ 2],m[ 3],
129 :     m[ 4],m[ 6],m[ 7],
130 :     m[ 8],m[10],m[11])/det;
131 :    
132 :     i[ 8] = det3x3(m[ 4],m[ 5],m[ 7],
133 :     m[ 8],m[ 9],m[11],
134 :     m[12],m[13],m[15])/det;
135 :    
136 :     i[ 9] = -det3x3(m[ 0],m[ 1],m[ 3],
137 :     m[ 8],m[ 9],m[11],
138 :     m[12],m[13],m[15])/det;
139 :    
140 :     i[10] = det3x3(m[ 0],m[ 1],m[ 3],
141 :     m[ 4],m[ 5],m[ 7],
142 :     m[12],m[13],m[15])/det;
143 :    
144 :     i[11] = -det3x3(m[ 0],m[ 1],m[ 3],
145 :     m[ 4],m[ 5],m[ 7],
146 :     m[ 8],m[ 9],m[11])/det;
147 :    
148 :     i[12] = -det3x3(m[ 4],m[ 5],m[ 6],
149 :     m[ 8],m[ 9],m[10],
150 :     m[12],m[13],m[14])/det;
151 :    
152 :     i[13] = det3x3(m[ 0],m[ 1],m[ 2],
153 :     m[ 8],m[ 9],m[10],
154 :     m[12],m[13],m[14])/det;
155 :    
156 :     i[14] = -det3x3(m[ 0],m[ 1],m[ 2],
157 :     m[ 4],m[ 5],m[ 6],
158 :     m[12],m[13],m[14])/det;
159 :    
160 :     i[15] = det3x3(m[ 0],m[ 1],m[ 2],
161 :     m[ 4],m[ 5],m[ 6],
162 :     m[ 8],m[ 9],m[10])/det;
163 :     }
164 :     void printMatrix(float * matrix, int rowSize)
165 :     {
166 :     int index = 0, end = 1, arraySize = rowSize * rowSize;
167 :    
168 : jhr 218 for(index = 0; index < arraySize; index++)
169 : lamonts 177 {
170 : jhr 218 if(end == rowSize)
171 : lamonts 177 {
172 :     printf(" %.2f\n",matrix[index]);
173 :     end = 1;
174 :     }
175 :     else
176 :     {
177 :     printf(" %.2f ",matrix[index]);
178 :     end++;
179 :     }
180 :     }
181 :     printf("\n");
182 :     }
183 :     void loadTransformMatrix(Nrrd * nin, cl_float16 transformMatrix)
184 :     {
185 :     int i,j, size = nin->spaceDim;
186 :     NrrdAxisInfo axisInfo;
187 :    
188 :     //Image axis Scaling and Rotation
189 :     for(i = 0; i < size; i++)
190 :     {
191 :     axisInfo = nin->axis[i];
192 :     for(j = 0; j < size; j++)
193 :     {
194 :     transformMatrix[ (size+ 1) * j + i] = axisInfo.spaceDirection[j];
195 :     }
196 :    
197 :     //Image Location
198 :     transformMatrix[ (i * (size + 1)) + size] = nin->spaceOrigin[i];
199 :    
200 :     //Bottom row of the Transform Matrix
201 :     transformMatrix[((size + 1) * (size)) + i ] = 0;
202 :     }
203 : jhr 218 transformMatrix[((size + 1) * (size)) + size ] = 1;
204 :    
205 : lamonts 177 }
206 :     Nrrd * loadNrrdFile(char * filename)
207 :     {
208 :     /* create a nrrd; at this point this is just an empty container */
209 :     Nrrd * nin;
210 :    
211 :     nin = nrrdNew();
212 :     char *err;
213 :    
214 :     /* read in the nrrd from file */
215 :     if (nrrdLoad(nin, filename, NULL)) {
216 :     err = biffGetDone(NRRD);
217 :     fprintf(stderr, "Mip: trouble reading \"%s\":\n%s", filename, err);
218 :     free(err);
219 :     return NULL;
220 :     }
221 :    
222 :     /* say something about the array
223 :     printf("Mip: \"%s\" is a %d-dimensional nrrd of type %d (%s)\n",
224 :     filename, nin->dim, nin->type,
225 :     airEnumStr(nrrdType, nin->type));
226 :     printf("Mip: the array contains %d elements, each %d bytes in size\n",
227 :     (int)nrrdElementNumber(nin), (int)nrrdElementSize(nin));*/
228 :    
229 :     return nin;
230 :    
231 :     }
232 : lamonts 302 int exe_MIP_Kernel(float * img, int imageDataSize, float inverseMatrix[16], int sAxis[3], float * out)
233 : lamonts 177 {
234 :    
235 :     cl_program program;
236 :     cl_kernel kernel;
237 :    
238 :     cl_command_queue queue;
239 :     cl_context context;
240 :    
241 :     cl_device_id cpu = NULL, device = NULL;
242 :    
243 :     cl_int err = 0;
244 : lamonts 302
245 : lamonts 177
246 : lamonts 302 cl_mem imageData_mem, out_mem,sAxis_mem;
247 : lamonts 177
248 :     /** Setup Device **/
249 :     err = clGetDeviceIDs(NULL,CL_DEVICE_TYPE_CPU,1,&cpu,NULL);
250 :     assert(err==CL_SUCCESS);
251 :    
252 :     err = clGetDeviceIDs(NULL,CL_DEVICE_TYPE_GPU,1,&device,NULL);
253 :     //if(err != CL_SUCCESS)
254 :     device = cpu;
255 :    
256 :     assert(device);
257 :    
258 :     /* Setup Context and Command Queue */
259 :     context = clCreateContext(0,1,&device,NULL,NULL,&err);
260 :     assert(err == CL_SUCCESS);
261 :    
262 :     queue = clCreateCommandQueue(context,device,0,NULL);
263 :    
264 :     /** Load the Kernel and Program **/
265 :     const char * filename = "mip.cl";
266 :     char * kernel_source = loadKernel(filename);
267 :    
268 :     assert(kernel_source != 0);
269 :    
270 :     program = clCreateProgramWithSource(context,1,(const char **)&kernel_source,NULL,&err);
271 :    
272 :     assert(err == CL_SUCCESS);
273 :    
274 :     err = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
275 :    
276 :    
277 :     /** Retrieve information about the program build to check for any possible errors **/
278 :     char * build_log;
279 :     size_t log_size;
280 :    
281 :     // First call to know the proper size
282 :     clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
283 :     build_log = (char *) malloc(log_size+1);
284 :     // Second call to get the log
285 :     clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, log_size, build_log, NULL);
286 :     build_log[log_size] = '\0';
287 :     printf("\nBuild Log:\n%s\n",build_log);
288 :     free(build_log);
289 :    
290 :     assert(err == CL_SUCCESS);
291 :    
292 : lamonts 302 kernel = clCreateKernel(program,"mip",&err);
293 : lamonts 177
294 :     assert(err == CL_SUCCESS);
295 :    
296 :     /** Memory Allocation for the Matrices **/
297 :    
298 : lamonts 191 imageData_mem = clCreateBuffer(context,CL_MEM_READ_ONLY,sizeof(float) * imageDataSize,NULL,NULL);
299 :     err |= clEnqueueWriteBuffer(queue,imageData_mem,CL_TRUE,0,sizeof(float) * imageDataSize,
300 : lamonts 302 (void *)img ,0,NULL,NULL);
301 : lamonts 177
302 : lamonts 302 sAxis_mem = clCreateBuffer(context,CL_MEM_READ_ONLY,sizeof(int) * 3,NULL,NULL);
303 :     err |= clEnqueueWriteBuffer(queue,sAxis_mem,CL_TRUE,0,sizeof(int) * 3,
304 :     (void *)sAxis ,0,NULL,NULL);
305 : lamonts 177
306 :     assert(err == CL_SUCCESS);
307 :    
308 :     out_mem = clCreateBuffer(context,CL_MEM_READ_WRITE, sizeof(float) * (SIZE *SIZE),NULL,NULL);
309 :    
310 :     clFinish(queue);
311 :    
312 :     size_t global_work_size[2], local_work_size[2];
313 :    
314 : jhr 310 global_work_size[0] = SIZE;
315 :     global_work_size[1] = SIZE;
316 : lamonts 177
317 :     local_work_size[0] = 1;
318 :     local_work_size[1] = 1;
319 :    
320 : lamonts 302 cl_int2 workDim = {SIZE,SIZE};
321 : lamonts 177
322 : lamonts 302 int index = 0;
323 :    
324 :     err =clSetKernelArg(kernel,index++,sizeof(cl_mem), &imageData_mem);
325 :     err |=clSetKernelArg(kernel,index++,sizeof(cl_mem), &out_mem);
326 :     err |=clSetKernelArg(kernel,index++,sizeof(cl_float16), inverseMatrix);
327 :     err |=clSetKernelArg(kernel,index++,sizeof(cl_int2), &workDim);
328 :     err |=clSetKernelArg(kernel,index++,sizeof(cl_mem), &sAxis_mem);
329 :    
330 : lamonts 177 assert(err == CL_SUCCESS);
331 : lamonts 302
332 : lamonts 177
333 :     err = clEnqueueNDRangeKernel(queue,kernel,2,NULL,global_work_size,
334 :     local_work_size,0,NULL,NULL);
335 :    
336 :     assert(err == CL_SUCCESS);
337 :    
338 :     clFinish(queue);
339 :    
340 :     err = clEnqueueReadBuffer(queue,out_mem,CL_TRUE,0, sizeof(float) * (SIZE *SIZE),out,0,NULL,NULL);
341 : lamonts 191
342 : lamonts 203 saveResults(out,SIZE * SIZE);
343 : lamonts 177
344 :     clReleaseKernel(kernel);
345 :     clReleaseProgram(program);
346 :     clReleaseCommandQueue(queue);
347 :     clReleaseContext(context);
348 :     clReleaseMemObject(imageData_mem);
349 :     clReleaseMemObject(out_mem);
350 :    
351 :     return CL_SUCCESS;
352 : jhr 218 }
353 : jhr 310
354 : lamonts 177 int main (int argc, char ** argv)
355 :     {
356 : jhr 310 char * dataFile = "../../data/txs.nrrd";
357 : lamonts 302 float transformMatrix[16];
358 :     float inverseMatrix[16];
359 :    
360 : jhr 310 if (argc == 2) {
361 :     dataFile = argv[1];
362 :     }
363 :    
364 : lamonts 177 float * out;
365 :     out = (float *) malloc(sizeof(float) * (SIZE * SIZE));
366 : lamonts 302 Nrrd * nin = loadNrrdFile(dataFile);
367 : lamonts 177
368 : lamonts 302 int sAxis[] = {nin->axis[0].size, nin->axis[1].size,nin->axis[2].size};
369 :    
370 :     loadTransformMatrix(nin,transformMatrix);
371 :     invMatrix(transformMatrix,inverseMatrix);
372 :    
373 : lamonts 177
374 : lamonts 302 exe_MIP_Kernel((float *)nin->data, (int)nrrdElementNumber(nin),inverseMatrix, sAxis, out);
375 : lamonts 177
376 :     return 0;
377 : jhr 435 }

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0