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 191, Mon Aug 2 14:05:11 2010 UTC revision 310, Tue Aug 17 18:52:00 2010 UTC
# Line 1  Line 1 
1  __kernel void raycast (__global float * img,  /* mip.cl
2                                             __global float * h1,   *
3                                             __global float * h2,   * COPYRIGHT (c) 2010 The Diderot Project (http://diderot.cs.uchicago.edu)
4                                             __global float * out,   * All rights reserved.
5                                              float4  orig,   *
6                                              float4  eyeVector,   * An OpenCL file containing a probe kernel and mip kernel implementation
7                                              float4  cVec,   */
8                                              float4  rVec,  
9                                              float16  transformMatrix,  __constant float h[4][4] = {            // bspln3
10                                              float stepSize,          { 1.33333,   2.0,  1.0,  0.166667 },   // -2 .. -1
11                                              int length,          { 0.666667,  0.0, -1.0, -0.5 },        // -1 .. 0
12                                              int width)          { 0.666667,  0.0, -1.0,  0.5 },        //  0 .. 1
13  {          { 1.33333,  -2.0,  1.0, -0.166667 },   //  1 .. 2
14          int row = get_global_id(0), col = get_global_id(1);  };
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          if(row < 200 &&  col < 200)  
28          {  __kernel float probe(float * img, int * sAxis, float4 imgPos) {
29                  int i;  
30                  float t,x,y,z,probedVal, maxValue = -INFINITY;         float probedVal;
31                  float4 t_i, t_j, t_k, value, f, imgPt, pt, v;         float4 f, nf, t, hx, hy, hz;
                 int4 n1;  
32                  int4 n;                  int4 n;
33    
34            float4 d = (float4) (h[3][0],h[2][0],h[1][0],h[0][0]);
35            float4 c = (float4) (h[3][1],h[2][1],h[1][1],h[0][1]);
36            float4 b = (float4) (h[3][2],h[2][2],h[1][2],h[0][2]);
37            float4 a = (float4) (h[3][3],h[2][3],h[1][3],h[0][3]);
38    
39            f = modf(imgPos,&nf);
40            n = convert_int4(nf);
41    
42    
43            t = (float4) (f.x + 1, f.x, f.x -1, f.x - 2);
44            hx = d + t * (c + t * (b + t*a));
45            t = (float4) (f.y + 1, f.y, f.y - 1, f.y - 2);
46            hy = d + t * (c + t * (b + t*a));
47            t = (float4) (f.z + 1, f.z, f.z - 1, f.z -2);
48            hz = d + t * (c + t * (b + t*a));
49    
50            float vx[4];
51            float vy[4];
52    
53            for(int k = 1-s; k <= s; k++) {         // z is the slowest dimension
54              for(int j = 1-s;  j <= s; j++) {      // y is the medium dimension
55                 int index = sAxis[0]*sAxis[1]*(n.z+k) + sAxis[0]*(n.y+j) + (n.x-1);
56                 float4 v = (float4)(img[index], img[index+1], img[index+2], img[index+3]);
57                 vx[j+s-1] = dot(v,hx);
58              }
59              vy[k+s-1] = dot ((float4) (vx[0],vx[1],vx[2],vx[3]),hy);
60            }
61    
62            probedVal = dot((float4) (vy[0],vy[1],vy[2],vy[3]),hz);
63    
64                  float4 d = (float4) (h2[0],h1[0],h1[0],h2[0]);          return probedVal;
                 float4 c = (float4) (h2[1],h1[1],h1[1],h2[1]);  
                 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 =  (pos - eyeVector) / fabs(pos - eyeVector);  
   
                 pt.w = 1.0f;  
                 pos.w = 1.0f;  
                 dir.w = 1.0f;  
                 value.w = 1.0f;  
65    
66                  for(t = 0.0;  t < 20; t+= stepSize)  }
67    __kernel void mip ( float * img,
68                                            float * out,
69                                            float16  transformMatrix,
70                                            int2 workDim,
71                                            int * sAxis){
72    
73            int row = get_global_id(0), col = get_global_id(1);
74    
75            if(row <= workDim.x &&  col <= workDim.y)
76                  {                  {
77                    float t, probedVal, maxValue = -INFINITY;
78                    float4 imgPt;
79    
80                          pos = pos + stepSize * dir;                  float4 pos = origVec + (float)row * rVec + (float)col * cVec;
81                    float4 dir = normalize(pos - eyeVec);
82    
83                    pos.w = 0.0;
84                    dir.w = 0.0;
85    
86                          // Begin Probe Operation                  for(t = 0.0;  t < 200; t+= stepSize)
87                    {
88                            pos = pos + stepSize * dir;
89    
90                          // Transform the value to image space position.                          // Transform the value to image space position.
91                          imgPt = (float4) (dot(pos,transformMatrix.s0123),                          imgPt = (float4) (dot(pos,transformMatrix.s0123),
# Line 50  Line 94 
94                                                            dot(pos,transformMatrix.scdef));                                                            dot(pos,transformMatrix.scdef));
95    
96    
97                            if( (imgPt.x > 1 && imgPt.x < (sAxis[0] - 2)) &&
98                          f.xyzw = (float4) (modf(imgPt.x,&x),modf(imgPt.y,&y),modf(imgPt.z,&z),1.0f);                                  (imgPt.y > 1 && imgPt.y < (sAxis[1]- 2)) &&
99                          n1 = (int4) ((int)x,(int)y,(int)z,1.0f);                                  (imgPt.z > 1 && imgPt.z < (sAxis[2] - 2)))
   
                         if( (n1.x > 1 && n1.x < 11) &&  
                                 (n1.y > 1 && n1.y < 12) &&  
                                 (n1.z > 1 && n1.z < 14))  
100                          {                          {
101                                    probedVal = probe(img,sAxis,imgPt);
   
                                 // the t value for h(fx - i)  
                                 t_i = (float4) (-1.0-f.x, -f.x,f.x - 1.0, f.x - 2.0);  
   
                                 // the t value for h(fy - j)  
                                 t_j = (float4) (-1.0-f.y, -f.y, f.y - 1.0,f.y - 2.0);  
   
                                 // the t value for h(fx - k)  
                                 t_k = (float4) (-1.0-f.z, -f.z, f.z - 1.0,f.z - 2.0);  
   
   
                                 value = ( (d + t_i * (c + t_i * (b + t_i * a)))  *  // h(fx - i) *  
                                                   (d + t_j * (c + t_j * (b + t_j * a)))  *  // h(fy - j) *  
                                                   (d + t_k * (c + t_k * (b + t_k * a))));   // h(fz - k)  
   
   
                                 n = (int4) (((n1.x-1) * length * width + n1.y-1 * width + n1.z-1),  
                                             (n1.x * length * width +  n1.y * width + n1.z),  
                                             ((n1.x+1) * length * width +  n1.y+1 * width + n1.z+1),  
                                             ((n1.x+2) * length * width + n1.y+2 * width + n1.z+2));  
   
   
                                         v = (float4)(img[n.x],  
                                                                  img[n.y],  
                                                                  img[n.z],  
                                                                  img[n.w]);  
   
   
                            probedVal  = dot(v,value);  // V(n + <i,j,k>) * summations of the h(x) components  
   
                            // End Probe Operation  
   
102                            if(maxValue < probedVal)                            if(maxValue < probedVal)
103                                          maxValue = probedVal;                                          maxValue = probedVal;
104                          }                          }
                         if(row == 0 && col == 0)  
                                 printf("IMG PT (%f,%f,%f)\n", imgPt.x,imgPt.y,imgPt.z );  
105                  }                  }
106                  out[row * 200 + col] = maxValue;                  out[row * workDim.x + col] = maxValue;
107          }          }
108  }  }

Legend:
Removed from v.191  
changed lines
  Added in v.310

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