Update of /cvsroot/teem/teem/src/coil In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv3054 Added Files: coil.h GNUmakefile coreCoil.c defaultsCoil.c enumsCoil.c methodsCoil.c realmethods.c scalarCoil.c tensorCoil.c Log Message: this is the last teem library that Gordon will write while he does not have a PhD --- NEW FILE: defaultsCoil.c --- /* teem: Gordon Kindlmann's research software Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include "coil.h" const char * coilBiffKey = "coil"; int coilDefaultRadius = 1; --- NEW FILE: GNUmakefile --- # # teem: Gordon Kindlmann's research software # Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah # # This library is free software; you can redistribute it and/or # modify it under the terms of the GNU Lesser General Public # License as published by the Free Software Foundation; either # version 2.1 of the License, or (at your option) any later version. # # This library is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public # License along with this library; if not, write to the Free Software # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # # #### Library name #### #### L := coil #### #### #### # boilerplate: default targets and include tricks TEEM_ROOT ?= ../.. TEEM_SRC ?= .. ifeq (,$(DEF_TARGETS)) DEF_TARGETS = true dev : $(L)/dev install : $(L)/install clean : $(L)/clean clobber : $(L)/clobber include ../GNUmakefile endif ifeq (,$($(L).SEEN)) $(L).SEEN := true #### Describe library here #### (ell for macros) #### $(L).NEED = ten gage ell nrrd biff air $(L).PUBLIC_HEADERS = coil.h $(L).PRIVATE_HEADERS = $(L).OBJS = defaultsCoil.o enumsCoil.o scalarCoil.o tensorCoil.o \ realmethods.o methodsCoil.o coreCoil.o $(L).TESTS = test/coiler #### #### #### # boilerplate: declare rules for this library include $(TEEM_SRC)/make/template.mk endif ifeq (,$(INCLUDED)) include $(TEEM_SRC)/bin/GNUmakefile endif --- NEW FILE: coil.h --- /* teem: Gordon Kindlmann's research software Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #ifndef COIL_HAS_BEEN_INCLUDED #define COIL_HAS_BEEN_INCLUDED #include <teem/air.h> #include <teem/biff.h> #include <teem/ell.h> #include <teem/nrrd.h> #include <teem/gage.h> #include <teem/ten.h> #ifdef __cplusplus extern "C" { #endif #define COIL coilBiffKey /* ******** coil_t ** ** this is the very crude means by which you can control the type ** of values that coil works with: "float" or "double". It is an ** unfortunate but greatly simplifying restriction that this type ** is used for all kinds of volumes, and all methods of filtering ** ** So: choose double by defining TEEM_COIL_TYPE_DOUBLE and float ** otherwise. */ #ifdef TEEM_COIL_TYPE_DOUBLE typedef double coil_t; #define coil_nrrdType nrrdTypeDouble #define COIL_TYPE_FLOAT 0 #else typedef float coil_t; #define coil_nrrdType nrrdTypeFloat #define COIL_TYPE_FLOAT 1 #endif /* ******** #define COIL_PARMS_NUM ** ** maximum number of parameters that may be needed by any coil-driven ** filtering method */ #define COIL_PARMS_NUM 5 /* ******** coilMethodType* enum ** ** enumerated possibilities for different filtering methods */ enum { coilMethodTypeUnknown, /* 0 */ coilMethodTypeTesting, /* 1: basically a no-op */ coilMethodTypeIsotropic, /* 2 */ coilMethodTypeAnisotropic, /* 3 */ coilMethodTypeModifiedCurvature, /* 4 */ coilMethodTypeCurvatureFlow, /* 5 */ coilMethodTypeLast }; #define COIL_METHOD_TYPE_MAX 5 /* ******** coilMethod struct ** ** someday, there will be total orthogonality between kind and method. ** until then, this will have only a portion of the things relating to ** running one method, regardless of kind */ typedef struct { char name[AIR_STRLEN_SMALL]; int type; /* from coilMethodType* enum */ int numParm; /* number of parameters we need */ } coilMethod; /* ******** coilKindType* enum ** ** enumerated possibilities for different kinds */ enum { coilKindTypeUnknown, /* 0 */ coilKindTypeScalar, /* 1 */ coilKindType3Color, /* 2 */ coilKindType7Tensor, /* 3 */ coilKindTypeLast }; #define COIL_KIND_TYPE_MAX 3 /* ******** coilKind struct ** ** yes, there is some redunancy with the gageKind, but the the main ** reason is that its significantly easier to implement meaningful ** per-sample filtering of various quantities, than it is to do ** gage-style convolution-based measurements at arbitrary locations. ** So, there will probably be significantly more coilKinds than ** there are gageKinds. The two kind systems may play closer together ** at some point in the future where appropriate. */ typedef struct { char name[AIR_STRLEN_SMALL]; /* short identifying string for kind */ int valLen; /* number of scalars per data point 1 for plain scalars (baseDim=0), or something else (baseDim=1) */ /* all the available methods */ void (*filter[COIL_METHOD_TYPE_MAX+1])(coil_t *delta, coil_t *iv3, double spacing[3], double parm[COIL_PARMS_NUM]); void (*update)(coil_t *val, coil_t *delta); /* how to apply update */ } coilKind; struct coilContext_t; /* ******** coilTask ** ** passed to all worker threads */ typedef struct { struct coilContext_t *cctx; /* parent's context */ airThread *thread; /* my thread */ int threadIdx, /* which thread am I */ startZ, endZ; /* my index range (inclusive) of slices */ coil_t *iv3; /* cache of local values, ordering same as in gage: any multiple values per voxel are the slowest axis, not the fastest */ void *returnPtr; /* for airThreadJoin */ } coilTask; /* ******** coilContext struct ** ** bag of stuff relating to filtering one volume */ typedef struct coilContext_t { /* ---------- input */ const Nrrd *nin; /* input volume, converted to coil_t nvol */ const coilKind *kind; /* what kind of volume is nin */ const coilMethod *method; /* what method of filtering to use */ int radius, /* how big a neighborhood to look at when doing filtering (use 1 for 3x3x3 size) */ numThreads, /* number of threads to enlist */ verbose; /* blah blah blah */ double parm[COIL_PARMS_NUM]; /* all the parameters used to control the action of the filtering. The timestep is probably the first value. */ /* ---------- internal */ int size[3]; /* size of volume */ double spacing[3]; /* sample spacings we'll use- we perhaps should be using a gageShape, but this is actually all we really need... */ Nrrd *nvol; /* an interleaved volume of (1st) the last filtering result, and (2nd) the update values from the current iteration */ int finished; /* used to signal all threads to return */ coilTask **task; /* dynamically allocated array of tasks */ airThreadBarrier *filterBarrier, /* for when filtering is done, and the "update" values have been set */ *updateBarrier, /* after the update values have been applied to current values */ *finishBarrier; /* so that thread 0 can see if filtering should go onward, and set "finished" */ } coilContext; /* defaultsCoil.c */ TEEM_API const char *coilBiffKey; TEEM_API int coilDefaultRadius; /* enumsCoil.c */ TEEM_API airEnum *coilMethodType; TEEM_API airEnum *coilKindType; /* scalarCoil.c */ TEEM_API const coilKind *coilKindScalar; /* tensorCoil.c */ TEEM_API const coilKind *coilKindTensor; /* realmethods.c */ TEEM_API const coilMethod *coilMethodTesting; TEEM_API const coilMethod *coilMethodIsotropic; TEEM_API const coilMethod *coilMethodArray[COIL_METHOD_TYPE_MAX+1]; /* methodsCoil.c (sorry, confusing name!) */ TEEM_API coilContext *coilContextNew(); TEEM_API int coilVolumeCheck(const Nrrd *nin, const coilKind *kind); TEEM_API int coilContextAllSet(coilContext *cctx, const Nrrd *nin, const coilKind *kind, const coilMethod *method, int radius, int numThreads, int verbose, double parm[COIL_PARMS_NUM]); TEEM_API int coilOutputGet(Nrrd *nout, coilContext *cctx); TEEM_API coilContext *coilContextNix(coilContext *cctx); /* coreCoil.c */ TEEM_API int coilStart(coilContext *cctx); TEEM_API int coilIterate(coilContext *cctx, int numIterations); TEEM_API int coilFinish(coilContext *cctx); #ifdef __cplusplus } #endif #endif /* COIL_HAS_BEEN_INCLUDED */ --- NEW FILE: tensorCoil.c --- /* teem: Gordon Kindlmann's research software Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include "coil.h" void _coilKind7TensorFilterTesting(coil_t *delta, coil_t *iv3, double spacing[3], double parm[COIL_PARMS_NUM]) { delta[0] = 0; delta[1] = 0; delta[2] = 0; delta[3] = 0; delta[4] = 0; delta[5] = 0; delta[6] = 0; } void _coilKind7TensorFilterIsotropic(coil_t *delta, coil_t *iv3, double spacing[3], double parm[COIL_PARMS_NUM]) { } void _coilKindTensorUpdate(coil_t *val, coil_t *delta) { val[0] += delta[0]; /* WARNING: this could change confidence! */ val[1] += delta[1]; val[2] += delta[2]; val[3] += delta[3]; val[4] += delta[4]; val[5] += delta[5]; val[6] += delta[6]; } const coilKind _coilKindTensor = { "tensor", 7, {NULL, _coilKind7TensorFilterTesting, _coilKind7TensorFilterIsotropic, NULL, NULL, NULL}, _coilKindTensorUpdate }; const coilKind * coilKindTensor = &_coilKindTensor; --- NEW FILE: scalarCoil.c --- /* teem: Gordon Kindlmann's research software Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include "coil.h" /* ** x ----> X ** \ 0 1 2 ** | \ 3 4 5 ** | Y 6 7 8 ** | ** | 9 10 11 ** Z 12 13 14 ** 15 16 17 ** ** 18 19 20 ** 21 22 23 ** 24 25 26 */ coil_t _coilLaplacian3(coil_t *iv3, double spacing[3]) { return ((iv3[12] - 2*iv3[13] + iv3[14])/(spacing[0]*spacing[0]) + (iv3[10] - 2*iv3[13] + iv3[16])/(spacing[1]*spacing[1]) + (iv3[ 4] - 2*iv3[13] + iv3[22])/(spacing[2]*spacing[2])); } void _coilKindScalarFilterTesting(coil_t *delta, coil_t *iv3, double spacing[3], double parm[COIL_PARMS_NUM]) { delta[0] = 0; } void _coilKindScalarFilterIsotropic(coil_t *delta, coil_t *iv3, double spacing[3], double parm[COIL_PARMS_NUM]) { coil_t lapl; lapl = _coilLaplacian3(iv3, spacing); delta[0] = parm[0]*lapl; } void _coilKindScalarUpdate(coil_t *val, coil_t *delta) { val[0] += delta[0]; } const coilKind _coilKindScalar = { "scalar", 1, {NULL, _coilKindScalarFilterTesting, _coilKindScalarFilterIsotropic, NULL, NULL, NULL}, _coilKindScalarUpdate }; const coilKind * coilKindScalar = &_coilKindScalar; --- NEW FILE: realmethods.c --- /* teem: Gordon Kindlmann's research software Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include "coil.h" const coilMethod _coilMethodTesting = { "testing", coilMethodTypeTesting, 0 }; const coilMethod* coilMethodTesting = &_coilMethodTesting; const coilMethod _coilMethodIsotropic = { "isotropic", coilMethodTypeIsotropic, 0 }; const coilMethod* coilMethodIsotropic = &_coilMethodIsotropic; const coilMethod* coilMethodArray[COIL_METHOD_TYPE_MAX+1] = { NULL, &_coilMethodTesting, &_coilMethodIsotropic, NULL, NULL, NULL }; --- NEW FILE: coreCoil.c --- /* teem: Gordon Kindlmann's research software Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include "coil.h" /* ** ** iv3 is: diam x diam x diam x valLen */ void _coilIv3Fill(coil_t *iv3, coil_t *here, int radius, int valLen, int x0, int y0, int z0, int sizeX, int sizeY, int sizeZ) { int diam, vi, /* value index */ xni, yni, zni, /* neighborhood (iv3) indices */ xvi, yvi, zvi; /* volume indices */ /* this should be parameterized on both radius and valLen */ /* this should shuffle values within iv3 when possible */ diam = 1 + 2*radius; for (zni=0; zni<diam; zni++) { zvi = AIR_CLAMP(0, zni-radius+z0, sizeZ-1); for (yni=0; yni<diam; yni++) { yvi = AIR_CLAMP(0, yni-radius+y0, sizeY-1); for (xni=0; xni<diam; xni++) { xvi = AIR_CLAMP(0, xni-radius+x0, sizeX-1); for (vi=0; vi<valLen; vi++) { iv3[xni + diam*(yni + diam*(zni + diam*vi))] = here[vi + valLen*(xvi + sizeX*(yvi + sizeY*zvi))]; } } } } return; } void _coilProcess(coilTask *task, int filter) { int xi, yi, zi, sizeX, sizeY, sizeZ, valLen, radius; coil_t *here; sizeX = task->cctx->size[0]; sizeY = task->cctx->size[1]; sizeZ = task->cctx->size[2]; valLen = task->cctx->kind->valLen; radius = task->cctx->radius; here = (coil_t*)(task->cctx->nvol->data); here += 2*valLen*sizeX*sizeY*task->startZ; if (filter) { for (zi=task->startZ; zi<=task->endZ; zi++) { for (yi=0; yi<sizeY; yi++) { for (xi=0; xi<sizeX; xi++) { _coilIv3Fill(task->iv3, here, radius, valLen, xi, yi, zi, sizeX, sizeY, sizeZ); /* filter */ here += 2*valLen; } } } } else { for (zi=task->startZ; zi<=task->endZ; zi++) { for (yi=0; yi<sizeY; yi++) { for (xi=0; xi<sizeX; xi++) { task->cctx->kind->update(here + 0*valLen, here + 1*valLen); here += 2*valLen; } } } } return; } coilTask * _coilTaskNew(coilContext *cctx, int threadIdx, int sizeZ) { coilTask *task; int len, diam; len = cctx->kind->valLen; diam = 1 + 2*cctx->radius; task = (coilTask *)calloc(1, sizeof(coilTask)); if (task) { task->cctx = cctx; task->thread = airThreadNew(); task->startZ = threadIdx*sizeZ/cctx->numThreads; task->endZ = (threadIdx+1)*sizeZ/cctx->numThreads - 1; task->threadIdx = threadIdx; task->iv3 = (coil_t*)calloc(len*diam*diam*diam, sizeof(coil_t)); task->returnPtr = NULL; } return task; } coilTask * _coilTaskNix(coilTask *task) { if (task) { task->thread = airThreadNix(task->thread); task->iv3 = airFree(task->iv3); free(task); } return NULL; } void * _coilWorker(void *_task) { char me[]="_coilWorker"; coilTask *task; task = (coilTask *)_task; while (1) { /* wait until parent has set cctx->finished */ if (task->cctx->verbose) { fprintf(stderr, "%s(%d): waiting to check finished\n", me, task->threadIdx); } airThreadBarrierWait(task->cctx->finishBarrier); if (task->cctx->finished) { if (task->cctx->verbose) { fprintf(stderr, "%s(%d): done!\n", me, task->threadIdx); } break; } /* else there's work */ /* first: filter */ if (task->cctx->verbose) { fprintf(stderr, "%s(%d): filtering ... \n", me, task->threadIdx); } _coilProcess(task, AIR_TRUE); airThreadBarrierWait(task->cctx->filterBarrier); /* second: update */ if (task->cctx->verbose) { fprintf(stderr, "%s(%d): updating ... \n", me, task->threadIdx); } _coilProcess(task, AIR_FALSE); airThreadBarrierWait(task->cctx->updateBarrier); } return _task; } int coilStart(coilContext *cctx) { char me[]="coilStart", err[AIR_STRLEN_MED]; int tidx, elIdx, valIdx, valLen; coil_t (*lup)(const void*, size_t), *val; if (!cctx) { sprintf(err, "%s: got NULL pointer", me); biffAdd(COIL, err); return 1; } cctx->task = (coilTask **)calloc(cctx->numThreads, sizeof(coilTask *)); if (!(cctx->task)) { sprintf(err, "%s: couldn't allocate array of tasks", me); biffAdd(COIL, err); return 1; } /* we create tasks for ALL threads, including me, thread 0 */ cctx->task[0] = NULL; for (tidx=0; tidx<cctx->numThreads; tidx++) { cctx->task[tidx] = _coilTaskNew(cctx, tidx, cctx->size[2]); if (!(cctx->task[tidx])) { sprintf(err, "%s: couldn't allocate task %d", me, tidx); biffAdd(COIL, err); return 1; } } cctx->finished = AIR_FALSE; cctx->finishBarrier = airThreadBarrierNew(cctx->numThreads); cctx->filterBarrier = airThreadBarrierNew(cctx->numThreads); cctx->updateBarrier = airThreadBarrierNew(cctx->numThreads); /* start threads 1 and up running (they won't get far) */ if (cctx->verbose) { fprintf(stderr, "%s: parent doing slices Z=[%d,%d]\n", me, cctx->task[0]->startZ, cctx->task[0]->endZ); } for (tidx=1; tidx<cctx->numThreads; tidx++) { if (cctx->verbose) { fprintf(stderr, "%s: spawning thread %d (Z=[%d,%d])\n", me, tidx, cctx->task[tidx]->startZ, cctx->task[tidx]->endZ); } airThreadStart(cctx->task[tidx]->thread, _coilWorker, (void *)(cctx->task[tidx])); } /* initialize the values in cctx->nvol */ val = (coil_t*)(cctx->nvol->data); valLen = cctx->kind->valLen; #if COIL_TYPE_FLOAT lup = nrrdFLookup[cctx->nin->type]; #else lup = nrrdDLookup[cctx->nin->type]; #endif for (elIdx=0; elIdx<cctx->size[0]*cctx->size[1]*cctx->size[2]; elIdx++) { for (valIdx=0; valIdx<valLen; valIdx++) { val[valIdx + 0*valLen] = lup(cctx->nin->data, valIdx + valLen*elIdx); val[valIdx + 1*valLen] = 0; } val += 2*valLen; } return 0; } int _coilPause(int huge) { int i, j, num; num = 0; for (i=1; i<huge; i++) { for (j=0; j<huge; j++) { num *= i; } } return num; } /* ******** coilIterate ** ** (documentation) ** ** NB: this implements the body of thread 0 */ int coilIterate(coilContext *cctx, int numIterations) { char me[]="coilIterate", err[AIR_STRLEN_MED]; int iter; if (!cctx) { sprintf(err, "%s: got NULL pointer", me); biffAdd(COIL, err); return 1; } for (iter=0; iter<numIterations; iter++) { if (cctx->verbose) { fprintf(stderr, "%s: starting iter %d\n", me, iter); } cctx->finished = AIR_FALSE; airThreadBarrierWait(cctx->finishBarrier); /* first: filter */ if (cctx->verbose) { fprintf(stderr, "%s: filtering ... \n", me); } _coilProcess(cctx->task[0], AIR_TRUE); airThreadBarrierWait(cctx->filterBarrier); /* second: update */ if (cctx->verbose) { fprintf(stderr, "%s: updating ... \n", me); } _coilProcess(cctx->task[0], AIR_FALSE); airThreadBarrierWait(cctx->updateBarrier); } return 0; } int coilFinish(coilContext *cctx) { char me[]="coilFinish", err[AIR_STRLEN_MED]; int tidx; if (!cctx) { sprintf(err, "%s: got NULL pointer", me); biffAdd(COIL, err); return 1; } if (cctx->verbose) { fprintf(stderr, "%s: finishing workers\n", me); } cctx->finished = AIR_TRUE; airThreadBarrierWait(cctx->finishBarrier); for (tidx=1; tidx<cctx->numThreads; tidx++) { airThreadJoin(cctx->task[tidx]->thread, &(cctx->task[tidx]->returnPtr)); cctx->task[tidx]->thread = airThreadNix(cctx->task[tidx]->thread); } cctx->task[0]->thread = airThreadNix(cctx->task[0]->thread); return 0; } --- NEW FILE: enumsCoil.c --- /* teem: Gordon Kindlmann's research software Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include "coil.h" char _coilMethodTypeStr[COIL_METHOD_TYPE_MAX+1][AIR_STRLEN_SMALL] = { "(unknown_method)", "testing", "isotropic", "anisotropic", "modified curvature", "curvature flow" }; char _coilMethodTypeDesc[COIL_METHOD_TYPE_MAX+1][AIR_STRLEN_MED] = { "unknown_method", "nothing, actually, just here for testing", "isotropic diffusion (Gaussian blurring)", "anisotropic diffusion (Perona-Malik)", "modified curvature diffusion", "curvature flow" }; char _coilMethodTypeStrEqv[][AIR_STRLEN_SMALL] = { "test", "testing", "iso", "isotropic", "aniso", "anisotropic", "mcde", "flow", "" }; int _coilMethodTypeValEqv[] = { coilMethodTypeTesting, coilMethodTypeTesting, coilMethodTypeIsotropic, coilMethodTypeIsotropic, coilMethodTypeAnisotropic, coilMethodTypeAnisotropic, coilMethodTypeModifiedCurvature, coilMethodTypeCurvatureFlow }; airEnum _coilMethodType = { "method", COIL_METHOD_TYPE_MAX, _coilMethodTypeStr, NULL, _coilMethodTypeDesc, _coilMethodTypeStrEqv, _coilMethodTypeValEqv, AIR_FALSE }; airEnum * coilMethodType = &_coilMethodType; /* -------------------------------------------------- */ char _coilKindTypeStr[COIL_KIND_TYPE_MAX+1][AIR_STRLEN_SMALL] = { "(unknown_kind)", "scalar", "3color", "7tensor" }; char _coilKindTypeDesc[COIL_KIND_TYPE_MAX+1][AIR_STRLEN_MED] = { "unknown_kind", "plain old scalar quantities", "3-component color", "ten-style 7-valued tensor" }; char _coilKindTypeStrEqv[][AIR_STRLEN_SMALL] = { "scalar", "3color", "7tensor", "" }; int _coilKindTypeValEqv[] = { coilKindTypeScalar, coilKindType3Color, coilKindType7Tensor }; airEnum _coilKindType = { "kind", COIL_KIND_TYPE_MAX, _coilKindTypeStr, NULL, _coilKindTypeDesc, _coilKindTypeStrEqv, _coilKindTypeValEqv, AIR_FALSE }; airEnum * coilKindType = &_coilKindType; --- NEW FILE: methodsCoil.c --- /* teem: Gordon Kindlmann's research software Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include "coil.h" int coilVolumeCheck(const Nrrd *nin, const coilKind *kind) { char me[]="coilVolumeCheck", err[AIR_STRLEN_MED]; int baseDim; if (!(nin && kind)) { sprintf(err, "%s: got NULL pointer", me); biffAdd(COIL, err); return 1; } if (nrrdTypeBlock == nin->type) { sprintf(err, "%s: can only operate on scalar types, not %s", me, airEnumStr(nrrdType, nrrdTypeBlock)); biffAdd(COIL, err); return 1; } baseDim = (1 == kind->valLen ? 0 : 1); if (3 + baseDim != nin->dim) { sprintf(err, "%s: dim of input must be 3+%d (3 + baseDim), not %d", me, baseDim, nin->dim); biffAdd(COIL, err); return 1; } return 0; } coilContext * coilContextNew() { coilContext *cctx; cctx = (coilContext *)calloc(1, sizeof(coilContext)); if (cctx) { cctx->nin = NULL; cctx->radius = coilDefaultRadius; cctx->numThreads = 1; ELL_3V_SET(cctx->spacing, AIR_NAN, AIR_NAN, AIR_NAN); cctx->nvol = NULL; cctx->finished = AIR_FALSE; cctx->task = NULL; cctx->finishBarrier = NULL; cctx->filterBarrier = NULL; cctx->updateBarrier = NULL; } return cctx; } int coilContextAllSet(coilContext *cctx, const Nrrd *nin, const coilKind *kind, const coilMethod *method, int radius, int numThreads, int verbose, double parm[COIL_PARMS_NUM]) { char me[]="coilContextAllSet", err[AIR_STRLEN_MED]; int size[NRRD_DIM_MAX], sx, sy, sz, someExist, allExist, baseDim, pi; double xsp, ysp, zsp; airArray *mop; if (!( cctx && nin && kind && method )) { sprintf(err, "%s: got NULL pointer", me); biffAdd(COIL, err); return 1; } if (coilVolumeCheck(nin, kind)) { sprintf(err, "%s: input volume not usable as %s", me, kind->name); biffAdd(COIL, err); return 1; } if (!( radius >= 1 && numThreads >= 1 )) { sprintf(err, "%s: radius (%d) not >= 1 or numThreads (%d) not >= 1", me, radius, numThreads); biffAdd(COIL, err); return 1; } if (!( AIR_IN_OP(coilMethodTypeUnknown, method->type, coilMethodTypeLast) )) { sprintf(err, "%s: method->type %d not valid", me, method->type); biffAdd(COIL, err); return 1; } if (!kind->filter[method->type]) { sprintf(err, "%s: sorry, %s filtering not available on %s kind", me, method->name, kind->name); biffAdd(COIL, err); return 1; } /* warn if we can't do the multiple threads user wants */ if (numThreads > 1 && !airThreadCapable && airThreadNoopWarning) { fprintf(stderr, "%s: WARNING: this teem not thread capable: using 1 " "thread, not %d\n", me, numThreads); numThreads = 1; } mop = airMopNew(); /* set parms */ for (pi=0; pi<method->numParm; pi++) { if (!AIR_EXISTS(parm[pi])) { sprintf(err, "%s: parm[%d] (need %d) doesn't exist", me, pi, method->numParm); biffAdd(COIL, err); airMopError(mop); return 1; } cctx->parm[pi] = parm[pi]; } /* set sizes and spacings */ baseDim = (1 == kind->valLen ? 0 : 1); sx = nin->axis[0 + baseDim].size; sy = nin->axis[1 + baseDim].size; sz = nin->axis[2 + baseDim].size; if (sz < numThreads) { fprintf(stderr, "%s: wanted %d threads but volume only has %d slices, " "using %d threads instead\n", me, numThreads, sz, sz); numThreads = sz; } ELL_3V_SET(cctx->size, sx, sy, sz); xsp = nin->axis[0 + baseDim].spacing; ysp = nin->axis[1 + baseDim].spacing; zsp = nin->axis[2 + baseDim].spacing; someExist = AIR_EXISTS(xsp) || AIR_EXISTS(ysp) || AIR_EXISTS(zsp); allExist = AIR_EXISTS(xsp) && AIR_EXISTS(ysp) && AIR_EXISTS(zsp); if (!( someExist )) { fprintf(stderr, "%s: WARNING: assuming unit spacing for all axes\n", me); xsp = 1; ysp = 1; zsp = 1; } else { if ( !allExist ) { sprintf(err, "%s: spacings (%g,%g,%g) not uniformly existant", me, xsp, ysp, zsp); biffAdd(COIL, err); airMopError(mop); return 1; } } ELL_3V_SET(cctx->spacing, xsp, ysp, zsp); /* allocate nvol */ if (0 == baseDim) { ELL_4V_SET(size, 2, sx, sy, sz); } else { ELL_5V_SET(size, kind->valLen, 2, sx, sy, sz); } cctx->nvol = nrrdNew(); if (nrrdMaybeAlloc_nva(cctx->nvol, coil_nrrdType, 4 + baseDim, size)) { sprintf(err, "%s: couldn't allocate internal processing volume", me); biffMove(COIL, err, NRRD); airMopError(mop); return 1; } airMopAdd(mop, cctx->nvol, (airMopper)nrrdNuke, airMopOnError); cctx->nin = nin; cctx->kind = kind; cctx->method = method; cctx->radius = radius; cctx->numThreads = numThreads; cctx->verbose = verbose; airMopOkay(mop); return 0; } /* ******** coilOutputGet ** ** slice the present intermediate volume to get an output. ** ** No, this does not do quantization or rounding to match the input ** type (of cctx->nin). The reason is that after filtering, it is often ** the case that subtle differences in values emerge, and it may be ** reckless to dump them back into the limited type or value range ** that they started with. That sort of operation should be under ** explicit user control. */ int coilOutputGet(Nrrd *nout, coilContext *cctx) { char me[]="coilOutputGet", err[AIR_STRLEN_MED]; if (!(nout && cctx)) { sprintf(err, "%s: got NULL pointer", me); biffAdd(COIL, err); return 1; } if (nrrdSlice(nout, cctx->nvol, 0, 0)) { sprintf(err, "%s: trouble slicing to get output", me); biffMove(COIL, err, NRRD); return 1; } return 0; } coilContext * coilContextNix(coilContext *cctx) { if (cctx) { cctx->nvol = nrrdNuke(cctx->nvol); cctx = airFree(cctx); } return NULL; } |