__kernel void raycast ( float * img,  
					    float * h1,
					    float * h2, 
					    float * out,
					    float4  orig, 
					    float4  eyeVector, 
					    float4  cVec, 
					    float4  rVec, 
					    float16  transformMatrix, 
					    float stepSize, 
					    int length,
					    int width,
					    int height)
{  
	int row = get_global_id(0), col = get_global_id(1);

	if(row <= 200 &&  col <= 200) 
	{ 	  	
		float t,x,y,z,probedVal, maxValue = -INFINITY; 
		float4 t_i, t_j, t_k, value, f, imgPt, pt, v, nf;
		int4 n1; 
		int4 n;
	
		float4 d = (float4) (h2[0],h1[0],h1[0],h2[0]); 
		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 = normalize(pos - eyeVector);
		float4 NSize = (float4) ((float)height,(float)length,(float)width,1.0f);
				
		for(t = 0.0;  t < 40; t+= stepSize) 
		{ 
		
			pos = pos + stepSize * dir; 

			// Begin Probe Operation	
			
			// 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));

			nf = floor(imgPt);
			f = imgPt - nf;
			n1 = convert_int4(nf);

			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 
				if(maxValue < probedVal) 
					maxValue = probedVal;
			}
		}		
		if(row == 0 && col == 74) 
			printf("Max Value: %f\n", maxValue ); 
		out[row * 201 + 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:  0 &amp;&amp; col == 74) 
			printf(&quot;Max Value: %f\n&quot;, maxValue ); 
		out[row * 201 + col] = maxValue; 
	}
}</div>
</div></div>
</body></html>