#include "PmeBase.h"#include <math.h>Go to the source code of this file.
Functions | |
| void | compute_b_spline (double *frac, double *M, double *dM, int order) |
|
||||||||||||||||||||
|
Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved. Definition at line 10 of file PmeBase.C. References order. 00010 {
00011 int j, n;
00012 double x,y,z,x1,y1,z1, div;
00013 double *Mx, *My, *Mz, *dMx, *dMy, *dMz;
00014 Mx=M-1; My=M+order-1; Mz=M+2*order-1;
00015 dMx=dM-1; dMy =dM+order-1; dMz=dM+2*order-1;
00016 x=frac[0];
00017 y=frac[1];
00018 z=frac[2];
00019 x1=1.0-x; y1=1.0-y; z1=1.0-z;
00020 /* Do n=3 case first */
00021 Mx[1]=.5*x1*x1;
00022 Mx[2]=x1*x + .5;
00023 Mx[3]=0.5*x*x;
00024 Mx[order]=0.0;
00025 My[1]=.5*y1*y1;
00026 My[2]=y1*y + .5;
00027 My[3]=0.5*y*y;
00028 My[order]=0.0;
00029 Mz[1]=.5*z1*z1;
00030 Mz[2]=z1*z + .5;
00031 Mz[3]=0.5*z*z;
00032 Mz[order]=0.0;
00033
00034 /* Recursively fill in the rest. */
00035 for (n=4; n<=order-1; n++) {
00036 double div=1.0/(n-1);
00037 int j;
00038 Mx[n] = x*div*Mx[n-1];
00039 My[n] = y*div*My[n-1];
00040 Mz[n] = z*div*Mz[n-1];
00041 for (j=1; j<=n-2; j++) {
00042 Mx[n-j] = ((x+j)*Mx[n-j-1] + (n-x-j)*Mx[n-j])*div;
00043 My[n-j] = ((y+j)*My[n-j-1] + (n-y-j)*My[n-j])*div;
00044 Mz[n-j] = ((z+j)*Mz[n-j-1] + (n-z-j)*Mz[n-j])*div;
00045 }
00046 Mx[1] *= (1.0-x)*div;
00047 My[1] *= (1.0-y)*div;
00048 Mz[1] *= (1.0-z)*div;
00049 }
00050 /* Now do the derivatives */
00051 dMx[1]=-Mx[1]; dMy[1]=-My[1]; dMz[1]=-Mz[1];
00052 for (j=2; j <= order; j++) {
00053 dMx[j] = Mx[j-1] - Mx[j];
00054 dMy[j] = My[j-1] - My[j];
00055 dMz[j] = Mz[j-1] - Mz[j];
00056 }
00057 /* Now finish the job! */
00058 div=1.0/(order-1);
00059 Mx[order] = x*div*Mx[order-1];
00060 My[order] = y*div*My[order-1];
00061 Mz[order] = z*div*Mz[order-1];
00062 for (j=1; j<=order-2; j++) {
00063 Mx[order-j] = ((x+j)*Mx[order-j-1] + (order-x-j)*Mx[order-j])*div;
00064 My[order-j] = ((y+j)*My[order-j-1] + (order-y-j)*My[order-j])*div;
00065 Mz[order-j] = ((z+j)*Mz[order-j-1] + (order-z-j)*Mz[order-j])*div;
00066 }
00067 Mx[1] *= (1.0-x)*div;
00068 My[1] *= (1.0-y)*div;
00069 Mz[1] *= (1.0-z)*div;
00070 }
|
1.3.9.1