1 |
|
|
/* |
2 |
|
|
Teem: Tools to process and visualize scientific data and images . |
3 |
|
|
Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago |
4 |
|
|
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann |
5 |
|
|
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah |
6 |
|
|
|
7 |
|
|
This library is free software; you can redistribute it and/or |
8 |
|
|
modify it under the terms of the GNU Lesser General Public License |
9 |
|
|
(LGPL) as published by the Free Software Foundation; either |
10 |
|
|
version 2.1 of the License, or (at your option) any later version. |
11 |
|
|
The terms of redistributing and/or modifying this software also |
12 |
|
|
include exceptions to the LGPL that facilitate static linking. |
13 |
|
|
|
14 |
|
|
This library is distributed in the hope that it will be useful, |
15 |
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of |
16 |
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
17 |
|
|
Lesser General Public License for more details. |
18 |
|
|
|
19 |
|
|
You should have received a copy of the GNU Lesser General Public License |
20 |
|
|
along with this library; if not, write to Free Software Foundation, Inc., |
21 |
|
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
22 |
|
|
*/ |
23 |
|
|
|
24 |
|
|
#include <teem/air.h> |
25 |
|
|
#include <teem/hest.h> |
26 |
|
|
#include <teem/biff.h> |
27 |
|
|
#include <teem/nrrd.h> |
28 |
|
|
#include <teem/gage.h> |
29 |
|
|
#include <teem/limn.h> |
30 |
|
|
#include <teem/hoover.h> |
31 |
|
|
#include <teem/ten.h> |
32 |
|
|
#include <teem/meet.h> |
33 |
|
|
|
34 |
|
|
#define MREND "mrender" |
35 |
|
|
|
36 |
|
|
static const char *info = |
37 |
|
|
("A demonstration of hoover, gage, and nrrd measures. " |
38 |
|
|
"Uses hoover to cast rays through a volume (scalar, vector, or " |
39 |
|
|
"tensor), gage to " |
40 |
|
|
"measure one of various quantities along the rays, and a " |
41 |
|
|
"specified nrrd measure to reduce all the values along a ray " |
42 |
|
|
"down to one scalar, which is saved in the output (double) " |
43 |
|
|
"image."); |
44 |
|
|
|
45 |
|
|
/* -------------------------------------------------------------- */ |
46 |
|
|
|
47 |
|
|
/* Even though the gageContext is really thread-specific, and |
48 |
|
|
therefore doesn't really belong in mrendUser, the first context |
49 |
|
|
from which all others is copied is logically shared across threads, |
50 |
|
|
as are the input parameter it contains. There is a per-thread |
51 |
|
|
gageContext pointer in mrendThread */ |
52 |
|
|
|
53 |
|
|
typedef struct { |
54 |
|
|
Nrrd *nin; /* input volume to render */ |
55 |
|
|
gageKind *kind; /* the kind of volume it is */ |
56 |
|
|
int verbPixel[2]; /* which pixel to do verbose stuff on */ |
57 |
|
|
double rayStep, /* distance between sampling planes */ |
58 |
|
|
fromNaN; /* what to convert non-existent value to */ |
59 |
|
|
int whatq, /* what to measure along the ray */ |
60 |
|
|
measr; /* how to reduce the ray samples to a scalar */ |
61 |
|
|
/* we have a separate copy of the kernel specs so that hest can |
62 |
|
|
set these, and then we'll gageKernelSet() them in the context |
63 |
|
|
in order to do the proper error checking- hest can't do the |
64 |
|
|
error checking that we need. */ |
65 |
|
|
NrrdKernelSpec *ksp[GAGE_KERNEL_MAX+1]; |
66 |
|
|
gageShape *shape; /* used to document ItoW transform */ |
67 |
|
|
gageContext *gctx0; /* gage input and parent thread state */ |
68 |
|
|
hooverContext *hctx; /* hoover input and state */ |
69 |
|
|
char *outS; /* (managed by hest) output filename */ |
70 |
|
|
double imgOrig[3], /* output: set as extra info about camera */ |
71 |
|
|
imgU[3], imgV[3]; |
72 |
|
|
|
73 |
|
|
airArray *mrmop; |
74 |
|
|
} mrendUser; |
75 |
|
|
|
76 |
|
|
mrendUser * |
77 |
|
|
mrendUserNew(void) { |
78 |
|
|
mrendUser *uu; |
79 |
|
|
int i; |
80 |
|
|
|
81 |
|
|
uu = AIR_CALLOC(1, mrendUser); |
82 |
|
|
uu->nin = NULL; |
83 |
|
|
uu->kind = NULL; |
84 |
|
|
uu->rayStep = 0.0; |
85 |
|
|
uu->whatq = gageSclUnknown; |
86 |
|
|
uu->measr = nrrdMeasureUnknown; |
87 |
|
|
for (i=gageKernelUnknown+1; i<gageKernelLast; i++) { |
88 |
|
|
uu->ksp[i] = NULL; |
89 |
|
|
} |
90 |
|
|
uu->shape = gageShapeNew(); |
91 |
|
|
uu->gctx0 = gageContextNew(); |
92 |
|
|
uu->hctx = hooverContextNew(); |
93 |
|
|
uu->outS = NULL; |
94 |
|
|
uu->mrmop = airMopNew(); |
95 |
|
|
airMopAdd(uu->mrmop, uu->shape, (airMopper)gageShapeNix, airMopAlways); |
96 |
|
|
airMopAdd(uu->mrmop, uu->gctx0, (airMopper)gageContextNix, airMopAlways); |
97 |
|
|
airMopAdd(uu->mrmop, uu->hctx, (airMopper)hooverContextNix, airMopAlways); |
98 |
|
|
return uu; |
99 |
|
|
} |
100 |
|
|
|
101 |
|
|
mrendUser * |
102 |
|
|
mrendUserNix(mrendUser *uu) { |
103 |
|
|
|
104 |
|
|
if (uu) { |
105 |
|
|
airMopOkay(uu->mrmop); |
106 |
|
|
airFree(uu); |
107 |
|
|
} |
108 |
|
|
return NULL; |
109 |
|
|
} |
110 |
|
|
|
111 |
|
|
int |
112 |
|
|
mrendUserCheck(mrendUser *uu) { |
113 |
|
|
static const char me[]="mrendUserCheck"; |
114 |
|
|
|
115 |
|
|
if (3 + uu->kind->baseDim != uu->nin->dim) { |
116 |
|
|
biffAddf(MREND, "%s: input nrrd needs %d dimensions, not %d", |
117 |
|
|
me, + uu->kind->baseDim, uu->nin->dim); |
118 |
|
|
return 1; |
119 |
|
|
} |
120 |
|
|
if (!( uu->nin->axis[0].center == uu->nin->axis[1].center && |
121 |
|
|
uu->nin->axis[0].center == uu->nin->axis[2].center )) { |
122 |
|
|
biffAddf(MREND, "%s: axes 0,1,2 centerings (%s,%s,%s) not equal", me, |
123 |
|
|
airEnumStr(nrrdCenter, uu->nin->axis[0].center), |
124 |
|
|
airEnumStr(nrrdCenter, uu->nin->axis[1].center), |
125 |
|
|
airEnumStr(nrrdCenter, uu->nin->axis[2].center)); |
126 |
|
|
return 1; |
127 |
|
|
} |
128 |
|
|
if (1 != uu->kind->table[uu->whatq].answerLength) { |
129 |
|
|
biffAddf(MREND, "%s: quantity %s (in %s volumes) isn't a scalar; " |
130 |
|
|
"can't render it", |
131 |
|
|
me, airEnumStr(uu->kind->enm, uu->whatq), uu->kind->name); |
132 |
|
|
return 1; |
133 |
|
|
} |
134 |
|
|
|
135 |
|
|
return 0; |
136 |
|
|
} |
137 |
|
|
|
138 |
|
|
/* -------------------------------------------------------------- */ |
139 |
|
|
|
140 |
|
|
typedef struct mrendRender_t { |
141 |
|
|
double time0, time1; /* render start and end times */ |
142 |
|
|
Nrrd *nout; /* output image: always 2D array of doubles */ |
143 |
|
|
double *imgData; /* output image data */ |
144 |
|
|
int sx, sy, /* image dimensions */ |
145 |
|
|
totalSamples; /* total number of samples used for all rays */ |
146 |
|
|
struct mrendThread_t *tinfo[HOOVER_THREAD_MAX]; |
147 |
|
|
} mrendRender; |
148 |
|
|
|
149 |
|
|
typedef struct mrendThread_t { |
150 |
|
|
double *val, /* array of ray samples */ |
151 |
|
|
rayLen, /* length of ray segment between near and far */ |
152 |
|
|
rayStep; /* ray step needed FOR THIS RAY, to acheive sampling on |
153 |
|
|
planes (same scaling of uu->rayStep) */ |
154 |
|
|
int thrid, /* thread ID */ |
155 |
|
|
valLen, /* allocated size of val */ |
156 |
|
|
valNum, /* number of values set in val (index of next value) */ |
157 |
|
|
ui, vi, /* image coords */ |
158 |
|
|
numSamples, /* total number of samples this thread has done */ |
159 |
|
|
verbose; /* blah blah blah blah */ |
160 |
|
|
gageContext *gctx; /* thread-specific gage context (or copy of uu->gctx0 |
161 |
|
|
for the first thread) */ |
162 |
|
|
const double *answer; /* pointer to the SINGLE answer we care about */ |
163 |
|
|
} mrendThread; |
164 |
|
|
|
165 |
|
|
int |
166 |
|
|
mrendRenderBegin(mrendRender **rrP, mrendUser *uu) { |
167 |
|
|
static const char me[]="mrendRenderBegin"; |
168 |
|
|
gagePerVolume *pvl; |
169 |
|
|
int E; |
170 |
|
|
unsigned int thr; |
171 |
|
|
|
172 |
|
|
/* this assumes that mrendUserCheck(uu) has passed */ |
173 |
|
|
|
174 |
|
|
*rrP = AIR_CALLOC(1, mrendRender); |
175 |
|
|
airMopAdd(uu->mrmop, *rrP, airFree, airMopAlways); |
176 |
|
|
/* pvl managed via parent context */ |
177 |
|
|
|
178 |
|
|
(*rrP)->time0 = airTime(); |
179 |
|
|
|
180 |
|
|
E = 0; |
181 |
|
|
if (!E) E |= !(pvl = gagePerVolumeNew(uu->gctx0, uu->nin, uu->kind)); |
182 |
|
|
if (!E) E |= gagePerVolumeAttach(uu->gctx0, pvl); |
183 |
|
|
if (!E) E |= gageKernelSet(uu->gctx0, gageKernel00, |
184 |
|
|
uu->ksp[gageKernel00]->kernel, |
185 |
|
|
uu->ksp[gageKernel00]->parm); |
186 |
|
|
if (!E) E |= gageKernelSet(uu->gctx0, gageKernel11, |
187 |
|
|
uu->ksp[gageKernel11]->kernel, |
188 |
|
|
uu->ksp[gageKernel11]->parm); |
189 |
|
|
if (!E) E |= gageKernelSet(uu->gctx0, gageKernel22, |
190 |
|
|
uu->ksp[gageKernel22]->kernel, |
191 |
|
|
uu->ksp[gageKernel22]->parm); |
192 |
|
|
if (!E) E |= gageQueryItemOn(uu->gctx0, pvl, uu->whatq); |
193 |
|
|
if (!E) E |= gageUpdate(uu->gctx0); |
194 |
|
|
if (E) { |
195 |
|
|
biffMovef(MREND, GAGE, "%s: gage trouble", me); |
196 |
|
|
return 1; |
197 |
|
|
} |
198 |
|
|
fprintf(stderr, "%s: kernel support = %d^3 samples\n", me, |
199 |
|
|
2*uu->gctx0->radius); |
200 |
|
|
|
201 |
|
|
if (nrrdMaybeAlloc_va((*rrP)->nout=nrrdNew(), nrrdTypeDouble, 2, |
202 |
|
|
AIR_CAST(size_t, uu->hctx->imgSize[0]), |
203 |
|
|
AIR_CAST(size_t, uu->hctx->imgSize[1]))) { |
204 |
|
|
biffMovef(MREND, NRRD, "%s: nrrd trouble", me); |
205 |
|
|
return 1; |
206 |
|
|
} |
207 |
|
|
(*rrP)->nout->axis[0].min = uu->hctx->cam->uRange[0]; |
208 |
|
|
(*rrP)->nout->axis[0].max = uu->hctx->cam->uRange[1]; |
209 |
|
|
(*rrP)->nout->axis[1].min = uu->hctx->cam->vRange[0]; |
210 |
|
|
(*rrP)->nout->axis[1].max = uu->hctx->cam->vRange[1]; |
211 |
|
|
airMopAdd(uu->mrmop, (*rrP)->nout, (airMopper)nrrdNuke, airMopAlways); |
212 |
|
|
(*rrP)->imgData = AIR_CAST(double*, (*rrP)->nout->data); |
213 |
|
|
(*rrP)->sx = uu->hctx->imgSize[0]; |
214 |
|
|
(*rrP)->sy = uu->hctx->imgSize[1]; |
215 |
|
|
|
216 |
|
|
for (thr=0; thr<uu->hctx->numThreads; thr++) { |
217 |
|
|
(*rrP)->tinfo[thr] = AIR_CALLOC(1, mrendThread); |
218 |
|
|
airMopAdd(uu->mrmop, (*rrP)->tinfo[thr], airFree, airMopAlways); |
219 |
|
|
} |
220 |
|
|
|
221 |
|
|
return 0; |
222 |
|
|
} |
223 |
|
|
|
224 |
|
|
int |
225 |
|
|
mrendRenderEnd(mrendRender *rr, mrendUser *uu) { |
226 |
|
|
static const char me[]="mrendRenderEnd"; |
227 |
|
|
unsigned int thr; |
228 |
|
|
|
229 |
|
|
/* add up # samples from all threads */ |
230 |
|
|
rr->totalSamples = 0; |
231 |
|
|
for (thr=0; thr<uu->hctx->numThreads; thr++) { |
232 |
|
|
rr->totalSamples += rr->tinfo[thr]->numSamples; |
233 |
|
|
} |
234 |
|
|
|
235 |
|
|
rr->time1 = airTime(); |
236 |
|
|
fprintf(stderr, "\n"); |
237 |
|
|
fprintf(stderr, "%s: rendering time = %g secs\n", me, |
238 |
|
|
rr->time1 - rr->time0); |
239 |
|
|
fprintf(stderr, "%s: sampling rate = %g KHz\n", me, |
240 |
|
|
rr->totalSamples/(1000.0*(rr->time1 - rr->time0))); |
241 |
|
|
if (nrrdSave(uu->outS, rr->nout, NULL)) { |
242 |
|
|
biffMovef(MREND, NRRD, "%s: trouble saving image", me); |
243 |
|
|
return 1; |
244 |
|
|
} |
245 |
|
|
|
246 |
|
|
return 0; |
247 |
|
|
} |
248 |
|
|
|
249 |
|
|
/* -------------------------------------------------------------- */ |
250 |
|
|
|
251 |
|
|
int |
252 |
|
|
mrendThreadBegin(mrendThread **ttP, |
253 |
|
|
mrendRender *rr, mrendUser *uu, int whichThread) { |
254 |
|
|
|
255 |
|
|
/* allocating the mrendThreads should be part of the thread body, |
256 |
|
|
but as long as there isn't a mutex around registering them with |
257 |
|
|
the airMop in the mrendRender, then all that needs to be done |
258 |
|
|
as part of mrendRenderBegin (see above) */ |
259 |
|
|
(*ttP) = rr->tinfo[whichThread]; |
260 |
|
|
if (!whichThread) { |
261 |
|
|
/* this is the first thread- it just points to the parent gageContext */ |
262 |
|
|
(*ttP)->gctx = uu->gctx0; |
263 |
|
|
} else { |
264 |
|
|
/* we have to generate a new gageContext */ |
265 |
|
|
(*ttP)->gctx = gageContextCopy(uu->gctx0); |
266 |
|
|
} |
267 |
|
|
(*ttP)->answer = gageAnswerPointer((*ttP)->gctx, |
268 |
|
|
(*ttP)->gctx->pvl[0], uu->whatq); |
269 |
|
|
(*ttP)->val = NULL; |
270 |
|
|
(*ttP)->valLen = 0; |
271 |
|
|
(*ttP)->valNum = 0; |
272 |
|
|
(*ttP)->rayLen = 0; |
273 |
|
|
(*ttP)->thrid = whichThread; |
274 |
|
|
(*ttP)->numSamples = 0; |
275 |
|
|
(*ttP)->verbose = 0; |
276 |
|
|
return 0; |
277 |
|
|
} |
278 |
|
|
|
279 |
|
|
int |
280 |
|
|
mrendThreadEnd(mrendThread *tt, mrendRender *rr, mrendUser *uu) { |
281 |
|
|
|
282 |
|
|
AIR_UNUSED(rr); |
283 |
|
|
AIR_UNUSED(uu); |
284 |
|
|
if (tt->thrid) { |
285 |
|
|
tt->gctx = gageContextNix(tt->gctx); |
286 |
|
|
} |
287 |
|
|
tt->val = AIR_CAST(double*, airFree(tt->val)); |
288 |
|
|
|
289 |
|
|
return 0; |
290 |
|
|
} |
291 |
|
|
|
292 |
|
|
/* -------------------------------------------------------------- */ |
293 |
|
|
|
294 |
|
|
int |
295 |
|
|
mrendRayBegin(mrendThread *tt, mrendRender *rr, mrendUser *uu, |
296 |
|
|
int uIndex, |
297 |
|
|
int vIndex, |
298 |
|
|
double rayLen, |
299 |
|
|
double rayStartWorld[3], |
300 |
|
|
double rayStartIndex[3], |
301 |
|
|
double rayDirWorld[3], |
302 |
|
|
double rayDirIndex[3]) { |
303 |
|
|
static const char me[]="mrendRayBegin"; |
304 |
|
|
int newLen; |
305 |
|
|
|
306 |
|
|
AIR_UNUSED(rr); |
307 |
|
|
AIR_UNUSED(rayStartWorld); |
308 |
|
|
AIR_UNUSED(rayStartIndex); |
309 |
|
|
AIR_UNUSED(rayDirWorld); |
310 |
|
|
AIR_UNUSED(rayDirIndex); |
311 |
|
|
tt->ui = uIndex; |
312 |
|
|
tt->vi = vIndex; |
313 |
|
|
if (!( -1 == uu->verbPixel[0] && -1 == uu->verbPixel[1] )) { |
314 |
|
|
if (uIndex == uu->verbPixel[0] && vIndex == uu->verbPixel[1]) { |
315 |
|
|
fprintf(stderr, "\n%s: verbose for pixel (%d,%d)\n", me, |
316 |
|
|
uu->verbPixel[0], uu->verbPixel[1]); |
317 |
|
|
gageParmSet(tt->gctx, gageParmVerbose, 6); |
318 |
|
|
tt->verbose = 6; |
319 |
|
|
} else { |
320 |
|
|
gageParmSet(tt->gctx, gageParmVerbose, AIR_FALSE); |
321 |
|
|
tt->verbose = 0; |
322 |
|
|
} |
323 |
|
|
} |
324 |
|
|
tt->rayLen = rayLen; |
325 |
|
|
tt->rayStep = (uu->rayStep*tt->rayLen / |
326 |
|
|
(uu->hctx->cam->vspFaar - uu->hctx->cam->vspNeer)); |
327 |
|
|
newLen = AIR_ROUNDUP(rayLen/tt->rayStep) + 1; |
328 |
|
|
if (!tt->val || newLen > tt->valLen) { |
329 |
|
|
tt->val = AIR_CAST(double*, airFree(tt->val)); |
330 |
|
|
tt->valLen = newLen; |
331 |
|
|
tt->val = AIR_CALLOC(newLen, double); |
332 |
|
|
} |
333 |
|
|
tt->valNum = 0; |
334 |
|
|
if (!uIndex) { |
335 |
|
|
fprintf(stderr, "%d/%d ", vIndex, uu->hctx->imgSize[1]); |
336 |
|
|
fflush(stderr); |
337 |
|
|
} |
338 |
|
|
|
339 |
|
|
if (0 == uIndex && 0 == vIndex) { |
340 |
|
|
ELL_3V_COPY(uu->imgOrig, rayStartWorld); |
341 |
|
|
} else if (1 == uIndex && 0 == vIndex) { |
342 |
|
|
ELL_3V_COPY(uu->imgU, rayStartWorld); /* will fix later */ |
343 |
|
|
} else if (0 == uIndex && 1 == vIndex) { |
344 |
|
|
ELL_3V_COPY(uu->imgV, rayStartWorld); /* will fix later */ |
345 |
|
|
} |
346 |
|
|
|
347 |
|
|
fflush(stderr); |
348 |
|
|
return 0; |
349 |
|
|
} |
350 |
|
|
|
351 |
|
|
int |
352 |
|
|
mrendRayEnd(mrendThread *tt, mrendRender *rr, mrendUser *uu) { |
353 |
|
|
double answer; |
354 |
|
|
|
355 |
|
|
if (tt->valNum) { |
356 |
|
|
nrrdMeasureLine[uu->measr](&answer, |
357 |
|
|
nrrdTypeDouble, |
358 |
|
|
tt->val, nrrdTypeDouble, |
359 |
|
|
tt->valNum, |
360 |
|
|
0, tt->rayLen); |
361 |
|
|
answer = AIR_EXISTS(answer) ? answer : uu->fromNaN; |
362 |
|
|
rr->imgData[(tt->ui) + (rr->sx)*(tt->vi)] = answer; |
363 |
|
|
} else { |
364 |
|
|
rr->imgData[(tt->ui) + (rr->sx)*(tt->vi)] = 0.0; |
365 |
|
|
} |
366 |
|
|
|
367 |
|
|
return 0; |
368 |
|
|
} |
369 |
|
|
|
370 |
|
|
/* -------------------------------------------------------------- */ |
371 |
|
|
|
372 |
|
|
double |
373 |
|
|
mrendSample(mrendThread *tt, mrendRender *rr, mrendUser *uu, |
374 |
|
|
int num, double rayT, |
375 |
|
|
int inside, |
376 |
|
|
double samplePosWorld[3], |
377 |
|
|
double samplePosIndex[3]) { |
378 |
|
|
static const char me[]="mrendSample"; |
379 |
|
|
|
380 |
|
|
AIR_UNUSED(rr); |
381 |
|
|
AIR_UNUSED(uu); |
382 |
|
|
AIR_UNUSED(num); |
383 |
|
|
AIR_UNUSED(rayT); |
384 |
|
|
AIR_UNUSED(samplePosWorld); |
385 |
|
|
|
386 |
|
|
if (tt->verbose) { |
387 |
|
|
fprintf(stderr, "%s: wrld(%g,%g,%g) -> indx(%g,%g,%g) -> %s\n", me, |
388 |
|
|
samplePosWorld[0], samplePosWorld[1], samplePosWorld[2], |
389 |
|
|
samplePosIndex[0], samplePosIndex[1], samplePosIndex[2], |
390 |
|
|
inside ? "INSIDE" : "(outside)"); |
391 |
|
|
} |
392 |
|
|
if (inside) { |
393 |
|
|
if (gageProbe(tt->gctx, |
394 |
|
|
samplePosIndex[0], |
395 |
|
|
samplePosIndex[1], |
396 |
|
|
samplePosIndex[2])) { |
397 |
|
|
biffAddf(MREND, "%s: gage trouble: %s (%d)", me, |
398 |
|
|
tt->gctx->errStr, tt->gctx->errNum); |
399 |
|
|
return AIR_NAN; |
400 |
|
|
} |
401 |
|
|
if (tt->verbose) { |
402 |
|
|
fprintf(stderr, "%s: val[%d] = %g\n", me, |
403 |
|
|
tt->valNum, *(tt->answer)); |
404 |
|
|
} |
405 |
|
|
tt->val[tt->valNum++] = *(tt->answer); |
406 |
|
|
if (tt->verbose) { |
407 |
|
|
fprintf(stderr, " ........ %g\n", tt->val[tt->valNum-1]); |
408 |
|
|
} |
409 |
|
|
tt->numSamples++; |
410 |
|
|
} |
411 |
|
|
|
412 |
|
|
return tt->rayStep; |
413 |
|
|
} |
414 |
|
|
|
415 |
|
|
/* -------------------------------------------------------------- */ |
416 |
|
|
|
417 |
|
|
#if 0 |
418 |
|
|
|
419 |
|
|
this was nixed once mrender learned to handle volume of general |
420 |
|
|
kind, instead of being restricted to scalars |
421 |
|
|
|
422 |
|
|
|
423 |
|
|
/* |
424 |
|
|
** learned: if you're playing games with strings with two passes, where |
425 |
|
|
** you first generate the set of strings in order to calculate their |
426 |
|
|
** cumulative length, and then (2nd pass) concatenate the strings |
427 |
|
|
** together, be very sure that the generation of the strings on the |
428 |
|
|
** two passes is identical. Had a very troublesome memory error because |
429 |
|
|
** I was using short version of the description string to determine |
430 |
|
|
** allocation, and then the long version in the second pass |
431 |
|
|
*/ |
432 |
|
|
char * |
433 |
|
|
mrendGage(char *prefix) { |
434 |
|
|
char *line, *ret; |
435 |
|
|
int i, len; |
436 |
|
|
|
437 |
|
|
/* 1st pass through- determine needed buffer size */ |
438 |
|
|
len = 0; |
439 |
|
|
for (i=airEnumUnknown(gageScl)+1; !airEnumValCheck(gageScl, i); i++) { |
440 |
|
|
if (1 == gageKindScl->table[i].answerLength) { |
441 |
|
|
line = airEnumFmtDesc(gageScl, i, AIR_FALSE, "\n \b\bo \"%s\": %s"); |
442 |
|
|
len += strlen(line); |
443 |
|
|
free(line); |
444 |
|
|
} |
445 |
|
|
} |
446 |
|
|
ret = AIR_CALLOC(strlen(prefix) + len + 1, char); |
447 |
|
|
if (ret) { |
448 |
|
|
strcpy(ret, prefix); |
449 |
|
|
/* 2nd pass through: create output */ |
450 |
|
|
for (i=airEnumUnknown(gageScl)+1; !airEnumValCheck(gageScl, i); i++) { |
451 |
|
|
if (1 == gageKindScl->table[i].answerLength) { |
452 |
|
|
line = airEnumFmtDesc(gageScl, i, AIR_FALSE, "\n \b\bo \"%s\": %s"); |
453 |
|
|
strcat(ret, line); |
454 |
|
|
free(line); |
455 |
|
|
} |
456 |
|
|
} |
457 |
|
|
} |
458 |
|
|
return ret; |
459 |
|
|
} |
460 |
|
|
|
461 |
|
|
#endif |
462 |
|
|
|
463 |
|
|
int |
464 |
|
|
main(int argc, const char *argv[]) { |
465 |
|
|
hestOpt *hopt=NULL; |
466 |
|
|
hestParm *hparm; |
467 |
|
|
int E, Ecode, Ethread, renorm, offfr; |
468 |
|
|
const char *me; |
469 |
|
|
char *errS, *whatS; |
470 |
|
|
mrendUser *uu; |
471 |
|
|
float isScale; |
472 |
|
|
airArray *mop; |
473 |
|
|
double gmc, turn, eye[3], eyedist; |
474 |
|
|
|
475 |
|
|
me = argv[0]; |
476 |
|
|
mop = airMopNew(); |
477 |
|
|
hparm = hestParmNew(); |
478 |
|
|
hparm->respFileEnable = AIR_TRUE; |
479 |
|
|
uu = mrendUserNew(); |
480 |
|
|
|
481 |
|
|
airMopAdd(mop, hparm, (airMopper)hestParmFree, airMopAlways); |
482 |
|
|
airMopAdd(mop, uu, (airMopper)mrendUserNix, airMopAlways); |
483 |
|
|
|
484 |
|
|
hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &(uu->nin), NULL, |
485 |
|
|
"input nrrd to render", NULL, NULL, nrrdHestNrrd); |
486 |
|
|
hestOptAdd(&hopt, "k", "kind", airTypeOther, 1, 1, &(uu->kind), NULL, |
487 |
|
|
"\"kind\" of volume (\"scalar\", \"vector\", or \"tensor\")", |
488 |
|
|
NULL, NULL, meetHestGageKind); |
489 |
|
|
limnHestCameraOptAdd(&hopt, uu->hctx->cam, |
490 |
|
|
NULL, "0 0 0", "0 0 1", |
491 |
|
|
NULL, NULL, NULL, |
492 |
|
|
"nan nan", "nan nan", "20"); |
493 |
|
|
hestOptAdd(&hopt, "offfr", NULL, airTypeInt, 0, 0, &offfr, NULL, |
494 |
|
|
"the given eye point (\"-fr\") is to be interpreted " |
495 |
|
|
"as an offset from the at point."); |
496 |
|
|
hestOptAdd(&hopt, "turn", "angle", airTypeDouble, 1, 1, &turn, "0.0", |
497 |
|
|
"angle (degrees) by which to rotate the from point around " |
498 |
|
|
"true up, for making stereo pairs. Positive means move " |
499 |
|
|
"towards positive U (the right)"); |
500 |
|
|
hestOptAdd(&hopt, "is", "image size", airTypeInt, 2, 2, uu->hctx->imgSize, |
501 |
|
|
"256 256", "image dimensions"); |
502 |
|
|
hestOptAdd(&hopt, "iss", "scale", airTypeFloat, 1, 1, &isScale, "1.0", |
503 |
|
|
"scaling of image size (from \"is\")"); |
504 |
|
|
hestOptAdd(&hopt, "k00", "kernel", airTypeOther, 1, 1, |
505 |
|
|
&(uu->ksp[gageKernel00]), "tent", |
506 |
|
|
"value reconstruction kernel", |
507 |
|
|
NULL, NULL, nrrdHestKernelSpec); |
508 |
|
|
hestOptAdd(&hopt, "k11", "kernel", airTypeOther, 1, 1, |
509 |
|
|
&(uu->ksp[gageKernel11]), "cubicd:1,0", |
510 |
|
|
"first derivative kernel", |
511 |
|
|
NULL, NULL, nrrdHestKernelSpec); |
512 |
|
|
hestOptAdd(&hopt, "k22", "kernel", airTypeOther, 1, 1, |
513 |
|
|
&(uu->ksp[gageKernel22]), "cubicdd:1,0", |
514 |
|
|
"second derivative kernel", |
515 |
|
|
NULL, NULL, nrrdHestKernelSpec); |
516 |
|
|
hestOptAdd(&hopt, "rn", NULL, airTypeBool, 0, 0, &renorm, NULL, |
517 |
|
|
"renormalize kernel weights at each new sample location. " |
518 |
|
|
"\"Accurate\" kernels don't need this; doing it always " |
519 |
|
|
"makes things go slower"); |
520 |
|
|
hestOptAdd(&hopt, "q", "query", airTypeString, 1, 1, &whatS, NULL, |
521 |
|
|
"the quantity (scalar, vector, or matrix) to learn by probing"); |
522 |
|
|
hestOptAdd(&hopt, "m", "measure", airTypeEnum, 1, 1, &(uu->measr), NULL, |
523 |
|
|
"how to collapse list of ray samples into one scalar. " |
524 |
|
|
NRRD_MEASURE_DESC, |
525 |
|
|
NULL, nrrdMeasure); |
526 |
|
|
hestOptAdd(&hopt, "gmc", "min gradmag", airTypeDouble, 1, 1, &gmc, "0.0", |
527 |
|
|
"For curvature-related queries, set answer to zero when " |
528 |
|
|
"gradient magnitude is below this"); |
529 |
|
|
hestOptAdd(&hopt, "fn", "from nan", airTypeDouble, 1, 1, &(uu->fromNaN), |
530 |
|
|
"nan", "When histo-based measures generate NaN answers, the " |
531 |
|
|
"value that should be substituted for NaN."); |
532 |
|
|
hestOptAdd(&hopt, "step", "size", airTypeDouble, 1, 1, &(uu->rayStep), |
533 |
|
|
"0.01", "step size along ray in world space"); |
534 |
|
|
hestOptAdd(&hopt, "nt", "# threads", airTypeInt, 1, 1, |
535 |
|
|
&(uu->hctx->numThreads), |
536 |
|
|
"1", "number of threads hoover should use"); |
537 |
|
|
hestOptAdd(&hopt, "vp", "img coords", airTypeInt, 2, 2, &(uu->verbPixel), |
538 |
|
|
"-1 -1", "pixel coordinates for which to turn on all verbose " |
539 |
|
|
"debugging messages, or \"-1 -1\" to disable this."); |
540 |
|
|
hestOptAdd(&hopt, "o", "filename", airTypeString, 1, 1, &(uu->outS), "-", |
541 |
|
|
"file to write output nrrd to. Defaults to stdout (\"-\")."); |
542 |
|
|
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways); |
543 |
|
|
|
544 |
|
|
hestParseOrDie(hopt, argc-1, argv+1, hparm, |
545 |
|
|
me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE); |
546 |
|
|
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways); |
547 |
|
|
|
548 |
|
|
uu->hctx->imgSize[0] = AIR_CAST(int, isScale*uu->hctx->imgSize[0]); |
549 |
|
|
uu->hctx->imgSize[1] = AIR_CAST(int, isScale*uu->hctx->imgSize[1]); |
550 |
|
|
|
551 |
|
|
uu->whatq = airEnumVal(uu->kind->enm, whatS); |
552 |
|
|
if (-1 == uu->whatq) { |
553 |
|
|
/* -1 indeed always means "unknown" for any gageKind */ |
554 |
|
|
fprintf(stderr, "%s: couldn't parse \"%s\" as measure of \"%s\" volume\n", |
555 |
|
|
me, whatS, uu->kind->name); |
556 |
|
|
hestUsage(stderr, hopt, me, hparm); |
557 |
|
|
hestGlossary(stderr, hopt, hparm); |
558 |
|
|
airMopError(mop); |
559 |
|
|
return 1; |
560 |
|
|
} |
561 |
|
|
if (mrendUserCheck(uu)) { |
562 |
|
|
fprintf(stderr, "%s: problem with input parameters:\n%s\n", |
563 |
|
|
me, errS = biffGetDone(MREND)); free(errS); |
564 |
|
|
airMopError(mop); |
565 |
|
|
return 1; |
566 |
|
|
} |
567 |
|
|
if (gageShapeSet(uu->shape, uu->nin, uu->kind->baseDim)) { |
568 |
|
|
fprintf(stderr, "%s: problem with shape:\n%s\n", |
569 |
|
|
me, errS = biffGetDone(GAGE)); free(errS); |
570 |
|
|
airMopError(mop); |
571 |
|
|
return 1; |
572 |
|
|
} |
573 |
|
|
|
574 |
|
|
gageParmSet(uu->gctx0, gageParmGradMagCurvMin, gmc); |
575 |
|
|
/* gageParmSet(uu->gctx0, gageParmRequireAllSpacings, AIR_FALSE); */ |
576 |
|
|
gageParmSet(uu->gctx0, gageParmRenormalize, renorm); |
577 |
|
|
fprintf(stderr, "%s: will render %s of %s in %s volume\n", me, |
578 |
|
|
airEnumStr(nrrdMeasure, uu->measr), |
579 |
|
|
airEnumStr(uu->kind->enm, uu->whatq), uu->kind->name); |
580 |
|
|
|
581 |
|
|
if (offfr) { |
582 |
|
|
ELL_3V_INCR(uu->hctx->cam->from, uu->hctx->cam->at); |
583 |
|
|
} |
584 |
|
|
if (limnCameraAspectSet(uu->hctx->cam, |
585 |
|
|
uu->hctx->imgSize[0], uu->hctx->imgSize[1], |
586 |
|
|
nrrdCenterCell) |
587 |
|
|
|| limnCameraUpdate(uu->hctx->cam)) { |
588 |
|
|
airMopAdd(mop, errS = biffGetDone(LIMN), airFree, airMopAlways); |
589 |
|
|
fprintf(stderr, "%s: trouble setting camera:\n%s\n", me, errS); |
590 |
|
|
airMopError(mop); |
591 |
|
|
return 1; |
592 |
|
|
} |
593 |
|
|
if (turn) { |
594 |
|
|
turn *= AIR_PI/180; |
595 |
|
|
ELL_3V_SUB(eye, uu->hctx->cam->from, uu->hctx->cam->at); |
596 |
|
|
ELL_3V_NORM(eye, eye, eyedist); |
597 |
|
|
ELL_3V_SCALE_ADD2(uu->hctx->cam->from, |
598 |
|
|
cos(turn), eye, |
599 |
|
|
sin(turn), uu->hctx->cam->U); |
600 |
|
|
ELL_3V_SCALE(uu->hctx->cam->from, eyedist, uu->hctx->cam->from); |
601 |
|
|
if (limnCameraUpdate(uu->hctx->cam)) { |
602 |
|
|
airMopAdd(mop, errS = biffGetDone(LIMN), airFree, airMopAlways); |
603 |
|
|
fprintf(stderr, "%s: trouble setting camera (again):\n%s\n", me, errS); |
604 |
|
|
airMopError(mop); |
605 |
|
|
return 1; |
606 |
|
|
} |
607 |
|
|
} |
608 |
|
|
/* |
609 |
|
|
fprintf(stderr, "%s: camera info\n", me); |
610 |
|
|
fprintf(stderr, " U = {%g,%g,%g}\n", |
611 |
|
|
uu->hctx->cam->U[0], uu->hctx->cam->U[1], uu->hctx->cam->U[2]); |
612 |
|
|
fprintf(stderr, " V = {%g,%g,%g}\n", |
613 |
|
|
uu->hctx->cam->V[0], uu->hctx->cam->V[1], uu->hctx->cam->V[2]); |
614 |
|
|
fprintf(stderr, " N = {%g,%g,%g}\n", |
615 |
|
|
uu->hctx->cam->N[0], uu->hctx->cam->N[1], uu->hctx->cam->N[2]); |
616 |
|
|
*/ |
617 |
|
|
|
618 |
|
|
/* set remaining fields of hoover context */ |
619 |
|
|
if (0) { |
620 |
|
|
/* this is the old code that assumed everything about orientation |
621 |
|
|
was determined by per-axis spacing (pre-gageShape) */ |
622 |
|
|
unsigned int base; |
623 |
|
|
base = uu->kind->baseDim; |
624 |
|
|
uu->hctx->volSize[0] = uu->nin->axis[base+0].size; |
625 |
|
|
uu->hctx->volSize[1] = uu->nin->axis[base+1].size; |
626 |
|
|
uu->hctx->volSize[2] = uu->nin->axis[base+2].size; |
627 |
|
|
uu->hctx->volSpacing[0] = uu->nin->axis[base+0].spacing; |
628 |
|
|
uu->hctx->volSpacing[1] = uu->nin->axis[base+1].spacing; |
629 |
|
|
uu->hctx->volSpacing[2] = uu->nin->axis[base+2].spacing; |
630 |
|
|
if (nrrdCenterUnknown != uu->nin->axis[base].center) { |
631 |
|
|
uu->hctx->volCentering = uu->nin->axis[base].center; |
632 |
|
|
fprintf(stderr, "%s: setting volCentering to %s\n", me, |
633 |
|
|
airEnumStr(nrrdCenter, uu->nin->axis[base].center)); |
634 |
|
|
} |
635 |
|
|
fprintf(stderr, "!%s: uu->hctx->volCentering = %d\n", |
636 |
|
|
me, uu->hctx->volCentering); |
637 |
|
|
} else { |
638 |
|
|
uu->hctx->shape = uu->shape; |
639 |
|
|
} |
640 |
|
|
/* this is reasonable for now */ |
641 |
|
|
uu->hctx->imgCentering = nrrdCenterCell; |
642 |
|
|
uu->hctx->user = uu; |
643 |
|
|
uu->hctx->renderBegin = (hooverRenderBegin_t *)mrendRenderBegin; |
644 |
|
|
uu->hctx->threadBegin = (hooverThreadBegin_t *)mrendThreadBegin; |
645 |
|
|
uu->hctx->rayBegin = (hooverRayBegin_t *)mrendRayBegin; |
646 |
|
|
uu->hctx->sample = (hooverSample_t *)mrendSample; |
647 |
|
|
uu->hctx->rayEnd = (hooverRayEnd_t *)mrendRayEnd; |
648 |
|
|
uu->hctx->threadEnd = (hooverThreadEnd_t *)mrendThreadEnd; |
649 |
|
|
uu->hctx->renderEnd = (hooverRenderEnd_t *)mrendRenderEnd; |
650 |
|
|
|
651 |
|
|
if (!airThreadCapable && 1 != uu->hctx->numThreads) { |
652 |
|
|
fprintf(stderr, "%s: This Teem not compiled with " |
653 |
|
|
"multi-threading support.\n", me); |
654 |
|
|
fprintf(stderr, "%s: --> can't use %d threads; only using 1\n", |
655 |
|
|
me, uu->hctx->numThreads); |
656 |
|
|
uu->hctx->numThreads = 1; |
657 |
|
|
} |
658 |
|
|
|
659 |
|
|
E = hooverRender(uu->hctx, &Ecode, &Ethread); |
660 |
|
|
if (E) { |
661 |
|
|
if (hooverErrInit == E) { |
662 |
|
|
fprintf(stderr, "%s: ERROR (code %d, thread %d):\n%s\n", |
663 |
|
|
me, Ecode, Ethread, errS = biffGetDone(HOOVER)); |
664 |
|
|
free(errS); |
665 |
|
|
} else { |
666 |
|
|
fprintf(stderr, "%s: ERROR (code %d, thread %d):\n%s\n", |
667 |
|
|
me, Ecode, Ethread, errS = biffGetDone(MREND)); |
668 |
|
|
free(errS); |
669 |
|
|
} |
670 |
|
|
airMopError(mop); |
671 |
|
|
return 1; |
672 |
|
|
} |
673 |
|
|
|
674 |
|
|
if (1) { |
675 |
|
|
ELL_3V_SUB(uu->imgU, uu->imgU, uu->imgOrig); |
676 |
|
|
ELL_3V_SUB(uu->imgV, uu->imgV, uu->imgOrig); |
677 |
|
|
fprintf(stderr, "%s: loc(0,0) = (%g,%g,%g)\n", me, |
678 |
|
|
uu->imgOrig[0], uu->imgOrig[1], uu->imgOrig[2]); |
679 |
|
|
fprintf(stderr, "%s: loc(1,0) - loc(0,0) = (%g,%g,%g)\n", me, |
680 |
|
|
uu->imgU[0], uu->imgU[1], uu->imgU[2]); |
681 |
|
|
fprintf(stderr, "%s: loc(0,1) - loc(0,0) = (%g,%g,%g)\n", me, |
682 |
|
|
uu->imgV[0], uu->imgV[1], uu->imgV[2]); |
683 |
|
|
} |
684 |
|
|
|
685 |
|
|
airMopOkay(mop); |
686 |
|
|
return 0; |
687 |
|
|
} |