From: <kin...@us...> - 2008-05-19 21:47:17
|
Revision: 3842 http://teem.svn.sourceforge.net/teem/?rev=3842&view=rev Author: kindlmann Date: 2008-05-19 14:47:13 -0700 (Mon, 19 May 2008) Log Message: ----------- less broken than before Modified Paths: -------------- teem/trunk/src/pull/actionPull.c teem/trunk/src/pull/constraints.c teem/trunk/src/pull/contextPull.c teem/trunk/src/pull/pointPull.c teem/trunk/src/pull/privatePull.h teem/trunk/src/pull/pull.h Modified: teem/trunk/src/pull/actionPull.c =================================================================== --- teem/trunk/src/pull/actionPull.c 2008-05-15 13:34:21 UTC (rev 3841) +++ teem/trunk/src/pull/actionPull.c 2008-05-19 21:47:13 UTC (rev 3842) @@ -28,13 +28,12 @@ ** issues: ** does everything work on the first iteration ** how to handle the needed extra probe for d strength / d scale -** does it eventually catch non-existant energy or force ** how are force/energy along scale handled differently than in space? */ -int praying = 0; -double prayCorner[2][2][3]; -size_t prayRes[2] = {400,400}; +int _pullPraying = 0; +double _pullPrayCorner[2][2][3]; +size_t _pullPrayRes[2] = {60,20}; unsigned int _neighBinPoints(pullTask *task, pullBin *bin, pullPoint *point) { @@ -113,29 +112,20 @@ return enr; } -/* -** meanNeighDist is allowed to be NULL. If it is non-NULL, -** it must be set, and must be < 0 if there are no neighbors -*/ double _energyFromPoints(pullTask *task, pullBin *bin, pullPoint *point, /* output */ - double egradSum[4], double *meanNeighDist) { + double egradSum[4]) { /* char me[]="_energyFromPoints"; */ - double energySum, distSqSum, spaDistSqMax; + double energySum, distSqSum, spaDistSqMax, wghtSum; int nopt, /* optimiziation: we sometimes re-use neighbor lists */ ntrue; /* we search all possible neighbors, stored in the bins (either because !nopt, or, this iter we learn true subset of interacting neighbors). This could also be called "dontreuse" or something like that */ - unsigned int nidx, neighCount, + unsigned int nidx, nnum; /* how much of task->neigh[] we use */ - if (meanNeighDist && !egradSum) { - /* shouldn't happen */ - return AIR_NAN; - } - /* set nopt and ntrue */ if (task->pctx->neighborTrueProb < 1) { nopt = AIR_TRUE; @@ -178,9 +168,12 @@ fprintf(stderr, "%s: radiusSpace = %g -> spaDistSqMax = %g\n", me, task->pctx->radiusSpace, spaDistSqMax); */ + wghtSum = 0; energySum = 0; distSqSum = 0; - neighCount = 0; + point->neighInterNum = 0; + point->neighDist = 0.0; + point->neighMode = 0.0; if (egradSum) { ELL_4V_SET(egradSum, 0, 0, 0, 0); } @@ -212,9 +205,15 @@ if (egradSum) { ELL_4V_INCR(egradSum, egrad); if (ELL_4V_DOT(egrad, egrad)) { - if (meanNeighDist) { - neighCount++; - distSqSum += spaDistSq; + point->neighInterNum++; + point->neighDist = spaDistSq; + if (task->pctx->ispec[pullInfoTangentMode]) { + double w, m; + m = _pullPointScalar(task->pctx, herPoint, pullInfoTangentMode, + NULL, NULL); + w = 1.0/spaDistSq; + point->neighMode += w*m; + wghtSum += w; } if (nopt && ntrue) { unsigned int ii; @@ -225,13 +224,13 @@ } } - /* finish computing mean distance to neighbors */ - if (meanNeighDist) { - if (neighCount) { - *meanNeighDist = sqrt(distSqSum/neighCount); - } else { - *meanNeighDist = -1; - } + /* finish computing things averaged over neighbors */ + if (point->neighInterNum) { + point->neighDist = sqrt(point->neighDist/point->neighInterNum); + point->neighMode /= wghtSum; + } else { + point->neighDist = -1; + point->neighMode = AIR_NAN; } return energySum; @@ -266,7 +265,8 @@ } if (task->pctx->ispec[pullInfoHeight] && !task->pctx->ispec[pullInfoHeight]->constraint - && !task->pctx->ispec[pullInfoHeightLaplacian]->constraint) { + && (!task->pctx->ispec[pullInfoHeightLaplacian] + || !task->pctx->ispec[pullInfoHeightLaplacian]->constraint)) { energy += _pullPointScalar(task->pctx, point, pullInfoHeight, grad3, NULL); if (egradSum) { @@ -294,8 +294,8 @@ double _pullPointEnergyTotal(pullTask *task, pullBin *bin, pullPoint *point, /* output */ - double force[4], double *neighDist) { - /* char me[]="_pullPointEnergyTotal"; */ + double force[4]) { + char me[]="_pullPointEnergyTotal"; double enrIm, enrPt, egradIm[4], egradPt[4], energy; ELL_4V_SET(egradIm, 0, 0, 0, 0); /* sssh */ @@ -303,7 +303,7 @@ enrIm = _energyFromImage(task, point, force ? egradIm : NULL); enrPt = _energyFromPoints(task, bin, point, - force ? egradPt : NULL, neighDist); + force ? egradPt : NULL); energy = AIR_LERP(task->pctx->alpha, enrIm, enrPt); /* fprintf(stderr, "!%s(%u): energy = lerp(%g, im %g, pt %g) = %g\n", me, @@ -340,9 +340,27 @@ } } } + if (!AIR_EXISTS(energy)) { + fprintf(stderr, "!%s(%u): HEY! non-exist energy %g\n", me, point->idtag, + energy); + } return energy; } +double +_pullDistLimit(pullTask *task, pullPoint *point) { + double ret; + + if (point->neighDist < 0 /* no neighbors! */ + || pullEnergyZero == task->pctx->energySpec->energy) { + ret = task->pctx->radiusSpace; + } else { + ret = point->neighDist; + } + /* task->pctx->constraintVoxelSize might be considered here */ + return ret; +} + int _pullPointProcess(pullTask *task, pullBin *bin, pullPoint *point) { char me[]="pullPointProcess", err[BIFF_STRLEN]; @@ -351,28 +369,42 @@ in a per-iteration way */ int stepBad, giveUp; + if (point->debug) { + return 0; + } if (!point->stepEnergy) { sprintf(err, "%s: whoa, point %u step is zero!", me, point->idtag); biffAdd(PULL, err); return 1; } - if (0 && 541 == point->idtag && 10 == task->pctx->iter) { + if (0 && 162 == point->idtag) { fprintf(stderr, "!%s: !!!!!!!!!!!! praying ON!\n", me); - praying = AIR_TRUE; - ELL_3V_COPY(prayCorner[0][0], point->pos); + _pullPraying = AIR_TRUE; + ELL_3V_COPY(_pullPrayCorner[0][0], point->pos); + _pullPrayCorner[0][0][2] -= 1; + ELL_3V_COPY(_pullPrayCorner[0][1], point->pos); + _pullPrayCorner[0][1][2] += 1; + fprintf(stderr, "!%s: corner[0][0] = %g %g %g\n", me, + _pullPrayCorner[0][0][0], + _pullPrayCorner[0][0][1], + _pullPrayCorner[0][0][2]); + fprintf(stderr, "!%s: corner[0][1] = %g %g %g\n", me, + _pullPrayCorner[0][1][0], + _pullPrayCorner[0][1][1], + _pullPrayCorner[0][1][2]); } else { - praying = AIR_FALSE; + _pullPraying = AIR_FALSE; } - if (praying) { + if (_pullPraying) { fprintf(stderr, "%s: =============================== (%u) hi @ %g %g %g\n", me, point->idtag, point->pos[0], point->pos[1], point->pos[2]); } - energyOld = _pullPointEnergyTotal(task, bin, point, force, &distLimit); - if (praying) { + energyOld = _pullPointEnergyTotal(task, bin, point, force); + if (_pullPraying) { fprintf(stderr, "!%s: =================== point %u has:\n " " energy = %g ; ndist = %g, force %g %g %g %g\n", me, - point->idtag, energyOld, distLimit, + point->idtag, energyOld, point->neighDist, force[0], force[1], force[2], force[3]); } if (!( AIR_EXISTS(energyOld) && ELL_4V_EXISTS(force) )) { @@ -383,29 +415,10 @@ if (task->pctx->constraint) { /* we have a constraint, so do something to get the force more tangential to the constriant surface */ - double vec[4], proj[9], nvec[3], outer[9], pfrc[3], len; - - switch (task->pctx->constraint) { - case pullInfoHeight: - /* HEY: need an approximate constraint surface tangent */ - ELL_3V_SET(vec, 0, 0, 0); - break; - case pullInfoHeightLaplacian: - _pullPointScalar(task->pctx, point, pullInfoHeight, vec, NULL); - break; - case pullInfoIsovalue: - _pullPointScalar(task->pctx, point, pullInfoIsovalue, vec, NULL); - break; - } - ELL_3V_NORM(nvec, vec, len); - if (len) { - ELL_3MV_OUTER(outer, nvec, nvec); - ELL_3M_IDENTITY_SET(proj); - ELL_3M_SUB(proj, proj, outer); - ELL_3MV_MUL(pfrc, proj, force); - ELL_3V_COPY(force, pfrc); - } - /* if !len, gradient is zero, and no change is made */ + double proj[9], pfrc[3]; + _pullConstraintTangent(task, point, proj); + ELL_3MV_MUL(pfrc, proj, force); + ELL_3V_COPY(force, pfrc); } point->status = 0; @@ -418,11 +431,7 @@ point->energy = energyOld; return 0; } - if (distLimit < 0 /* no neighbors! */ - || pullEnergyZero == task->pctx->energySpec->energy) { - distLimit = task->pctx->radiusSpace; - } - /* maybe voxel size should also be considered for finding distLimit */ + distLimit = _pullDistLimit(task, point); /* find capscl */ ELL_3V_SCALE(capvec, point->stepEnergy, force); @@ -432,10 +441,75 @@ } else { capscl = 1; } - if (praying) { - fprintf(stderr, "%s: ======= (%u) capscl = %g\n", me, point->idtag, capscl); + if (_pullPraying) { + fprintf(stderr, "%s: ======= (%u) capscl = %g\n", me, + point->idtag, capscl); } + if (_pullPraying) { + double nfrc[3], len, phold[4], ee; + int cfail; + ELL_4V_COPY(phold, point->pos); + + fprintf(stderr, "!%s(%u,%u): energy(%g,%g,%g) = %f (original)\n", + me, task->pctx->iter, point->idtag, + point->pos[0], point->pos[1], point->pos[2], energyOld); + + ELL_4V_SCALE_ADD2(point->pos, 1.0, posOld, + capscl*point->stepEnergy, force); + ELL_3V_COPY(_pullPrayCorner[1][0], point->pos); + _pullPrayCorner[1][0][2] -= 1; + ELL_3V_COPY(_pullPrayCorner[1][1], point->pos); + _pullPrayCorner[1][1][2] += 1; + fprintf(stderr, "!%s: corner[1][0] = %g %g %g\n", me, + _pullPrayCorner[1][0][0], + _pullPrayCorner[1][0][1], + _pullPrayCorner[1][0][2]); + fprintf(stderr, "!%s: corner[1][1] = %g %g %g\n", me, + _pullPrayCorner[1][1][0], + _pullPrayCorner[1][1][1], + _pullPrayCorner[1][1][2]); + +#define PROBE(l) if (_pullProbe(task, point)) { \ + sprintf(err, "%s: while praying", me); \ + biffAdd(PULL, err); return 1; \ + } \ + (l) = _pullPointScalar(task->pctx, point, \ + pullInfoHeight, NULL, NULL); + if (1) { + double *enr, *lpl, uu, vv, vpos[2][3], ll, mid[3]; + unsigned int ui, vi; + Nrrd *nenr, *nlpl; + nenr = nrrdNew(); + nlpl = nrrdNew(); + nrrdMaybeAlloc_nva(nenr, nrrdTypeDouble, 2, _pullPrayRes); + enr = AIR_CAST(double *, nenr->data); + nrrdMaybeAlloc_nva(nlpl, nrrdTypeDouble, 2, _pullPrayRes); + lpl = AIR_CAST(double *, nlpl->data); + for (vi=0; vi<_pullPrayRes[1]; vi++) { + vv = AIR_AFFINE(0, vi, _pullPrayRes[1]-1, 0, 1); + ELL_3V_LERP(vpos[0], vv, _pullPrayCorner[0][0], _pullPrayCorner[0][1]); + ELL_3V_LERP(vpos[1], vv, _pullPrayCorner[1][0], _pullPrayCorner[1][1]); + for (ui=0; ui<_pullPrayRes[0]; ui++) { + uu = AIR_AFFINE(0, ui, _pullPrayRes[0]-1, 0, 1); + ELL_3V_LERP(point->pos, uu, vpos[0], vpos[1]); + PROBE(ll); + lpl[ui + _pullPrayRes[0]*vi] = ll; + enr[ui + _pullPrayRes[0]*vi] = _pullPointEnergyTotal(task, bin, + point, NULL); + } + } + nrrdSave("nenr.nrrd", nenr, NULL); + nrrdSave("nlpl.nrrd", nlpl, NULL); + nenr = nrrdNuke(nenr); + nlpl = nrrdNuke(nlpl); + } +#undef PROBE + + ELL_4V_COPY(point->pos, phold); + _pullPointEnergyTotal(task, bin, point, NULL); + } + do { int constrFail; @@ -444,7 +518,7 @@ capscl*point->stepEnergy, force); _pullPointHistAdd(point, pullCondEnergyTry); - if (praying) { + if (_pullPraying) { fprintf(stderr, "!%s: ======= (%u) try step %g to pos %g %g %g %g\n", me, point->idtag, capscl*point->stepEnergy, point->pos[0], point->pos[1], point->pos[2], point->pos[3]); @@ -459,31 +533,40 @@ } if (constrFail) { energyNew = AIR_NAN; + if (_pullPraying) { + fprintf(stderr, "!%s: ======= constraint fail\n", me); + } } else { - energyNew = _pullPointEnergyTotal(task, bin, point, NULL, NULL); - if (praying) { - fprintf(stderr, "!%s: ======= e new = %g %s old %g %s\n", me, energyNew, - energyNew > energyOld ? ">" : "<=", energyOld, + energyNew = _pullPointEnergyTotal(task, bin, point, NULL); + if (_pullPraying) { + fprintf(stderr, "!%s: ======= e new = %g %s old %g %s\n", me, + energyNew, energyNew > energyOld ? ">" : "<=", energyOld, energyNew > energyOld ? "!! BADSTEP !!" : "ok"); } } stepBad = constrFail || (energyNew > energyOld); if (stepBad) { point->stepEnergy *= task->pctx->stepScale; - _pullPointHistAdd(point, pullCondEnergyBad); + if (constrFail) { + _pullPointHistAdd(point, pullCondConstraintFail); + } else { + _pullPointHistAdd(point, pullCondEnergyBad); + } /* you have a problem if you had a non-trivial force, but you can't ever seem to take a small enough step to reduce energy */ if (point->stepEnergy < 0.000000000000001) { - fprintf(stderr, "\n%s: %u (%g,%g,%g,%g) stepEnr %g stuck; " - "capscl %g <<<<\n\n\n", - me, point->idtag, - point->pos[0], point->pos[1], point->pos[2], point->pos[3], - point->stepEnergy, capscl); + if (task->pctx->verbose > 1) { + fprintf(stderr, "\n%s: %u (%g,%g,%g,%g) stepEnr %g stuck; " + "capscl %g <<<<\n\n\n", me, point->idtag, + point->pos[0], point->pos[1], point->pos[2], point->pos[3], + point->stepEnergy, capscl); + } /* This point is fuct, may as well reset its step, maybe things will go better next time. Without this resetting, it will stay effectively frozen */ ELL_4V_COPY(point->pos, posOld); - point->stepEnergy = task->pctx->stepInitial; + energyNew = energyOld; /* to be copied into point->energy below */ + point->stepEnergy = task->pctx->stepInitial/100; point->status = 1; giveUp = AIR_TRUE; } @@ -496,6 +579,11 @@ /* not recorded for the sake of this function, but for system accounting */ point->energy = energyNew; + if (!AIR_EXISTS(energyNew)) { + sprintf(err, "%s: point %u has non-exist final energy %g\n", + me, point->idtag, energyNew); + biffAdd(PULL, err); return 1; + } return 0; } Modified: teem/trunk/src/pull/constraints.c =================================================================== --- teem/trunk/src/pull/constraints.c 2008-05-15 13:34:21 UTC (rev 3841) +++ teem/trunk/src/pull/constraints.c 2008-05-19 21:47:13 UTC (rev 3842) @@ -47,13 +47,13 @@ (v) = _pullPointScalar(task->pctx, point, \ pullInfoIsovalue, (g), NULL); \ (av) = AIR_ABS(v) -#define SAVE(state, aval, val, grad, pos) \ +#define SAVE(state, aval, val, grad, pos) \ state[0] = aval; \ state[1] = val; \ ELL_3V_COPY(state + 1 + 1, grad); \ ELL_3V_COPY(state + 1 + 1 + 3, pos) -#define RESTORE(aval, val, grad, pos, state) \ - aval = state[0]; \ +#define RESTORE(aval, val, grad, pos, state) \ + aval = state[0]; \ val = state[1]; \ ELL_3V_COPY(grad, state + 1 + 1); \ ELL_3V_COPY(pos, state + 1 + 1 + 3) @@ -71,7 +71,7 @@ for-loop, instead of a nested do-while loop */ grad[4], dir[3], len, state[1 + 1 + 3 + 3]; unsigned int iter = 0; /* 0: initial probe, 1..iterMax: probes in loop */ - + PROBE(val, aval, grad); SAVE(state, aval, val, grad, point->pos); hack = 1; @@ -226,45 +226,66 @@ ? 1 - airIntPow(1-m, p) \ : airIntPow(m+1, p) - 1) +#define MODENAVG(pnt, m) \ + ( ((pnt)->neighInterNum*(pnt)->neighMode + (m))/ \ + (1 + (pnt)->neighInterNum) ) + + static int -probeHeight(pullTask *task, pullPoint *point, int tang2Use, int useMode, +probeHeight(pullTask *task, pullPoint *point, /* output */ - double *heightP, double grad[3], double hess[9], double proj[9]) { + double *heightP, double grad[3], double hess[9]) { char me[]="probeHeight", err[BIFF_STRLEN]; - double *tng; if (_pullProbe(task, point)) { sprintf(err, "%s: trouble", me); biffAdd(PULL, err); return 1; } *heightP = _pullPointScalar(task->pctx, point, pullInfoHeight, grad, hess); + return 0; +} + +static void +creaseProj(pullTask *task, pullPoint *point, int tang2Use, int modeUse, + /* output */ + double proj[9]) { + char me[]="creaseProj"; + double *tng; + tng = point->info + task->pctx->infoIdx[pullInfoTangent1]; + if (_pullPraying) { + fprintf(stderr, "!%s: tng1 = %g %g %g\n", me, tng[0], tng[1], tng[2]); + } ELL_3MV_OUTER(proj, tng, tng); if (tang2Use) { double proj2[9]; tng = point->info + task->pctx->infoIdx[pullInfoTangent2]; ELL_3MV_OUTER(proj2, tng, tng); - if (useMode) { + if (modeUse) { double mode; ELL_3M_ADD2(proj2, proj, proj2); mode = _pullPointScalar(task->pctx, point, pullInfoTangentMode, NULL, NULL); - mode = MODEWARP(mode, 8); + if (point->neighInterNum) { + mode = MODENAVG(point, mode); + } + mode = MODEWARP(mode, 4); mode = AIR_AFFINE(-1, mode, 1, 0, 1); ELL_3M_LERP(proj, mode, proj, proj2); } else { ELL_3M_ADD2(proj, proj, proj2); } } - return 0; + return; } #define PROBE(height, grad, hess, proj) \ - if (probeHeight(task, point, tang2Use, modeUse, \ - &(height), (grad), (hess), (proj))) { \ + if (probeHeight(task, point, \ + &(height), (grad), (hess))) { \ sprintf(err, "%s: trouble on iter %u", me, iter); \ biffAdd(PULL, err); return 1; \ - } + } \ + creaseProj(task, point, tang2Use, modeUse, proj) #define SAVE(state, height, grad, hess, proj, pos) \ state[0] = height; \ ELL_3V_COPY(state + 1, grad); \ @@ -289,19 +310,31 @@ d1 = ELL_3V_DOT(grad, pdir); \ d2 = ELL_3MV_CONTR(hess, pdir) -double _tmpv[3]; +/* double _tmpv[3]; */ static int constraintSatHght(pullTask *task, pullPoint *point, int tang2Use, int modeUse, double stepMax, unsigned int iterMax, int *constrFailP) { char me[]="constraintSatHght", err[BIFF_STRLEN]; - double val, grad[3], hess[9], proj[9], + double val, grad[3], hess[9], proj[9], _tmpv[3]={0,0,0}, state[1+3+9+9+3], hack, step, d1, d2, pdir[3], plen, pgrad[3]; unsigned int iter = 0; /* 0: initial probe, 1..iterMax: probes in loop */ /* http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization */ + /* + if (_pullPraying) { + pullPoint *npnt; + fprintf(stderr, "%s: starting at %g %g %g\n", me, + point->pos[0], point->pos[1], point->pos[2]); + npnt = pullPointNew(task->pctx); + npnt->debug = AIR_TRUE; + ELL_4V_COPY(npnt->pos, point->pos); + _pullProbe(task, npnt); + pullBinsPointAdd(task->pctx, npnt); + } + */ _pullPointHistAdd(point, pullCondOld); PROBE(val, grad, hess, proj); SAVE(state, val, grad, hess, proj, point->pos); @@ -309,21 +342,31 @@ for (iter=1; iter<=iterMax; iter++) { NORM(d1, d2, pdir, plen, pgrad, grad, hess, proj); step = (d2 <= 0 ? -plen : -d1/d2); - /* - fprintf(stderr, "!%s: iter %u step = (%g <= 0 ? %g : %g) --> %g\n", me, - iter, d2, -plen, -d1/d2, step); - */ + if (_pullPraying) { + fprintf(stderr, "!%s: iter %u step = (%g <= 0 ? %g : %g) --> %g\n", me, + iter, d2, -plen, -d1/d2, step); + } step = step > 0 ? AIR_MIN(stepMax, step) : AIR_MAX(-stepMax, step); - /* - fprintf(stderr, " -> %g, |pdir| = %g\n", step, ELL_3V_LEN(pdir)); - */ - ELL_3V_COPY(_tmpv, point->pos); + if (_pullPraying) { + fprintf(stderr, " -> %g, |pdir| = %g\n", step, ELL_3V_LEN(pdir)); + ELL_3V_COPY(_tmpv, point->pos); + } ELL_3V_SCALE_INCR(point->pos, hack*step, pdir); - ELL_3V_SUB(_tmpv, _tmpv, point->pos); + if (_pullPraying) { + ELL_3V_SUB(_tmpv, _tmpv, point->pos); + fprintf(stderr, " -> moved %g\n", ELL_3V_LEN(_tmpv)); + } + _pullPointHistAdd(point, pullCondConstraintSatA); /* - fprintf(stderr, " -> moved %g\n", ELL_3V_LEN(_tmpv)); + if (_pullPraying) { + pullPoint *npnt; + npnt = pullPointNew(task->pctx); + npnt->debug = AIR_TRUE; + ELL_4V_COPY(npnt->pos, point->pos); + _pullProbe(task, npnt); + pullBinsPointAdd(task->pctx, npnt); + } */ - _pullPointHistAdd(point, pullCondConstraintSatA); PROBE(val, grad, hess, proj); if (val <= state[0]) { if (AIR_ABS(step) < stepMax*task->pctx->constraintStepMin) { @@ -366,10 +409,19 @@ char me[]="_pullConstraintSatisfy", err[BIFF_STRLEN]; double stepMax; unsigned int iterMax; + int wantLine, wantSurf, failLine, failSurf; + double mode, oldPos[4], posLine[4], posSurf[4], + modeLine = AIR_NAN, modeSurf = AIR_NAN; stepMax = task->pctx->constraintVoxelSize; iterMax = task->pctx->constraintIterMax; /* + dlim = _pullDistLimit(task, point); + if (iterMax*stepMax > dlim) { + stepMax = dlim/iterMax; + } + */ + /* fprintf(stderr, "!%s(%d): hi %g %g %g, stepMax = %g, iterMax = %u\n", me, point->idtag, point->pos[0], point->pos[1], point->pos[2], stepMax, iterMax); @@ -388,6 +440,7 @@ } break; case pullInfoHeight: + /* if (constraintSatHght(task, point, (task->pctx->ispec[pullInfoTangent2] ? AIR_TRUE @@ -395,22 +448,152 @@ (task->pctx->ispec[pullInfoTangentMode] ? AIR_TRUE : AIR_FALSE), - stepMax, iterMax, constrFailP)) { + stepMax/2, 2*iterMax, constrFailP)) { sprintf(err, "%s: trouble", me); biffAdd(PULL, err); return 1; } + */ + + wantSurf = wantLine = AIR_FALSE; + mode = AIR_NAN; + if (task->pctx->ispec[pullInfoTangentMode]) { + double md; + if (_pullProbe(task, point)) { + sprintf(err, "%s: trouble", me); + biffAdd(PULL, err); return 1; + } + mode = _pullPointScalar(task->pctx, point, pullInfoTangentMode, + NULL, NULL); + if (task->pctx->iter < 2) { + wantSurf = (mode < 0); + wantLine = !wantSurf; + } else { + if (point->neighInterNum) { + mode = MODENAVG(point, mode); + } + md = MODEWARP(mode, 5); + md = (1 + md)/2; + if (md < 0.01) { + wantSurf = AIR_TRUE; + } else if (md > 0.99) { + wantLine = AIR_TRUE; + } else { + wantSurf = wantLine = AIR_TRUE; + } + } + } else { + if (task->pctx->ispec[pullInfoTangent2]) { + wantLine = AIR_TRUE; + } else if (task->pctx->ispec[pullInfoTangent1]) { + wantSurf = AIR_TRUE; + } + } + if (!( wantLine || wantSurf )) { + sprintf(err, "%s: confusion, should want line (%d) or surf (%d)", me, + wantLine, wantSurf); + biffAdd(PULL, err); return 1; + } + ELL_4V_SET(posSurf, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN); + ELL_4V_SET(posLine, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN); + ELL_4V_COPY(oldPos, point->pos); + if (wantSurf) { + if (constraintSatHght(task, point, AIR_FALSE, + AIR_FALSE, stepMax, iterMax, &failSurf)) { + sprintf(err, "%s: surf trouble", me); + biffAdd(PULL, err); return 1; + } + if (wantLine) { + modeSurf = _pullPointScalar(task->pctx, point, pullInfoTangentMode, + NULL, NULL); + } + ELL_4V_COPY(posSurf, point->pos); + } + /* yes, the starting point of the line constraint is the + already-determined surface constraint */ + if (wantLine) { + if (constraintSatHght(task, point, AIR_TRUE, + AIR_FALSE, stepMax, iterMax, &failLine)) { + sprintf(err, "%s: line trouble", me); + biffAdd(PULL, err); return 1; + } + if (wantSurf) { + modeLine = _pullPointScalar(task->pctx, point, pullInfoTangentMode, + NULL, NULL); + } + ELL_4V_COPY(posLine, point->pos); + } + if (wantSurf && wantLine) { + if (!failSurf && !failLine) { + ELL_4V_LERP(point->pos, (1+mode)/2, posLine, posSurf); + *constrFailP = AIR_FALSE; + } else if (!failSurf) { + ELL_4V_COPY(point->pos, posSurf); + /* ELL_4V_LERP(point->pos, (1+mode)/2, oldPos, posSurf); */ + *constrFailP = AIR_FALSE; + } else { + ELL_4V_COPY(point->pos, oldPos); + *constrFailP = AIR_TRUE; + } + } else if (wantSurf) { + ELL_4V_COPY(point->pos, posSurf); + *constrFailP = failSurf; + } else { + ELL_4V_COPY(point->pos, posLine); + *constrFailP = failLine; + } + break; default: fprintf(stderr, "%s: constraint on %s (%d) unimplemented!!\n", me, airEnumStr(pullInfo, task->pctx->constraint), task->pctx->constraint); } - if (*constrFailP) { - fprintf(stderr, "!%s(%u) bye, fail = %d, %g %g %g\n", me, - point->idtag, *constrFailP, - point->pos[0], point->pos[1], point->pos[2]); - } + /* + fprintf(stderr, "!%s(%u) bye, fail = %d, %g %g %g\n", me, + point->idtag, *constrFailP, + point->pos[0], point->pos[1], point->pos[2]); + */ return 0; } #undef NORMALIZE + +/* +** this can assume that probe() has just been called +*/ +void +_pullConstraintTangent(pullTask *task, pullPoint *point, + /* output */ + double proj[9]) { + double vec[4], nvec[3], outer[9], len, aproj[9]; + + ELL_3M_IDENTITY_SET(proj); + switch (task->pctx->constraint) { + case pullInfoHeight: + creaseProj(task, point, + (task->pctx->ispec[pullInfoTangent2] + ? AIR_TRUE + : AIR_FALSE), + (task->pctx->ispec[pullInfoTangentMode] + ? AIR_TRUE + : AIR_FALSE), + aproj); + ELL_3M_SUB(proj, proj, aproj); + break; + case pullInfoHeightLaplacian: + case pullInfoIsovalue: + if (pullInfoHeightLaplacian == task->pctx->constraint) { + /* using gradient of height as approx normal to laplacian 0-crossing */ + _pullPointScalar(task->pctx, point, pullInfoHeight, vec, NULL); + } else { + _pullPointScalar(task->pctx, point, pullInfoIsovalue, vec, NULL); + } + ELL_3V_NORM(nvec, vec, len); + if (len) { + ELL_3MV_OUTER(outer, nvec, nvec); + ELL_3M_SUB(proj, proj, outer); + } + break; + } + return; +} Modified: teem/trunk/src/pull/contextPull.c =================================================================== --- teem/trunk/src/pull/contextPull.c 2008-05-15 13:34:21 UTC (rev 3841) +++ teem/trunk/src/pull/contextPull.c 2008-05-19 21:47:13 UTC (rev 3842) @@ -54,7 +54,7 @@ pctx->opporStepScale = 1.0; pctx->stepScale = 0.5; pctx->energyImprovMin = 0.01; - pctx->constraintStepMin = 0.00001; + pctx->constraintStepMin = 0.0001; pctx->wall = 1; pctx->seedRNG = 42; @@ -306,7 +306,7 @@ CHECK(radiusSpace, 0.000001, 15.0); CHECK(neighborTrueProb, 0.02, 1.0); CHECK(probeProb, 0.02, 1.0); - CHECK(opporStepScale, 1.0, 2.0); + CHECK(opporStepScale, 1.0, 10.0); CHECK(stepScale, 0.01, 0.99); CHECK(energyImprovMin, -0.2, 1.0); CHECK(constraintStepMin, 0.00000000000000001, 0.1); @@ -334,6 +334,7 @@ pullPoint *point; pointNum = _pullPointNumber(pctx); + fprintf(stderr, "!%s: pointNum = %u\n", me, pointNum); if (pointNum != pctx->idtagNext) { /* HEY: should really be checking that all IDs between 0 and pointNum-1 contiguous are used. Will have to do something @@ -383,8 +384,9 @@ eval[2] = AIR_MIN(eceil, 1.0/eval[2]); ELL_3V_NORM(eval, eval, len); tenMakeSingle_d(tout, 1, eval, evec); - } else if (pctx->ispec[pullInfoHeightGradient] - || pctx->ispec[pullInfoIsovalueGradient]) { + } else if (pctx->constraint + && (pctx->ispec[pullInfoHeightGradient] + || pctx->ispec[pullInfoIsovalueGradient])) { double *grad, norm[3], len, mat[9], out[9]; grad = point->info + (pctx->ispec[pullInfoHeightGradient] ? pctx->infoIdx[pullInfoHeightGradient] @@ -392,7 +394,7 @@ ELL_3V_NORM(norm, grad, len); ELL_3MV_OUTER(out, norm, norm); ELL_3M_IDENTITY_SET(mat); - ELL_3M_SCALE_INCR(mat, -0.92, out); + ELL_3M_SCALE_INCR(mat, -0.95, out); TEN_M2T(tout, mat); tout[0] = 1; } else { @@ -550,6 +552,9 @@ case pullCondEnergyBad: ELL_3V_SET(rgb, 255, 0, 0); break; + case pullCondConstraintFail: + ELL_3V_SET(rgb, 255, 0, 255); + break; case pullCondNew: ELL_3V_SET(rgb, 128, 255, 128); break; Modified: teem/trunk/src/pull/pointPull.c =================================================================== --- teem/trunk/src/pull/pointPull.c 2008-05-15 13:34:21 UTC (rev 3841) +++ teem/trunk/src/pull/pointPull.c 2008-05-19 21:47:13 UTC (rev 3842) @@ -55,11 +55,14 @@ pnt->neighArr = airArrayNew((void**)&(pnt->neighPoint), &(pnt->neighNum), sizeof(pullPoint *), PULL_POINT_NEIGH_INCR); pnt->neighArr->noReallocWhenSmaller = AIR_TRUE; + pnt->neighDist = pnt->neighMode = AIR_NAN; + pnt->neighInterNum = 0; pnt->phist = NULL; pnt->phistNum = 0; pnt->phistArr = airArrayNew((void**)&(pnt->phist), &(pnt->phistNum), 5*sizeof(double), 32); pnt->status = 0; + pnt->debug = AIR_FALSE; ELL_4V_SET(pnt->pos, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN); pnt->energy = AIR_NAN; ELL_4V_SET(pnt->force, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN); @@ -86,10 +89,14 @@ dst->neighPoint = src->neighPoint; dst->neighNum = src->neighNum; dst->neighArr = src->neighArr; + dst->neighDist = src->neighDist; + dst->neighMode = src->neighMode; + dst->neighInterNum = src->neighInterNum; dst->phist = src->phist; dst->phistNum = src->phistNum; dst->phistArr = src->phistArr; dst->status = src->status; + dst->debug = src->debug; ELL_4V_COPY(dst->pos, src->pos); dst->energy = src->energy; ELL_4V_COPY(dst->force, src->force); @@ -154,6 +161,9 @@ bin = pctx->bin + binIdx; for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) { point = bin->point[pointIdx]; + if (point->debug) { + continue; + } sum += point->energy; } } @@ -170,6 +180,9 @@ bin = pctx->bin + binIdx; for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) { point = bin->point[pointIdx]; + if (point->debug) { + continue; + } point->stepEnergy *= scale; } } @@ -190,6 +203,9 @@ pointNum += bin->pointNum; for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) { point = bin->point[pointIdx]; + if (point->debug) { + continue; + } sum += point->stepEnergy; } } @@ -211,6 +227,9 @@ pointNum += bin->pointNum; for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) { point = bin->point[pointIdx]; + if (point->debug) { + continue; + } sum += point->stepConstr; } } @@ -269,15 +288,22 @@ infoIdx = pctx->infoIdx; ispec = pctx->ispec[sclInfo]; scl = point->info[infoIdx[sclInfo]]; + scl = (scl - ispec->zero)*ispec->scale; + /* + learned: this wasn't thought through: the idea was that the height + *laplacian* answer should be transformed by the *height* zero and + scale. scale might make sense, but not zero. This cost a few + hours of tracking down the fact that the first zero-crossing + detection phase of the lapl constraint was failing because the + laplacian was facilating around hspec->zero, not 0.0 ... if (pullInfoHeightLaplacian == sclInfo) { - /* NOTE: _pullContextCheck makes sure that Height is set - if the Laplacian is requested */ const pullInfoSpec *hspec; hspec = pctx->ispec[pullInfoHeight]; scl = (scl - hspec->zero)*hspec->scale; } else { scl = (scl - ispec->zero)*ispec->scale; } + */ /* fprintf(stderr, "%s = (%g - %g)*%g = %g*%g = %g = %g\n", airEnumStr(pullInfo, sclInfo), @@ -314,6 +340,11 @@ continue; } if (task->vol[ii]->ninSingle) { + /* + fprintf(stderr, "!%s: probing vol[%u] \"%s\" @ wsp %g %g %g\n", me, + ii, task->vol[ii]->name, + point->pos[0], point->pos[1], point->pos[2]); + */ gret = gageProbeSpace(task->vol[ii]->gctx, point->pos[0], point->pos[1], point->pos[2], AIR_FALSE, AIR_TRUE); @@ -343,6 +374,16 @@ alen = _pullInfoAnswerLen[ii]; aidx = task->pctx->infoIdx[ii]; _pullInfoAnswerCopy[alen](point->info + aidx, task->ans[ii]); + if (0 && 1 == alen) { + pullVolume *vol; + pullInfoSpec *isp; + isp = task->pctx->ispec[ii]; + vol = task->pctx->vol[isp->volIdx]; + fprintf(stderr, "!%s: info[%u] %s: %s \"%s\" = %g\n", me, + ii, airEnumStr(pullInfo, ii), + airEnumStr(vol->kind->enm, isp->item), + vol->name, task->ans[ii][0]); + } } } return 0; @@ -365,7 +406,7 @@ double *posData; airRandMTState *rng; pullBin *bin; - int reject; + int reject, threshFail, constrFail; airArray *mop; Nrrd *npos; @@ -408,6 +449,7 @@ biffAdd(PULL, err); return 1; } } else { + unsigned int threshFailCount = 0, constrFailCount = 0; do { _pullPointHistInit(point); ELL_3V_SET(point->pos, @@ -428,20 +470,39 @@ sprintf(err, "%s: probing pointIdx %u of world", me, pointIdx); biffAdd(PULL, err); return 1; } - reject = AIR_FALSE; if (pctx->ispec[pullInfoSeedThresh]) { double val; val = _pullPointScalar(pctx, point, pullInfoSeedThresh, NULL, NULL); - reject |= val < 0; + /* + fprintf(stderr, "!%s: val(%g %g %g) = %g -> %g\n", me, + point->pos[0], point->pos[1], point->pos[2], + point->info[pctx->infoIdx[pullInfoSeedThresh]], val); + */ + threshFailCount += (threshFail = val < 0); + } else { + threshFail = AIR_FALSE; } - if (!reject && pctx->constraint) { - int constrFail; + if (!threshFail && pctx->constraint) { if (_pullConstraintSatisfy(pctx->task[0], point, &constrFail)) { sprintf(err, "%s: trying constraint on point %u", me, pointIdx); biffAdd(PULL, err); return 1; } - reject = constrFail; + constrFailCount += constrFail; + } else { + constrFail = AIR_FALSE; } + reject = threshFail || constrFail; + if (reject) { + if (threshFailCount + constrFailCount > 1000) { + sprintf(err, "%s: failed too often placing %u " + "(thresh %d/%u, constr %d/%u)", + me, pointIdx, threshFail, threshFailCount, + constrFail, constrFailCount); + biffAdd(PULL, err); return 1; + } + } else { + threshFailCount = constrFailCount = 0; + } } while (reject); } if (pullBinsPointAdd(pctx, point)) { @@ -469,7 +530,7 @@ for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) { point = bin->point[pointIdx]; point->energy = _pullPointEnergyTotal(pctx->task[0], bin, point, - point->force, NULL); + point->force); } } Modified: teem/trunk/src/pull/privatePull.h =================================================================== --- teem/trunk/src/pull/privatePull.h 2008-05-15 13:34:21 UTC (rev 3841) +++ teem/trunk/src/pull/privatePull.h 2008-05-19 21:47:13 UTC (rev 3842) @@ -48,10 +48,16 @@ extern void _pullTaskFinish(pullContext *pctx); /* actionPull.c */ +extern int _pullPraying; +extern double _pullPrayCorner[2][2][3]; +extern size_t _pullPrayRes[2]; +extern double _pullDistLimit(pullTask *task, pullPoint *point); /* constraints.c */ extern int _pullConstraintSatisfy(pullTask *task, pullPoint *point, int *constrFailP); +extern void _pullConstraintTangent(pullTask *task, pullPoint *point, + double proj[9]); /* pointPull.c */ extern double _pullPointScalar(const pullContext *pctx, @@ -67,7 +73,7 @@ extern double _pullEnergyTotal(const pullContext *pctx); extern double _pullPointEnergyTotal(pullTask *task, pullBin *bin, pullPoint *point, - double force[4], double *neighDist); + double force[4]); extern int _pullProbe(pullTask *task, pullPoint *point); extern void _pullPointStepEnergyScale(pullContext *pctx, double scale); extern int _pullPointSetup(pullContext *pctx); Modified: teem/trunk/src/pull/pull.h =================================================================== --- teem/trunk/src/pull/pull.h 2008-05-15 13:34:21 UTC (rev 3841) +++ teem/trunk/src/pull/pull.h 2008-05-19 21:47:13 UTC (rev 3842) @@ -120,8 +120,9 @@ pullCondConstraintSatA, /* 2 */ pullCondConstraintSatB, /* 3 */ pullCondEnergyTry, /* 4 */ - pullCondEnergyBad, /* 5 */ - pullCondNew, /* 6 */ + pullCondConstraintFail, /* 5 */ + pullCondEnergyBad, /* 6 */ + pullCondNew, /* 7 */ pullCondLast }; @@ -154,12 +155,15 @@ unsigned int neighNum; airArray *neighArr; /* airArray around neigh and neigNum (no callbacks used here) */ + double neighDist, neighMode; + unsigned int neighInterNum; double *phist; /* history of positions tried in the last iter, in sets of 5 doubles: (x,y,z,t,info) */ unsigned int phistNum; /* number of positions stored */ airArray *phistArr; /* airArray around phist */ unsigned int status; /* bit-flag of status info, though right now its just a boolean for having gotten stuck */ + int debug; /* this point is only here for debugging */ double pos[4], /* position in space and scale */ energy, /* energy accumulator for this iteration */ force[4], /* force accumulator for this iteration */ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |