1 /* 2 Teem: Tools to process and visualize scientific data and images . 3 Copyright (C) 2012, 2011, 2010, 2009 University of Chicago 4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann 5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah 6 7 This library is free software; you can redistribute it and/or 8 modify it under the terms of the GNU Lesser General Public License 9 (LGPL) as published by the Free Software Foundation; either 10 version 2.1 of the License, or (at your option) any later version. 11 The terms of redistributing and/or modifying this software also 12 include exceptions to the LGPL that facilitate static linking. 13 14 This library is distributed in the hope that it will be useful, 15 but WITHOUT ANY WARRANTY; without even the implied warranty of 16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 Lesser General Public License for more details. 18 19 You should have received a copy of the GNU Lesser General Public License 20 along with this library; if not, write to Free Software Foundation, Inc., 21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 22 */ 23 24 #ifndef PUSH_HAS_BEEN_INCLUDED 25 #define PUSH_HAS_BEEN_INCLUDED 26 27 #include <math.h> 28 29 #include <teem/air.h> 30 #include <teem/biff.h> 31 #include <teem/nrrd.h> 32 #include <teem/ell.h> 33 #include <teem/gage.h> 34 #include <teem/ten.h> 35 36 #if defined(_WIN32) && !defined(__CYGWIN__) && !defined(TEEM_STATIC) 37 # if defined(TEEM_BUILD) || defined(push_EXPORTS) || defined(teem_EXPORTS) 38 # define PUSH_EXPORT extern __declspec(dllexport) 39 # else 40 # define PUSH_EXPORT extern __declspec(dllimport) 41 # endif 42 #else /* TEEM_STATIC || UNIX */ 43 # define PUSH_EXPORT extern 44 #endif 45 46 #ifdef __cplusplus 47 extern "C" { 48 #endif 49 50 #define PUSH pushBiffKey 51 #define PUSH_THREAD_MAXNUM 512 52 53 /* 54 ******** pushPoint 55 ** 56 ** information about a point in the simulation. There are really two 57 ** kinds of information here: "pos", "enr", "frc" pertain to the 58 ** simulation of point dynamics, while "ten", "inv", "cnt", "grav", 59 ** "gravGrad", "seedThresh" are properties of the field sampled at the 60 ** point. 61 */ 62 typedef struct pushPoint_t { 63 unsigned int ttaagg; 64 double pos[3], /* position in world space */ 65 enr, /* energy accumulator (current iteration) */ 66 frc[3], /* force accumulator (current iteration) */ 67 ten[7], /* tensor here */ 68 inv[7], /* inverse of tensor */ 69 cnt[3], /* mask's containment gradient */ 70 grav, gravGrad[3], /* gravity stuff */ 71 seedThresh; /* seed thresh */ 72 /* per-point list of active neighbors- which is updated only periodically. 73 In addition to spatial binning, this greatly reduces the number of 74 pair-wise interactions computed (based on idea from Meyer et al.) */ 75 struct pushPoint_t **neigh; 76 unsigned int neighNum; 77 airArray *neighArr; 78 } pushPoint; 79 80 /* 81 ******** pushBin 82 ** 83 ** the data structure for doing spatial binning. 84 ** 85 ** in tractlet-less push, bins do own the points they contain 86 */ 87 typedef struct pushBin_t { 88 unsigned int pointNum; /* # of points in this bin */ 89 pushPoint **point; /* dyn. alloc. array of point pointers */ 90 airArray *pointArr; /* airArray around point and pointNum */ 91 struct pushBin_t **neighbor; /* pre-computed NULL-terminated list of all 92 neighboring bins, including myself */ 93 } pushBin; 94 95 /* 96 ******** pushTask 97 ** 98 ** The information specific for a thread. 99 */ 100 typedef struct pushTask_t { 101 struct pushContext_t *pctx; /* parent's context */ 102 gageContext *gctx; /* result of gageContextCopy(pctx->gctx) */ 103 const double *tenAns, /* results of gage probing */ 104 *invAns, *cntAns, 105 *gravAns, *gravGradAns, 106 *seedThreshAns; 107 airThread *thread; /* my thread */ 108 unsigned int threadIdx, /* which thread am I */ 109 pointNum; /* # points I let live this iteration */ 110 double energySum, /* sum of energies of points I processed */ 111 deltaFracSum; /* contribution to pctx->deltaFrac */ 112 airRandMTState *rng; /* state for my RNG */ 113 void *returnPtr; /* for airThreadJoin */ 114 } pushTask; 115 116 /* 117 ******** pushEnergyType* enum 118 ** 119 ** the different shapes of potential energy profiles that can be used 120 */ 121 enum { 122 pushEnergyTypeUnknown, /* 0 */ 123 pushEnergyTypeSpring, /* 1 */ 124 pushEnergyTypeGauss, /* 2 */ 125 pushEnergyTypeCoulomb, /* 3 */ 126 pushEnergyTypeCotan, /* 4 */ 127 pushEnergyTypeZero, /* 5 */ 128 pushEnergyTypeLast 129 }; 130 #define PUSH_ENERGY_TYPE_MAX 5 131 #define PUSH_ENERGY_PARM_NUM 3 132 133 /* 134 ******** pushEnergy 135 ** 136 ** the functions which determine inter-point forces 137 ** 138 ** NOTE: the eval() function probably does NOT check to see it was passed 139 ** non-NULL pointers into which to store energy and force 140 */ 141 typedef struct { 142 char name[AIR_STRLEN_SMALL]; 143 unsigned int parmNum; 144 void (*eval)(double *energy, double *force, 145 double dist, const double parm[PUSH_ENERGY_PARM_NUM]); 146 double (*support)(const double parm[PUSH_ENERGY_PARM_NUM]); 147 } pushEnergy; 148 149 typedef struct { 150 const pushEnergy *energy; 151 double parm[PUSH_ENERGY_PARM_NUM]; 152 } pushEnergySpec; 153 154 /* 155 ******** pushContext 156 ** 157 ** everything for doing one simulation computation 158 ** 159 */ 160 typedef struct pushContext_t { 161 /* INPUT ----------------------------- */ 162 unsigned int pointNum; /* number points to start simulation w/ */ 163 Nrrd *nin, /* 3D image of 3D masked tensors, though 164 it may only be a single slice */ 165 *npos; /* positions to start with 166 (overrides pointNum) */ 167 double stepInitial, /* initial time step in integration 168 (which will be reduced as the system 169 converges) */ 170 scale, /* scaling from tensor to glyph size */ 171 wall, /* spring constant of walls */ 172 cntScl, /* magnitude of containment gradient */ 173 deltaLimit, /* speed limit on particles' motion, as a 174 fraction of glyph radius along 175 direction of motion */ 176 deltaFracMin, /* lowest value of deltaFrac (see below) 177 that is allowed without decreasing 178 step size */ 179 energyStepFrac, /* when energy goes up instead of down, the 180 fraction by which to scale step size */ 181 deltaFracStepFrac, /* when deltaFrac goes below deltaFracMin, 182 fraction by which to scale step size */ 183 neighborTrueProb, /* probability that we find the true 184 neighbors of the particle, as opposed to 185 using a cached list */ 186 probeProb, /* probability that we gageProbe() to find 187 the local tensor value, instead of 188 re-using last value */ 189 energyImprovMin; /* convergence threshold: stop when 190 fracional improvement (decrease) in 191 energy dips below this */ 192 int detReject, /* determinant-based rejection at init */ 193 midPntSmp, /* sample midpoint btw part.s for physics */ 194 verbose; /* blah blah blah */ 195 unsigned int seedRNG, /* seed value for random number generator */ 196 threadNum, /* number of threads to use */ 197 maxIter, /* if non-zero, max number of iterations */ 198 snap; /* if non-zero, interval between iterations 199 at which output snapshots are saved */ 200 int gravItem, /* tenGage item (scalar) for "height" 201 potential energy associated w/ gravity */ 202 gravGradItem; /* tenGage item (vector) for gravity */ 203 double gravScl, /* sign and magnitude of gravity's effect: 204 when this is positive, higher values of 205 gravItem have higher potential energy */ 206 gravZero; /* the height that corresponds to zero 207 potential energy from gravity */ 208 209 int seedThreshItem, /* item for constraining random seeding */ 210 seedThreshSign; /* +1: need val > thresh; -1: opposite */ 211 double seedThresh; /* threshold for seed constraint */ 212 213 pushEnergySpec *ensp; /* potential energy function to use */ 214 215 int binSingle; /* disable binning (for debugging) */ 216 unsigned int binIncr; /* increment for per-bin airArray */ 217 218 NrrdKernelSpec *ksp00, /* for sampling tensor field */ 219 *ksp11, /* for gradient of mask, other 1st derivs */ 220 *ksp22; /* for 2nd derivatives */ 221 222 /* INTERNAL -------------------------- */ 223 224 unsigned int ttaagg; /* next value for per-point ID */ 225 Nrrd *nten, /* 3D image of 3D masked tensors */ 226 *ninv, /* pre-computed inverse of nten */ 227 *nmask; /* mask image from nten */ 228 gageContext *gctx; /* gage context around nten, ninv, nmask */ 229 gagePerVolume *tpvl, *ipvl; /* gage pervolumes around nten and ninv */ 230 int finished; /* used to signal all threads to return */ 231 unsigned int dimIn, /* dim (2 or 3) of input, meaning whether 232 it was a single slice or a full volume */ 233 sliceAxis; /* got a single 3-D slice, which axis had 234 only a single sample */ 235 236 pushBin *bin; /* volume of bins (see binsEdge, binNum) */ 237 unsigned int binsEdge[3], /* # bins along each volume edge, 238 determined by maxEval and scale */ 239 binNum, /* total # bins in grid */ 240 binIdx; /* *next* bin of points needing to be 241 processed. Stage is done when 242 binIdx == binNum */ 243 airThreadMutex *binMutex; /* mutex around bin */ 244 245 double step, /* current working step size */ 246 maxDist, /* max distance btween interacting points */ 247 maxEval, meanEval, /* max and mean principal eval in field */ 248 maxDet, 249 energySum; /* potential energy of entire particles */ 250 pushTask **task; /* dynamically allocated array of tasks */ 251 airThreadBarrier *iterBarrierA; /* barriers between iterations */ 252 airThreadBarrier *iterBarrierB; /* barriers between iterations */ 253 double deltaFrac; /* mean (over all particles in last 254 iteration) of fraction of distance 255 actually travelled to distance that it 256 wanted to travel (due to speed limit) */ 257 258 /* OUTPUT ---------------------------- */ 259 260 double timeIteration, /* time needed for last (single) iter */ 261 timeRun; /* total time spent in computation */ 262 unsigned int iter; /* how many iterations were needed */ 263 Nrrd *noutPos, /* list of 2D or 3D positions */ 264 *noutTen; /* list of 2D or 3D masked tensors */ 265 } pushContext; 266 267 typedef union { 268 pushPoint ***point; 269 void **v; 270 } pushPtrPtrUnion; 271 272 /* defaultsPush.c */ 273 PUSH_EXPORT const int pushPresent; 274 PUSH_EXPORT const char *pushBiffKey; 275 276 /* methodsPush.c */ 277 PUSH_EXPORT pushPoint *pushPointNew(pushContext *pctx); 278 PUSH_EXPORT pushPoint *pushPointNix(pushPoint *pnt); 279 PUSH_EXPORT pushContext *pushContextNew(void); 280 PUSH_EXPORT pushContext *pushContextNix(pushContext *pctx); 281 282 /* forces.c (legacy name for info about (derivatives of) energy functions) */ 283 PUSH_EXPORT const airEnum *const pushEnergyType; 284 PUSH_EXPORT const pushEnergy *const pushEnergyUnknown; 285 PUSH_EXPORT const pushEnergy *const pushEnergySpring; 286 PUSH_EXPORT const pushEnergy *const pushEnergyGauss; 287 PUSH_EXPORT const pushEnergy *const pushEnergyCoulomb; 288 PUSH_EXPORT const pushEnergy *const pushEnergyCotan; 289 PUSH_EXPORT const pushEnergy *const pushEnergyZero; 290 PUSH_EXPORT const pushEnergy *const pushEnergyAll[PUSH_ENERGY_TYPE_MAX+1]; 291 PUSH_EXPORT pushEnergySpec *pushEnergySpecNew(void); 292 PUSH_EXPORT void pushEnergySpecSet(pushEnergySpec *ensp, 293 const pushEnergy *energy, 294 const double parm[PUSH_ENERGY_PARM_NUM]); 295 PUSH_EXPORT pushEnergySpec *pushEnergySpecNix(pushEnergySpec *ensp); 296 PUSH_EXPORT int pushEnergySpecParse(pushEnergySpec *ensp, const char *str); 297 PUSH_EXPORT hestCB *pushHestEnergySpec; 298 299 /* corePush.c */ 300 PUSH_EXPORT int pushStart(pushContext *pctx); 301 PUSH_EXPORT int pushIterate(pushContext *pctx); 302 PUSH_EXPORT int pushRun(pushContext *pctx); 303 PUSH_EXPORT int pushFinish(pushContext *pctx); 304 305 /* binning.c */ 306 PUSH_EXPORT void pushBinInit(pushBin *bin, unsigned int incr); 307 PUSH_EXPORT void pushBinDone(pushBin *bin); 308 PUSH_EXPORT int pushBinPointAdd(pushContext *pctx, pushPoint *point); 309 PUSH_EXPORT void pushBinAllNeighborSet(pushContext *pctx); 310 PUSH_EXPORT int pushRebin(pushContext *pctx); 311 312 /* action.c */ 313 PUSH_EXPORT int pushBinProcess(pushTask *task, unsigned int myBinIdx); 314 PUSH_EXPORT int pushOutputGet(Nrrd *nPos, Nrrd *nTen, Nrrd *nEnr, 315 pushContext *pctx); 316 317 #ifdef __cplusplus 318 } 319 #endif 320 321 #endif /* PUSH_HAS_BEEN_INCLUDED */ 322 323