Menu
▾
▴
[Teem-commit] teem/src/push action.c, 1.62, 1.63 binning.c, 1.16, 1.17 corePush.c, 1.22, 1.23 forces.c, 1.19, 1.20 methodsPush.c, 1.23, 1.24 push.h, 1.41, 1.42 setup.c, 1.13, 1.14
|
From: Gordon K. <kin...@us...> - 2006-07-02 23:30:37
|
Update of /cvsroot/teem/teem/src/push In directory sc8-pr-cvs10.sourceforge.net:/tmp/cvs-serv21437 Modified Files: action.c binning.c corePush.c forces.c methodsPush.c push.h setup.c Log Message: less broken than before Index: push.h =================================================================== RCS file: /cvsroot/teem/teem/src/push/push.h,v retrieving revision 1.41 retrieving revision 1.42 diff -C2 -d -r1.41 -r1.42 *** push.h 21 Apr 2006 17:52:29 -0000 1.41 --- push.h 2 Jul 2006 23:30:30 -0000 1.42 *************** *** 165,171 **** */ typedef struct { char name[AIR_STRLEN_SMALL]; double (*func)(double dist, const double parm[PUSH_FORCE_PARM_MAXNUM]); ! double (*maxDist)(const double parm[PUSH_FORCE_PARM_MAXNUM]); double parm[PUSH_FORCE_PARM_MAXNUM]; } pushForce; --- 165,172 ---- */ typedef struct { + int noop; char name[AIR_STRLEN_SMALL]; double (*func)(double dist, const double parm[PUSH_FORCE_PARM_MAXNUM]); ! double (*extent)(const double parm[PUSH_FORCE_PARM_MAXNUM]); double parm[PUSH_FORCE_PARM_MAXNUM]; } pushForce; *************** *** 193,196 **** --- 194,198 ---- double drag, /* to slow fast things down */ preDrag, /* different drag pre-min-iter */ + velWarp, /* velocity warping stability hack */ step, /* time step in integration */ mass, /* mass of particles */ *************** *** 203,207 **** tlThresh, tlSoft, tlStep, /* tractlet formation parameters */ minMeanVel; /* stop if mean velocity drops below this */ ! int tlFrenet, /* use Frenet frames for tractlet forces */ singleBin, /* disable binning (for debugging) */ driftCorrect, /* prevent sliding near anisotropy edges */ --- 205,210 ---- tlThresh, tlSoft, tlStep, /* tractlet formation parameters */ minMeanVel; /* stop if mean velocity drops below this */ ! int tlUse, /* enable tractlets */ ! tlFrenet, /* use Frenet frames for tractlet forces */ singleBin, /* disable binning (for debugging) */ driftCorrect, /* prevent sliding near anisotropy edges */ *************** *** 234,238 **** *ninv, /* pre-computed inverse of nten */ *nmask; /* mask image from nten */ ! gageContext *gctx; /* gage context around nten and nmask */ gagePerVolume *tpvl, *ipvl; /* gage pervolumes around nten and ninv */ tenFiberContext *fctx; /* tenFiber context around nten */ --- 237,241 ---- *ninv, /* pre-computed inverse of nten */ *nmask; /* mask image from nten */ ! gageContext *gctx; /* gage context around nten, ninv, nmask */ gagePerVolume *tpvl, *ipvl; /* gage pervolumes around nten and ninv */ tenFiberContext *fctx; /* tenFiber context around nten */ *************** *** 260,264 **** *stageBarrierB; /* OUTPUT ---------------------------- */ ! double time; /* how long it took to run */ unsigned int iter; /* how many iterations were needed */ Nrrd *noutPos, /* list of 2D or 3D positions */ --- 263,267 ---- *stageBarrierB; /* OUTPUT ---------------------------- */ ! double time; /* total time spent in computation */ unsigned int iter; /* how many iterations were needed */ Nrrd *noutPos, /* list of 2D or 3D positions */ *************** *** 283,286 **** --- 286,290 ---- PUSH_EXPORT int pushStart(pushContext *pctx); PUSH_EXPORT int pushIterate(pushContext *pctx); + PUSH_EXPORT int pushFreeze(pushContext *pctx); PUSH_EXPORT int pushRun(pushContext *pctx); PUSH_EXPORT int pushFinish(pushContext *pctx); Index: setup.c =================================================================== RCS file: /cvsroot/teem/teem/src/push/setup.c,v retrieving revision 1.13 retrieving revision 1.14 diff -C2 -d -r1.13 -r1.14 *** setup.c 21 Apr 2006 17:52:29 -0000 1.13 --- setup.c 2 Jul 2006 23:30:30 -0000 1.14 *************** *** 84,87 **** --- 84,91 ---- TEN_T_COPY(ten, _ten); TEN_T_INV(inv, ten, det); + if (!det || !AIR_EXISTS(det)) { + fprintf(stderr, "!%s: tensor %u/%u has determinant %g\n", me, + AIR_CAST(unsigned int, ii), AIR_CAST(unsigned int, NN), det); + } TEN_T_COPY(_inv, inv); _ten += 7; *************** *** 333,341 **** /* ------------------------ find maxEval, maxDet, and set up binning */ nn = nrrdElementNumber(pctx->nten)/7; - tdata = (float*)pctx->nten->data; pctx->maxEval = 0; pctx->maxDet = 0; pctx->meanEval = 0; count = 0; for (ii=0; ii<nn; ii++) { tenEigensolve_f(eval, NULL, tdata); --- 337,345 ---- /* ------------------------ find maxEval, maxDet, and set up binning */ nn = nrrdElementNumber(pctx->nten)/7; pctx->maxEval = 0; pctx->maxDet = 0; pctx->meanEval = 0; count = 0; + tdata = (float*)pctx->nten->data; for (ii=0; ii<nn; ii++) { tenEigensolve_f(eval, NULL, tdata); *************** *** 345,349 **** pctx->meanEval += eval[0]; pctx->maxEval = AIR_MAX(pctx->maxEval, eval[0]); ! pctx->maxDet = AIR_MAX(pctx->maxDet, eval[0] + eval[1] + eval[2]); } tdata += 7; --- 349,353 ---- pctx->meanEval += eval[0]; pctx->maxEval = AIR_MAX(pctx->maxEval, eval[0]); ! pctx->maxDet = AIR_MAX(pctx->maxDet, eval[0]*eval[1]*eval[2]); } tdata += 7; *************** *** 351,355 **** pctx->meanEval /= count; pctx->maxDist = (2*pctx->scale*pctx->maxEval ! *pctx->force->maxDist(pctx->force->parm)); if (pctx->singleBin) { pctx->binsEdge[0] = 1; --- 355,359 ---- pctx->meanEval /= count; pctx->maxDist = (2*pctx->scale*pctx->maxEval ! *pctx->force->extent(pctx->force->parm)); if (pctx->singleBin) { pctx->binsEdge[0] = 1; *************** *** 361,379 **** ELL_4MV_COL1_GET(col[1], pctx->gctx->shape->ItoW); col[1][3] = 0.0; ELL_4MV_COL2_GET(col[2], pctx->gctx->shape->ItoW); col[2][3] = 0.0; ! volEdge[0] = ELL_3V_LEN(col[0])*(pctx->gctx->shape->size[0]-1); ! volEdge[1] = ELL_3V_LEN(col[1])*(pctx->gctx->shape->size[1]-1); ! volEdge[2] = ELL_3V_LEN(col[2])*(pctx->gctx->shape->size[2]-1); fprintf(stderr, "!%s: volEdge = %g %g %g\n", me, volEdge[0], volEdge[1], volEdge[2]); ! pctx->binsEdge[0] = (int)floor(volEdge[0]/pctx->maxDist); pctx->binsEdge[0] = pctx->binsEdge[0] ? pctx->binsEdge[0] : 1; ! pctx->binsEdge[1] = (int)floor(volEdge[1]/pctx->maxDist); pctx->binsEdge[1] = pctx->binsEdge[1] ? pctx->binsEdge[1] : 1; ! pctx->binsEdge[2] = (int)floor(volEdge[2]/pctx->maxDist); pctx->binsEdge[2] = pctx->binsEdge[2] ? pctx->binsEdge[2] : 1; if (2 == pctx->dimIn) { pctx->binsEdge[pctx->sliceAxis] = 1; } ! fprintf(stderr, "!%s: maxEval=%g -> maxDist=%g -> binsEdge=%u,%u,%u\n", me, pctx->maxEval, pctx->maxDist, pctx->binsEdge[0], pctx->binsEdge[1], pctx->binsEdge[2]); --- 365,386 ---- ELL_4MV_COL1_GET(col[1], pctx->gctx->shape->ItoW); col[1][3] = 0.0; ELL_4MV_COL2_GET(col[2], pctx->gctx->shape->ItoW); col[2][3] = 0.0; ! volEdge[0] = ELL_3V_LEN(col[0])*pctx->gctx->shape->size[0]; ! volEdge[1] = ELL_3V_LEN(col[1])*pctx->gctx->shape->size[1]; ! volEdge[2] = ELL_3V_LEN(col[2])*pctx->gctx->shape->size[2]; fprintf(stderr, "!%s: volEdge = %g %g %g\n", me, volEdge[0], volEdge[1], volEdge[2]); ! pctx->binsEdge[0] = AIR_CAST(unsigned int, ! floor(volEdge[0]/pctx->maxDist)); pctx->binsEdge[0] = pctx->binsEdge[0] ? pctx->binsEdge[0] : 1; ! pctx->binsEdge[1] = AIR_CAST(unsigned int, ! floor(volEdge[1]/pctx->maxDist)); pctx->binsEdge[1] = pctx->binsEdge[1] ? pctx->binsEdge[1] : 1; ! pctx->binsEdge[2] = AIR_CAST(unsigned int, ! floor(volEdge[2]/pctx->maxDist)); pctx->binsEdge[2] = pctx->binsEdge[2] ? pctx->binsEdge[2] : 1; if (2 == pctx->dimIn) { pctx->binsEdge[pctx->sliceAxis] = 1; } ! fprintf(stderr, "!%s: maxEval=%g -> maxDist=%g -> binsEdge=(%u,%u,%u)\n", me, pctx->maxEval, pctx->maxDist, pctx->binsEdge[0], pctx->binsEdge[1], pctx->binsEdge[2]); *************** *** 459,463 **** _pushThingSetup(pushContext *pctx) { char me[]="_pushThingSetup", err[BIFF_STRLEN]; ! double (*lup)(const void *v, size_t I); unsigned int *stn, pointIdx, baseIdx, thingIdx; pushThing *thing; --- 466,470 ---- _pushThingSetup(pushContext *pctx) { char me[]="_pushThingSetup", err[BIFF_STRLEN]; ! double (*lup)(const void *v, size_t I), maxDet; unsigned int *stn, pointIdx, baseIdx, thingIdx; pushThing *thing; *************** *** 470,473 **** --- 477,486 ---- lup = pctx->npos ? nrrdDLookup[pctx->npos->type] : NULL; stn = pctx->nstn ? (unsigned int*)pctx->nstn->data : NULL; + fprintf(stderr, "!%s: initilizing/seeding ... \n", me); + /* HEY: we end up keeping a local copy of maxDet because convolution + can produce a tensor with higher determinant than that of any + original sample. However, if this is going into effect, + detReject should probably *not* be enabled... */ + maxDet = pctx->maxDet; for (thingIdx=0; thingIdx<pctx->thingNum; thingIdx++) { double detProbe; *************** *** 532,536 **** */ _pushProbe(pctx->task[0], thing->vert + 0); ! detProbe = TEN_T_TRACE(thing->vert[0].ten); /* assuming that we're not using some very blurring kernel, this will eventually succeed, because we previously checked --- 545,550 ---- */ _pushProbe(pctx->task[0], thing->vert + 0); ! detProbe = TEN_T_DET(thing->vert[0].ten); ! maxDet = AIR_MAX(maxDet, detProbe); /* assuming that we're not using some very blurring kernel, this will eventually succeed, because we previously checked *************** *** 548,552 **** *pctx->seedThreshSign > 0) ) ! || (pctx->detReject && (airDrandMT() < detProbe/pctx->maxDet)) ); } --- 562,566 ---- *pctx->seedThreshSign > 0) ) ! || (pctx->detReject && (airDrandMT() < detProbe/maxDet)) ); } *************** *** 565,568 **** --- 579,583 ---- } } + fprintf(stderr, "!%s: ... seeding DONE\n", me); /* { Index: forces.c =================================================================== RCS file: /cvsroot/teem/teem/src/push/forces.c,v retrieving revision 1.19 retrieving revision 1.20 diff -C2 -d -r1.19 -r1.20 *** forces.c 21 Apr 2006 17:52:29 -0000 1.19 --- forces.c 2 Jul 2006 23:30:30 -0000 1.20 *************** *** 82,87 **** double ! _pushForceUnknownMaxDist(const double *parm) { ! char me[]="_pushForceUnknownMaxDist"; AIR_UNUSED(parm); --- 82,87 ---- double ! _pushForceUnknownExtent(const double *parm) { ! char me[]="_pushForceUnknownExtent"; AIR_UNUSED(parm); *************** *** 98,102 **** double _pushForceSpringFunc(double dist, const double *parm) { ! char me[]="_pushForceSpringFunc"; double xx, ret, pull; --- 98,102 ---- double _pushForceSpringFunc(double dist, const double *parm) { ! /* char me[]="_pushForceSpringFunc"; */ double xx, ret, pull; *************** *** 120,124 **** double ! _pushForceSpringMaxDist(const double *parm) { return 1.0 + parm[0]; --- 120,124 ---- double ! _pushForceSpringExtent(const double *parm) { return 1.0 + parm[0]; *************** *** 147,151 **** double ! _pushForceGaussMaxDist(const double *parm) { return (1.0/SQRTTHREE)*parm[0]; --- 147,151 ---- double ! _pushForceGaussExtent(const double *parm) { return (1.0/SQRTTHREE)*parm[0]; *************** *** 166,170 **** double ! _pushForceChargeMaxDist(const double *parm) { return parm[0]; --- 166,170 ---- double ! _pushForceChargeExtent(const double *parm) { return parm[0]; *************** *** 186,190 **** double ! _pushForceCotanMaxDist(const double *parm) { AIR_UNUSED(parm); --- 186,190 ---- double ! _pushForceCotanExtent(const double *parm) { AIR_UNUSED(parm); *************** *** 206,210 **** double ! _pushForceNoneMaxDist(const double *parm) { AIR_UNUSED(parm); --- 206,210 ---- double ! _pushForceNoneExtent(const double *parm) { AIR_UNUSED(parm); *************** *** 238,248 **** double ! (*_pushForceMaxDist[PUSH_FORCE_MAX+1])(const double *parm) = { ! _pushForceUnknownMaxDist, ! _pushForceSpringMaxDist, ! _pushForceGaussMaxDist, ! _pushForceChargeMaxDist, ! _pushForceCotanMaxDist, ! _pushForceNoneMaxDist, }; --- 238,248 ---- double ! (*_pushForceExtent[PUSH_FORCE_MAX+1])(const double *parm) = { ! _pushForceUnknownExtent, ! _pushForceSpringExtent, ! _pushForceGaussExtent, ! _pushForceChargeExtent, ! _pushForceCotanExtent, ! _pushForceNoneExtent, }; *************** *** 254,259 **** force = (pushForce *)calloc(1, sizeof(pushForce)); if (force) { force->func = NULL; ! force->maxDist = NULL; for (pi=0; pi<PUSH_FORCE_PARM_MAXNUM; pi++) { force->parm[pi] = AIR_NAN; --- 254,260 ---- force = (pushForce *)calloc(1, sizeof(pushForce)); if (force) { + force->noop = AIR_FALSE; /* by default, not a no-op */ force->func = NULL; ! force->extent = NULL; for (pi=0; pi<PUSH_FORCE_PARM_MAXNUM; pi++) { force->parm[pi] = AIR_NAN; *************** *** 292,298 **** if (!strcmp(airEnumStr(pushForceEnum, pushForceNone), str)) { /* special case: no parameters */ strcpy(force->name, _pushForceStr[pushForceNone]); force->func = _pushForceFunc[pushForceNone]; ! force->maxDist = _pushForceMaxDist[pushForceNone]; airMopOkay(mop); return force; --- 293,300 ---- if (!strcmp(airEnumStr(pushForceEnum, pushForceNone), str)) { /* special case: no parameters */ + force->noop = AIR_TRUE; strcpy(force->name, _pushForceStr[pushForceNone]); force->func = _pushForceFunc[pushForceNone]; ! force->extent = _pushForceExtent[pushForceNone]; airMopOkay(mop); return force; *************** *** 303,307 **** strcpy(force->name, _pushForceStr[pushForceCotan]); force->func = _pushForceFunc[pushForceCotan]; ! force->maxDist = _pushForceMaxDist[pushForceCotan]; airMopOkay(mop); return force; --- 305,309 ---- strcpy(force->name, _pushForceStr[pushForceCotan]); force->func = _pushForceFunc[pushForceCotan]; ! force->extent = _pushForceExtent[pushForceCotan]; airMopOkay(mop); return force; *************** *** 358,362 **** strcpy(force->name, _pushForceStr[fri]); force->func = _pushForceFunc[fri]; ! force->maxDist = _pushForceMaxDist[fri]; airMopOkay(mop); --- 360,364 ---- strcpy(force->name, _pushForceStr[fri]); force->func = _pushForceFunc[fri]; ! force->extent = _pushForceExtent[fri]; airMopOkay(mop); Index: binning.c =================================================================== RCS file: /cvsroot/teem/teem/src/push/binning.c,v retrieving revision 1.16 retrieving revision 1.17 diff -C2 -d -r1.16 -r1.17 *** binning.c 28 Mar 2006 11:48:58 -0000 1.16 --- binning.c 2 Jul 2006 23:30:30 -0000 1.17 *************** *** 77,81 **** ELL_34V_HOMOG(posIdx, posIdx); for (axi=0; axi<3; axi++) { ! eidx[axi] = airIndexClamp(0, posIdx[axi], pctx->gctx->shape->size[axi], pctx->binsEdge[axi]); } --- 77,83 ---- ELL_34V_HOMOG(posIdx, posIdx); for (axi=0; axi<3; axi++) { ! eidx[axi] = airIndexClamp(-0.5, ! posIdx[axi], ! pctx->gctx->shape->size[axi]-0.5, pctx->binsEdge[axi]); } *************** *** 183,187 **** void pushBinAllNeighborSet(pushContext *pctx) { ! char me[]="pushBinAllNeighborSet"; pushBin *nei[3*3*3]; unsigned int neiNum, xi, yi, zi, xx, yy, zz, xmax, ymax, zmax, binIdx; --- 185,189 ---- void pushBinAllNeighborSet(pushContext *pctx) { ! /* char me[]="pushBinAllNeighborSet"; */ pushBin *nei[3*3*3]; unsigned int neiNum, xi, yi, zi, xx, yy, zz, xmax, ymax, zmax, binIdx; Index: corePush.c =================================================================== RCS file: /cvsroot/teem/teem/src/push/corePush.c,v retrieving revision 1.22 retrieving revision 1.23 diff -C2 -d -r1.22 -r1.23 *** corePush.c 21 Apr 2006 17:52:29 -0000 1.22 --- corePush.c 2 Jul 2006 23:30:30 -0000 1.23 *************** *** 305,308 **** --- 305,309 ---- char me[]="pushIterate", *_err, err[BIFF_STRLEN]; unsigned int ti, thingNum; + double time0, time1; if (!pctx) { *************** *** 314,317 **** --- 315,321 ---- fprintf(stderr, "%s: starting iteration\n", me); } + + time0 = airTime(); + /* the _pushWorker checks finished after the barriers */ pctx->finished = AIR_FALSE; *************** *** 362,369 **** --- 366,408 ---- _pushForceSample(pctx, 300, 300); } + + time1 = airTime(); + pctx->time += time1 - time0; return 0; } + /* + ******** pushFreeze + ** + ** zeros out the velocity of all points + */ + int + pushFreeze(pushContext *pctx) { + char me[]="pushFreeze", err[BIFF_STRLEN]; + unsigned int binIdx, thingIdx, pointNum, pointIdx; + pushBin *bin; + pushThing *thing; + pushPoint *point; + + if (!pctx) { + sprintf(err, "%s: got NULL pointer", me); + biffAdd(PUSH, err); return 1; + } + + pointNum = 0; + for (binIdx=0; binIdx<pctx->binNum; binIdx++) { + bin = pctx->bin + binIdx; + for (thingIdx=0; thingIdx<bin->thingNum; thingIdx++) { + thing = bin->thing[thingIdx]; + for (pointIdx=0; pointIdx<thing->vertNum; pointIdx++) { + point = thing->vert + pointIdx; + ELL_3V_SET(point->vel, 0.0, 0.0, 0.0); + } + } + } + return 0; + } + int pushRun(pushContext *pctx) { *************** *** 437,440 **** --- 476,480 ---- || pctx->iter < pctx->maxIter)) ); pctx->time1 = airTime(); + /* this over-writes the value incremented within pushIterate() */ pctx->time = pctx->time1 - pctx->time0; Index: action.c =================================================================== RCS file: /cvsroot/teem/teem/src/push/action.c,v retrieving revision 1.62 retrieving revision 1.63 diff -C2 -d -r1.62 -r1.63 *** action.c 21 Apr 2006 17:52:29 -0000 1.62 --- action.c 2 Jul 2006 23:30:30 -0000 1.63 *************** *** 31,36 **** int _pushProbe(pushTask *task, pushPoint *point) { ! char me[]="_pushProbe"; ! double eval[3], sum, posWorld[4], posIdx[4], det; int inside, ret; --- 31,36 ---- int _pushProbe(pushTask *task, pushPoint *point) { ! /* char me[]="_pushProbe"; */ ! double eval[3], sum, posWorld[4], posIdx[4]; int inside, ret; *************** *** 49,61 **** TEN_T_COPY(point->ten, task->tenAns); ! /* TEN_T_INV(point->inv, point->ten, det); */ TEN_T_COPY(point->inv, task->invAns); ! tenEigensolve_d(eval, NULL, point->ten); ! /* sadly, the fact that tenAnisoCalc_f exists only for floats is part ! of the motivation for hard-wiring the aniso measure to Cl1 */ ! sum = eval[0] + eval[1] + eval[2]; ! point->aniso = (eval[0] - eval[1])/(sum + FLT_EPSILON); ELL_3V_COPY(point->cnt, task->cntAns); if (tenGageUnknown != task->pctx->gravItem) { --- 49,66 ---- TEN_T_COPY(point->ten, task->tenAns); ! /* now we're probing the inverse from the pre-computed ninv TEN_T_INV(point->inv, point->ten, det); */ TEN_T_COPY(point->inv, task->invAns); ! if (task->pctx->tlUse) { ! tenEigensolve_d(eval, NULL, point->ten); ! /* sadly, the fact that tenAnisoCalc_f exists only for floats is part ! of the motivation for hard-wiring the aniso measure to Cl1 */ ! /* HEY: with _tenAnisoEval_f[](), that's no longer true!!! */ ! sum = eval[0] + eval[1] + eval[2]; ! point->aniso = (eval[0] - eval[1])/(sum + FLT_EPSILON); ! } else { ! point->aniso = 0.0; ! } ELL_3V_COPY(point->cnt, task->cntAns); if (tenGageUnknown != task->pctx->gravItem) { *************** *** 76,80 **** } ! int _pushThingTotal(pushContext *pctx) { unsigned int binIdx, thingNum; --- 81,85 ---- } ! unsigned int _pushThingTotal(pushContext *pctx) { unsigned int binIdx, thingNum; *************** *** 87,91 **** } ! int _pushPointTotal(pushContext *pctx) { unsigned int binIdx, thingIdx, pointNum; --- 92,96 ---- } ! unsigned int _pushPointTotal(pushContext *pctx) { unsigned int binIdx, thingIdx, pointNum; *************** *** 171,189 **** _pushPairwiseForce(pushContext *pctx, double fvec[3], pushForce *force, pushPoint *myPoint, pushPoint *herPoint) { ! char me[]="_pushPairwiseForce", err[BIFF_STRLEN]; ! double inv[7], dist, mag, scl, ! U[3], nU[3], lenU, lenUsqr, ! myV[3], mylenV, herV[3], herlenV, V[3], lenV; ! ! /* in case we return early */ ! ELL_3V_SET(fvec, 0, 0, 0); ! if (!strcmp("none", force->name)) { /* no actual force */ return 0; } ELL_3V_SUB(U, herPoint->pos, myPoint->pos); lenUsqr = ELL_3V_DOT(U, U); if (lenUsqr < FLT_EPSILON) { /* myPoint and herPoint are overlapping */ --- 176,195 ---- _pushPairwiseForce(pushContext *pctx, double fvec[3], pushForce *force, pushPoint *myPoint, pushPoint *herPoint) { ! /* char me[]="_pushPairwiseForce", err[BIFF_STRLEN] */; ! double inv[7], dist, mag, ! U[3], lenUsqr, ! V[3], nV[3], W[3], lenV; ! #if 0 ! if (force->noop) { /* no actual force */ + ELL_3V_SET(fvec, 0, 0, 0); return 0; } + #endif ELL_3V_SUB(U, herPoint->pos, myPoint->pos); lenUsqr = ELL_3V_DOT(U, U); + #if 0 if (lenUsqr < FLT_EPSILON) { /* myPoint and herPoint are overlapping */ *************** *** 192,236 **** myPoint->pos[0], myPoint->pos[1], myPoint->pos[2]); */ return 0; } if (lenUsqr >= (pctx->maxDist)*(pctx->maxDist)) { /* too far away to influence each other */ return 0; } - lenU = sqrt(lenUsqr); - ELL_3V_SCALE(nU, 1.0/lenU, U); - TEN_T_SCALE_ADD2(inv, 0.5, myPoint->inv, 0.5, herPoint->inv); - - TEN_TV_MUL(myV, myPoint->inv, U); - mylenV = ELL_3V_LEN(myV); - TEN_TV_MUL(herV, herPoint->inv, U); - herlenV = ELL_3V_LEN(herV); - TEN_TV_MUL(V, inv, U); ! lenV = ELL_3V_LEN(V); ! dist = lenV/pctx->scale; ! scl = lenU/(DBL_EPSILON + lenV); ! if (pctx->driftCorrect) { ! scl *= (herlenV/(DBL_EPSILON + mylenV))*(herlenV/(DBL_EPSILON + mylenV)); ! /* ! scl *= (herlenV - mylenV)/(lenU*mylenV); ! */ ! } ! mag = force->func(dist/2, force->parm)*scl; ! if (!( AIR_EXISTS(mag) && ELL_3V_EXISTS(nU) )) { ! fprintf(stderr, "!%s: mag=%g, nU=(%g,%g,%g), dist=|%g,%g,%g|/%g=%g\n", ! me, mag, ! nU[0], nU[1], nU[2], ! V[0], V[1], V[2], pctx->scale, dist); ELL_3V_SET(fvec, 0, 0, 0); } else { ! ELL_3V_SCALE(fvec, mag, nU); } return 0; --- 198,237 ---- myPoint->pos[0], myPoint->pos[1], myPoint->pos[2]); */ + ELL_3V_SET(fvec, 0, 0, 0); return 0; } + #endif if (lenUsqr >= (pctx->maxDist)*(pctx->maxDist)) { /* too far away to influence each other */ + ELL_3V_SET(fvec, 0, 0, 0); return 0; } TEN_T_SCALE_ADD2(inv, 0.5, myPoint->inv, 0.5, herPoint->inv); TEN_TV_MUL(V, inv, U); ! ELL_3V_NORM(nV, V, lenV); ! dist = lenV/(2*pctx->scale); ! mag = force->func(dist, force->parm); ! TEN_TV_MUL(W, inv, nV); ! ELL_3V_SCALE(fvec, mag, W); ! #if 0 ! if (!( AIR_EXISTS(mag) )) { ! fprintf(stderr, "!%s: myPos = %g %g %g\n", me, ! myPoint->pos[0], myPoint->pos[1], myPoint->pos[2]); ! fprintf(stderr, "!%s: mag=%g, dist=%g=(|V=%g,%g,%g|=%g)/%g\n", ! me, mag, dist, ! V[0], V[1], V[2], lenV, ! pctx->scale); ! fprintf(stderr, "!%s: U=%g,%g,%g\n", me, U[0], U[1], U[2]); ELL_3V_SET(fvec, 0, 0, 0); } else { ! TEN_TV_MUL(W, inv, nV); ! ELL_3V_SCALE(fvec, mag, W); } + #endif return 0; *************** *** 317,321 **** AIR_UNUSED(parm); - /* fprintf(stderr, "!%s: bingo 0 %d\n", me, myBinIdx); */ myBin = task->pctx->bin + myBinIdx; --- 318,321 ---- *************** *** 329,334 **** } - /* fprintf(stderr, "!%s: bingo 1 %d\n", me, myBinIdx); */ - /* compute pair-wise forces between all points in this bin, and all points in all neighboring bins */ --- 329,332 ---- *************** *** 377,381 **** if (task->pctx->bigTrace) { double scl, trc; ! trc = myPoint->ten[1] + myPoint->ten[4] + myPoint->ten[6]; scl = task->pctx->bigTrace*trc/(task->pctx->bigTrace + trc); ELL_3V_SCALE(myPoint->frc, scl, myPoint->frc); --- 375,379 ---- if (task->pctx->bigTrace) { double scl, trc; ! trc = TEN_T_TRACE(myPoint->ten); scl = task->pctx->bigTrace*trc/(task->pctx->bigTrace + trc); ELL_3V_SCALE(myPoint->frc, scl, myPoint->frc); *************** *** 401,407 **** } - /* fprintf(stderr, "!%s: bingo 1.3 %d\n", me, myBinIdx); */ /* each point in this thing also potentially experiences wall forces */ if (task->pctx->wall) { double posWorld[4], posIdx[4], len, frcIdx[4], frcWorld[4]; ELL_3V_COPY(posWorld, myPoint->pos); posWorld[3] = 1.0; --- 399,407 ---- } /* each point in this thing also potentially experiences wall forces */ if (task->pctx->wall) { + /* these wall forces are NOT the same as springs; the scheme + ensures that the forces are C1, not just C0, which seems to + help with stability */ double posWorld[4], posIdx[4], len, frcIdx[4], frcWorld[4]; ELL_3V_COPY(posWorld, myPoint->pos); posWorld[3] = 1.0; *************** *** 412,422 **** frcIdx[ci] = 0; } else { ! len = posIdx[ci] - (-0.5); if (len < 0) { ! frcIdx[ci] = -task->pctx->wall*len; } else { len = posIdx[ci] - (task->pctx->gctx->shape->size[ci] - 0.5); if (len > 0) { ! frcIdx[ci] = -task->pctx->wall*len; } else { frcIdx[ci] = 0; --- 412,423 ---- frcIdx[ci] = 0; } else { ! len = posIdx[ci] - 0.5; if (len < 0) { ! len *= -1; ! frcIdx[ci] = task->pctx->wall*len*len; } else { len = posIdx[ci] - (task->pctx->gctx->shape->size[ci] - 0.5); if (len > 0) { ! frcIdx[ci] = -task->pctx->wall*len*len; } else { frcIdx[ci] = 0; *************** *** 431,458 **** } /* for myPointIdx ... */ ! /* fprintf(stderr, "!%s: bingo 2 %d\n", me, myBinIdx); */ ! ! /* drag and nudging are computed per-thing, not per-point */ ! for (myThingIdx=0; myThingIdx<myBin->thingNum; myThingIdx++) { ! double len; ! myThing = myBin->thing[myThingIdx]; ! if (task->pctx->minIter ! && task->pctx->iter < task->pctx->minIter) { ! drag = AIR_CAST(double, ! AIR_AFFINE(0, task->pctx->iter, task->pctx->minIter, ! task->pctx->preDrag, task->pctx->drag)); ! } else { ! drag = AIR_CAST(double, task->pctx->drag); ! } ! ELL_3V_SCALE_INCR(myThing->point.frc, -drag, myThing->point.vel); ! if (task->pctx->nudge) { ! ELL_3V_NORM(dir, myThing->point.pos, len); ! if (len) { ! ELL_3V_SCALE_INCR(myThing->point.frc, -task->pctx->nudge*len, dir); } } } - /* fprintf(stderr, "!%s: bingo 100\n", me); */ return 0; } --- 432,460 ---- } /* for myPointIdx ... */ ! if (task->pctx->preDrag ! || task->pctx->drag ! || task->pctx->nudge) { ! /* drag and nudging are computed per-thing, not per-point */ ! for (myThingIdx=0; myThingIdx<myBin->thingNum; myThingIdx++) { ! double len; ! myThing = myBin->thing[myThingIdx]; ! if (task->pctx->minIter ! && task->pctx->iter < task->pctx->minIter) { ! drag = AIR_CAST(double, ! AIR_AFFINE(0, task->pctx->iter, task->pctx->minIter, ! task->pctx->preDrag, task->pctx->drag)); ! } else { ! drag = AIR_CAST(double, task->pctx->drag); ! } ! ELL_3V_SCALE_INCR(myThing->point.frc, -drag, myThing->point.vel); ! if (task->pctx->nudge) { ! ELL_3V_NORM(dir, myThing->point.pos, len); ! if (len) { ! ELL_3V_SCALE_INCR(myThing->point.frc, -task->pctx->nudge*len, dir); ! } } } } return 0; } *************** *** 479,483 **** /* the now-single point does have to be binned */ _pushBinPointAdd(task->pctx, oldBin, &(thing->point)); ! thing->point.charge = _pushThingPointCharge(task->pctx, thing); /* free vertex info */ airFree(thing->vert); --- 481,487 ---- /* the now-single point does have to be binned */ _pushBinPointAdd(task->pctx, oldBin, &(thing->point)); ! thing->point.charge = (task->pctx->tlUse ! ? _pushThingPointCharge(task->pctx, thing) ! : 1.0); /* free vertex info */ airFree(thing->vert); *************** *** 701,766 **** /* update dynamics: applies equally to points and tractlets */ ! mass = _pushThingMass(task->pctx, thing); ! if (mass) { ! if (task->pctx->verbose) { ! fprintf(stderr, "vel(%f,%f) * step(%f) -+-> pos(%f,%f)\n", ! thing->point.vel[0], thing->point.vel[0], step, ! thing->point.pos[0], thing->point.pos[1]); ! fprintf(stderr, "frc(%f,%f)*step(%f)/mass(%f) (%f) -+-> vel(%f,%f)\n", ! thing->point.frc[0], thing->point.frc[0], ! step, mass, step/mass, ! thing->point.vel[0], thing->point.vel[1]); ! } ! if (!ELL_3V_EXISTS(thing->point.vel)) { ! ELL_3V_SET(thing->point.vel, 0, 0, 0); ! } ! ELL_3V_SCALE_INCR(thing->point.pos, step, thing->point.vel); ! if (2 == task->pctx->dimIn) { ! double posIdx[4], posWorld[4]; ! ELL_3V_COPY(posWorld, thing->point.pos); posWorld[3] = 1.0; ! ELL_4MV_MUL(posIdx, task->pctx->gctx->shape->WtoI, posWorld); ! ELL_4V_HOMOG(posIdx, posIdx); ! posIdx[task->pctx->sliceAxis] = 0.0; ! ELL_4MV_MUL(posWorld, task->pctx->gctx->shape->ItoW, posIdx); ! ELL_34V_HOMOG(thing->point.pos, posWorld); ! } ! if (!ELL_3V_EXISTS(thing->point.frc)) { ! ELL_3V_SET(thing->point.frc, 0, 0, 0); ! } ! ELL_3V_SCALE_INCR(thing->point.vel, step/mass, thing->point.frc); ! if (task->pctx->verbose) { ! fprintf(stderr, "thing %d: pos(%f,%f); vel(%f,%f)\n", ! thing->ttaagg, ! thing->point.pos[0], thing->point.pos[1], ! thing->point.vel[0], thing->point.vel[0]); ! fprintf(stderr, "sumVel = %f ---> ", task->sumVel); ! } ! task->sumVel += ELL_3V_LEN(thing->point.vel); ! if (task->pctx->verbose) { ! fprintf(stderr, "%f (exists %d)\n", task->sumVel, ! AIR_EXISTS(task->sumVel)); ! } ! if (!AIR_EXISTS(task->sumVel)) { ! sprintf(err, "%s(%d): sumVel went NaN (from vel (%g,%g,%g), force " ! "(%g,%g,%g)) on thing %d (%d verts) %p of bin %d", ! me, task->threadIdx, ! thing->point.vel[0], ! thing->point.vel[1], ! thing->point.vel[2], ! thing->point.frc[0], ! thing->point.frc[1], ! thing->point.frc[2], ! thing->ttaagg, thing->vertNum, AIR_CAST(void*, thing), binIdx); ! biffAdd(PUSH, err); return 1; } } else { ! /* no mass */ ! ELL_3V_SCALE_INCR(thing->point.pos, step, thing->point.frc); } task->thingNum += 1; /* while _pushProbe clamps positions to inside domain before calling gageProbe, we can exert no such control over the gageProbe called within tenFiberTraceSet. So for now, things turn to points ! as soon as they leave the domain. Tough fucking luck. */ /* sample field at new point location */ inside = _pushProbe(task, &(thing->point)); --- 705,787 ---- /* update dynamics: applies equally to points and tractlets */ ! mass = task->pctx->tlUse ? _pushThingMass(task->pctx, thing) : 1.0; ! if (task->pctx->verbose) { ! fprintf(stderr, "vel(%f,%f) * step(%f) -+-> pos(%f,%f)\n", ! thing->point.vel[0], thing->point.vel[0], step, ! thing->point.pos[0], thing->point.pos[1]); ! fprintf(stderr, "frc(%f,%f)*step(%f)/mass(%f) (%f) -+-> vel(%f,%f)\n", ! thing->point.frc[0], thing->point.frc[0], ! step, mass, step/mass, ! thing->point.vel[0], thing->point.vel[1]); ! } ! if (!ELL_3V_EXISTS(thing->point.vel)) { ! ELL_3V_SET(thing->point.vel, 0, 0, 0); ! } ! if (2 == task->pctx->dimIn) { ! double posIdx[4], posWorld[4]; ! ELL_3V_COPY(posWorld, thing->point.pos); posWorld[3] = 1.0; ! ELL_4MV_MUL(posIdx, task->pctx->gctx->shape->WtoI, posWorld); ! ELL_4V_HOMOG(posIdx, posIdx); ! posIdx[task->pctx->sliceAxis] = 0.0; ! ELL_4MV_MUL(posWorld, task->pctx->gctx->shape->ItoW, posIdx); ! ELL_34V_HOMOG(thing->point.pos, posWorld); ! } ! if (!ELL_3V_EXISTS(thing->point.frc)) { ! ELL_3V_SET(thing->point.frc, 0, 0, 0); ! } ! if (task->pctx->velWarp) { ! /* HEY: try to optimize out some sqrt()s */ ! double dpos[3], ndpos[3], dposLen, tmp[3], rad, speed, wspeed; ! /* dpos = step that we would take, without velocity warping */ ! ELL_3V_SCALE(dpos, step*step/mass, thing->point.frc); ! dposLen = ELL_3V_LEN(dpos); ! if (dposLen) { ! ELL_3V_SCALE(ndpos, 1.0/dposLen, dpos); ! } else { ! ELL_3V_SET(ndpos, 1, 0, 0); } + TEN_T3V_MUL(tmp, thing->point.inv, ndpos); + /* speed = magnitude of velocity in glyphs per time step */ + rad = task->pctx->scale/ELL_3V_LEN(tmp); + speed = dposLen / rad; + wspeed = task->pctx->velWarp*speed/(task->pctx->velWarp + speed); + /* fix thing->point.vel accordingly */ + ELL_3V_SCALE_INCR(thing->point.vel, + task->pctx->scale*wspeed/step, ndpos); } else { ! ELL_3V_SCALE_INCR(thing->point.vel, step/mass, thing->point.frc); ! } ! if (task->pctx->verbose) { ! fprintf(stderr, "thing %d: pos(%f,%f); vel(%f,%f)\n", ! thing->ttaagg, ! thing->point.pos[0], thing->point.pos[1], ! thing->point.vel[0], thing->point.vel[0]); ! fprintf(stderr, "sumVel = %f ---> ", task->sumVel); ! } ! task->sumVel += ELL_3V_LEN(thing->point.vel); ! if (task->pctx->verbose) { ! fprintf(stderr, "%f (exists %d)\n", task->sumVel, ! AIR_EXISTS(task->sumVel)); } + if (!AIR_EXISTS(task->sumVel)) { + sprintf(err, "%s(%d): sumVel went NaN (from vel (%g,%g,%g), force " + "(%g,%g,%g)) on thing %d (%d verts) %p of bin %d", + me, task->threadIdx, + thing->point.vel[0], + thing->point.vel[1], + thing->point.vel[2], + thing->point.frc[0], + thing->point.frc[1], + thing->point.frc[2], + thing->ttaagg, thing->vertNum, AIR_CAST(void*, thing), binIdx); + biffAdd(PUSH, err); return 1; + } + ELL_3V_SCALE_INCR(thing->point.pos, step, thing->point.vel); + task->thingNum += 1; /* while _pushProbe clamps positions to inside domain before calling gageProbe, we can exert no such control over the gageProbe called within tenFiberTraceSet. So for now, things turn to points ! as soon as they leave the domain. Tough luck. */ /* sample field at new point location */ inside = _pushProbe(task, &(thing->point)); Index: methodsPush.c =================================================================== RCS file: /cvsroot/teem/teem/src/push/methodsPush.c,v retrieving revision 1.23 retrieving revision 1.24 diff -C2 -d -r1.23 -r1.24 *** methodsPush.c 21 Apr 2006 17:52:29 -0000 1.23 --- methodsPush.c 2 Jul 2006 23:30:30 -0000 1.24 *************** *** 40,44 **** thg->vertNum = vertNum; if (1 == vertNum) { ! thg->vert = &(thg->point); } else { thg->vert = (pushPoint *)calloc(vertNum, sizeof(pushPoint)); --- 40,44 ---- thg->vertNum = vertNum; if (1 == vertNum) { ! thg->vert = &(thg->point); /* HEY: is this just inviting confusion? */ } else { thg->vert = (pushPoint *)calloc(vertNum, sizeof(pushPoint)); *************** *** 58,62 **** if (thg) { ! if (thg->vert != &(thg->point)) { thg->vert = (pushPoint *)airFree(thg->vert); } --- 58,62 ---- if (thg) { ! if (1 < thg->vertNum) { thg->vert = (pushPoint *)airFree(thg->vert); } *************** *** 78,81 **** --- 78,82 ---- pctx->drag = 0.1; pctx->preDrag = 1.0; + pctx->velWarp = 0.0; pctx->step = 0.01; pctx->mass = 1.0; *************** *** 89,93 **** --- 90,101 ---- pctx->tlSoft = 0.0; pctx->minMeanVel = 0.0; + pctx->tlUse = AIR_FALSE; + pctx->tlFrenet = AIR_FALSE; + pctx->singleBin = AIR_FALSE; + pctx->driftCorrect = AIR_FALSE; + pctx->detReject = AIR_FALSE; + pctx->verbose = 0; pctx->seed = 42; + pctx->tlStepNum = 5; pctx->binIncr = 128; pctx->thingNum = 0; *************** *** 104,111 **** pctx->gravScl = 0.0; pctx->seedThresh = 0.0; - pctx->singleBin = AIR_FALSE; - pctx->driftCorrect = AIR_TRUE; - pctx->detReject = AIR_FALSE; - pctx->verbose = 0; pctx->force = NULL; pctx->ksp00 = nrrdKernelSpecNew(); --- 112,115 ---- *************** *** 123,136 **** pctx->gctx = NULL; pctx->tpvl = NULL; pctx->fctx = NULL; pctx->dimIn = 0; ! pctx->sliceAxis = 5280; ! /* binsEdge and binNum are found later */ ELL_3V_SET(pctx->binsEdge, 0, 0, 0); pctx->binNum = 0; - pctx->finished = AIR_FALSE; pctx->stageIdx = pctx->binIdx = 0; pctx->bin = NULL; pctx->maxDist = AIR_NAN; pctx->meanVel = 0; pctx->time0 = pctx->time1 = 0; --- 127,144 ---- pctx->gctx = NULL; pctx->tpvl = NULL; + pctx->ipvl = NULL; pctx->fctx = NULL; + pctx->finished = AIR_FALSE; pctx->dimIn = 0; ! pctx->sliceAxis = 5280; /* an invalid value */ ! /* binsEdge and binNum are set later */ ELL_3V_SET(pctx->binsEdge, 0, 0, 0); pctx->binNum = 0; pctx->stageIdx = pctx->binIdx = 0; pctx->bin = NULL; pctx->maxDist = AIR_NAN; + pctx->maxEval = AIR_NAN; + pctx->meanEval = AIR_NAN; + pctx->maxDet = AIR_NAN; pctx->meanVel = 0; pctx->time0 = pctx->time1 = 0; |