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

Parent Directory | Revision Log | Patch

	revision 301, Sun Aug 15 23:13:20 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	}	}

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

Legend: Removed from v.301   changed lines   Added in v.302

---

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