itkfemutility.cxx

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMUtility.cxx,v $
  Language:  C++
  Date:      $Date: 2003-09-10 14:29:44 $
  Version:   $Revision: 1.8 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

// disable debug warnings in MS compiler
#ifdef _MSC_VER
#pragma warning(disable: 4786)
#endif

#include "itkFEMUtility.h"

namespace itk {
namespace fem {




/**
 * Numerical integration (Gauss-Legendre formula).
 * Integrates function f(x) from x=a to x=b in n points.
 */
double GaussIntegrate::Integrate(double (*f)(double), double a, double b, int n)
{

/*
  This subprogram produces the gauss-legendre numerical
  integral for
  
        b
    int   f(x)*dx
        a
  
  The number of nodes n must satisfy
      2 <= n <= 20
  The following data statements contain the gauss-legendre
  nodes for the interval [-1,1].
*/

  double scale, t, tl, tu, sum;
  int i, m, ibase;


  /*  Begin integration  */
      
  scale = (b - a)/two;
  if ( (n&1)==0 )  
  {
    m = n/2;
    ibase = m*m;
    sum = zero;
  }
  else  
  {
    m = (n - 1)/2;
    ibase = (n*n - 1)/4;
    sum = w[ibase+m]*(*f)((a + b)/two);
  }
        
  for (i=1; i <= m; i++)  
  {
    t = z[ibase + i - 1];
    tl = (a*(one + t) + (one - t)*b)/two;
    tu = (a*(one - t) + (one + t)*b)/two;
    sum = sum + w[ibase + i - 1]*((*f)(tl)  +(*f)(tu));
  }

  return scale*sum;

}


const double GaussIntegrate::zero = 0.0;
const double GaussIntegrate::one = 1.0;
const double GaussIntegrate::two = 2.0;
const double GaussIntegrate::z[110] = {
  0.0, 0.577350269189626, 0.774596669241483, 0.0,
  0.861136311594053, 0.339981043584856, 0.906179845938664,
  0.538469310105683, 0.0, 0.932469514203152,
  0.661209386466265, 0.238619186083197, 0.949107912342759,
  0.741531185599394, 0.405845151377397, 0.0,
  0.960289856497536, 0.796666477413627, 0.525532409916329,
  0.183434642495650, 0.968160239507626, 0.836031107326636,
  0.613371432700590, 0.324253423403809, 0.0,
  0.973906528517172, 0.865063366688985, 0.679409568299024,
  0.433395394129247, 0.148874338981631, 0.978228658146057,
  0.887062599768095, 0.730152005574049, 0.519096129206812,
  0.269543155952345, 0.0, 0.981560634246719,
  0.904117256370475, 0.769902674194305, 0.587317954286617,
  0.367831498998180, 0.125233408511469, 0.984183054718588,
  0.917598399222978, 0.801578090733310, 0.642349339440340,
  0.448492751036447, 0.230458315955135, 0.0, 
  0.986283808696812, 0.928434883663574, 0.827201315069765,
  0.687292904811685, 0.515248636358154, 0.319112368927890,
  0.108054948707344, 0.987992518020485, 0.937273392400706,
  0.848206583410427, 0.724417731360170, 0.570972172608539,
  0.394151347077563, 0.201194093997435, 0.0, 
  0.989400934991650, 0.944575023073233, 0.865631202387832,
  0.755404408355003, 0.617876244402644, 0.458016777657227,
  0.281603550779259, 0.0950125098376374,0.990575475314417,
  0.950675521768768, 0.880239153726986, 0.781514003896801,
  0.657671159216691, 0.512690537086477, 0.351231763453876,
  0.178484181495848, 0.0, 0.991565168420931, 
  0.955823949571398, 0.892602466497556, 0.803704958972523,
  0.691687043060353, 0.559770831073948, 0.411751161462843,
  0.251886225691506, 0.0847750130417353,0.992406843843584,
  0.960208152134830, 0.903155903614818, 0.822714656537143,
  0.720966177335229, 0.600545304661681, 0.464570741375961,
  0.316564099963630, 0.160358645640225, 0.0, 
  0.993128599185095, 0.963971927277914, 0.912234428251326,
  0.839116971822219, 0.746331906460151, 0.636053680726515,
  0.510867001950827, 0.373706088715420, 0.227785851141645,
  0.0765265211334973
};

const double GaussIntegrate::w[110] = {
  0.0, 1.0, 0.555555555555556, 0.888888888888889,
  0.347854845137454, 0.652145154862546, 0.236926885056189,
  0.478628670499366, 0.568888888888889, 0.171324492379170,
  0.360761573048139, 0.467913934572691, 0.129484966168870,
  0.279705391489277, 0.381830050505119, 0.417959183673469,
  0.101228536290376, 0.222381034453374, 0.313706645877887,
  0.362683783378362, 0.0812743883615744, 0.180648160694857,
  0.260610696402935, 0.312347077040003, 0.330239355001260,
  0.0666713443086881, 0.149451349150581, 0.219086362515982,
  0.269266719309996, 0.295524224714753, 0.0556685671161737,
  0.125580369464905, 0.186290210927734, 0.233193764591990,
  0.262804544510247, 0.272925086777901, 0.0471753363865118,
  0.106939325995318, 0.160078328543346, 0.203167426723066,
  0.233492536538355, 0.249147045813403, 0.0404840047653159,
  0.0921214998377284, 0.138873510219787, 0.178145980761946,
  0.207816047536889, 0.226283180262897, 0.232551553230874,
  0.0351194603317519,0.0801580871597602, 0.121518570687903,
   0.157203167158194, 0.185538397477938, 0.205198463721296,
   0.215263853463158, 0.0307532419961173,0.0703660474881081,
   0.107159220467172, 0.139570677926154, 0.166269205816994,
  0.186161000015562, 0.198431485327112, 0.202578241925561,
   0.0271524594117541,0.0622535239386478,0.0951585116824928,
   0.124628971255534, 0.149595988816577, 0.169156519395003,
  0.182603415044924, 0.189450610455068, 0.0241483028685479,
   0.0554595293739872, 0.0850361483171792,0.111883847193404,
   0.135136368468525, 0.154045761076810, 0.168004102156450,
   0.176562705366993, 0.179446470356207, 0.0216160135264833,
   0.0497145488949698, 0.0764257302548891,0.100942044106287,
   0.122555206711478, 0.140642914670651, 0.154684675126265,
   0.164276483745833, 0.169142382963144, 0.0194617882297265,
   0.0448142267656996,0.0690445427376412,0.0914900216224500,
   0.111566645547334, 0.128753962539336, 0.142606702173607,
   0.152766042065860, 0.158968843393954, 0.161054449848784,
   0.0176140071391521,0.0406014298003869,0.0626720483341091,
   0.0832767415767047, 0.101930119817240, 0.118194531961518,
   0.131688638449177, 0.142096109318382, 0.149172986472604,
    0.152753387130726
};

}} // end namespace itk::fem