| File: | src/echo/sqd.c |
| Location: | line 368, column 9 |
| Description: | Value stored to 'Vmax' is never read |
| 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 "echo.h" |
| 25 | #include "privateEcho.h" |
| 26 | |
| 27 | echoPos_t |
| 28 | _echo_SuperquadX_v(echoPos_t x, echoPos_t y, echoPos_t z, |
| 29 | echoPos_t A, echoPos_t B) { |
| 30 | echoPos_t xxb, yya, zza; |
| 31 | |
| 32 | xxb = pow(x*x, 1/B); yya = pow(y*y, 1/A); zza = pow(z*z, 1/A); |
| 33 | return pow(yya + zza, A/B) + xxb - 1; |
| 34 | } |
| 35 | |
| 36 | echoPos_t |
| 37 | _echo_SuperquadY_v(echoPos_t x, echoPos_t y, echoPos_t z, |
| 38 | echoPos_t A, echoPos_t B) { |
| 39 | echoPos_t xxa, yyb, zza; |
| 40 | |
| 41 | xxa = pow(x*x, 1/A); yyb = pow(y*y, 1/B); zza = pow(z*z, 1/A); |
| 42 | return pow(xxa + zza, A/B) + yyb - 1; |
| 43 | } |
| 44 | |
| 45 | echoPos_t |
| 46 | _echo_SuperquadZ_v(echoPos_t x, echoPos_t y, echoPos_t z, |
| 47 | echoPos_t A, echoPos_t B) { |
| 48 | echoPos_t xxa, yya, zzb; |
| 49 | |
| 50 | xxa = pow(x*x, 1/A); yya = pow(y*y, 1/A); zzb = pow(z*z, 1/B); |
| 51 | return pow(xxa + yya, A/B) + zzb - 1; |
| 52 | } |
| 53 | |
| 54 | /* -------------------------------------------------------- */ |
| 55 | |
| 56 | echoPos_t |
| 57 | _echo_SuperquadX_vg(echoPos_t grad[3], |
| 58 | echoPos_t x, echoPos_t y, echoPos_t z, |
| 59 | echoPos_t A, echoPos_t B) { |
| 60 | echoPos_t R, xxb, yya, zza; |
| 61 | |
| 62 | xxb = pow(x*x, 1/B); yya = pow(y*y, 1/A); zza = pow(z*z, 1/A); |
| 63 | R = pow(yya + zza, (A/B)-1); |
| 64 | ELL_3V_SET(grad, 2*xxb/(B*x), 2*R*yya/(B*y), 2*R*zza/(B*z))((grad)[0] = (2*xxb/(B*x)), (grad)[1] = (2*R*yya/(B*y)), (grad )[2] = (2*R*zza/(B*z))); |
| 65 | return pow(yya + zza, A/B) + xxb - 1; |
| 66 | } |
| 67 | |
| 68 | echoPos_t |
| 69 | _echo_SuperquadY_vg(echoPos_t grad[3], |
| 70 | echoPos_t x, echoPos_t y, echoPos_t z, |
| 71 | echoPos_t A, echoPos_t B) { |
| 72 | echoPos_t R, xxa, yyb, zza; |
| 73 | |
| 74 | xxa = pow(x*x, 1/A); yyb = pow(y*y, 1/B); zza = pow(z*z, 1/A); |
| 75 | R = pow(xxa + zza, (A/B)-1); |
| 76 | ELL_3V_SET(grad, 2*R*xxa/(B*x), 2*yyb/(B*y), 2*R*zza/(B*z))((grad)[0] = (2*R*xxa/(B*x)), (grad)[1] = (2*yyb/(B*y)), (grad )[2] = (2*R*zza/(B*z))); |
| 77 | return pow(xxa + zza, A/B) + yyb - 1; |
| 78 | } |
| 79 | |
| 80 | echoPos_t |
| 81 | _echo_SuperquadZ_vg(echoPos_t grad[3], |
| 82 | echoPos_t x, echoPos_t y, echoPos_t z, |
| 83 | echoPos_t A, echoPos_t B) { |
| 84 | echoPos_t R, xxa, yya, zzb; |
| 85 | |
| 86 | xxa = pow(x*x, 1/A); yya = pow(y*y, 1/A); zzb = pow(z*z, 1/B); |
| 87 | R = pow(xxa + yya, (A/B)-1); |
| 88 | ELL_3V_SET(grad, 2*R*xxa/(B*x), 2*R*yya/(B*y), 2*zzb/(B*z))((grad)[0] = (2*R*xxa/(B*x)), (grad)[1] = (2*R*yya/(B*y)), (grad )[2] = (2*zzb/(B*z))); |
| 89 | return pow(xxa + yya, A/B) + zzb - 1; |
| 90 | } |
| 91 | |
| 92 | /* -------------------------------------------------------- */ |
| 93 | |
| 94 | echoPos_t |
| 95 | _echo_SuperquadX_lvg(echoPos_t grad[3], |
| 96 | echoPos_t x, echoPos_t y, echoPos_t z, |
| 97 | echoPos_t A, echoPos_t B) { |
| 98 | echoPos_t R, xxb, yya, zza, larg; |
| 99 | |
| 100 | echoPos_t ret; |
| 101 | |
| 102 | xxb = pow(x*x, 1/B); yya = pow(y*y, 1/A); zza = pow(z*z, 1/A); |
| 103 | R = pow(yya + zza, 1-(A/B))*xxb; |
| 104 | ELL_3V_SET(grad,((grad)[0] = (2/(B*x*(1 + pow(yya + zza, A/B)/xxb))), (grad)[ 1] = (2*yya/(B*y*(yya + zza + R))), (grad)[2] = (2*zza/(B*z*( yya + zza + R)))) |
| 105 | 2/(B*x*(1 + pow(yya + zza, A/B)/xxb)),((grad)[0] = (2/(B*x*(1 + pow(yya + zza, A/B)/xxb))), (grad)[ 1] = (2*yya/(B*y*(yya + zza + R))), (grad)[2] = (2*zza/(B*z*( yya + zza + R)))) |
| 106 | 2*yya/(B*y*(yya + zza + R)),((grad)[0] = (2/(B*x*(1 + pow(yya + zza, A/B)/xxb))), (grad)[ 1] = (2*yya/(B*y*(yya + zza + R))), (grad)[2] = (2*zza/(B*z*( yya + zza + R)))) |
| 107 | 2*zza/(B*z*(yya + zza + R)))((grad)[0] = (2/(B*x*(1 + pow(yya + zza, A/B)/xxb))), (grad)[ 1] = (2*yya/(B*y*(yya + zza + R))), (grad)[2] = (2*zza/(B*z*( yya + zza + R)))); |
| 108 | larg = pow(yya + zza, A/B) + xxb; |
| 109 | ret= larg > 0 ? log(larg) : ECHO_POS_MIN(-1.7976931348623157e+308); |
| 110 | /* |
| 111 | if (!( AIR_EXISTS(grad[0]) && AIR_EXISTS(grad[1]) && AIR_EXISTS(grad[2]) )) { |
| 112 | fprintf(stderr, "_echo_SuperquadX_lvg: PANIC\n"); |
| 113 | fprintf(stderr, "x = %g, y = %g, z = %g, A = %g, B = %g\n", |
| 114 | x, y, z, A, B); |
| 115 | fprintf(stderr, "pow(%g * %g = %g, 1/%g = %g) = %g\n", |
| 116 | x, x, x*x, B, 1/B, pow(x*x, 1/B)); |
| 117 | fprintf(stderr, "xxb = %g, yya = %g, zza = %g\n", |
| 118 | xxb, yya, zza); |
| 119 | fprintf(stderr, "R: pow(%g + %g = %g, 1-(%g/%g = %g) = %g) = %g*%g = %g\n", |
| 120 | yya, zza, yya + zza, |
| 121 | A, B, A/B, 1-(A/B), |
| 122 | pow(yya + zza, 1-(A/B)), xxb, |
| 123 | pow(yya + zza, 1-(A/B))*xxb); |
| 124 | fprintf(stderr, "grad[0]: 2/(%g * %g * (1 + pow(%g + %g = %g, %g/%g = %g)/%g = %g)) = %g\n", |
| 125 | B, x, yya, zza, yya+zza, A, B, A/B, xxb, |
| 126 | pow(yya + zza, A/B)/xxb, grad[0]); |
| 127 | fprintf(stderr, "grad[1]: 2*%g/(%g*%g*(%g + %g + %g = %g) = %g) = %g\n", |
| 128 | yya, B, y, yya, zza, R, yya + zza + R, |
| 129 | B*y*(yya + zza + R), |
| 130 | 2*yya/(B*y*(yya + zza + R))); |
| 131 | |
| 132 | fprintf(stderr, "log(pow(%g + %g = %g, %g) = %g + %g) = %g\n", |
| 133 | yya, zza, yya+zza, A/B, pow(yya + zza, A/B), xxb, ret); |
| 134 | fprintf(stderr, "grad = %g %g %g\n", grad[0], grad[1], grad[2]); |
| 135 | fprintf(stderr, "\n----------\n\n"); |
| 136 | } |
| 137 | */ |
| 138 | return ret; |
| 139 | } |
| 140 | |
| 141 | echoPos_t |
| 142 | _echo_SuperquadY_lvg(echoPos_t grad[3], |
| 143 | echoPos_t x, echoPos_t y, echoPos_t z, |
| 144 | echoPos_t A, echoPos_t B) { |
| 145 | echoPos_t R, xxa, yyb, zza, larg; |
| 146 | |
| 147 | xxa = pow(x*x, 1/A); yyb = pow(y*y, 1/B); zza = pow(z*z, 1/A); |
| 148 | R = pow(xxa + zza, 1-(A/B))*yyb; |
| 149 | ELL_3V_SET(grad,((grad)[0] = (2*xxa/(B*x*(xxa + zza + R))), (grad)[1] = (2/(B *y*(1 + pow(xxa + zza, A/B)/yyb))), (grad)[2] = (2*zza/(B*z*( xxa + zza + R)))) |
| 150 | 2*xxa/(B*x*(xxa + zza + R)),((grad)[0] = (2*xxa/(B*x*(xxa + zza + R))), (grad)[1] = (2/(B *y*(1 + pow(xxa + zza, A/B)/yyb))), (grad)[2] = (2*zza/(B*z*( xxa + zza + R)))) |
| 151 | 2/(B*y*(1 + pow(xxa + zza, A/B)/yyb)),((grad)[0] = (2*xxa/(B*x*(xxa + zza + R))), (grad)[1] = (2/(B *y*(1 + pow(xxa + zza, A/B)/yyb))), (grad)[2] = (2*zza/(B*z*( xxa + zza + R)))) |
| 152 | 2*zza/(B*z*(xxa + zza + R)))((grad)[0] = (2*xxa/(B*x*(xxa + zza + R))), (grad)[1] = (2/(B *y*(1 + pow(xxa + zza, A/B)/yyb))), (grad)[2] = (2*zza/(B*z*( xxa + zza + R)))); |
| 153 | larg = pow(xxa + zza, A/B) + yyb; |
| 154 | return larg > 0 ? log(larg) : ECHO_POS_MIN(-1.7976931348623157e+308); |
| 155 | } |
| 156 | |
| 157 | echoPos_t |
| 158 | _echo_SuperquadZ_lvg(echoPos_t grad[3], |
| 159 | echoPos_t x, echoPos_t y, echoPos_t z, |
| 160 | echoPos_t A, echoPos_t B) { |
| 161 | echoPos_t R, xxa, yya, zzb, larg; |
| 162 | |
| 163 | echoPos_t ret; |
| 164 | |
| 165 | xxa = pow(x*x, 1/A); yya = pow(y*y, 1/A); zzb = pow(z*z, 1/B); |
| 166 | R = pow(xxa + yya, 1-(A/B))*zzb; |
| 167 | ELL_3V_SET(grad,((grad)[0] = (2*xxa/(B*x*(xxa + yya + R))), (grad)[1] = (2*yya /(B*y*(xxa + yya + R))), (grad)[2] = (2/(B*z*(1 + pow(xxa + yya , A/B)/zzb)))) |
| 168 | 2*xxa/(B*x*(xxa + yya + R)),((grad)[0] = (2*xxa/(B*x*(xxa + yya + R))), (grad)[1] = (2*yya /(B*y*(xxa + yya + R))), (grad)[2] = (2/(B*z*(1 + pow(xxa + yya , A/B)/zzb)))) |
| 169 | 2*yya/(B*y*(xxa + yya + R)),((grad)[0] = (2*xxa/(B*x*(xxa + yya + R))), (grad)[1] = (2*yya /(B*y*(xxa + yya + R))), (grad)[2] = (2/(B*z*(1 + pow(xxa + yya , A/B)/zzb)))) |
| 170 | 2/(B*z*(1 + pow(xxa + yya, A/B)/zzb)))((grad)[0] = (2*xxa/(B*x*(xxa + yya + R))), (grad)[1] = (2*yya /(B*y*(xxa + yya + R))), (grad)[2] = (2/(B*z*(1 + pow(xxa + yya , A/B)/zzb)))); |
| 171 | larg = pow(xxa + yya, A/B) + zzb; |
| 172 | |
| 173 | ret = larg > 0 ? log(larg) : ECHO_POS_MIN(-1.7976931348623157e+308); |
| 174 | /* |
| 175 | if (!AIR_EXISTS(ret)) { |
| 176 | fprintf(stderr, "_echo_SuperquadZ_lvg: PANIC\n"); |
| 177 | fprintf(stderr, "x = %g, y = %g, z = %g, A = %g, B = %g\n", |
| 178 | x, y, z, A, B); |
| 179 | fprintf(stderr, "pow(%g*%g = %g, %g) = %g\n", |
| 180 | x, x, x*x, 1/A, xxa); |
| 181 | fprintf(stderr, "pow(%g*%g = %g, %g) = %g\n", |
| 182 | y, y, y*y, 1/A, yya); |
| 183 | fprintf(stderr, "log(pow(%g, %g) = %g + %g) = %g\n", |
| 184 | xxa + yya, A/B, pow(xxa + yya, A/B), zzb, ret); |
| 185 | exit(0); |
| 186 | } |
| 187 | */ |
| 188 | return ret; |
| 189 | } |
| 190 | |
| 191 | /* -------------------------------------------------------- */ |
| 192 | |
| 193 | int |
| 194 | _echoRayIntx_Superquad(RAYINTX_ARGS(Superquad)echoIntx *intx, echoRay *ray, echoSuperquad *obj, echoRTParm * parm, echoThreadState *tstate) { |
| 195 | char me[]="_echoRayIntx_Superquad"; |
| 196 | echoPos_t TT=0, Tmin, Tmax, t0, t1, t2, t3, v1, v2, diff, tol, |
| 197 | saveTmin, Vmin, Vmax, VV=0, dV, dVmin, dVmax, tmp, |
| 198 | (*v)(echoPos_t, echoPos_t, echoPos_t, |
| 199 | echoPos_t, echoPos_t), |
| 200 | (*vg)(echoPos_t[3], |
| 201 | echoPos_t, echoPos_t, echoPos_t, |
| 202 | echoPos_t, echoPos_t), |
| 203 | (*lvg)(echoPos_t[3], |
| 204 | echoPos_t, echoPos_t, echoPos_t, |
| 205 | echoPos_t, echoPos_t), |
| 206 | from[3], grad[3], pos[3]; /* these two used only by macros */ |
| 207 | int iter; |
| 208 | |
| 209 | if (!_echoRayIntx_CubeSolid(&Tmin, &Tmax, |
| 210 | -1-2*ECHO_EPSILON0.00005, 1+2*ECHO_EPSILON0.00005, |
| 211 | -1-2*ECHO_EPSILON0.00005, 1+2*ECHO_EPSILON0.00005, |
| 212 | -1-2*ECHO_EPSILON0.00005, 1+2*ECHO_EPSILON0.00005, ray)) { |
| 213 | return AIR_FALSE0; |
| 214 | } |
| 215 | switch(obj->axis) { |
| 216 | case 0: |
| 217 | v = _echo_SuperquadX_v; |
| 218 | vg = _echo_SuperquadX_vg; |
| 219 | lvg = _echo_SuperquadX_lvg; |
| 220 | break; |
| 221 | case 1: |
| 222 | v = _echo_SuperquadY_v; |
| 223 | vg = _echo_SuperquadY_vg; |
| 224 | lvg = _echo_SuperquadY_lvg; |
| 225 | break; |
| 226 | case 2: default: |
| 227 | v = _echo_SuperquadZ_v; |
| 228 | vg = _echo_SuperquadZ_vg; |
| 229 | lvg = _echo_SuperquadZ_lvg; |
| 230 | break; |
| 231 | } |
| 232 | if (tstate->verbose) { |
| 233 | fprintf(stderr__stderrp, "%s%s: Tmin, Tmax = %g, %g, ax = %d\n", |
| 234 | _echoDot(tstate->depth), me, Tmin, Tmax, obj->axis); |
| 235 | } |
| 236 | |
| 237 | #define VAL(TT)(((pos)[0] = (1)*(from)[0] + ((TT))*(ray->dir)[0], (pos)[1 ] = (1)*(from)[1] + ((TT))*(ray->dir)[1], (pos)[2] = (1)*( from)[2] + ((TT))*(ray->dir)[2]), v(pos[0], pos[1], pos[2] , obj->A, obj->B)) \ |
| 238 | (ELL_3V_SCALE_ADD2(pos, 1, from, (TT), ray->dir)((pos)[0] = (1)*(from)[0] + ((TT))*(ray->dir)[0], (pos)[1] = (1)*(from)[1] + ((TT))*(ray->dir)[1], (pos)[2] = (1)*(from )[2] + ((TT))*(ray->dir)[2]), \ |
| 239 | v(pos[0], pos[1], pos[2], obj->A, obj->B)) |
| 240 | |
| 241 | #define VALGRAD(VV, DV, TT)((pos)[0] = (1)*(from)[0] + ((TT))*(ray->dir)[0], (pos)[1] = (1)*(from)[1] + ((TT))*(ray->dir)[1], (pos)[2] = (1)*(from )[2] + ((TT))*(ray->dir)[2]); (VV) = vg(grad, pos[0], pos[ 1], pos[2], obj->A, obj->B); (DV) = ((grad)[0]*(ray-> dir)[0] + (grad)[1]*(ray->dir)[1] + (grad)[2]*(ray->dir )[2]) \ |
| 242 | ELL_3V_SCALE_ADD2(pos, 1, from, (TT), ray->dir)((pos)[0] = (1)*(from)[0] + ((TT))*(ray->dir)[0], (pos)[1] = (1)*(from)[1] + ((TT))*(ray->dir)[1], (pos)[2] = (1)*(from )[2] + ((TT))*(ray->dir)[2]); \ |
| 243 | (VV) = vg(grad, pos[0], pos[1], pos[2], obj->A, obj->B); \ |
| 244 | (DV) = ELL_3V_DOT(grad, ray->dir)((grad)[0]*(ray->dir)[0] + (grad)[1]*(ray->dir)[1] + (grad )[2]*(ray->dir)[2]) |
| 245 | |
| 246 | #define LVALGRAD(VV, DV, TT)((pos)[0] = (1)*(from)[0] + ((TT))*(ray->dir)[0], (pos)[1] = (1)*(from)[1] + ((TT))*(ray->dir)[1], (pos)[2] = (1)*(from )[2] + ((TT))*(ray->dir)[2]); (VV) = lvg(grad, pos[0], pos [1], pos[2], obj->A, obj->B); (DV) = ((grad)[0]*(ray-> dir)[0] + (grad)[1]*(ray->dir)[1] + (grad)[2]*(ray->dir )[2]) \ |
| 247 | ELL_3V_SCALE_ADD2(pos, 1, from, (TT), ray->dir)((pos)[0] = (1)*(from)[0] + ((TT))*(ray->dir)[0], (pos)[1] = (1)*(from)[1] + ((TT))*(ray->dir)[1], (pos)[2] = (1)*(from )[2] + ((TT))*(ray->dir)[2]); \ |
| 248 | (VV) = lvg(grad, pos[0], pos[1], pos[2], obj->A, obj->B); \ |
| 249 | (DV) = ELL_3V_DOT(grad, ray->dir)((grad)[0]*(ray->dir)[0] + (grad)[1]*(ray->dir)[1] + (grad )[2]*(ray->dir)[2]) |
| 250 | |
| 251 | #define RR0.61803399 0.61803399 |
| 252 | #define CC(1.0-0.61803399) (1.0-RR0.61803399) |
| 253 | #define SHIFT3(a,b,c,d)(a)=(b); (b)=(c); (c)=(d) (a)=(b); (b)=(c); (c)=(d) |
| 254 | #define SHIFT2(a,b,c)(a)=(b); (b)=(c) (a)=(b); (b)=(c) |
| 255 | |
| 256 | /* testing */ |
| 257 | ELL_3V_SCALE_ADD2(from, 1, ray->from, Tmin, ray->dir)((from)[0] = (1)*(ray->from)[0] + (Tmin)*(ray->dir)[0], (from)[1] = (1)*(ray->from)[1] + (Tmin)*(ray->dir)[1], (from)[2] = (1)*(ray->from)[2] + (Tmin)*(ray->dir)[2]); |
| 258 | saveTmin = Tmin; |
| 259 | Tmin = 0; |
| 260 | Tmax -= saveTmin; |
| 261 | /* |
| 262 | ELL_3V_COPY(from, ray->from); |
| 263 | saveTmin = 0; |
| 264 | */ |
| 265 | |
| 266 | /* evaluate value and derivatives at Tmin and Tmax */ |
| 267 | VALGRAD(Vmin, dVmin, Tmin)((pos)[0] = (1)*(from)[0] + ((Tmin))*(ray->dir)[0], (pos)[ 1] = (1)*(from)[1] + ((Tmin))*(ray->dir)[1], (pos)[2] = (1 )*(from)[2] + ((Tmin))*(ray->dir)[2]); (Vmin) = vg(grad, pos [0], pos[1], pos[2], obj->A, obj->B); (dVmin) = ((grad) [0]*(ray->dir)[0] + (grad)[1]*(ray->dir)[1] + (grad)[2] *(ray->dir)[2]); |
| 268 | VALGRAD(Vmax, dVmax, Tmax)((pos)[0] = (1)*(from)[0] + ((Tmax))*(ray->dir)[0], (pos)[ 1] = (1)*(from)[1] + ((Tmax))*(ray->dir)[1], (pos)[2] = (1 )*(from)[2] + ((Tmax))*(ray->dir)[2]); (Vmax) = vg(grad, pos [0], pos[1], pos[2], obj->A, obj->B); (dVmax) = ((grad) [0]*(ray->dir)[0] + (grad)[1]*(ray->dir)[1] + (grad)[2] *(ray->dir)[2]); |
| 269 | |
| 270 | /* if the segment start and end are both positive or both negative, |
| 271 | and if the derivatives also don't change sign, there's no root. |
| 272 | Also, due to convexity, if values at start and end are negative, |
| 273 | then there is no root */ |
| 274 | if ( (Vmin*Vmax >= 0 && dVmin*dVmax >= 0) |
| 275 | || (Vmin <= 0 && Vmax <= 0) ) { |
| 276 | return AIR_FALSE0; |
| 277 | } |
| 278 | if (tstate->verbose) { |
| 279 | fprintf(stderr__stderrp, "%s%s: dVmin = %g, dVmax = %g, Vmin = %g, Vmax = %g\n", |
| 280 | _echoDot(tstate->depth), me, dVmin, dVmax, Vmin, Vmax); |
| 281 | } |
| 282 | |
| 283 | /* either the value changed sign, or the derivative changed sign, or |
| 284 | both. If, as is common, the derivatives changed sign, but the |
| 285 | values didn't (which means they are both positive, due to a test |
| 286 | above), we need to limit the interval by minimizing the value |
| 287 | until either we see a negative value, or until the minimization |
| 288 | converged. Based on Golden Section Search, NR pg.401 */ |
| 289 | if (dVmin*dVmax < 0 && Vmin*Vmax >= 0) { |
| 290 | t0 = Tmin; |
| 291 | t1 = RR0.61803399*Tmin + CC(1.0-0.61803399)*Tmax; |
| 292 | t2 = CC(1.0-0.61803399)*Tmin + RR0.61803399*Tmax; |
| 293 | t3 = Tmax; |
| 294 | v1 = VAL(t1)(((pos)[0] = (1)*(from)[0] + ((t1))*(ray->dir)[0], (pos)[1 ] = (1)*(from)[1] + ((t1))*(ray->dir)[1], (pos)[2] = (1)*( from)[2] + ((t1))*(ray->dir)[2]), v(pos[0], pos[1], pos[2] , obj->A, obj->B)); |
| 295 | v2 = VAL(t2)(((pos)[0] = (1)*(from)[0] + ((t2))*(ray->dir)[0], (pos)[1 ] = (1)*(from)[1] + ((t2))*(ray->dir)[1], (pos)[2] = (1)*( from)[2] + ((t2))*(ray->dir)[2]), v(pos[0], pos[1], pos[2] , obj->A, obj->B)); |
| 296 | if (tstate->verbose) { |
| 297 | fprintf(stderr__stderrp, "%s%s: \n" |
| 298 | " t0 = % 31.15f\n" |
| 299 | " t1 = % 31.15f -> v1 = % 31.15f\n" |
| 300 | " t2 = % 31.15f -> v2 = % 31.15f\n" |
| 301 | " t3 = % 31.15f\n", |
| 302 | _echoDot(tstate->depth), me, |
| 303 | t0, t1, v1, t2, v2, t3); |
| 304 | } |
| 305 | tol = sqrt(ECHO_POS_EPS2.2204460492503131e-16); |
| 306 | while ( (t3-t0 > tol*(t1+t2)) /* still haven't converged */ |
| 307 | && (v1 > 0 && v2 > 0) ) { /* v1 and v2 both positive */ |
| 308 | diff = v2 - v1; |
| 309 | if (v1 < v2) { |
| 310 | SHIFT3(t3, t2, t1, CC*t0 + RR*t2)(t3)=(t2); (t2)=(t1); (t1)=((1.0-0.61803399)*t0 + 0.61803399* t2); |
| 311 | SHIFT2(v2, v1, VAL(t1))(v2)=(v1); (v1)=((((pos)[0] = (1)*(from)[0] + ((t1))*(ray-> dir)[0], (pos)[1] = (1)*(from)[1] + ((t1))*(ray->dir)[1], ( pos)[2] = (1)*(from)[2] + ((t1))*(ray->dir)[2]), v(pos[0], pos[1], pos[2], obj->A, obj->B))); |
| 312 | } else { |
| 313 | SHIFT3(t0, t1, t2, RR*t1 + CC*t3)(t0)=(t1); (t1)=(t2); (t2)=(0.61803399*t1 + (1.0-0.61803399)* t3); |
| 314 | SHIFT2(v1, v2, VAL(t2))(v1)=(v2); (v2)=((((pos)[0] = (1)*(from)[0] + ((t2))*(ray-> dir)[0], (pos)[1] = (1)*(from)[1] + ((t2))*(ray->dir)[1], ( pos)[2] = (1)*(from)[2] + ((t2))*(ray->dir)[2]), v(pos[0], pos[1], pos[2], obj->A, obj->B))); |
| 315 | } |
| 316 | if (tstate->verbose) { |
| 317 | fprintf(stderr__stderrp, "%s%s: %s ---> \n" |
| 318 | " t0 = % 31.15f\n" |
| 319 | " t1 = % 31.15f -> v1 = % 31.15f\n" |
| 320 | " t2 = % 31.15f -> v2 = % 31.15f\n" |
| 321 | " t3 = % 31.15f\n", |
| 322 | _echoDot(tstate->depth), me, |
| 323 | diff > 0 ? "v1 < v2" : "v1 > v2", |
| 324 | t0, t1, v1, t2, v2, t3); |
| 325 | } |
| 326 | } |
| 327 | if (v1 > 0 && v2 > 0) { |
| 328 | /* the minimization stopped, yet both v1 and v2 are still positive, |
| 329 | so there's no way we can have a root */ |
| 330 | if (tstate->verbose) { |
| 331 | fprintf(stderr__stderrp, "%s%s: minimization found no root\n", |
| 332 | _echoDot(tstate->depth), me); |
| 333 | } |
| 334 | return AIR_FALSE0; |
| 335 | } |
| 336 | /* else either v1 or v2 <= 0, so there is a root (actually two). |
| 337 | By construction, f(t0) is positive, so we can now bracket the |
| 338 | root between t0 and t1 or t2, whichever one is associated with |
| 339 | a smaller value */ |
| 340 | Tmin = t0; |
| 341 | /* HEY: shouldn't I just be using whichever one is closer? */ |
| 342 | Tmax = v1 < v2 ? t1 : t2; |
| 343 | Vmin = VAL(Tmin)(((pos)[0] = (1)*(from)[0] + ((Tmin))*(ray->dir)[0], (pos) [1] = (1)*(from)[1] + ((Tmin))*(ray->dir)[1], (pos)[2] = ( 1)*(from)[2] + ((Tmin))*(ray->dir)[2]), v(pos[0], pos[1], pos [2], obj->A, obj->B)); |
| 344 | Vmax = VAL(Tmax)(((pos)[0] = (1)*(from)[0] + ((Tmax))*(ray->dir)[0], (pos) [1] = (1)*(from)[1] + ((Tmax))*(ray->dir)[1], (pos)[2] = ( 1)*(from)[2] + ((Tmax))*(ray->dir)[2]), v(pos[0], pos[1], pos [2], obj->A, obj->B)); |
| 345 | } |
| 346 | |
| 347 | /* the value isn't necessarily monotonic between Tmin and Tmax, but |
| 348 | we know that there is only one root. Find it with newton-raphson, |
| 349 | using the log of function, both for values and for derivatives */ |
| 350 | iter = 0; |
| 351 | TT = (Tmin + Tmax)/2; |
| 352 | LVALGRAD(VV, dV, TT)((pos)[0] = (1)*(from)[0] + ((TT))*(ray->dir)[0], (pos)[1] = (1)*(from)[1] + ((TT))*(ray->dir)[1], (pos)[2] = (1)*(from )[2] + ((TT))*(ray->dir)[2]); (VV) = lvg(grad, pos[0], pos [1], pos[2], obj->A, obj->B); (dV) = ((grad)[0]*(ray-> dir)[0] + (grad)[1]*(ray->dir)[1] + (grad)[2]*(ray->dir )[2]); |
| 353 | |
| 354 | while (iter < parm->sqNRI && AIR_ABS(VV)((VV) > 0.0f ? (VV) : -(VV)) > parm->sqTol |
| 355 | && AIR_EXISTS(VV)(((int)(!((VV) - (VV))))) && AIR_EXISTS(dV)(((int)(!((dV) - (dV)))))) { |
| 356 | if (tstate->verbose) { |
| 357 | fprintf(stderr__stderrp, "%s%s: iter = %d: TT = %g, VV = %g, dV = %g\n", |
| 358 | _echoDot(tstate->depth), me, iter, TT, VV, dV); |
| 359 | } |
| 360 | TT -= VV/dV; |
| 361 | if (!AIR_IN_OP(Tmin, TT, Tmax)((Tmin) < (TT) && (TT) < (Tmax))) { |
| 362 | /* newton-raphson sent us flying out of bounds; find a tighter |
| 363 | [Tmin,Tmax] bracket with bisection and try again */ |
| 364 | TT = (Tmin + Tmax)/2; |
| 365 | VV = VAL(TT)(((pos)[0] = (1)*(from)[0] + ((TT))*(ray->dir)[0], (pos)[1 ] = (1)*(from)[1] + ((TT))*(ray->dir)[1], (pos)[2] = (1)*( from)[2] + ((TT))*(ray->dir)[2]), v(pos[0], pos[1], pos[2] , obj->A, obj->B)); |
| 366 | if (Vmin*VV < 0) { |
| 367 | Tmax = TT; |
| 368 | Vmax = VV; |
Value stored to 'Vmax' is never read | |
| 369 | } else { |
| 370 | Tmin = TT; |
| 371 | Vmin = VV; |
| 372 | } |
| 373 | TT = (Tmin + Tmax)/2; |
| 374 | } |
| 375 | LVALGRAD(VV, dV, TT)((pos)[0] = (1)*(from)[0] + ((TT))*(ray->dir)[0], (pos)[1] = (1)*(from)[1] + ((TT))*(ray->dir)[1], (pos)[2] = (1)*(from )[2] + ((TT))*(ray->dir)[2]); (VV) = lvg(grad, pos[0], pos [1], pos[2], obj->A, obj->B); (dV) = ((grad)[0]*(ray-> dir)[0] + (grad)[1]*(ray->dir)[1] + (grad)[2]*(ray->dir )[2]); |
| 376 | iter++; |
| 377 | } |
| 378 | |
| 379 | if (!( AIR_EXISTS(VV)(((int)(!((VV) - (VV))))) && AIR_EXISTS(dV)(((int)(!((dV) - (dV))))) )) { |
| 380 | /* we bailed out of the loop above because |
| 381 | we got screwed by numerical errors |
| 382 | --> pretend that there was no intersection, |
| 383 | and HEY this will have to be debugged later */ |
| 384 | if (tstate->verbose) { |
| 385 | fprintf(stderr__stderrp, "%s%s: SORRY, numerical problems!\n", |
| 386 | _echoDot(tstate->depth), me); |
| 387 | } |
| 388 | return AIR_FALSE0; |
| 389 | } |
| 390 | |
| 391 | /* else we succeeded in finding the intersection */ |
| 392 | intx->t = TT + saveTmin; |
| 393 | VALGRAD(VV, dV, TT)((pos)[0] = (1)*(from)[0] + ((TT))*(ray->dir)[0], (pos)[1] = (1)*(from)[1] + ((TT))*(ray->dir)[1], (pos)[2] = (1)*(from )[2] + ((TT))*(ray->dir)[2]); (VV) = vg(grad, pos[0], pos[ 1], pos[2], obj->A, obj->B); (dV) = ((grad)[0]*(ray-> dir)[0] + (grad)[1]*(ray->dir)[1] + (grad)[2]*(ray->dir )[2]); /* puts gradient into grad */ |
| 394 | ELL_3V_NORM(intx->norm, grad, tmp)(tmp = (sqrt((((grad))[0]*((grad))[0] + ((grad))[1]*((grad))[ 1] + ((grad))[2]*((grad))[2]))), ((intx->norm)[0] = (1.0/tmp )*(grad)[0], (intx->norm)[1] = (1.0/tmp)*(grad)[1], (intx-> norm)[2] = (1.0/tmp)*(grad)[2])); |
| 395 | intx->obj = OBJECT(obj)((echoObject*)obj); |
| 396 | /* set in intx: |
| 397 | yes: t, norm |
| 398 | no: u, v |
| 399 | */ |
| 400 | return AIR_TRUE1; |
| 401 | } |