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 3349 - (view) (download) (as text)

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

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