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

SCM Repository

[diderot] Diff of /branches/pure-cfg/test/MIP/mip.cl
ViewVC logotype

Diff of /branches/pure-cfg/test/MIP/mip.cl

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 218, Tue Aug 3 19:40:56 2010 UTC revision 441, Wed Oct 20 19:21:32 2010 UTC
# Line 1  Line 1 
1  __kernel void raycast ( float * img,  /* mip.cl
2                                              float * h1,   *
3                                              float * h2,   * COPYRIGHT (c) 2010 The Diderot Project (http://diderot-language.cs.uchicago.edu)
4     * All rights reserved.
5     *
6     * An OpenCL file containing a probe kernel and mip kernel implementation
7     */
8    
9    __constant float h[4][4] = {            // bspln3
10            { 1.33333,   2.0,  1.0,  0.166667 },   // -2 .. -1
11            { 0.666667,  0.0, -1.0, -0.5 },        // -1 .. 0
12            { 0.666667,  0.0, -1.0,  0.5 },        //  0 .. 1
13            { 1.33333,  -2.0,  1.0, -0.166667 },   //  1 .. 2
14    };
15    
16    //__constant float4 eyeVec = (float4)(25,15,10,1);
17    //__constant float4 origVec = (float4)(8.83877,2.5911,7.65275,0);
18    //__constant float4 cVec= (float4)(-0.0151831,0.0278357,0,0);
19    //__constant float4 rVec = (float4)(0.0074887,0.00408474,-0.0305383,0);
20    __constant float4 eyeVec = (float4)(127.331, -1322.05, 272.53, 0);
21    __constant float4 origVec = (float4)(122.835,17.7112,188.044, 0);
22    __constant float4 cVec= (float4)(-0.00403611,-0.029826,-0.244066, 0);
23    __constant float4 rVec = (float4)(-0.245595,-0.0112916,0.00544129, 0);
24    __constant float stepSize = 0.5;
25    __constant int s = 2;
26    
27    
28    __kernel float probe(float * img, int * sAxis, float4 imgPos)
29    {
30    
31        float probedVal;
32        float4 f, nf, t, hx, hy, hz;
33        int4 n;
34    
35        float4 d = (float4) (h[3][0],h[2][0],h[1][0],h[0][0]);
36        float4 c = (float4) (h[3][1],h[2][1],h[1][1],h[0][1]);
37        float4 b = (float4) (h[3][2],h[2][2],h[1][2],h[0][2]);
38        float4 a = (float4) (h[3][3],h[2][3],h[1][3],h[0][3]);
39    
40        f = modf(imgPos,&nf);
41        n = convert_int4(nf);
42    
43    
44        t = (float4) (f.x + 1, f.x, f.x -1, f.x - 2);
45        hx = d + t * (c + t * (b + t*a));
46        t = (float4) (f.y + 1, f.y, f.y - 1, f.y - 2);
47        hy = d + t * (c + t * (b + t*a));
48        t = (float4) (f.z + 1, f.z, f.z - 1, f.z -2);
49        hz = d + t * (c + t * (b + t*a));
50    
51        float vx[4];
52        float vy[4];
53    
54        for(int k = 1-s; k <= s; k++) {             // z is the slowest dimension
55            for(int j = 1-s;  j <= s; j++) {        // y is the medium dimension
56                int index = sAxis[0]*sAxis[1]*(n.z+k) + sAxis[0]*(n.y+j) + (n.x-1);
57                float4 v = (float4)(img[index], img[index+1], img[index+2], img[index+3]);
58                vx[j+s-1] = dot(v,hx);
59            }
60            vy[k+s-1] = dot ((float4) (vx[0],vx[1],vx[2],vx[3]),hy);
61        }
62    
63        probedVal = dot((float4) (vy[0],vy[1],vy[2],vy[3]),hz);
64    
65        return probedVal;
66    
67    }
68    
69    __kernel void mip ( float * img,
70                                              float * out,                                              float * out,
                                             float4  orig,  
                                             float4  eyeVector,  
                                             float4  cVec,  
                                             float4  rVec,  
71                                              float16  transformMatrix,                                              float16  transformMatrix,
72                                              float stepSize,                      int2 workDim,
73                                              int length,                      int * sAxis)
                                             int width,  
                                             int height)  
