<commit_before>/*
  Teem: Tools to process and visualize scientific data and images              
  Copyright (C) 2011, 2010, 2009, 2008 <NAME>

  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 "elf.h"

#if TEEM_LEVMAR
#include <lm.h>
#endif

/* Routines for estimating the ball-and-multi-stick model from single-shell
 * DWI data. Implements the ideas from
 *
 * <NAME>, <NAME>, <NAME>: Multi-Diffusion-Tensor
 * Fitting via Spherical Deconvolution: A Unifying Framework. MICCAI
 * 2010, Part I, LNCS 6361, pp. 673-680. Springer, 2010
 */

/* elfKernelStick:
 *
 * Computes the deconvolution kernel corresponding to a single stick, for
 * desired ESH order, bd (b value times diffusivity) and b0 (non-DWI signal)
 * If delta!=0, deconvolves to delta peak, else to rank-1 peak
 *
 * returns 0 on success, 1 if order is not supported
 */
int elfKernelStick_f(float *kernel, unsigned int order, float bd,
                     float b0, int delta) {
  float ebd=exp(bd);
  float embd=exp(-bd);
  float sbd=sqrt(bd);
  float erfsbd=airErf(sbd);
  float spi=sqrt(AIR_PI);
  kernel[0]=b0*AIR_PI*erfsbd/sbd;
  if (order>=2) {
    kernel[1]=-b0/(4.0*bd*sbd)*embd*sqrt(5*AIR_PI)*
      (6*sbd+(-3+2*bd)*ebd*spi*erfsbd);
    if (order>=4) {
      kernel[2]=b0/(32.0*bd*bd*sbd)*embd*spi*
	(-30*sbd*(21+2*bd)+9*(35+4*bd*(-5+bd))*ebd*spi*erfsbd);
      if (order>=6) { /* At order 6, noise starts to take over! */
	kernel[3]=b0/(128*bd*bd*bd*sbd)*embd*sqrt(13.0)*
	  (-42*sbd*(165+4*bd*(5+bd))*spi-
	   5*(-693+378*bd-84*bd*bd+8*bd*bd*bd)*
	   ebd*AIR_PI*erfsbd);
	if (order>=8) {
	  kernel[4]=b0/(2048*bd*bd*bd*bd*sbd)*embd*sqrt(17.0)*
	    (-6*sbd*(225225+2*bd*(15015+2*bd*(1925+62*bd)))*spi+
	     35*(19305+8*bd*(-1287+bd*(297+2*(-18+bd)*bd)))*
	     ebd*AIR_PI*erfsbd);
          if (order>8)
            return 1;
	}
      }
    }
  }
  if (delta) {
    return tijk_esh_make_kernel_delta_f(kernel,kernel,order);
  }
  return tijk_esh_make_kernel_rank1_f(kernel,kernel,order);
}

/* elfBallStickODF:
 *
 * Deconvolves single-shell DWI measurements to an ODF, using a kernel derived
 * from the ball-and-stick model
 *
 * Output:
 *  odf   - ESH coefficients of computed ODF
 *  fiso  - if non-NULL, estimated isotropic volume fraction
 *  d     - if non-NULL, estimated ADC
 * Input:
 *  dwi   - DWI measurement and parameters
 *  T     - matrix computed by elfESHEstimMatrix
 *  order - desired ESH order of odf
 *  delta - whether to use delta peak (set to 0 when using elfBallStickPredict)
 *
 * Returns 0 on success, 1 if order is not supported
 */
