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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 302 - (download) (annotate)
Tue Aug 17 03:52:21 2010 UTC (9 years, 9 months ago) by lamonts
File size: 3439 byte(s)
Updated verision of the mip code that appears to be working fine.
/* 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 float stepSize = 0.1;
__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 vy[4] = {0.0f,0.0f,0.0f,0.0f}; 
        float vz[4] = {0.0f,0.0f,0.0f,0.0f};
        
        for(int i = 1-s; i <= s; i++) { 
         for(int j = 1-s;  j <= s; j++) { 
            float4 v = (float4) (
                     img[((int)n.x+i) + sAxis[0]*(((int)n.y+j) + sAxis[1]*((int)n.z-1))], 
                     img[((int)n.x+i) + sAxis[0]*(((int)n.y+j) + sAxis[1]*((int)n.z))],
                     img[((int)n.x+i) + sAxis[0]*(((int)n.y+j) + sAxis[1]*((int)n.z+1))],
                     img[((int)n.x+i) + sAxis[0]*(((int)n.y+j) + sAxis[1]*((int)n.z+2))]);
            vz[j+s-1] += dot(v,hz); 
         }
         	vy[i+s-1] = dot ((float4) (vz[0],vz[1],vz[2],vz[3]),hy); 
        } 
       
        probedVal = dot((float4) (vy[0],vy[1],vy[2],vy[3]),hx); 
                
        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 < 40; 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; 
	}
} 

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