/*

  Teem: Tools to process and visualize scientific data and images             .

  Copyright (C) 2013, 2012, 2011, 2010, 2009  University of Chicago

  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann

  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

  (LGPL) as published by the Free Software Foundation; either

  version 2.1 of the License, or (at your option) any later version.

  The terms of redistributing and/or modifying this software also

  include exceptions to the LGPL that facilitate static linking.

  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 Free Software Foundation, Inc.,

  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

*/

#include "gage.h"

#include "privateGage.h"

double

gageStackWtoI(gageContext *ctx, double swrl, int *outside) {

  double si;

    unsigned int sidx;

      /* we'll extrapolate from stackPos[0] and [1] */

      sidx = 0;

      /* extrapolate from stackPos[ctx->pvlNum-3] and [ctx->pvlNum-2];

         gageStackPerVolumeAttach ensures that we there are at least two

         blurrings pvls & one base pvl ==> pvlNum >= 3 ==> pvlNum-3 >= 0 */

      /* HEY: stupid linear search */

          break;

        }

      }

        /* search failure */

      }

      *outside = AIR_FALSE;

    }

                    sidx, sidx+1);

  }

}

double

gageStackItoW(gageContext *ctx, double si, int *outside) {

  unsigned int sidx;

  double swrl, sfrac;

      sidx = 0;

    }

    /*

    fprintf(stderr, "!%s: si %g (%u) --> %u + %g --> [%g,%g] -> %g\n", me,

            si, ctx->pvlNum, sidx, sfrac,

            ctx->stackPos[sidx], ctx->stackPos[sidx+1], swrl);

    */

  }

}

/*

** this is a little messy: the given pvlStack array has to be allocated

** by the caller to hold blNum gagePerVolume pointers, BUT, the values

** of pvlStack[i] shouldn't be set to anything: as with gagePerVolumeNew(),

** gage allocates the pervolume itself.

*/

int

gageStackPerVolumeNew(gageContext *ctx,

                      gagePerVolume **pvlStack,

                      const Nrrd *const *nblur, unsigned int blNum,

                      const gageKind *kind) {

  static const char me[]="gageStackPerVolumeNew";

  unsigned int blIdx;

  }

  }

    }

  }

}

/*

** the "base" pvl is the LAST pvl, ctx->pvl[pvlNum-1]

*/

int

gageStackPerVolumeAttach(gageContext *ctx, gagePerVolume *pvlBase,

                         gagePerVolume **pvlStack, const double *stackPos,

                         unsigned int blNum) {

  static const char me[]="gageStackPerVolumeAttach";

  unsigned int blIdx;

  }

    /* this constraint is important for the logic of stack reconstruction:

       minimum number of node-centered samples is 2, and the number of

       pvls has to be at least 3 (two blurrings + one base pvl) */

  }

             "prior to stack attachment", me, ctx->pvlNum);

  }

               stackPos[blIdx]);

    }

                 blIdx, stackPos[blIdx], blIdx+1, stackPos[blIdx+1]);

      }

    }

  }

  /* the base volume is LAST, after all the stack samples */

    }

  }

  }

  }

  }

}

/*

** _gageStackBaseIv3Fill

**

** after the individual iv3's in the stack have been filled, this does

** the across-stack filtering to fill the iv3 of pvl[pvlNum-1] (the

** "base" pvl)

*/

int

