GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/mite/ray.c Lines: 0 160 0.0 %
Date: 2017-05-26 Branches: 0 84 0.0 %

Line Branch Exec Source
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 "mite.h"
25
#include "privateMite.h"
26
27
int
28
miteRayBegin(miteThread *mtt, miteRender *mrr, miteUser *muu,
29
             int uIndex, int vIndex,
30
             double rayLen,
31
             double rayStartWorld[3], double rayStartIndex[3],
32
             double rayDirWorld[3], double rayDirIndex[3]) {
33
  airPtrPtrUnion appu;
34
  AIR_UNUSED(mrr);
35
  AIR_UNUSED(rayStartWorld);
36
  AIR_UNUSED(rayStartIndex);
37
  AIR_UNUSED(rayDirIndex);
38
39
  mtt->ui = uIndex;
40
  mtt->vi = vIndex;
41
  mtt->rayStep = (muu->rayStep*rayLen /
42
                  (muu->hctx->cam->vspFaar - muu->hctx->cam->vspNeer));
43
  if (!uIndex) {
44
    fprintf(stderr, "%d/%d ", vIndex, muu->hctx->imgSize[1]);
45
    fflush(stderr);
46
  }
47
  mtt->verbose = (uIndex == muu->verbUi && vIndex == muu->verbVi);
48
  mtt->skip = (muu->verbUi >= 0 && muu->verbVi >= 0
49
               && !mtt->verbose);
50
  if (mtt->verbose) {
51
    /* create muu->ndebug */
52
    muu->ndebug = nrrdNew();
53
    /* we want to store the value and index for each txf domain variable,
54
       plus the RGBAZ computed for that sample */
55
    muu->ndebug->axis[0].size = 2*mtt->stageNum + 5;
56
    /* we really do want to associate ndebug with the miteUser's mop,
57
       because the information stored in it has to persist for as long as
58
       the user wants: mite itself doesn't call miteUserNix */
59
    airMopAdd(muu->umop, muu->ndebug, (airMopper)nrrdNuke, airMopAlways);
60
    /* but the scope of the debug array allocation is within this ray */
61
    muu->debugArr = airArrayNew((appu.d = &(muu->debug), appu.v),
62
                                NULL, sizeof(double), 128);
63
  }
64
  mtt->raySample = 0;
65
  mtt->RR = mtt->GG = mtt->BB = 0.0;
66
  mtt->TT = 1.0;
67
  mtt->ZZ = AIR_NAN;
68
  ELL_3V_SCALE(mtt->V, -1, rayDirWorld);
69
70
  return 0;
71
}
72
73
void
74
_miteRGBACalc(mite_t *R, mite_t *G, mite_t *B, mite_t *A,
75
              miteThread *mtt, miteRender *mrr, miteUser *muu) {
76
  static const char me[]="_miteRGBACalc";
77
  mite_t tmp,
78
    ad[3],                          /* ambient+diffuse light contribution */
79
    s[3] = {0,0,0},                 /* specular light contribution */
80
    col[3], E, ka, kd, ks, sp,      /* txf-determined rendering variables */
81
    LdotN=0, HdotN, H[3], N[3];     /* for lighting calculation */
82
83
  col[0] = mtt->range[miteRangeRed];
84
  col[1] = mtt->range[miteRangeGreen];
85
  col[2] = mtt->range[miteRangeBlue];
86
  E = mtt->range[miteRangeEmissivity];
87
  ka = mtt->range[miteRangeKa];
88
  kd = mtt->range[miteRangeKd];
89
  ks = mtt->range[miteRangeKs];
90
  ELL_3V_SCALE(ad, ka, muu->lit->amb);
91
  switch (mrr->shadeSpec->method) {
92
  case miteShadeMethodNone:
93
    /* nothing to do */
94
    break;
95
  case miteShadeMethodPhong:
96
    if (kd || ks) {
97
      ELL_3V_NORM(N, mtt->shadeVec0, tmp);
98
      if (1 == muu->normalSide) {
99
        ELL_3V_SCALE(N, -1, N);
100
      }
101
      /* else -1==side --> N = -1*-1*N = N
102
         or 0==side --> N = N, so there's nothing to do */
103
      if (kd) {
104
        LdotN = ELL_3V_DOT(muu->lit->dir[0], N);
105
        if (!muu->normalSide) {
106
          LdotN = AIR_ABS(LdotN);
107
        }
108
        if (LdotN > 0) {
109
          ELL_3V_SCALE_INCR(ad, LdotN*kd, muu->lit->col[0]);
110
        }
111
      }
112
      if (ks) {
113
        sp = mtt->range[miteRangeSP];
114
        ELL_3V_ADD2(H, muu->lit->dir[0], mtt->V);
115
        ELL_3V_NORM(H, H, tmp);
116
        HdotN = ELL_3V_DOT(H, N);
117
        if (!muu->normalSide) {
118
          HdotN = AIR_ABS(HdotN);
119
        }
120
        if (HdotN > 0) {
121
          HdotN = pow(HdotN, sp);
122
          ELL_3V_SCALE(s, HdotN*ks, muu->lit->col[0]);
123
        }
124
      }
125
    }
126
    break;
127
  case miteShadeMethodLitTen:
128
    fprintf(stderr, "!%s: lit-tensor not yet implemented\n", me);
129
    break;
130
  default:
131
    fprintf(stderr, "!%s: PANIC, shadeMethod %d unimplemented\n",
132
            me, mrr->shadeSpec->method);
133
    exit(1);
134
    break;
135
  }
136
  *R = (E - 1 + ad[0])*col[0] + s[0];
137
  *G = (E - 1 + ad[1])*col[1] + s[1];
138
  *B = (E - 1 + ad[2])*col[2] + s[2];
139
  *A = mtt->range[miteRangeAlpha];
140
  *A = AIR_CLAMP(0.0, *A, 1.0);
141
  /*
142
  if (mtt->verbose) {
143
    fprintf(stderr, "%s: col[] = %g,%g,%g; A,E = %g,%g; Kads = %g,%g,%g\n", me,
144
            col[0], col[1], col[2], mtt->range[miteRangeAlpha], E, ka, kd, ks);
145
    fprintf(stderr, "%s: N = (%g,%g,%g), L = (%g,%g,%g) ---> LdotN = %g\n",
146
            me, N[0], N[1], N[2], muu->lit->dir[0][0], muu->lit->dir[0][1],
147
            muu->lit->dir[0][2], LdotN);
148
    fprintf(stderr, "%s: ad[] = %g,%g,%g\n", me, ad[0], ad[1], ad[2]);
149
    fprintf(stderr, "%s:  --> R,G,B,A = %g,%g,%g,%g\n", me, *R, *G, *B, *A);
150
  }
151
  */
152
  return;
153
}
154
155
double
156
miteSample(miteThread *mtt, miteRender *mrr, miteUser *muu,
157
           int num, double rayT, int inside,
158
           double samplePosWorld[3],
159
           double samplePosIndex[3]) {
160
  static const char me[]="miteSample";
161
  mite_t R, G, B, A;
162
  double *NN;
163
  double NdotV, kn[3], knd[3], ref[3], len, *dbg=NULL;
164
165
  if (!inside) {
166
    return mtt->rayStep;
167
  }
168
169
  if (mtt->skip) {
170
    /* we have one verbose pixel, but we're not on it */
171
    return 0.0;
172
  }
173
174
  /* early ray termination */
175
  if (1-mtt->TT >= muu->opacNear1) {
176
    mtt->TT = 0.0;
177
    return 0.0;
178
  }
179
180
  /* set (fake) view based on fake from */
181
  if (AIR_EXISTS(muu->fakeFrom[0])) {
182
    ELL_3V_SUB(mtt->V, samplePosWorld, muu->fakeFrom);
183
    ELL_3V_NORM(mtt->V, mtt->V, len);
184
  }
185
186
  /* do probing at this location to determine values of everything
187
     that might appear in the txf domain */
188
  if (gageProbe(mtt->gctx,
189
                samplePosIndex[0],
190
                samplePosIndex[1],
191
                samplePosIndex[2])) {
192
    biffAddf(MITE, "%s: gage trouble: %s (%d)", me,
193
             mtt->gctx->errStr, mtt->gctx->errNum);
194
    return AIR_NAN;
195
  }
196
197
  if (mrr->queryMiteNonzero) {
198
    /* There is some optimal trade-off between slowing things down
199
       with too many branches on all possible checks of queryMite,
200
       and slowing things down with doing the work of setting them all.
201
       This code has not been profiled whatsoever */
202
    mtt->directAnsMiteVal[miteValXw][0] = samplePosWorld[0];
203
    mtt->directAnsMiteVal[miteValXi][0] = samplePosIndex[0];
204
    mtt->directAnsMiteVal[miteValYw][0] = samplePosWorld[1];
205
    mtt->directAnsMiteVal[miteValYi][0] = samplePosIndex[1];
206
    mtt->directAnsMiteVal[miteValZw][0] = samplePosWorld[2];
207
    mtt->directAnsMiteVal[miteValZi][0] = samplePosIndex[2];
208
    mtt->directAnsMiteVal[miteValRw][0] = ELL_3V_LEN(samplePosWorld);
209
    mtt->directAnsMiteVal[miteValRi][0] = ELL_3V_LEN(samplePosIndex);
210
    mtt->directAnsMiteVal[miteValTw][0] = rayT;
211
    mtt->directAnsMiteVal[miteValTi][0] = num;
212
    ELL_3V_COPY(mtt->directAnsMiteVal[miteValView], mtt->V);
213
    NN = mtt->directAnsMiteVal[miteValNormal];
214
    if (mtt->_normal) {
215
      if (1 == muu->normalSide) {
216
        ELL_3V_SCALE(NN, -1, mtt->_normal);
217
      } else {
218
        ELL_3V_COPY(NN, mtt->_normal);
219
      }
220
    }
221
222
    if ((GAGE_QUERY_ITEM_TEST(mrr->queryMite, miteValNdotV)
223
         || GAGE_QUERY_ITEM_TEST(mrr->queryMite, miteValNdotL)
224
         || GAGE_QUERY_ITEM_TEST(mrr->queryMite, miteValVrefN))) {
225
      mtt->directAnsMiteVal[miteValNdotV][0] = ELL_3V_DOT(NN, mtt->V);
226
      mtt->directAnsMiteVal[miteValNdotL][0] =
227
        ELL_3V_DOT(NN, muu->lit->dir[0]);
228
      if (!muu->normalSide) {
229
        mtt->directAnsMiteVal[miteValNdotV][0] =
230
          AIR_ABS(mtt->directAnsMiteVal[miteValNdotV][0]);
231
        mtt->directAnsMiteVal[miteValNdotL][0] =
232
          AIR_ABS(mtt->directAnsMiteVal[miteValNdotL][0]);
233
      }
234
      NdotV = mtt->directAnsMiteVal[miteValNdotV][0];
235
      ELL_3V_SCALE_ADD2(ref, 2*NdotV, NN, -1, mtt->V);
236
      ELL_3V_NORM(mtt->directAnsMiteVal[miteValVrefN], ref, len);
237
    }
238
239
    if (GAGE_QUERY_ITEM_TEST(mrr->queryMite, miteValGTdotV)) {
240
      ELL_3MV_MUL(kn, mtt->nPerp, mtt->V);
241
      ELL_3V_NORM(kn, kn, len);
242
      ELL_3MV_MUL(knd, mtt->geomTens, kn);
243
      mtt->directAnsMiteVal[miteValGTdotV][0] = ELL_3V_DOT(knd, kn);
244
    }
245
  }
246
247
248
  /* initialize txf range quantities, and apply all txfs */
249
  if (mtt->verbose) {
250
    muu->debugIdx = airArrayLenIncr(muu->debugArr, muu->ndebug->axis[0].size);
251
  }
252
253
  memcpy(mtt->range, muu->rangeInit, MITE_RANGE_NUM*sizeof(mite_t));
254
  _miteStageRun(mtt, muu);
255
256
  /* if there's opacity, do shading and compositing */
257
  if (mtt->range[miteRangeAlpha]) {
258
    /* fprintf(stderr, "%s: mtt->TT = %g\n", me, mtt->TT); */
259
    /*
260
    if (mtt->verbose) {
261
      fprintf(stderr, "%s: before compositing: RGBT = %g,%g,%g,%g\n",
262
              me, mtt->RR, mtt->GG, mtt->BB, mtt->TT);
263
    }
264
    */
265
    _miteRGBACalc(&R, &G, &B, &A, mtt, mrr, muu);
266
    mtt->RR += mtt->TT*A*R;
267
    mtt->GG += mtt->TT*A*G;
268
    mtt->BB += mtt->TT*A*B;
269
    mtt->TT *= 1-A;
270
    /*
271
    if (mtt->verbose) {
272
      fprintf(stderr, "%s: after compositing: RGBT = %g,%g,%g,%g\n",
273
              me, mtt->RR, mtt->GG, mtt->BB, mtt->TT);
274
    }
275
    */
276
    /* fprintf(stderr, "%s: mtt->TT = %g\n", me, mtt->TT); */
277
  } else {
278
    R = G = B = A = 0;
279
  }
280
  if (mtt->verbose) {
281
    dbg = muu->debug + muu->debugIdx;
282
    dbg[0 + 2*mtt->stageNum] = R;
283
    dbg[1 + 2*mtt->stageNum] = G;
284
    dbg[2 + 2*mtt->stageNum] = B;
285
    dbg[3 + 2*mtt->stageNum] = A;
286
    dbg[4 + 2*mtt->stageNum] = rayT;
287
  }
288
289
  /* set Z if it hasn't been set already */
290
  if (1-mtt->TT >= muu->opacMatters && !AIR_EXISTS(mtt->ZZ)) {
291
    mtt->ZZ = rayT;
292
  }
293
294
  /* this is used to index mtt->debug */
295
  mtt->raySample += 1;
296
297
  return mtt->rayStep;
298
}
299
300
int
301
miteRayEnd(miteThread *mtt, miteRender *mrr, miteUser *muu) {
302
  int idx, slen, stageIdx;
303
  mite_t *imgData;
304
  double A;
305
306
  AIR_UNUSED(mrr);
307
  mtt->samples += mtt->raySample;
308
  idx = mtt->ui + (muu->nout->axis[1].size)*mtt->vi;
309
  imgData = (mite_t*)muu->nout->data;
310
  A = 1 - mtt->TT;
311
  if (A) {
312
    ELL_5V_SET(imgData + 5*idx, mtt->RR/A, mtt->GG/A, mtt->BB/A,
313
               A, mtt->ZZ);
314
  } else {
315
    ELL_5V_SET(imgData + 5*idx, 0, 0, 0, 0, AIR_NAN);
316
  }
317
  if (mtt->verbose) {
318
    /* muu->debug may be over-allocated, but that's harmless */
319
    muu->ndebug->axis[1].size = mtt->raySample;
320
    nrrdWrap_va(muu->ndebug, muu->debug, nrrdTypeDouble, 2,
321
                AIR_CAST(size_t, muu->ndebug->axis[0].size),
322
                AIR_CAST(size_t, mtt->raySample));
323
    airArrayNix(muu->debugArr);
324
    slen = 0;
325
    for (stageIdx=0; stageIdx<mtt->stageNum; stageIdx++) {
326
      slen += strlen(mtt->stage[stageIdx].label) + 2;
327
    }
328
    slen += strlen("R,G,B,A,Z") + 1;
329
    muu->ndebug->axis[0].label = (char *)calloc(slen, sizeof(char));
330
    for (stageIdx=0; stageIdx<mtt->stageNum; stageIdx++) {
331
      strcat(muu->ndebug->axis[0].label, mtt->stage[stageIdx].label);
332
      strcat(muu->ndebug->axis[0].label, ",,");
333
    }
334
    strcat(muu->ndebug->axis[0].label, "R,G,B,A,Z");
335
  }
336
  return 0;
337
}
338