/* mip.cl
 *
 * COPYRIGHT (c) 2010 The Diderot Project (http://diderot.cs.uchicago.edu)
 * All rights reserved.
 *
 * An OpenCL file containing a probe kernel and mip kernel implementation 
 */
 
__constant float h[4][4] = {		// bspln3
	{ 1.33333,   2.0,  1.0,  0.166667 },   // -2 .. -1
	{ 0.666667,  0.0, -1.0, -0.5 },        // -1 .. 0
	{ 0.666667,  0.0, -1.0,  0.5 },        //  0 .. 1
	{ 1.33333,  -2.0,  1.0, -0.166667 },   //  1 .. 2
};

//__constant float4 eyeVec = (float4)(25,15,10,1); 
//__constant float4 origVec = (float4)(8.83877,2.5911,7.65275,0); 
//__constant float4 cVec= (float4)(-0.0151831,0.0278357,0,0);  
//__constant float4 rVec = (float4)(0.0074887,0.00408474,-0.0305383,0);
__constant float4 eyeVec = (float4)(127.331, -1322.05, 272.53, 0); 
__constant float4 origVec = (float4)(122.835,17.7112,188.044, 0); 
__constant float4 cVec= (float4)(-0.00403611,-0.029826,-0.244066, 0);  
__constant float4 rVec = (float4)(-0.245595,-0.0112916,0.00544129, 0);
__constant float stepSize = 0.5;
__constant int s = 2;


__kernel float probe(float * img, int * sAxis, float4 imgPos) { 

       float probedVal;
       float4 f, nf, t, hx, hy, hz;
       int4 n; 
          
        float4 d = (float4) (h[3][0],h[2][0],h[1][0],h[0][0]); 
        float4 c = (float4) (h[3][1],h[2][1],h[1][1],h[0][1]); 
        float4 b = (float4) (h[3][2],h[2][2],h[1][2],h[0][2]); 
        float4 a = (float4) (h[3][3],h[2][3],h[1][3],h[0][3]);     
        
        f = modf(imgPos,&nf); 
        n = convert_int4(nf);
       

        t = (float4) (f.x + 1, f.x, f.x -1, f.x - 2); 
        hx = d + t * (c + t * (b + t*a)); 
        t = (float4) (f.y + 1, f.y, f.y - 1, f.y - 2); 
        hy = d + t * (c + t * (b + t*a)); 
        t = (float4) (f.z + 1, f.z, f.z - 1, f.z -2); 
        hz = d + t * (c + t * (b + t*a));  

        float vx[4]; 
        float vy[4];
        
        for(int k = 1-s; k <= s; k++) {		// z is the slowest dimension
	  for(int j = 1-s;  j <= s; j++) {	// y is the medium dimension
	     int index = sAxis[0]*sAxis[1]*(n.z+k) + sAxis[0]*(n.y+j) + (n.x-1);
	     float4 v = (float4)(img[index], img[index+1], img[index+2], img[index+3]);
	     vx[j+s-1] = dot(v,hx); 
	  }
	  vy[k+s-1] = dot ((float4) (vx[0],vx[1],vx[2],vx[3]),hy); 
        } 
       
        probedVal = dot((float4) (vy[0],vy[1],vy[2],vy[3]),hz); 
                
        return probedVal; 
       
}
__kernel void mip ( float * img,  
					float * out,
					float16  transformMatrix, 
					int2 workDim, 
					int * sAxis){
					   
	int row = get_global_id(0), col = get_global_id(1);
	
	if(row <= workDim.x &&  col <= workDim.y) 
	{ 	  	
		float t, probedVal, maxValue = -INFINITY; 
		float4 imgPt;

		float4 pos = origVec + (float)row * rVec + (float)col * cVec;
		float4 dir = normalize(pos - eyeVec);
				
		pos.w = 0.0; 
		dir.w = 0.0; 

		for(t = 0.0;  t < 200; t+= stepSize) 
		{ 		
			pos = pos + stepSize * dir; 
			
			// Transform the value to image space position.  
			imgPt = (float4) (dot(pos,transformMatrix.s0123),
					  dot(pos,transformMatrix.s4567),
					  dot(pos,transformMatrix.s89ab),
					  dot(pos,transformMatrix.scdef));

				
			if( (imgPt.x > 1 && imgPt.x < (sAxis[0] - 2)) && 
				(imgPt.y > 1 && imgPt.y < (sAxis[1]- 2)) && 
				(imgPt.z > 1 && imgPt.z < (sAxis[2] - 2))) 
			{
				probedVal = probe(img,sAxis,imgPt);  	
				if(maxValue < probedVal) 
					maxValue = probedVal;
			}
		}		
		out[row * workDim.x + col] = maxValue; 
	}
}

<script type="text/javascript">//<![CDATA[
		function toggle_ffErrors() {
			var errorsblock = document.getElementById("ffErrorsBlock");
			var errorsgroup = document.getElementById("ffErrors");
			if (errorsblock.style.display == "none") {
				errorsblock.style.display = "block";
				errorsgroup.style.right = "10px";
			} else {
				errorsblock.style.display = "none";
				errorsgroup.style.right = "300px";
			}
		}
//]]></script>
<div id="ffErrors">
<a href="javascript:toggle_ffErrors();">Click to toggle</a>
<div id="ffErrorsBlock">
	<div class="error">does not end with &lt;/html&gt; tag</div>
	<div class="error">does not end with &lt;/body&gt; tag</div>
	<div class="info">The output has ended thus: e &lt; probedVal) 
					maxValue = probedVal;
			}
		}		
		out[row * workDim.x + col] = maxValue; 
	}
}</div>
</div></div>
</body></html>