65 |
|
|
66 |
STATIC_INLINE bool inside(double pos[3]) |
STATIC_INLINE bool inside(double pos[3]) |
67 |
{ |
{ |
68 |
// XXX - Hack |
// XXX - Hack specific to ../../data/vfrhand-nohip.nhdr" |
69 |
if(pos[0] < -0.5 || pos[0] > 175.5) return false; |
// but for some reason still not matching what Diderot says |
70 |
if(pos[1] < -0.5 || pos[1] > 186.5) return false; |
if(pos[0] < 2*0.9375 || pos[0] > (176-3)*0.9375) return false; |
71 |
if(pos[2] < -0.5 || pos[2] > 189.5) return false; |
if(pos[1] < 2*0.9375 || pos[1] > (187-3)*0.9375) return false; |
72 |
|
if(pos[2] < 2.0 || pos[2] > (190-3)) return false; |
73 |
return true; |
return true; |
74 |
} |
} |
75 |
|
|
87 |
{ |
{ |
88 |
const char *me = argv[0]; |
const char *me = argv[0]; |
89 |
char *infile = "../../data/vfrhand-nohip.nhdr"; |
char *infile = "../../data/vfrhand-nohip.nhdr"; |
90 |
char *outfile = "vr-lite-cam.pgm"; |
char *outfile = "bmark-teem.nrrd"; |
91 |
double camEye[3] = {127.331, -1322.05, 272.53}; |
double camEye[3] = {127.331, -1322.05, 272.53}; |
92 |
double camAt[3] = {63.0, 82.6536, 98.0}; |
double camAt[3] = {63.0, 82.6536, 98.0}; |
93 |
double camUp[3] = {0.9987, 0.0459166, -0.0221267}; |
double camUp[3] = {0.9987, 0.0459166, -0.0221267}; |
96 |
double camFOV = 5.0; |
double camFOV = 5.0; |
97 |
int imgResU = 480; |
int imgResU = 480; |
98 |
int imgResV = 345; |
int imgResV = 345; |
99 |
double rayStep = 0.5; |
double rayStep = 0.1; |
100 |
double valOpacMin = 400.0; // 400.0 for skin, 1150.0 for bone |
double valOpacMin = 400.0; // 400.0 for skin, 1150.0 for bone |
101 |
double valOpacMax = 700.0; // 700.0 for skin, 1450.0 for bone |
double valOpacMax = 700.0; // 700.0 for skin, 1450.0 for bone |
102 |
|
|
103 |
|
|
104 |
double temp[3]; |
double temp[3]; |
105 |
double temp2[3]; |
double temp2[3]; |
106 |
double temp3[3]; |
double temp3[3]; |
140 |
nin = nrrdNew(); |
nin = nrrdNew(); |
141 |
airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways); |
airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways); |
142 |
|
|
143 |
|
// Get data for mip |
144 |
if (nrrdLoad(nin, infile, NULL)) { |
if (nrrdLoad(nin, infile, NULL)) { |
145 |
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
146 |
fprintf(stderr, "Trouble reading \"%s\":\n%s", infile, err); |
fprintf(stderr, "%s: Trouble reading \"%s\":\n%s", me, infile, err); |
147 |
airMopError(mop); |
airMopError(mop); |
148 |
return -1; |
return -1; |
149 |
} |
} |
150 |
|
|
|
// Get data for mip |
|
151 |
gageContext *ctx; |
gageContext *ctx; |
152 |
gagePerVolume *pvl; |
gagePerVolume *pvl; |
153 |
const double *val; |
const double *val; |
154 |
const double *grad; |
const double *grad; |
155 |
double kparm[NRRD_KERNEL_PARMS_NUM] = {0.0}; |
double kparm[NRRD_KERNEL_PARMS_NUM] = {0.0}; |
|
int E; |
|
156 |
|
|
157 |
ctx = gageContextNew(); |
ctx = gageContextNew(); |
158 |
airMopAdd(mop, ctx, (airMopper)gageContextNix, airMopAlways); |
airMopAdd(mop, ctx, (airMopper)gageContextNix, airMopAlways); |
159 |
pvl = gagePerVolumeNew(ctx, nin, gageKindScl); |
if (!( pvl = gagePerVolumeNew(ctx, nin, gageKindScl) )) { |
160 |
if (!( ctx && pvl )) { |
airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways); |
161 |
fprintf(stderr, "%s: couldn't allocate gage ctx or pvl", me); |
fprintf(stderr, "%s: couldn't create gage context:\n%s", me, err); |
162 |
airMopError(mop); |
airMopError(mop); |
163 |
return -1; |
return -1; |
164 |
} |
} |
169 |
|| gageQueryItemOn(ctx, pvl, gageSclGradVec) |
|| gageQueryItemOn(ctx, pvl, gageSclGradVec) |
170 |
|| gageUpdate(ctx)) { |
|| gageUpdate(ctx)) { |
171 |
airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways); |
airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways); |
172 |
fprintf(stderr, "Trouble setting up gage context:\n%s", err); |
fprintf(stderr, "%s: Trouble setting up gage context:\n%s", me, err); |
173 |
airMopError(mop); |
airMopError(mop); |
174 |
return -1; |
return -1; |
175 |
} |
} |
176 |
|
|
177 |
if (!( (val = gageAnswerPointer(ctx, pvl, gageSclValue)) |
if (!( (val = gageAnswerPointer(ctx, pvl, gageSclValue)) |
178 |
&& (grad = gageAnswerPointer(ctx, pvl, gageSclGradVec)) )) { |
&& (grad = gageAnswerPointer(ctx, pvl, gageSclGradVec)) )) { |
179 |
fprintf(stderr, "Trouble getting answer pointer\n"); |
fprintf(stderr, "%s: Trouble getting answer pointer\n", me); |
180 |
airMopError(mop); |
airMopError(mop); |
181 |
return -1; |
return -1; |
182 |
} |
} |
184 |
int ui, vi; |
int ui, vi; |
185 |
double rayU, rayV, rayN; |
double rayU, rayV, rayN; |
186 |
double rayVec[3]; |
double rayVec[3]; |
187 |
double rayDir[3]; |
double toEye[3]; |
188 |
double rayPos[3]; |
double rayPos[3]; |
189 |
double transp, gray; |
double transp, gray; |
|
double output[4]; |
|
190 |
double norm[3]; |
double norm[3]; |
191 |
double alpha; |
double alpha; |
192 |
double ld; |
double ld; |
193 |
double hd; |
double hd; |
194 |
double mat; |
double mat; |
|
bool set = false; |
|
|
double max_gray = 0; |
|
|
double min_gray = 0; |
|
195 |
|
|
196 |
// allocate image and output |
// allocate output image |
197 |
Nrrd *nimg = nrrdNew(); |
Nrrd *nimg = nrrdNew(); |
198 |
airMopAdd(mop, nimg, (airMopper)nrrdNuke, airMopAlways); |
airMopAdd(mop, nimg, (airMopper)nrrdNuke, airMopAlways); |
|
Nrrd *nout = nrrdNew(); |
|
|
airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways); |
|
199 |
if (nrrdAlloc_va(nimg, nrrdTypeDouble, 3, |
if (nrrdAlloc_va(nimg, nrrdTypeDouble, 3, |
200 |
AIR_CAST(size_t, 4), |
AIR_CAST(size_t, 4), |
201 |
AIR_CAST(size_t, imgResU), |
AIR_CAST(size_t, imgResU), |
|
AIR_CAST(size_t, imgResV)) |
|
|
|| nrrdAlloc_va(nout, nrrdTypeUChar, 2, |
|
|
AIR_CAST(size_t, imgResU), |
|
202 |
AIR_CAST(size_t, imgResV))) { |
AIR_CAST(size_t, imgResV))) { |
203 |
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
204 |
fprintf(stderr, "Trouble allocating image and output:\n%s", err); |
fprintf(stderr, "%s: Trouble allocating image output:\n%s", me, err); |
205 |
airMopError(mop); |
airMopError(mop); |
206 |
return -1; |
return -1; |
207 |
} |
} |
208 |
double *out_data = AIR_CAST(double *, nimg->data); |
double *out_data = AIR_CAST(double *, nimg->data); |
209 |
unsigned char *uc_out_data = AIR_CAST(unsigned char *, nout->data); |
|
210 |
|
fprintf(stderr, "hit return\n"); |
211 |
|
fgetc(stdin); |
212 |
|
fprintf(stderr, "here we go\n"); |
213 |
|
|
214 |
double t0 = airTime(); // start timing |
double t0 = airTime(); // start timing |
215 |
for(vi = 0; vi < imgResV; vi++) { |
for(vi = 0; vi < imgResV; vi++) { |
216 |
rayV = lerp(-camVmax, camVmax, -0.5, (double)(vi), (double)(imgResV) - 0.5); |
rayV = lerp(-camVmax, camVmax, -0.5, (double)(vi), (double)(imgResV) - 0.5); |
|
scale_vec3(rayV, camV, temp); |
|
217 |
for(ui = 0; ui < imgResU; ui++) { |
for(ui = 0; ui < imgResU; ui++) { |
218 |
|
double len; |
219 |
rayU = lerp(-camUmax, camUmax, -0.5, (double)(ui), (double)(imgResU) - 0.5); |
rayU = lerp(-camUmax, camUmax, -0.5, (double)(ui), (double)(imgResU) - 0.5); |
220 |
scale_vec3(rayU, camU, temp2); |
ELL_3V_SCALE_ADD3(rayVec, 1.0, camN, rayU/camDist, camU, rayV/camDist, camV); |
221 |
add3(temp, temp2, temp2); |
ELL_3V_SCALE(toEye, -1, rayVec); |
222 |
scale_vec3(camDist, camN, temp3); |
ELL_3V_NORM(toEye, toEye, len); |
|
add3(temp2, temp3, temp3); |
|
|
scale_vec3(1 / camDist, temp3, rayVec); |
|
|
norm3(rayVec, rayDir); |
|
223 |
transp = 1; |
transp = 1; |
224 |
gray = 0; |
gray = 0; |
225 |
output[0] = 0; |
|
|
output[1] = 0; |
|
|
output[2] = 0; |
|
|
output[3] = 0; |
|
226 |
for (rayN = camVspNear; rayN <= camVspFar; rayN += rayStep) { |
for (rayN = camVspNear; rayN <= camVspFar; rayN += rayStep) { |
227 |
scale_vec3(rayN, rayDir, temp2); |
ELL_3V_SCALE_ADD2(rayPos, 1.0, camEye, rayN, rayVec); |
228 |
add3(temp2, camEye, rayPos); |
|
229 |
if (inside(rayPos)) { |
if (inside(rayPos)) { |
230 |
if (gageProbe(ctx, rayPos[0], rayPos[1], rayPos[2])) { |
if (gageProbeSpace(ctx, rayPos[0], rayPos[1], rayPos[2], |
231 |
fprintf(stderr, "Trouble with gageProbe:\n%s", ctx->errStr); |
AIR_FALSE /* index space */, AIR_TRUE /* clamp */)) { |
232 |
|
fprintf(stderr, "%s: Trouble with gageProbe:\n%s", me, ctx->errStr); |
233 |
airMopError(mop); return -1; |
airMopError(mop); return -1; |
234 |
} |
} |
235 |
if (*val > valOpacMin) { |
if (*val > valOpacMin) { |
236 |
temp2[0] = -*(grad + 0); |
ELL_3V_SCALE(temp2, -1, grad); |
237 |
temp2[1] = -*(grad + 1); |
ELL_3V_NORM(norm, temp2, len); |
|
temp2[2] = -*(grad + 2); |
|
|
norm3(temp2, norm); |
|
238 |
alpha = lerp(0.0, 1.0, valOpacMin, *val, valOpacMax); |
alpha = lerp(0.0, 1.0, valOpacMin, *val, valOpacMax); |
239 |
alpha = (alpha < 1.0) ? alpha : 1; |
alpha = (alpha < 1.0) ? alpha : 1; |
240 |
ld = dot3(norm, lightDir); |
ld = dot3(norm, lightDir); |
241 |
ld = (ld < 0) ? 0 : ld; |
ld = (ld < 0) ? 0 : ld; |
242 |
sub3(camEye, rayPos, temp2); |
ELL_3V_ADD2(temp2, toEye, lightDir); |
243 |
norm3(temp2, temp2); |
ELL_3V_NORM(temp2, temp2, len); |
|
add3(temp2, lightDir, temp2); |
|
|
norm3(temp2, temp2); |
|
244 |
hd = dot3(norm, temp2); |
hd = dot3(norm, temp2); |
245 |
hd = (hd < 0) ? 0 : hd; |
hd = (hd < 0) ? 0 : hd; |
246 |
mat = phongKa + |
mat = (phongKa + |
247 |
phongKd * ((ld > 0) ? ld : 0) + |
phongKd * ld + |
248 |
phongKs * ((hd > 0) ? pow(hd, phongSp) : 0); |
phongKs * pow(hd, phongSp)); |
249 |
gray = gray + transp * alpha * mat; |
gray = gray + transp * alpha * mat; |
250 |
transp *= 1 - alpha; |
transp *= 1 - alpha; |
251 |
} |
} |
255 |
break; |
break; |
256 |
} |
} |
257 |
} |
} |
|
*(out_data + ui * imgResV * 4 + vi * 4 + 0) = gray; |
|
|
*(out_data + ui * imgResV * 4 + vi * 4 + 1) = gray; |
|
|
*(out_data + ui * imgResV * 4 + vi * 4 + 2) = gray; |
|
|
*(out_data + ui * imgResV * 4 + vi * 4 + 3) = 1.0f - transp; |
|
|
if(!set) { |
|
|
max_gray = gray; |
|
|
min_gray = gray; |
|
|
set = true; |
|
|
} |
|
|
max_gray = (gray > max_gray) ? gray : max_gray; |
|
|
min_gray = (gray < min_gray) ? gray : min_gray; |
|
|
} |
|
|
} |
|
258 |
|
|
259 |
for(vi = 0; vi < imgResV; vi++) { |
ELL_4V_SET(out_data + 4*(ui + imgResU*vi), gray, gray, gray, 1.0-transp); |
260 |
for(ui = 0; ui < imgResU; ui++) { |
|
|
*(uc_out_data + vi * imgResU + ui) = ulerp(min_gray, *(out_data + ui * imgResV * 4 + vi * 4 + 0), max_gray); |
|
261 |
} |
} |
262 |
} |
} |
263 |
|
|
264 |
double totalTime = airTime() - t0; |
double totalTime = airTime() - t0; |
265 |
printf("usr=%f\n", totalTime); |
printf("usr=%f\n", totalTime); |
266 |
|
|
267 |
if (nrrdSave(outfile, nout, NULL)) { |
if (nrrdSave(outfile, nimg, NULL)) { |
268 |
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
269 |
fprintf(stderr, "Trouble writing output:\n%s", err); |
fprintf(stderr, "%s: Trouble writing output:\n%s", me, err); |
270 |
airMopError(mop); |
airMopError(mop); |
271 |
return -1; |
return -1; |
272 |
} |
} |