74  {  {
75    
76          int row = get_global_id(0), col = get_global_id(1);          int row = get_global_id(0), col = get_global_id(1);
77    
78          if(row <= 200 &&  col <= 200)      if(row <= workDim.x &&  col <= workDim.y) {
79          {          float t, probedVal, maxValue = -INFINITY;
80                  float t,x,y,z,probedVal, maxValue = -INFINITY;          float4 imgPt;
                 float4 t_i, t_j, t_k, value, f, imgPt, pt, v, nf;  
                 int4 n1;  
                 int4 n;  
81    
82                  float4 d = (float4) (h2[0],h1[0],h1[0],h2[0]);          float4 pos = origVec + (float)row * rVec + (float)col * cVec;
83                  float4 c = (float4) (h2[1],h1[1],h1[1],h2[1]);          float4 dir = normalize(pos - eyeVec);
                 float4 b = (float4) (h2[2],h1[2],h1[2],h2[2]);  
                 float4 a = (float4) (h2[3],h1[3],h1[3],h2[3]);  
   
                 float4 pos = orig + (float)row * rVec + (float)col * cVec;  
                 float4 dir = normalize(pos - eyeVector);  
                 float4 NSize = (float4) ((float)height,(float)length,(float)width,1.0f);  
84    
85                  for(t = 0.0;  t < 40; t+= stepSize)          pos.w = 0.0;
86                  {          dir.w = 0.0;
87    
88            for(t = 0.0;  t < 200; t+= stepSize) {
89                          pos = pos + stepSize * dir;                          pos = pos + stepSize * dir;
90    
                         // Begin Probe Operation  
   
91                          // Transform the value to image space position.                          // Transform the value to image space position.
92                          imgPt = (float4) (dot(pos,transformMatrix.s0123),                          imgPt = (float4) (dot(pos,transformMatrix.s0123),
93                                            dot(pos,transformMatrix.s4567),                                            dot(pos,transformMatrix.s4567),
94                                            dot(pos,transformMatrix.s89ab),                                            dot(pos,transformMatrix.s89ab),
95                                            dot(pos,transformMatrix.scdef));                                            dot(pos,transformMatrix.scdef));
96    
97                          nf = floor(imgPt);              if ((imgPt.x > 1 && imgPt.x < (sAxis[0] - 2)
98                          f = imgPt - nf;              &&  (imgPt.y > 1 && imgPt.y < (sAxis[1] - 2))
99                          n1 = convert_int4(nf);              &&  (imgPt.z > 1 && imgPt.z < (sAxis[2] - 2))) {
100                    probedVal = probe(img,sAxis,imgPt);
                         if( (imgPt.x > 1 && imgPt.x < NSize.x - 2) &&  
                                 (imgPt.y > 1 && imgPt.y < NSize.y - 2) &&  
                                 (imgPt.z > 1 && imgPt.z < NSize.z - 2))  
                         {  
                                 // the t value for h(fx - i)  
                                 float4 t_i = (float4) (-1.0 - f.x, -f.x, f.x - 1.0, f.x - 2.0);  
                                 // the t value for h(fy - j)  
                                 float4 t_j = (float4) (-1.0 - f.y, -f.y, f.y - 1.0, f.y - 2.0);  
                                 float4 t_k = (float4) (-1.0 - f.z, -f.z, f.z - 1.0, f.z - 2.0);  
                                 float tx[4], ty[4];  
                                 for (int i = -1;  i <= 2;  i++) {  
                                         for (int j = -1;  j <= 2;  j++) {  
                                             // compute z axis using vectors  
                                                 v = (float4)(  
                                                         img[(n1.x+i) * length * width + (n1.y+j) * width + n1.z-1],  
                                                         img[(n1.x+i) * length * width + (n1.y+j) * width + n1.z],  
                                                         img[(n1.x+i) * length * width + (n1.y+j) * width + n1.z+1],  
                                                         img[(n1.x+i) * length * width + (n1.y+j) * width + n1.z+2]);  
                                                 ty[j+1] = dot(v, d + t_k * (c + t_k * (b + t_k * a)));  
                                         }  
                                         tx[i+1] = dot((float4)(ty[0], ty[1], ty[2], ty[3]), d + t_j * (c + t_j * (b + t_j * a)));  
                                 }  
                                 probedVal = dot((float4)(tx[0],tx[1],tx[2],tx[3]), d + t_i * (c + t_i * (b + t_i * a)));  
   
 if ((row == 100) && (col == 100))  
         printf("%f %f %f\n", t, f.z, probedVal);  
                           // End Probe Operation  
101                                  if(maxValue < probedVal)                                  if(maxValue < probedVal)
102                                          maxValue = probedVal;                                          maxValue = probedVal;
103                          }                          }
104                  }                  }
105                  if(row == 0 && col == 74)          out[row * workDim.x + col] = maxValue;
                         printf("Max Value: %f\n", maxValue );  
                 out[row * 201 + col] = maxValue;  
106          }          }
107  }  }

Legend:
Removed from v.218  
changed lines
  Added in v.441

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