_gageStackBaseIv3Fill(gageContext *ctx) {

  static const char me[]="_gageStackBaseIv3Fill";

  unsigned int fd, pvlIdx, cacheIdx, cacheLen, baseIdx, valLen;

  /* the "base" pvl is the LAST pvl */

  }

    unsigned int iii, xi, yi, zi, blurIdx, valIdx, fdd, sz, sy;

    double xx, *iv3, *iv30, *iv31, sigma0, sigma1,

      val0, val1, drv0, drv1, lapl0, lapl1;

    fdd = fd*fd;

    /* initialize the output iv3 to all zeros, since we won't be

       usefully setting the values on the boundary (the boundary which

       is required in the rest of the stack's iv3s in order to do the

       laplacian-based spline recon), and we can't have any

       non-existent values creeping in.  We shouldn't need to do any

       kind of nrrdBoundaryBleed thing here, because the kernel

       weights really should be zero on the boundary. */

    }

    /* find the interval in the pre-blurred volumes containing the

       desired scale location */

        /* has to be non-zero somewhere, since _gageLocationSet()

           gives an error if there aren't non-zero stackFw[i] */

        break;

      }

    }

    /* so no way that pvlIdx == pvlNum-1 */

      /* pvlNum-2 is pvl index of last pre-blurred volume */

      /* gageStackPerVolumeAttach() enforces getting at least two

         pre-blurred volumes --> pvlNum >= 3 --> blurIdx >= 0 */

      xx = 1;

      blurIdx = pvlIdx;

      /* by design, the hermite non-kernel generates the same values as

         the tent kernel (with scale forced == 1), so we can use those

         to control the interpolation */

    }

        /* (1-D data: HEY copy and paste; see 2D case below) */

          /* nixed "for zi" and "for yi" loop; zi==yi==1 */

            /* can do a 1D instead of 2D discrete laplacian */

            /*

            fprintf(stderr, "!%s: (%u): val %g %g, lapl %g %g, sigma %g %g, drv %g %g --> iv3[%u] = %g\n", me,

                    xi, val0, val1, lapl0, lapl1, sigma0, sigma1, drv0, drv1, iii, iv3[iii]);

            fprintf(stderr, "!%s: lapl0 %g = %g + %g + %g + %g - 4*%g\n", me, lapl0,

                    iv30[iii+1] , iv30[iii-1] , iv30[iii+fd] , iv30[iii-fd] , val0);

            fprintf(stderr, "!%s: lapl1 %g = %g + %g + %g + %g - 4*%g\n", me, lapl1,

                    iv31[iii+1] , iv31[iii-1] , iv31[iii+fd] , iv31[iii-fd] , val1);

            */

          }

              }

            }

          }

        }

      } else {

        /* as in gageIv3Fill; we do some special-case-ing for 2-D images;

           (HEY copy and paste; see sz > 1 below for explanatory comments) */

          /* nixed "for zi" loop; zi==1 */

              /* can do a 2D instead of 3D discrete laplacian */

              /*

              fprintf(stderr, "!%s: (%u,%u): val %g %g, lapl %g %g, sigma %g %g, drv %g %g --> iv3[%u] = %g\n", me,

                      xi, yi, val0, val1, lapl0, lapl1, sigma0, sigma1, drv0, drv1, iii, iv3[iii]);

              fprintf(stderr, "!%s: lapl0 %g = %g + %g + %g + %g - 4*%g\n", me, lapl0,

                      iv30[iii+1] , iv30[iii-1] , iv30[iii+fd] , iv30[iii-fd] , val0);

              fprintf(stderr, "!%s: lapl1 %g = %g + %g + %g + %g - 4*%g\n", me, lapl1,

                      iv31[iii+1] , iv31[iii-1] , iv31[iii+fd] , iv31[iii-fd] , val1);

              */

            }

          }

              }

            }

          }

        }

      }

    } else {

      /* sz > 1 */

              /* note that iv3 axis ordering is x, y, z, tuple */

              /* the (sigma1 - sigma0) factor is needed to convert the

                 derivative with respect to sigma (sigma*lapl) into the

                 derivative with respect to xx (ranges from 0 to 1) */

              /* This inner loop is the bottleneck for some uses of

                 scale-space; a re-arrangement of the Hermite spline

                 evaluation (thanks Mathematica) does save a little time */

            }

          }

        }

      }

    }

    /* we're doing simple convolution-based recon on the stack */

    /* NOTE we are treating the 4D fd*fd*fd*valLen iv3 as a big 1-D array */

    double wght, val;

      val = 0;

                : 0);

      }

    }

  }

}

/*

******** gageStackProbe()

*/

int

gageStackProbe(gageContext *ctx,

               double xi, double yi, double zi, double stackIdx) {

  static const char me[]="gageStackProbe";

  }

    }

  }

}

int

gageStackProbeSpace(gageContext *ctx,

                    double xx, double yy, double zz, double ss,

                    int indexSpace, int clamp) {

  static const char me[]="gageStackProbeSpace";

  int ret;

  /*

  fprintf(stderr, "!%s(%g,%g,%g,%g, %d,%d): hello\n", me,

          xx, yy, zz, ss, indexSpace, clamp);

  if (AIR_ABS(-98.2539 - xx) < 0.1 && AIR_ABS(yy) < 0.1 && AIR_ABS(zz) < 0.1 && AIR_ABS(495.853 - ss)) {

    ctx->verbose += 10;

  }

  */

  }

    }

  }

  /* fprintf(stderr, "!%s: returning %d\n", me, ret); */

}