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

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

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