int elfBallStickODF_f(float *odf, float *fiso, float *d,
                      const elfSingleShellDWI *dwi,
                      const float *T, unsigned int order, int delta)
{
  unsigned int C = tijk_esh_len[order/2], k, l;
  float mean=0, _d=0;
  float *odfp = odf;
  float isovf0, isovf1, _fiso;
  float kernel[10]; /* should be tijk_max_esh_order/2+1,
                     * but ANSI C doesn't allow that :/ */

  if (order>tijk_max_esh_order || order%2!=0)
    return 1;
<commit_msg>Working around a small numerical glitch in ball+stick code.<commit_after>/*
  Teem: Tools to process and visualize scientific data and images              
  Copyright (C) 2011, 2010, 2009, 2008 <NAME>

  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 "elf.h"

#if TEEM_LEVMAR
#include <lm.h>
#endif

/* Routines for estimating the ball-and-multi-stick model from single-shell
 * DWI data. Implements the ideas from
 *
 * <NAME>, <NAME>, <NAME>: Multi-Diffusion-Tensor
 * Fitting via Spherical Deconvolution: A Unifying Framework. MICCAI
 * 2010, Part I, LNCS 6361, pp. 673-680. Springer, 2010
 */

/* elfKernelStick:
 *
 * Computes the deconvolution kernel corresponding to a single stick, for
 * desired ESH order, bd (b value times diffusivity) and b0 (non-DWI signal)
 * If delta!=0, deconvolves to delta peak, else to rank-1 peak
 *
 * returns 0 on success, 1 if order is not supported
 */
int elfKernelStick_f(float *kernel, unsigned int order, float bd,
                     float b0, int delta) {
  double ebd=exp(bd);
  double embd=exp(-bd);
  double sbd=sqrt(bd);
  double erfsbd=airErf(sbd);
  double spi=sqrt(AIR_PI);
  kernel[0]=b0*AIR_PI*erfsbd/sbd;
  if (order>=2) {
    kernel[1]=-b0/(4.0*bd*sbd)*embd*sqrt(5*AIR_PI)*
      (6*sbd+(-3+2*bd)*ebd*spi*erfsbd);
    if (order>=4) {
      kernel[2]=b0/(32.0*bd*bd*sbd)*embd*spi*
	(-30*sbd*(21+2*bd)+9*(35+4*bd*(-5+bd))*ebd*spi*erfsbd);
      if (order>=6) { /* At order 6, noise starts to take over! */
	kernel[3]=b0/(128*bd*bd*bd*sbd)*embd*sqrt(13.0)*
	  (-42*sbd*(165+4*bd*(5+bd))*spi-
	   5*(-693+378*bd-84*bd*bd+8*bd*bd*bd)*
	   ebd*AIR_PI*erfsbd);
	if (order>=8) {
	  kernel[4]=b0/(2048*bd*bd*bd*bd*sbd)*embd*sqrt(17.0)*
	    (-6*sbd*(225225+2*bd*(15015+2*bd*(1925+62*bd)))*spi+
	     35*(19305+8*bd*(-1287+bd*(297+2*(-18+bd)*bd)))*
	     ebd*AIR_PI*erfsbd);
          if (order>8)
            return 1;
	}
      }
    }
  }
  if (delta) {
    return tijk_esh_make_kernel_delta_f(kernel,kernel,order);
  }
  return tijk_esh_make_kernel_rank1_f(kernel,kernel,order);
}

/* elfBallStickODF:
 *
 * Deconvolves single-shell DWI measurements to an ODF, using a kernel derived
 * from the ball-and-stick model
 *
 * Output:
 *  odf   - ESH coefficients of computed ODF
 *  fiso  - if non-NULL, estimated isotropic volume fraction
 *  d     - if non-NULL, estimated ADC
 * Input:
 *  dwi   - DWI measurement and parameters
 *  T     - matrix computed by elfESHEstimMatrix
 *  order - desired ESH order of odf
 *  delta - whether to use delta peak (set to 0 when using elfBallStickPredict)
 *
 * Returns 0 on success, 1 if order is not supported
 */
int elfBallStickODF_f(float *odf, float *fiso, float *d,
                      const elfSingleShellDWI *dwi,
                      const float *T, unsigned int order, int delta)
{
  unsigned int C = tijk_esh_len[order/2], k, l;
  float mean=0, _d=1e-20; /* small epsilon to avoid potential divide by zero */
  float *odfp = odf;
  float isovf0, isovf1, _fiso;
  float kernel[10]; /* should be tijk_max_esh_order/2+1,
                     * but ANSI C doesn't allow that :/ */

  if (order>tijk_max_esh_order || order%2!=0)
    return 1;

{'fraction_non_alphanumeric': 0.1191723080853238, 'fraction_numerical': 0.04445012134372206, 'mean_word_length': 3.3816452154448795, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3912465', 'n_tokens_mistral': 2976, 'n_tokens_neox': 2727, 'n_words': 860}