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

SCM Repository

[diderot] Diff of /trunk/test/MIP/mip.cl
ViewVC logotype

Diff of /trunk/test/MIP/mip.cl

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

revision 220, Tue Aug 3 20:08:09 2010 UTC revision 302, Tue Aug 17 03:52:21 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.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 float stepSize = 0.1;
21    __constant int s = 2;
22    
23    
24    __kernel float probe(float * img, int * sAxis, float4 imgPos) {
25    
26           float probedVal;
27           float4 f, nf, t, hx, hy, hz;
28           int4 n;
29    
30            float4 d = (float4) (h[3][0],h[2][0],h[1][0],h[0][0]);
31            float4 c = (float4) (h[3][1],h[2][1],h[1][1],h[0][1]);
32            float4 b = (float4) (h[3][2],h[2][2],h[1][2],h[0][2]);
33            float4 a = (float4) (h[3][3],h[2][3],h[1][3],h[0][3]);
34    
35            f = modf(imgPos,&nf);
36            n = convert_int4(nf);
37    
38    
39            t = (float4) (f.x + 1, f.x, f.x -1, f.x - 2);
40            hx = d + t * (c + t * (b + t*a));
41            t = (float4) (f.y + 1, f.y, f.y - 1, f.y - 2);
42            hy = d + t * (c + t * (b + t*a));
43            t = (float4) (f.z + 1, f.z, f.z - 1, f.z -2);
44            hz = d + t * (c + t * (b + t*a));
45    
46            float vy[4] = {0.0f,0.0f,0.0f,0.0f};
47            float vz[4] = {0.0f,0.0f,0.0f,0.0f};
48    
49            for(int i = 1-s; i <= s; i++) {
50             for(int j = 1-s;  j <= s; j++) {
51                float4 v = (float4) (
52                         img[((int)n.x+i) + sAxis[0]*(((int)n.y+j) + sAxis[1]*((int)n.z-1))],
53                         img[((int)n.x+i) + sAxis[0]*(((int)n.y+j) + sAxis[1]*((int)n.z))],
54                         img[((int)n.x+i) + sAxis[0]*(((int)n.y+j) + sAxis[1]*((int)n.z+1))],
55                         img[((int)n.x+i) + sAxis[0]*(((int)n.y+j) + sAxis[1]*((int)n.z+2))]);
56                vz[j+s-1] += dot(v,hz);
57             }
58                    vy[i+s-1] = dot ((float4) (vz[0],vz[1],vz[2],vz[3]),hy);
59            }
60    
61            probedVal = dot((float4) (vy[0],vy[1],vy[2],vy[3]),hx);
62    
63            return probedVal;
64    
65    }
66    __kernel void mip ( float * img,
67                                              float * out,                                              float * out,
                                             float4  orig,  
                                             float4  eyeVector,  
                                             float4  cVec,  
                                             float4  rVec,  
68                                              float16  transformMatrix,                                              float16  transformMatrix,
69                                              float stepSize,                                          int2 workDim,
70                                              int length,                                          int * sAxis){
71                                              int width,  
                                             int height)  
 {  
72          int row = get_global_id(0), col = get_global_id(1);          int row = get_global_id(0), col = get_global_id(1);
73    
74          if(row <= 200 &&  col <= 200)          if(row <= workDim.x &&  col <= workDim.y)
75          {          {
76                  float t,x,y,z,probedVal, maxValue = -INFINITY;                  float t, probedVal, maxValue = -INFINITY;
77                  float4 t_i, t_j, t_k, value, f, imgPt, pt, v, nf;                  float4 imgPt;
78                  int4 n1;  
79                  int4 n;                  float4 pos = origVec + (float)row * rVec + (float)col * cVec;
80                    float4 dir = normalize(pos - eyeVec);
81    
82                  float4 d = (float4) (h2[0],h1[0],h1[0],h2[0]);                  pos.w = 0.0;
83                  float4 c = (float4) (h2[1],h1[1],h1[1],h2[1]);                  dir.w = 0.0;
                 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)                  for(t = 0.0;  t < 40; t+= stepSize)
86                  {                  {
   
87                          pos = pos + stepSize * dir;                          pos = pos + stepSize * dir;
88    
                         // Begin Probe Operation  
   
89                          // Transform the value to image space position.                          // Transform the value to image space position.
90                          imgPt = (float4) (dot(pos,transformMatrix.s0123),                          imgPt = (float4) (dot(pos,transformMatrix.s0123),
91                                            dot(pos,transformMatrix.s4567),                                            dot(pos,transformMatrix.s4567),
92                                            dot(pos,transformMatrix.s89ab),                                            dot(pos,transformMatrix.s89ab),
93                                            dot(pos,transformMatrix.scdef));                                            dot(pos,transformMatrix.scdef));
94    
95                          nf = floor(imgPt);  
96                          f = imgPt - nf;                          if( (imgPt.x > 1 && imgPt.x < (sAxis[0] - 2)) &&
97                          n1 = convert_int4(nf);                                  (imgPt.y > 1 && imgPt.y < (sAxis[1]- 2)) &&
98                                    (imgPt.z > 1 && imgPt.z < (sAxis[2] - 2)))
                         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))  
99                          {                          {
100                                  // the t value for h(fx - i)                                  probedVal = probe(img,sAxis,imgPt);
                                 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  
                                                 int index = (n1.x+i) * length * width + (n1.y+j) * width + n1.z-1;  
                                                 v = (float4)(  
                                                         img[index],  
                                                         img[index+1],  
                                                         img[index+2],  
                                                         img[index+3]);  
                                                 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.220  
changed lines
  Added in v.302

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