[Teem-commit] teem/src/coil coil.h,NONE,1.1 GNUmakefile,NONE,1.1 coreCoil.c,NONE,1.1 defaultsCoil.c,NONE,1.1 enumsCoil.c,NONE,1.1 methodsCoil.c,NONE,1.1 realmethods.c,NONE,1.1 scalarCoil.c,NONE,1.1 tensorCoil.c,NONE,1.1

[Gordon K. <kin...@us...>]

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;
}