Menu
▾
▴
[Teem-commit] teem/src/push action.c, 1.66, 1.67 binning.c, 1.18, 1.19 corePush.c, 1.25, 1.26 forces.c, 1.21, 1.22 methodsPush.c, 1.26, 1.27 privatePush.h, 1.16, 1.17 push.h, 1.44, 1.45 setup.c, 1.15, 1.16
|
From: Gordon K. <kin...@us...> - 2006-10-25 12:19:26
|
Update of /cvsroot/teem/teem/src/push In directory sc8-pr-cvs10.sourceforge.net:/tmp/cvs-serv10340 Modified Files: action.c binning.c corePush.c forces.c methodsPush.c privatePush.h push.h setup.c Log Message: massively reworked: 1) tractlets gone, 2) asynchronous (Gauss-Seidel) updates, 3) potential energy gradient descent. Also two optimizations: only periodically gageProbe() field to get new values, and only periodically update list of active neighbors Index: push.h =================================================================== RCS file: /cvsroot/teem/teem/src/push/push.h,v retrieving revision 1.44 retrieving revision 1.45 diff -C2 -d -r1.44 -r1.45 *** push.h 22 Oct 2006 12:22:15 -0000 1.44 --- push.h 25 Oct 2006 12:19:19 -0000 1.45 *************** *** 51,125 **** /* ! ******** pointPoint ** ** information about a point in the simulation. There are really two ** kinds of information here: "pos", "vel", and "frc" pertain to the simulation ! ** of point dynamics, while "ten", "aniso", "inv", and "cnt" are properties of ! ** the field sampled at the point. "tan" and "nor" are only meaningful in ! ** tractlets. ! ** ! ** HEY: with the addition of grav, gravNot, seedThresh, GK wonders why so much ! ** information has to be stored per point, and not in the task ... ! ** ! ** NB: there is no constructor for this, nor does there really need to be */ typedef struct pushPoint_t { double pos[3], /* position in world space */ ! vel[3], /* velocity */ ! frc[3], /* force accumulator for current iteration */ ten[7], /* tensor here */ - aniso, /* value of Cl1 (Westin ISMRM '97) */ inv[7], /* inverse of tensor */ cnt[3], /* mask's containment gradient */ ! tan[3], /* tangent: unit direction through me */ ! nor[3], /* change in tangent, normalized */ ! grav[3], gravNot[2][3], /* gravity stuff */ seedThresh; /* seed thresh */ } pushPoint; /* - ******** pushThing struct - ** - ** Represents both single points, and tractlets, as follows: - ** - ** for single points: "point" tells the whole story of the point, - ** but point.{tan,nor} is meaningless. For the sake of easily computing all - ** pair-wise point interactions between things, "vertNum" is 1, and - ** "vert" points to "point". "len" is 0. - ** - ** for tractlets: the "pos", "vel", "frc" fields of "point" summarize the - ** the dynamics of the entire tractlet, while the field attributes - ** ("ten", "inv", "cnt") pertain exactly to the seed point. For example, - ** a particular tensor anisotropy in "point.ten" may have resulted in this - ** thing turning from a point into a tractlet, or vice versa. "vertNum" is - ** the number of tractlet vertices; "vert" is the array of them. The - ** only field of the vertex points that is not meaningful is "vel": the - ** tractlet velocity is "point.vel" - */ - typedef struct pushThing_t { - int ttaagg; - pushPoint point; /* information about single point, or a - seed point, hard to say exactly */ - unsigned int vertNum; /* 1 for single point, else length of vert[] */ - pushPoint *vert; /* dyn. alloc. array of tractlet vertices - (*not* pointers to pushPoints), or, just - the address of "point" for single point */ - unsigned int seedIdx; /* which of the vertices is the seed point */ - double len; /* 0 for point, else (world-space) length of - tractlet */ - } pushThing; - - /* ******** pushBin ** ! ** the data structure for doing spatial binning. This really serves two ! ** related purposes- to bin *all* the points in the simulation-- both ! ** single points as well as vertices of tractlets-- and also to bin the ! ** all "things". The binning of points is used for the force calculation ! ** stage, and the binning of things is used for the update stage. ** ! ** The bins are the *only* way to access the things in the simulation, so ! ** bins own the things they contain. However, because things point to ! ** points, bins do not own the points they contain. */ typedef struct pushBin_t { --- 51,82 ---- /* ! ******** pushPoint ** ** information about a point in the simulation. There are really two ** kinds of information here: "pos", "vel", and "frc" pertain to the simulation ! ** of point dynamics, while "ten", "inv", and "cnt" are properties of ! ** the field sampled at the point. */ typedef struct pushPoint_t { + unsigned int ttaagg; double pos[3], /* position in world space */ ! enr, /* energy accumulator (current iteration) */ ! frc[3], /* force accumulator (current iteration) */ ten[7], /* tensor here */ inv[7], /* inverse of tensor */ cnt[3], /* mask's containment gradient */ ! grav, gravGrad[3], /* gravity stuff */ seedThresh; /* seed thresh */ + struct pushPoint_t **neigh; + unsigned int neighNum; + airArray *neighArr; } pushPoint; /* ******** pushBin ** ! ** the data structure for doing spatial binning. ** ! ** in tractlet-less push, bins do own the points they contain */ typedef struct pushBin_t { *************** *** 127,133 **** pushPoint **point; /* dyn. alloc. array of point pointers */ airArray *pointArr; /* airArray around point and pointNum */ - unsigned int thingNum; /* # of things in this bin */ - pushThing **thing; /* dyn. alloc. array of thing pointers */ - airArray *thingArr; /* airArray around thing and thingNum */ struct pushBin_t **neighbor; /* pre-computed NULL-terminated list of all neighboring bins, including myself */ --- 84,87 ---- *************** *** 142,154 **** struct pushContext_t *pctx; /* parent's context */ gageContext *gctx; /* result of gageContextCopy(pctx->gctx) */ ! const double *tenAns, *invAns, *cntAns, ! *gravAns, *gravNotAns[2], /* results of gage probing */ ! *seedThreshAns; /* seed threshold answer */ ! tenFiberContext *fctx; /* result of tenFiberContextCopy(pctx->fctx) */ airThread *thread; /* my thread */ unsigned int threadIdx, /* which thread am I */ ! thingNum; /* # things I let live this iteration */ ! double sumVel, /* sum of velocities of things in my bins */ ! *vertBuff; /* buffer for tractlet vertices */ void *returnPtr; /* for airThreadJoin */ } pushTask; --- 96,107 ---- struct pushContext_t *pctx; /* parent's context */ gageContext *gctx; /* result of gageContextCopy(pctx->gctx) */ ! const double *tenAns, /* results of gage probing */ ! *invAns, *cntAns, ! *gravAns, *gravGradAns, ! *seedThreshAns; airThread *thread; /* my thread */ unsigned int threadIdx, /* which thread am I */ ! pointNum; /* # points I let live this iteration */ ! double energySum; /* sum of energies of points I processed */ void *returnPtr; /* for airThreadJoin */ } pushTask; *************** *** 175,184 **** ** ** the functions which determine inter-point forces */ typedef struct { char name[AIR_STRLEN_SMALL]; unsigned int parmNum; ! double (*energy)(double dist, const double parm[PUSH_ENERGY_PARM_NUM]); ! double (*force)(double dist, const double parm[PUSH_ENERGY_PARM_NUM]); double (*support)(const double parm[PUSH_ENERGY_PARM_NUM]); } pushEnergy; --- 128,140 ---- ** ** the functions which determine inter-point forces + ** + ** NOTE: the eval() function probably does NOT check to see it was passed + ** non-NULL pointers into which to store energy and force */ typedef struct { char name[AIR_STRLEN_SMALL]; unsigned int parmNum; ! void (*eval)(double *energy, double *force, ! double dist, const double parm[PUSH_ENERGY_PARM_NUM]); double (*support)(const double parm[PUSH_ENERGY_PARM_NUM]); } pushEnergy; *************** *** 196,209 **** typedef struct pushContext_t { /* INPUT ----------------------------- */ Nrrd *nin, /* 3D image of 3D masked tensors, though it may only be a single slice */ ! *npos, /* positions to start with ! (overrides thingNum) */ ! *nstn; /* start/nums for tractlets in npos */ double step, /* time step in integration */ scale, /* scaling from tensor to glyph size */ wall, /* spring constant of walls */ cntScl, /* magnitude of containment gradient */ - bigTrace, /* a last minute hack */ minMeanVel; /* stop if mean velocity drops below this */ int detReject, /* determinant-based rejection at init */ --- 152,164 ---- typedef struct pushContext_t { /* INPUT ----------------------------- */ + unsigned int pointNum; /* number points to start simulation w/ */ Nrrd *nin, /* 3D image of 3D masked tensors, though it may only be a single slice */ ! *npos; /* positions to start with ! (overrides pointNum) */ double step, /* time step in integration */ scale, /* scaling from tensor to glyph size */ wall, /* spring constant of walls */ cntScl, /* magnitude of containment gradient */ minMeanVel; /* stop if mean velocity drops below this */ int detReject, /* determinant-based rejection at init */ *************** *** 211,237 **** verbose; /* blah blah blah */ unsigned int seedRNG, /* seed value for random number generator */ - thingNum, /* number things to start simulation w/ */ threadNum, /* number of threads to use */ maxIter, /* if non-zero, max number of iterations */ snap; /* if non-zero, interval between iterations at which output snapshots are saved */ ! int gravItem, /* tenGage item (vector) for gravity */ ! gravNotItem[2]; /* for constraining gravity */ ! double gravScl; /* sign and magnitude of gravity's effect */ ! int seedThreshItem, /* item for constraining random seeding */ seedThreshSign; /* +1: need val > thresh; -1: opposite */ double seedThresh; /* threshold for seed constraint */ ! const pushEnergySpec *ensp; /* potential energy function to use */ ! ! int tltUse, /* enable tractlets */ ! tltFrenet; /* use Frenet frames for tractlet forces */ ! unsigned int tltStepNum; /* max # points on each tractlet half */ ! double tltThresh, ! tltSoft, tltStep; /* tractlet formation parameters */ int binSingle; /* disable binning (for debugging) */ ! unsigned int binIncr; /* increment for per-bin thing airArray */ NrrdKernelSpec *ksp00, /* for sampling tensor field */ --- 166,196 ---- verbose; /* blah blah blah */ unsigned int seedRNG, /* seed value for random number generator */ threadNum, /* number of threads to use */ + iterNeighbor, /* recompute list of active neighbors on + iterations that are a multiple of this, + or, 0 to disable this optimization */ + iterProbe, /* only do field probing (which is slow) on + iterations that are a multiple of this, + or, 0 to disable this optimization */ maxIter, /* if non-zero, max number of iterations */ snap; /* if non-zero, interval between iterations at which output snapshots are saved */ ! int gravItem, /* tenGage item (scalar) for "height" ! potential energy associated w/ gravity */ ! gravGradItem; /* tenGage item (vector) for gravity */ ! double gravScl, /* sign and magnitude of gravity's effect: ! when this is positive, higher values of ! gravItem have higher potential energy */ ! gravZero; /* the height that corresponds to zero ! potential energy from gravity */ ! int seedThreshItem, /* item for constraining random seeding */ seedThreshSign; /* +1: need val > thresh; -1: opposite */ double seedThresh; /* threshold for seed constraint */ ! pushEnergySpec *ensp; /* potential energy function to use */ int binSingle; /* disable binning (for debugging) */ ! unsigned int binIncr; /* increment for per-bin airArray */ NrrdKernelSpec *ksp00, /* for sampling tensor field */ *************** *** 241,244 **** --- 200,204 ---- /* INTERNAL -------------------------- */ + unsigned int ttaagg; /* next value for per-point ID */ Nrrd *nten, /* 3D image of 3D masked tensors */ *ninv, /* pre-computed inverse of nten */ *************** *** 246,250 **** gageContext *gctx; /* gage context around nten, ninv, nmask */ gagePerVolume *tpvl, *ipvl; /* gage pervolumes around nten and ninv */ - tenFiberContext *fctx; /* tenFiber context around nten */ int finished; /* used to signal all threads to return */ unsigned int dimIn, /* dim (2 or 3) of input, meaning whether --- 206,209 ---- *************** *** 265,269 **** maxEval, meanEval, /* max and mean principal eval in field */ maxDet, ! meanVel; /* latest mean velocity of particles */ pushTask **task; /* dynamically allocated array of tasks */ airThreadBarrier *iterBarrierA; /* barriers between iterations */ --- 224,228 ---- maxEval, meanEval, /* max and mean principal eval in field */ maxDet, ! energySum; /* potential energy of entire particles */ pushTask **task; /* dynamically allocated array of tasks */ airThreadBarrier *iterBarrierA; /* barriers between iterations */ *************** *** 283,288 **** /* methodsPush.c */ ! PUSH_EXPORT pushThing *pushThingNew(unsigned int vertNum); ! PUSH_EXPORT pushThing *pushThingNix(pushThing *thg); PUSH_EXPORT pushContext *pushContextNew(void); PUSH_EXPORT pushContext *pushContextNix(pushContext *pctx); --- 242,247 ---- /* methodsPush.c */ ! PUSH_EXPORT pushPoint *pushPointNew(pushContext *pctx); ! PUSH_EXPORT pushPoint *pushPointNix(pushPoint *pnt); PUSH_EXPORT pushContext *pushContextNew(void); PUSH_EXPORT pushContext *pushContextNix(pushContext *pctx); *************** *** 298,301 **** --- 257,263 ---- PUSH_EXPORT const pushEnergy *const pushEnergyAll[PUSH_ENERGY_TYPE_MAX+1]; PUSH_EXPORT pushEnergySpec *pushEnergySpecNew(); + PUSH_EXPORT void pushEnergySpecSet(pushEnergySpec *ensp, + const pushEnergy *energy, + const double parm[PUSH_ENERGY_PARM_NUM]); PUSH_EXPORT pushEnergySpec *pushEnergySpecNix(pushEnergySpec *ensp); PUSH_EXPORT int pushEnergySpecParse(pushEnergySpec *ensp, const char *str); *************** *** 311,329 **** PUSH_EXPORT void pushBinInit(pushBin *bin, unsigned int incr); PUSH_EXPORT void pushBinDone(pushBin *bin); - PUSH_EXPORT int pushBinThingAdd(pushContext *pctx, pushThing *thing); PUSH_EXPORT int pushBinPointAdd(pushContext *pctx, pushPoint *point); PUSH_EXPORT void pushBinAllNeighborSet(pushContext *pctx); PUSH_EXPORT int pushRebin(pushContext *pctx); - /* setup.c */ - PUSH_EXPORT int pushTaskFiberReSetup(pushContext *pctx, - double tltThresh, - double tltSoft, - double tltStep, - unsigned int tltStepNum); - /* action.c */ PUSH_EXPORT int pushBinProcess(pushTask *task, unsigned int myBinIdx); ! PUSH_EXPORT int pushOutputGet(Nrrd *nPosOut, Nrrd *nTenOut, Nrrd *nStnOut, pushContext *pctx); --- 273,283 ---- PUSH_EXPORT void pushBinInit(pushBin *bin, unsigned int incr); PUSH_EXPORT void pushBinDone(pushBin *bin); PUSH_EXPORT int pushBinPointAdd(pushContext *pctx, pushPoint *point); PUSH_EXPORT void pushBinAllNeighborSet(pushContext *pctx); PUSH_EXPORT int pushRebin(pushContext *pctx); /* action.c */ PUSH_EXPORT int pushBinProcess(pushTask *task, unsigned int myBinIdx); ! PUSH_EXPORT int pushOutputGet(Nrrd *nPosOut, Nrrd *nTenOut, pushContext *pctx); Index: privatePush.h =================================================================== RCS file: /cvsroot/teem/teem/src/push/privatePush.h,v retrieving revision 1.16 retrieving revision 1.17 diff -C2 -d -r1.16 -r1.17 *** privatePush.h 22 Oct 2006 12:22:15 -0000 1.16 --- privatePush.h 25 Oct 2006 12:19:19 -0000 1.17 *************** *** 30,35 **** /* binning.c */ extern pushBin *_pushBinLocate(pushContext *pctx, double *pos); - extern int _pushBinPointNullify(pushContext *pctx, - pushBin *oldBin, pushPoint *point); extern void _pushBinPointAdd(pushContext *pctx, pushBin *bin, pushPoint *point); --- 30,33 ---- *************** *** 40,57 **** extern int _pushTensorFieldSetup(pushContext *pctx); extern int _pushGageSetup(pushContext *pctx); - extern int _pushFiberSetup(pushContext *pctx); extern int _pushTaskSetup(pushContext *pctx); extern int _pushBinSetup(pushContext *pctx); ! extern int _pushThingSetup(pushContext *pctx); /* action.c */ ! extern double _pushThingPointCharge(pushContext *pctx, pushThing *thg); ! extern int _pushForceSample(pushContext *pctx, ! unsigned int sx, unsigned int sy); ! extern int _pushProbe(pushTask *task, pushPoint *point); ! extern int _pushInputProcess(pushContext *pctx); ! extern void _pushInitialize(pushContext *pctx); ! extern int _pushForce(pushTask *task, int bin, const double *parm); ! extern int _pushUpdate(pushTask *task, int bin, const double *parm); #ifdef __cplusplus --- 38,47 ---- extern int _pushTensorFieldSetup(pushContext *pctx); extern int _pushGageSetup(pushContext *pctx); extern int _pushTaskSetup(pushContext *pctx); extern int _pushBinSetup(pushContext *pctx); ! extern int _pushPointSetup(pushContext *pctx); /* action.c */ ! extern void _pushProbe(pushTask *task, pushPoint *point); #ifdef __cplusplus Index: setup.c =================================================================== RCS file: /cvsroot/teem/teem/src/push/setup.c,v retrieving revision 1.15 retrieving revision 1.16 diff -C2 -d -r1.15 -r1.16 *** setup.c 22 Oct 2006 12:22:15 -0000 1.15 --- setup.c 25 Oct 2006 12:19:19 -0000 1.16 *************** *** 138,147 **** if (tenGageUnknown != pctx->gravItem) { if (!E) E |= gageQueryItemOn(pctx->gctx, pctx->tpvl, pctx->gravItem); ! } ! if (tenGageUnknown != pctx->gravNotItem[0]) { ! if (!E) E |= gageQueryItemOn(pctx->gctx, pctx->tpvl, pctx->gravNotItem[0]); ! } ! if (tenGageUnknown != pctx->gravNotItem[1]) { ! if (!E) E |= gageQueryItemOn(pctx->gctx, pctx->tpvl, pctx->gravNotItem[1]); } /* set up tensor inverse probing */ --- 138,142 ---- if (tenGageUnknown != pctx->gravItem) { if (!E) E |= gageQueryItemOn(pctx->gctx, pctx->tpvl, pctx->gravItem); ! if (!E) E |= gageQueryItemOn(pctx->gctx, pctx->tpvl, pctx->gravGradItem); } /* set up tensor inverse probing */ *************** *** 175,218 **** } - /* - ** _pushFiberSetup sets: - **** pctx->fctx - */ - int - _pushFiberSetup(pushContext *pctx) { - char me[]="_pushFiberSetup", err[BIFF_STRLEN]; - int E; - - pctx->fctx = tenFiberContextNew(pctx->nten); - if (!pctx->fctx) { - sprintf(err, "%s: couldn't create fiber context", me); - biffMove(PUSH, err, TEN); return 1; - } - E = AIR_FALSE; - if (!E) E |= tenFiberStopSet(pctx->fctx, tenFiberStopNumSteps, - pctx->tltStepNum); - if (!E) E |= tenFiberStopSet(pctx->fctx, tenFiberStopAniso, - tenAniso_Cl1, - pctx->tltThresh - pctx->tltSoft); - if (!E) E |= tenFiberTypeSet(pctx->fctx, tenFiberTypeEvec1); - if (!E) E |= tenFiberKernelSet(pctx->fctx, - pctx->ksp00->kernel, pctx->ksp00->parm); - /* if (!E) E |= tenFiberIntgSet(pctx->fctx, tenFiberIntgRK4); */ - if (!E) E |= tenFiberIntgSet(pctx->fctx, tenFiberIntgMidpoint); - /* if (!E) E |= tenFiberIntgSet(pctx->fctx, tenFiberIntgEuler); */ - if (!E) E |= tenFiberParmSet(pctx->fctx, tenFiberParmStepSize, - pctx->tltStep/pctx->tltStepNum); - if (!E) E |= tenFiberAnisoSpeedSet(pctx->fctx, tenAniso_Cl1, - 1 /* lerp */ , - pctx->tltThresh /* thresh */, - pctx->tltSoft); - if (!E) E |= tenFiberUpdate(pctx->fctx); - if (E) { - sprintf(err, "%s: trouble setting up fiber context", me); - biffMove(PUSH, err, TEN); return 1; - } - return 0; - } - pushTask * _pushTaskNew(pushContext *pctx, int threadIdx) { --- 170,173 ---- *************** *** 223,227 **** task->pctx = pctx; task->gctx = gageContextCopy(pctx->gctx); - task->fctx = tenFiberContextCopy(pctx->fctx); /* ** HEY: its a limitation in gage that we have to know a priori --- 178,181 ---- *************** *** 234,261 **** task->cntAns = gageAnswerPointer(task->gctx, task->gctx->pvl[2], gageSclGradVec); ! task->gravAns = (tenGageUnknown != task->pctx->gravItem ! ? gageAnswerPointer(task->gctx, task->gctx->pvl[0], ! task->pctx->gravItem) ! : NULL); ! task->gravNotAns[0] = (tenGageUnknown != task->pctx->gravNotItem[0] ! ? gageAnswerPointer(task->gctx, task->gctx->pvl[0], ! task->pctx->gravNotItem[0]) ! : NULL); ! task->gravNotAns[1] = (tenGageUnknown != task->pctx->gravNotItem[1] ! ? gageAnswerPointer(task->gctx, task->gctx->pvl[0], ! task->pctx->gravNotItem[1]) ! : NULL); ! task->seedThreshAns = (tenGageUnknown != task->pctx->seedThreshItem ! ? gageAnswerPointer(task->gctx, task->gctx->pvl[0], ! task->pctx->seedThreshItem) ! : NULL); if (threadIdx) { task->thread = airThreadNew(); } task->threadIdx = threadIdx; ! task->thingNum = 0; ! task->sumVel = 0; ! task->vertBuff = (double*)calloc(3*(1 + 2*pctx->tltStepNum), ! sizeof(double)); task->returnPtr = NULL; } --- 188,212 ---- task->cntAns = gageAnswerPointer(task->gctx, task->gctx->pvl[2], gageSclGradVec); ! if (tenGageUnknown != task->pctx->gravItem) { ! task->gravAns = gageAnswerPointer(task->gctx, task->gctx->pvl[0], ! task->pctx->gravItem); ! task->gravGradAns = gageAnswerPointer(task->gctx, task->gctx->pvl[0], ! task->pctx->gravGradItem); ! } else { ! task->gravAns = NULL; ! task->gravGradAns = NULL; ! } ! if (tenGageUnknown != task->pctx->seedThreshItem) { ! task->seedThreshAns = gageAnswerPointer(task->gctx, task->gctx->pvl[0], ! task->pctx->seedThreshItem); ! } else { ! task->seedThreshAns = NULL; ! } if (threadIdx) { task->thread = airThreadNew(); } task->threadIdx = threadIdx; ! task->pointNum = 0; ! task->energySum = 0; task->returnPtr = NULL; } *************** *** 268,276 **** if (task) { task->gctx = gageContextNix(task->gctx); - task->fctx = tenFiberContextNix(task->fctx); if (task->threadIdx) { task->thread = airThreadNix(task->thread); } - task->vertBuff = (double *)airFree(task->vertBuff); airFree(task); } --- 219,225 ---- *************** *** 294,297 **** --- 243,249 ---- } for (tidx=0; tidx<pctx->threadNum; tidx++) { + if (pctx->verbose) { + fprintf(stderr, "%s: creating task %u/%u\n", me, tidx, pctx->threadNum); + } pctx->task[tidx] = _pushTaskNew(pctx, tidx); if (!(pctx->task[tidx])) { *************** *** 383,442 **** /* ! ** THIS IS A COMPLETE HACK!!! ! */ ! int ! pushTaskFiberReSetup(pushContext *pctx, ! double tltThresh, double tltSoft, double tltStep, ! unsigned int tltStepNum) { ! char me[]="pushTaskFiberReSetup", err[BIFF_STRLEN]; ! tenFiberContext *fctx; ! unsigned int taskIdx; ! int E = 0; ! ! fprintf(stderr, "!%s: %d\n", me, __LINE__); ! if (!pctx->task) { ! fprintf(stderr, "!%s: %d bailing\n", me, __LINE__); ! return 0; ! } ! pctx->tltStepNum = tltStepNum; ! pctx->tltThresh = tltThresh; ! pctx->tltSoft = tltSoft; ! pctx->tltStep = tltStep; ! ! fprintf(stderr, "!%s: %d\n", me, __LINE__); ! for (taskIdx=0; pctx->threadNum; taskIdx++) { ! fctx = pctx->task[taskIdx]->fctx; ! fprintf(stderr, "!%s: %d %p\n", me, __LINE__, fctx); ! if (!E) E |= tenFiberStopSet(fctx, tenFiberStopNumSteps, ! pctx->tltStepNum); ! fprintf(stderr, "!%s: %d %p\n", me, __LINE__, fctx); ! if (!E) E |= tenFiberStopSet(fctx, tenFiberStopAniso, ! tenAniso_Cl1, ! pctx->tltThresh - pctx->tltSoft); ! fprintf(stderr, "!%s: %d %p\n", me, __LINE__, fctx); ! if (!E) E |= tenFiberParmSet(fctx, tenFiberParmStepSize, ! pctx->tltStep/pctx->tltStepNum); ! fprintf(stderr, "!%s: %d %p\n", me, __LINE__, fctx); ! if (!E) E |= tenFiberAnisoSpeedSet(fctx, tenAniso_Cl1, ! 1 /* lerp */ , ! pctx->tltThresh /* thresh */, ! pctx->tltSoft); ! fprintf(stderr, "!%s: %d %p\n", me, __LINE__, fctx); ! if (!E) E |= tenFiberUpdate(fctx); ! fprintf(stderr, "!%s: %d %p\n", me, __LINE__, fctx); ! if (E) { ! sprintf(err, "%s: trouble resetting task %u fiber context", me, taskIdx); ! biffMove(PUSH, err, TEN); return 1; ! } ! fprintf(stderr, "!%s: %d %p\n", me, __LINE__, fctx); ! } ! fprintf(stderr, "!%s: %d\n", me, __LINE__); ! return 1; ! } ! ! ! /* ! ** _pushThingSetup sets: ! **** pctx->thingNum (in case pctx->nstn and/or pctx->npos) ** ** This is only called by the master thread --- 335,340 ---- /* ! ** _pushPointSetup sets: ! **** pctx->pointNum (in case pctx->npos) ** ** This is only called by the master thread *************** *** 446,462 **** */ int ! _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; ! pctx->thingNum = (pctx->nstn ! ? pctx->nstn->axis[1].size ! : (pctx->npos ! ? pctx->npos->axis[1].size ! : pctx->thingNum)); 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 --- 344,362 ---- */ int ! _pushPointSetup(pushContext *pctx) { ! char me[]="_pushPointSetup", err[BIFF_STRLEN]; double (*lup)(const void *v, size_t I), maxDet; ! unsigned int pointIdx; ! pushPoint *point; ! /* ! double posIdxHack[2][4] = { ! {49.99999, 50, 0, 1}, ! {50, 50, 0, 1}}; ! */ ! pctx->pointNum = (pctx->npos ! ? pctx->npos->axis[1].size ! : pctx->pointNum); lup = pctx->npos ? nrrdDLookup[pctx->npos->type] : NULL; fprintf(stderr, "!%s: initilizing/seeding ... \n", me); /* HEY: we end up keeping a local copy of maxDet because convolution *************** *** 465,505 **** detReject should probably *not* be enabled... */ maxDet = pctx->maxDet; ! for (thingIdx=0; thingIdx<pctx->thingNum; thingIdx++) { double detProbe; /* ! fprintf(stderr, "!%s: thingIdx = %u/%u\n", me, thingIdx, pctx->thingNum); */ ! if (pctx->nstn) { ! baseIdx = stn[0 + 3*thingIdx]; ! thing = pushThingNew(stn[1 + 3*thingIdx]); ! for (pointIdx=0; pointIdx<thing->vertNum; pointIdx++) { ! ELL_3V_SET(thing->vert[pointIdx].pos, ! lup(pctx->npos->data, 0 + 3*(pointIdx + baseIdx)), ! lup(pctx->npos->data, 1 + 3*(pointIdx + baseIdx)), ! lup(pctx->npos->data, 2 + 3*(pointIdx + baseIdx))); ! _pushProbe(pctx->task[0], thing->vert + pointIdx); ! } ! thing->seedIdx = stn[2 + 3*thingIdx]; ! if (1 < thing->vertNum) { ! /* info about seedpoint has to be set separately */ ! ELL_3V_SET(thing->point.pos, ! lup(pctx->npos->data, 0 + 3*(thing->seedIdx + baseIdx)), ! lup(pctx->npos->data, 1 + 3*(thing->seedIdx + baseIdx)), ! lup(pctx->npos->data, 2 + 3*(thing->seedIdx + baseIdx))); ! _pushProbe(pctx->task[0], &(thing->point)); ! } /* ! fprintf(stderr, "!%s: thingNum(%d) = %d\n", "_pushThingSetup", ! thingIdx, thing->vertNum); */ - } else if (pctx->npos) { - thing = pushThingNew(1); - ELL_3V_SET(thing->vert[0].pos, - lup(pctx->npos->data, 0 + 3*thingIdx), - lup(pctx->npos->data, 1 + 3*thingIdx), - lup(pctx->npos->data, 2 + 3*thingIdx)); - _pushProbe(pctx->task[0], thing->vert + 0); - } else { - thing = pushThingNew(1); do { double posIdx[4], posWorld[4]; --- 365,387 ---- detReject should probably *not* be enabled... */ maxDet = pctx->maxDet; ! for (pointIdx=0; pointIdx<pctx->pointNum; pointIdx++) { double detProbe; /* ! fprintf(stderr, "!%s: pointIdx = %u/%u\n", me, pointIdx, pctx->pointNum); */ ! point = pushPointNew(pctx); ! if (pctx->npos) { ! ELL_3V_SET(point->pos, ! lup(pctx->npos->data, 0 + 3*pointIdx), ! lup(pctx->npos->data, 1 + 3*pointIdx), ! lup(pctx->npos->data, 2 + 3*pointIdx)); ! _pushProbe(pctx->task[0], point); ! } else { /* ! double posWorld[4]; ! ELL_4MV_MUL(posWorld, pctx->gctx->shape->ItoW, posIdxHack[pointIdx]); ! ELL_34V_HOMOG(point->pos, posWorld); ! _pushProbe(pctx->task[0], point); */ do { double posIdx[4], posWorld[4]; *************** *** 515,519 **** } ELL_4MV_MUL(posWorld, pctx->gctx->shape->ItoW, posIdx); ! ELL_34V_HOMOG(thing->vert[0].pos, posWorld); /* fprintf(stderr, "%s: posIdx = %g %g %g --> posWorld = %g %g %g " --- 397,401 ---- } ELL_4MV_MUL(posWorld, pctx->gctx->shape->ItoW, posIdx); ! ELL_34V_HOMOG(point->pos, posWorld); /* fprintf(stderr, "%s: posIdx = %g %g %g --> posWorld = %g %g %g " *************** *** 521,529 **** posIdx[0], posIdx[1], posIdx[2], posWorld[0], posWorld[1], posWorld[2], ! thing->vert[0].pos[0], thing->vert[0].pos[1], ! thing->vert[0].pos[2]); */ ! _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, --- 403,410 ---- posIdx[0], posIdx[1], posIdx[2], posWorld[0], posWorld[1], posWorld[2], ! point->pos[0], point->pos[1], point->pos[2]); */ ! _pushProbe(pctx->task[0], point); ! detProbe = TEN_T_DET(point->ten); maxDet = AIR_MAX(maxDet, detProbe); /* assuming that we're not using some very blurring kernel, *************** *** 533,543 **** seedThresh enabled! */ /* ! fprintf(stderr, "!%s: ten[0] = %g\n", me, thing->vert[0].ten[0]); */ ! /* we OR together all the things that would make us REJECT this last sample */ ! } while (thing->vert[0].ten[0] < 0.5 || (tenGageUnknown != pctx->seedThreshItem ! && ((pctx->seedThresh - thing->vert[0].seedThresh) *pctx->seedThreshSign > 0) ) --- 414,424 ---- seedThresh enabled! */ /* ! fprintf(stderr, "!%s: ten[0] = %g\n", me, point->ten[0]); */ ! /* we OR together all the tests that would make us REJECT this last sample */ ! } while (point->ten[0] < 0.5 || (tenGageUnknown != pctx->seedThreshItem ! && ((pctx->seedThresh - point->seedThresh) *pctx->seedThreshSign > 0) ) *************** *** 545,588 **** ); } ! for (pointIdx=0; pointIdx<thing->vertNum; pointIdx++) { ! if (pushBinPointAdd(pctx, thing->vert + pointIdx)) { ! sprintf(err, "%s: trouble binning vert %d of thing %d", me, ! pointIdx, thingIdx); ! biffAdd(PUSH, err); return 1; ! } ! ELL_3V_SET(thing->vert[pointIdx].vel, 0, 0, 0); ! } ! if (pushBinThingAdd(pctx, thing)) { ! sprintf(err, "%s: trouble thing %d", me, thingIdx); biffAdd(PUSH, err); return 1; } } fprintf(stderr, "!%s: ... seeding DONE\n", me); - /* - { - Nrrd *nten, *npos, *nstn; - char me[]="dammit", err[BIFF_STRLEN], poutS[AIR_STRLEN_MED], - toutS[AIR_STRLEN_MED], soutS[AIR_STRLEN_MED]; - nten = nrrdNew(); - npos = nrrdNew(); - nstn = nrrdNew(); - sprintf(poutS, "snap-pre.%06d.pos.nrrd", -1); - sprintf(toutS, "snap-pre.%06d.ten.nrrd", -1); - sprintf(soutS, "snap-pre.%06d.stn.nrrd", -1); - if (pushOutputGet(npos, nten, nstn, pctx)) { - sprintf(err, "%s: couldn't get snapshot for iter %d", me, -1); - biffAdd(PUSH, err); return 1; - } - if (nrrdSave(poutS, npos, NULL) - || nrrdSave(toutS, nten, NULL) - || nrrdSave(soutS, nstn, NULL)) { - sprintf(err, "%s: couldn't save snapshot for iter %d", me, -1); - biffMove(PUSH, err, NRRD); return 1; - } - nten = nrrdNuke(nten); - npos = nrrdNuke(npos); - nstn = nrrdNuke(nstn); - } - */ return 0; } --- 426,435 ---- ); } ! if (pushBinPointAdd(pctx, point)) { ! sprintf(err, "%s: trouble binning point %u", me, point->ttaagg); biffAdd(PUSH, err); return 1; } } fprintf(stderr, "!%s: ... seeding DONE\n", me); return 0; } Index: forces.c =================================================================== RCS file: /cvsroot/teem/teem/src/push/forces.c,v retrieving revision 1.21 retrieving revision 1.22 diff -C2 -d -r1.21 -r1.22 *** forces.c 22 Oct 2006 12:22:15 -0000 1.21 --- forces.c 25 Oct 2006 12:19:19 -0000 1.22 *************** *** 66,87 **** ** ---------------------------------------------------------------- */ ! double ! _pushEnergyUnknownEnergy(double dist, const double *parm) { ! char me[]="_pushEnergyUnknownEnergy"; ! ! AIR_UNUSED(dist); ! AIR_UNUSED(parm); ! fprintf(stderr, "%s: this is not good.\n", me); ! return AIR_NAN; ! } ! ! double ! _pushEnergyUnknownForce(double dist, const double *parm) { ! char me[]="_pushEnergyUnknownForce"; AIR_UNUSED(dist); AIR_UNUSED(parm); ! fprintf(stderr, "%s: this is not good.\n", me); ! return AIR_NAN; } --- 66,80 ---- ** ---------------------------------------------------------------- */ ! void ! _pushEnergyUnknownEval(double *enr, double *frc, ! double dist, const double *parm) { ! char me[]="_pushEnergyUnknownEval"; AIR_UNUSED(dist); AIR_UNUSED(parm); ! *enr = AIR_NAN; ! *frc = AIR_NAN; ! fprintf(stderr, "%s: ERROR- using unknown energy.\n", me); ! return; } *************** *** 91,95 **** AIR_UNUSED(parm); ! fprintf(stderr, "%s: this is not good.\n", me); return AIR_NAN; } --- 84,88 ---- AIR_UNUSED(parm); ! fprintf(stderr, "%s: ERROR- using unknown energy.\n", me); return AIR_NAN; } *************** *** 99,104 **** "unknown", 0, ! _pushEnergyUnknownEnergy, ! _pushEnergyUnknownForce, _pushEnergyUnknownSupport }; --- 92,96 ---- "unknown", 0, ! _pushEnergyUnknownEval, _pushEnergyUnknownSupport }; *************** *** 110,151 **** ** ---------------------------------------------------------------- ** 1 parms: ! ** 0: pull distance */ ! double ! _pushEnergySpringEnergy(double dist, const double *parm) { ! /* char me[]="_pushEnergySpringEnergy"; */ ! double xx, ret, pull; ! ! pull = parm[0]; ! xx = dist - 1.0; ! if (xx > pull) { ! ret = 0; ! } else if (xx > 0) { ! ret = xx*xx*(xx*xx/(4*pull*pull) - 2*xx/(3*pull) + 1/2); ! } else { ! ret = xx*xx/2; ! } ! /* ! if (!AIR_EXISTS(ret)) { ! fprintf(stderr, "!%s: dist=%g, pull=%g, blah=%d --> ret=%g\n", ! me, dist, pull, blah, ret); ! } ! */ ! return ret; ! } ! ! double ! _pushEnergySpringForce(double dist, const double *parm) { ! /* char me[]="_pushEnergySpringForce"; */ ! double xx, ret, pull; pull = parm[0]; xx = dist - 1.0; if (xx > pull) { ! ret = 0; } else if (xx > 0) { ! ret = xx*(xx*xx/(pull*pull) - 2*xx/pull + 1); } else { ! ret = xx; } /* --- 102,126 ---- ** ---------------------------------------------------------------- ** 1 parms: ! ** parm[0]: width of pull region (beyond 1.0) ! ** ! ** learned: "1/2" is not 0.5 !!!!! */ ! void ! _pushEnergySpringEval(double *enr, double *frc, ! double dist, const double *parm) { ! /* char me[]="_pushEnergySpringEval"; */ ! double xx, pull; pull = parm[0]; xx = dist - 1.0; if (xx > pull) { ! *enr = 0; ! *frc = 0; } else if (xx > 0) { ! *enr = xx*xx*(xx*xx/(4*pull*pull) - 2*xx/(3*pull) + 1.0/2.0); ! *frc = xx*(xx*xx/(pull*pull) - 2*xx/pull + 1); } else { ! *enr = xx*xx/2; ! *frc = xx; } /* *************** *** 155,159 **** } */ ! return ret; } --- 130,134 ---- } */ ! return; } *************** *** 168,173 **** SPRING, 1, ! _pushEnergySpringEnergy, ! _pushEnergySpringForce, _pushEnergySpringSupport }; --- 143,147 ---- SPRING, 1, ! _pushEnergySpringEval, _pushEnergySpringSupport }; *************** *** 179,206 **** ** ---------------------------------------------------------------- ** 1 parms: ! ** (scale: distance to inflection point of force function) ! ** parm[0]: cut-off (as a multiple of standard dev) */ ! #define _DGAUSS(x, sig, cut) ( \ ! x >= sig*cut ? 0 \ ! : -exp(-x*x/(2.0*sig*sig))*x) ! #define SQRTTHREE 1.73205080756887729352 ! ! double ! _pushEnergyGaussEnergy(double dist, const double *parm) { ! double sig, cut; ! sig = 1.0/SQRTTHREE; ! cut = parm[0]; ! return AIR_CAST(double, _DGAUSS(dist, sig, cut)); ! } ! double ! _pushEnergyGaussForce(double dist, const double *parm) { ! double sig, cut; - sig = 1.0/SQRTTHREE; cut = parm[0]; ! return AIR_CAST(double, _DGAUSS(dist, sig, cut)); } --- 153,177 ---- ** ---------------------------------------------------------------- ** 1 parms: ! ** (distance to inflection point of force function is always 1.0) ! ** parm[0]: cut-off (as a multiple of standard dev (which is 1.0)) */ ! /* HEY: copied from teem/src/nrrd/kernel.c */ ! #define _GAUSS(x, sig, cut) ( \ ! x >= sig*cut ? 0 \ ! : exp(-x*x/(2.0*sig*sig))/(sig*2.50662827463100050241)) ! #define _DGAUSS(x, sig, cut) ( \ ! x >= sig*cut ? 0 \ ! : -exp(-x*x/(2.0*sig*sig))*x/(sig*sig*sig*2.50662827463100050241)) ! void ! _pushEnergyGaussEval(double *enr, double *frc, ! double dist, const double *parm) { ! double cut; cut = parm[0]; ! *enr = _GAUSS(dist, 1.0, cut); ! *frc = _DGAUSS(dist, 1.0, cut); ! return; } *************** *** 208,212 **** _pushEnergyGaussSupport(const double *parm) { ! return (1.0/SQRTTHREE)*parm[0]; } --- 179,183 ---- _pushEnergyGaussSupport(const double *parm) { ! return parm[0]; } *************** *** 215,220 **** GAUSS, 1, ! _pushEnergyGaussEnergy, ! _pushEnergyGaussForce, _pushEnergyGaussSupport }; --- 186,190 ---- GAUSS, 1, ! _pushEnergyGaussEval, _pushEnergyGaussSupport }; *************** *** 229,242 **** ** parm[0]: cut-off (as multiple of "1.0") */ ! double ! _pushEnergyCoulombEnergy(double dist, const double *parm) { ! ! return (dist > parm[0] ? 0 : 1.0/dist)); ! } ! ! double ! _pushEnergyCoulombForce(double dist, const double *parm) { ! return (dist > parm[0] ? 0 : -1.0/(dist*dist)); } --- 199,209 ---- ** parm[0]: cut-off (as multiple of "1.0") */ ! void ! _pushEnergyCoulombEval(double *enr, double *frc, ! double dist, const double *parm) { ! *enr = (dist > parm[0] ? 0 : 1.0/dist); ! *frc = (dist > parm[0] ? 0 : -1.0/(dist*dist)); ! return; } *************** *** 251,256 **** COULOMB, 1, ! _pushEnergyCoulombEnergy, ! _pushEnergyCoulombForce, _pushEnergyCoulombSupport }; --- 218,222 ---- COULOMB, 1, ! _pushEnergyCoulombEval, _pushEnergyCoulombSupport }; *************** *** 263,282 **** ** 0 parms! */ ! double ! _pushEnergyCotanEnergy(double dist, const double *parm) { ! double ss; ! ! AIR_UNUSED(parm); ! ss = sin(dist*AIR_PI/2.0); ! return (AIR_PI/2.0)*(dist > 1 ? 0 : 1.0 - 1.0/(ss*ss)); ! } ! ! double ! _pushEnergyCotanForce(double dist, const double *parm) { ! double ss; AIR_UNUSED(parm); ! ss = sin(dist*AIR_PI/2.0); ! return (AIR_PI/2.0)*(dist > 1 ? 0 : 1.0 - 1.0/(ss*ss)); } --- 229,243 ---- ** 0 parms! */ ! void ! _pushEnergyCotanEval(double *enr, double *frc, ! double dist, const double *parm) { ! double pot, cc; AIR_UNUSED(parm); ! pot = AIR_PI/2.0; ! cc = 1.0/(FLT_MIN + tan(dist*pot)); ! *enr = dist > 1 ? 0 : cc + dist*pot - pot; ! *frc = dist > 1 ? 0 : -cc*cc*pot; ! return; } *************** *** 292,297 **** COTAN, 0, ! _pushEnergyCotanEnergy, ! _pushEnergyCotanForce, _pushEnergyCotanSupport }; --- 253,257 ---- COTAN, 0, ! _pushEnergyCotanEval, _pushEnergyCotanSupport }; *************** *** 304,313 **** ** 0 parms: */ ! double ! _pushEnergyZeroFunc(double dist, const double *parm) { AIR_UNUSED(dist); AIR_UNUSED(parm); ! return 0.0; } --- 264,276 ---- ** 0 parms: */ ! void ! _pushEnergyZeroEval(double *enr, double *frc, ! double dist, const double *parm) { AIR_UNUSED(dist); AIR_UNUSED(parm); ! *enr = 0; ! *frc = 0; ! return; } *************** *** 323,328 **** ZERO, 0, ! _pushEnergyZeroFunc, ! _pushEnergyZeroFunc, _pushEnergyZeroSupport }; --- 286,290 ---- ZERO, 0, ! _pushEnergyZeroEval, _pushEnergyZeroSupport }; *************** *** 359,362 **** --- 321,338 ---- } + void + pushEnergySpecSet(pushEnergySpec *ensp, const pushEnergy *energy, + const double parm[PUSH_ENERGY_PARM_NUM]) { + unsigned int pi; + + if (ensp && energy && parm) { + ensp->energy = energy; + for (pi=0; pi<PUSH_ENERGY_PARM_NUM; pi++) { + ensp->parm[pi] = parm[pi]; + } + } + return; + } + pushEnergySpec * pushEnergySpecNix(pushEnergySpec *ensp) { Index: binning.c =================================================================== RCS file: /cvsroot/teem/teem/src/push/binning.c,v retrieving revision 1.18 retrieving revision 1.19 diff -C2 -d -r1.18 -r1.19 *** binning.c 22 Oct 2006 12:22:15 -0000 1.18 --- binning.c 25 Oct 2006 12:19:19 -0000 1.19 *************** *** 24,27 **** --- 24,31 ---- #include "privatePush.h" + /* + ** because the pushContext keeps an array of bins (not pointers to them) + ** we have Init and Done functions (not New and Nix) + */ void pushBinInit(pushBin *bin, unsigned int incr) { *************** *** 31,39 **** bin->pointArr = airArrayNew((void**)&(bin->point), &(bin->pointNum), sizeof(pushPoint *), incr); - bin->thingNum = 0; - bin->thing = NULL; - bin->thingArr = airArrayNew((void**)&(bin->thing), &(bin->thingNum), - sizeof(pushThing *), incr); - /* airArray callbacks are tempting but super confusing .... */ bin->neighbor = NULL; return; --- 35,38 ---- *************** *** 41,46 **** /* ! ** bins own the "thing" they contain, when you nix a bin, you nix the ! ** the things inside, but not the points (they belong to things) */ void --- 40,44 ---- /* ! ** bins own the points they contain- so this frees them */ void *************** *** 48,56 **** unsigned int idx; ! bin->pointArr = airArrayNuke(bin->pointArr); ! for (idx=0; idx<bin->thingNum; idx++) { ! bin->thing[idx] = pushThingNix(bin->thing[idx]); } ! bin->thingArr = airArrayNuke(bin->thingArr); bin->neighbor = (pushBin **)airFree(bin->neighbor); return; --- 46,53 ---- unsigned int idx; ! for (idx=0; idx<bin->pointNum; idx++) { ! bin->point[idx] = pushPointNix(bin->point[idx]); } ! bin->pointArr = airArrayNuke(bin->pointArr); bin->neighbor = (pushBin **)airFree(bin->neighbor); return; *************** *** 59,70 **** /* ! ** bins on boundary now extend to infinity; so this never returns NULL */ pushBin * _pushBinLocate(pushContext *pctx, double *_posWorld) { ! /* char me[]="_pushBinLocate"; */ double posWorld[4], posIdx[4]; unsigned int axi, eidx[3], binIdx; if (pctx->binSingle) { binIdx = 0; --- 56,74 ---- /* ! ** bins on boundary now extend to infinity; so the only time this ! ** returns NULL (indicating error) is for non-existant positions */ pushBin * _pushBinLocate(pushContext *pctx, double *_posWorld) { ! char me[]="_pushBinLocate", err[BIFF_STRLEN]; double posWorld[4], posIdx[4]; unsigned int axi, eidx[3], binIdx; + if (!ELL_3V_EXISTS(_posWorld)) { + sprintf(err, "%s: non-existant position (%g,%g,%g)", me, + _posWorld[0], _posWorld[1], _posWorld[2]); + biffAdd(PUSH, err); return NULL; + } + if (pctx->binSingle) { binIdx = 0; *************** *** 93,110 **** /* ! ** HEY: doesn't check to make sure thing isn't already in bin */ void - _pushBinThingAdd(pushContext *pctx, pushBin *bin, pushThing *thing) { - int thgI; - - AIR_UNUSED(pctx); - thgI = airArrayLenIncr(bin->thingArr, 1); - bin->thing[thgI] = thing; - - return; - } - - void _pushBinPointAdd(pushContext *pctx, pushBin *bin, pushPoint *point) { int pntI; --- 97,103 ---- /* ! ** this makes the bin the owner of the point */ void _pushBinPointAdd(pushContext *pctx, pushBin *bin, pushPoint *point) { int pntI; *************** *** 118,135 **** /* ! ** remove the pointer to the thing from the bin, don't kill the thing ! ** because we don't know here if this is part of a killing, or a rebinning */ void - _pushBinThingRemove(pushContext *pctx, pushBin *bin, int loseIdx) { - - AIR_UNUSED(pctx); - bin->thing[loseIdx] = bin->thing[bin->thingNum-1]; - airArrayLenIncr(bin->thingArr, -1); - - return; - } - - void _pushBinPointRemove(pushContext *pctx, pushBin *bin, int loseIdx) { --- 111,117 ---- /* ! ** the bin loses track of the point, caller responsible for ownership */ void _pushBinPointRemove(pushContext *pctx, pushBin *bin, int loseIdx) { *************** *** 141,171 **** } - int - _pushBinPointNullify(pushContext *pctx, pushBin *oldBin, pushPoint *point) { - char me[]="_pushBinPointNullify", err[BIFF_STRLEN]; - pushBin *bin; - unsigned int pointIdx; - - if (!( bin = oldBin ? oldBin : _pushBinLocate(pctx, point->pos) )) { - sprintf(err, "%s: NULL bin for point %p (%g,%g,%g)", me, - AIR_CAST(void*, point), - point->pos[0], point->pos[1], point->pos[2]); - biffAdd(PUSH, err); return 1; - } - for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) { - if (point == bin->point[pointIdx]) { - break; - } - } - if (pointIdx == bin->pointNum) { - sprintf(err, "%s: point %p (%g,%g,%g) wasn't in its bin", me, - AIR_CAST(void*, point), - point->pos[0], point->pos[1], point->pos[2]); - biffAdd(PUSH, err); return 1; - } - bin->point[pointIdx] = NULL; - return 0; - } - void _pushBinNeighborSet(pushBin *bin, pushBin **nei, unsigned int num) { --- 123,126 ---- *************** *** 225,242 **** int - pushBinThingAdd(pushContext *pctx, pushThing *thing) { - char me[]="pushBinThingAdd", err[BIFF_STRLEN]; - pushBin *bin; - - if (!( bin = _pushBinLocate(pctx, thing->point.pos) )) { - sprintf(err, "%s: can't locate point of thing %p", - me, AIR_CAST(void*, thing)); - biffAdd(PUSH, err); return 1; - } - _pushBinThingAdd(pctx, bin, thing); - return 0; - } - - int pushBinPointAdd(pushContext *pctx, pushPoint *point) { char me[]="pushBinPointAdd", err[BIFF_STRLEN]; --- 180,183 ---- *************** *** 244,248 **** if (!( bin = _pushBinLocate(pctx, point->pos) )) { ! sprintf(err, "%s: can't locate point %p", me, AIR_CAST(void*, point)); biffAdd(PUSH, err); return 1; } --- 185,190 ---- if (!( bin = _pushBinLocate(pctx, point->pos) )) { ! sprintf(err, "%s: can't locate point %p %u", ! me, AIR_CAST(void*, point), point->ttaagg); biffAdd(PUSH, err); return 1; } *************** *** 257,320 **** int pushRebin(pushContext *pctx) { ! char me[]="pushRebin"; ! unsigned int oldBinIdx, thingIdx, pointIdx; pushBin *oldBin, *newBin; - pushThing *thing; pushPoint *point; ! /* even if there is a single bin, we have to toss out-of-bounds ! things, and prune nullified points */ ! for (oldBinIdx=0; oldBinIdx<pctx->binNum; oldBinIdx++) { ! oldBin = pctx->bin + oldBinIdx; ! ! /* quietly clear out pointers to points that got nullified, ! or that went out of bounds (not an error here) */ ! for (pointIdx=0; pointIdx<oldBin->pointNum; /* nope! */) { ! point = oldBin->point[pointIdx]; ! if (!point) { ! /* this point got nullified */ ! _pushBinPointRemove(pctx, oldBin, pointIdx); ! } else { ! if (!( newBin = _pushBinLocate(pctx, point->pos) )) { ! /* this point is out of bounds */ ! _pushBinPointRemove(pctx, oldBin, pointIdx); ! } else { ! if (oldBin != newBin) { ! _pushBinPointRemove(pctx, oldBin, pointIdx); ! _pushBinPointAdd(pctx, newBin, point); ! } else { ! /* its in the right bin, move on */ ! pointIdx++; ! } } - } - } /* for pointIdx */ - - for (thingIdx=0; thingIdx<oldBin->thingNum; /* nope! */) { - thing = oldBin->thing[thingIdx]; - if (!( newBin = _pushBinLocate(pctx, thing->point.pos) )) { - /* bad thing! I kill you now */ - fprintf(stderr, "%s: killing thing at (%g,%g,%g)\n", me, - thing->point.pos[0], - thing->point.pos[1], - thing->point.pos[2]); - /* any real out-of-bounds points in the bins have already - been removed, so we don't need to nullify here */ - _pushBinThingRemove(pctx, oldBin, thingIdx); - pushThingNix(thing); - } else { if (oldBin != newBin) { ! _pushBinThingRemove(pctx, oldBin, thingIdx); ! _pushBinThingAdd(pctx, newBin, thing); ! /* don't increment thingIdx; the *next* thing index ! is now at thingIdx */ } else { ! /* this thing is already in the right bin, move on */ ! thingIdx++; } ! } ! } /* for thingIdx */ ! ! } /* for oldBinIdx */ return 0; --- 199,230 ---- int pushRebin(pushContext *pctx) { ! char me[]="pushRebin", err[BIFF_STRLEN]; ! unsigned int oldBinIdx, pointIdx; pushBin *oldBin, *newBin; pushPoint *point; ! if (!pctx->binSingle) { ! for (oldBinIdx=0; oldBinIdx<pctx->binNum; oldBinIdx++) { ! oldBin = pctx->bin + oldBinIdx; ! ! for (pointIdx=0; pointIdx<oldBin->pointNum; /* nope! */) { ! point = oldBin->point[pointIdx]; ! newBin = _pushBinLocate(pctx, point->pos); ! if (!newBin) { ! sprintf(err, "%s: can't locate point %p %u", ! me, AIR_CAST(void*, point), point->ttaagg); ! biffAdd(PUSH, err); return 1; } if (oldBin != newBin) { ! _pushBinPointRemove(pctx, oldBin, pointIdx); ! _pushBinPointAdd(pctx, newBin, point); } else { ! /* its in the right bin, move on */ ! pointIdx++; } ! } /* for pointIdx */ ! ! } /* for oldBinIdx */ ! } return 0; Index: corePush.c =================================================================== RCS file: /cvsroot/teem/teem/src/push/corePush.c,v retrieving revision 1.25 retrieving revision 1.26 diff -C2 -d -r1.25 -r1.26 *** corePush.c 22 Oct 2006 12:22:15 -0000 1.25 --- corePush.c 25 Oct 2006 12:19:19 -0000 1.26 *************** *** 25,29 **** /* ! ** this is run once per task (thread) */ int --- 25,30 ---- /* ! ** this is the core of the worker threads: as long as there are bins ! ** left to process, get the next one, and process it */ int *************** *** 43,47 **** } } while (binIdx < task->pctx->binNum - && 0 == task->pctx->bin[binIdx].thingNum && 0 == task->pctx->bin[binIdx].pointNum); if (task->pctx->threadNum > 1) { --- 44,47 ---- *************** *** 66,70 **** void * _pushWorker(void *_task) { ! char me[]="_pushWorker", *err; pushTask *task; --- 66,70 ---- void * _pushWorker(void *_task) { ! char me[]="_pushWorker", err[BIFF_STRLEN]; pushTask *task; *************** *** 73,77 **** while (1) { if (task->pctx->verbose > 1) { ! fprintf(stderr, "%s(%u): waiting to check finished\n", me, task->threadIdx); } --- 73,77 ---- while (1) { if (task->pctx->verbose > 1) { ! fprintf(stderr, "%s(%u): waiting on barrier A\n", me, task->threadIdx); } *************** *** 84,99 **** break; } ! /* else there's work to do ... */ ! if (task->pctx->verbose > 1) { fprintf(stderr, "%s(%u): starting to process\n", me, task->threadIdx); } if (_pushProcess(task)) { ! err = biffGetDone(PUSH); ! fprintf(stderr, "%s: thread %u trouble:\n%s", me, ! task->threadIdx, err); task->pctx->finished = AIR_TRUE; ! /* HEY: we should be using the "finished" mechanism to ! shut the whole production down */ } airThreadBarrierWait(task->pctx->iterBarrierB); --- 84,100 ---- break; } ! /* else there's work to do ... */ if (task->pctx->verbose > 1) { fprintf(stderr, "%s(%u): starting to process\n", me, task->threadIdx); } if (_pushProcess(task)) { ! /* HEY clearly not threadsafe ... */ ! sprintf(err, "%s: thread %u trouble", me, task->threadIdx); ! biffAdd(PUSH, err); task->pctx->finished = AIR_TRUE; ! } ! if (task->pctx->verbose > 1) { ! fprintf(stderr, "%s(%u): waiting on barrier B\n", ! me, task->threadIdx); } airThreadBarrierWait(task->pctx->iterBarrierB); *************** *** 112,117 **** biffAdd(PUSH, err); return 1; } ! if (!( pctx->thingNum >= 1 )) { ! sprintf(err, "%s: pctx->thingNum (%d) not >= 1\n", me, pctx->thingNum); biffAdd(PUSH, err); return 1; } --- 113,118 ---- biffAdd(PUSH, err); return 1; } ! if (!( pctx->pointNum >= 1 )) { ! sprintf(err, "%s: pctx->pointNum (%d) not >= 1\n", me, pctx->pointNum); biffAdd(PUSH, err); return 1; } *************** *** 150,201 **** } } - if (pctx->nstn) { - if (!pctx->npos) { - sprintf(err, "%s: can't have start/num nrrd w/out position nrrd", me); - biffAdd(PUSH, err); return 1; - } - if (nrrdCheck(pctx->nstn)) { - sprintf(err, "%s: got a broken start/num nrrd", me); - biffMove(PUSH, err, NRRD); return 1; - } - if (!( 2 == pctx->nstn->dim - && nrrdTypeUInt == pctx->nstn->type - && 3 == pctx->nstn->axis[0].size )) { - sprintf(err, "%s: start/num nrrd not 2-D 3-by-N array of %ss", me, - airEnumStr(nrrdType, nrrdTypeUInt)); - biffAdd(PUSH, err); return 1; - } - } if (tenGageUnknown != pctx->gravItem) { - unsigned int nii; if (airEnumValCheck(tenGage, pctx->gravItem)) { sprintf(err, "%s: gravity item %u invalid", me, pctx->gravItem); biffAdd(PUSH, err); return 1; } ! if (3 != tenGageKind->table[pctx->gravItem].answerLength) { ! sprintf(err, "%s: answer length of gravity item %s is %u, not 3", me, airEnumStr(tenGage, pctx->gravItem), tenGageKind->table[pctx->gravItem].answerLength); biffAdd(PUSH, err); return 1; } if (!AIR_EXISTS(pctx->gravScl)) { sprintf(err, "%s: gravity scaling doesn't exist", me); biffAdd(PUSH, err); return 1; } ! for (nii=0; nii<=1; nii++) { ! if (tenGageUnknown != pctx->gravNotItem[nii]) { ! if (airEnumValCheck(tenGage, pctx->gravNotItem[nii])) { ! sprintf(err, "%s: not gravity item[%u] %u invalid", ! me, nii, pctx->gravNotItem[nii]); ! biffAdd(PUSH, err); return 1; ! } ! if (3 != tenGageKind->table[pctx->gravItem].answerLength) { ! sprintf(err, "%s: answer length of not gravity item[%u] " ! "%s is %u, not 3", me, nii, ! air... [truncated message content] |