trunk/SHAPES/makeweights.c

/***************************************************************************
  **************************************************************************
  
                           S2kit 1.0

          A lite version of Spherical Harmonic Transform Kit

   Peter Kostelec, Dan Rockmore
   {geelong,rockmore}@cs.dartmouth.edu
  
   Contact: Peter Kostelec
            geelong@cs.dartmouth.edu
  
   Copyright 2004 Peter Kostelec, Dan Rockmore

   This file is part of S2kit.

   S2kit is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   S2kit 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 General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with S2kit; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

   See the accompanying LICENSE file for details.
  
  ************************************************************************
  ************************************************************************/


/*
  source file which contains the function that generates the
  weights for a bandwidth bw legendre transform. Basically,
  it contains the implementation of the formula as defined in
  the tensor paper, and also given in the so(3) paper. It's
  just mentioned in the s^2 paper!

  This formula is slightly different from the one given in the
  original DH paper because they were sampling at the poles in
  that paper, and now we're not.

  In pseudo-TeX, the formula for the bandwidth B weights is

  w_B(j) = 2/B sin((pi*(2j+1))/(4B)) *
           sum_{k=0}^{B-1} 1/(2k+1)*sin((2j+1)(2k+1)pi/(4B))

  where j = 0, 1, ..., 2 B - 1


  Note that if you want to use these weights for an *odd*
  order transform, given the way the code is set up, you
  have to MULTIPLY the j-th weight by sin(pi*(2j+1)/(4B))

*/

#include <math.h>

/*
  makeweights: given a bandwidth bw, make weights for
  both even *and* odd order Legendre transforms.

  bw -> bandwidth of transform
  weights -> pointer to double array of length 4*bw; this
             array will contain the even and odd weights;
             even weights start at weight[0], and odd weights
             start at weights[2*bw]

*/

void makeweights( int bw,
                  double *weights )
{
  int j, k ;
  double fudge ;
  double tmpsum ;

  fudge = M_PI/((double)(4*bw)) ;
  

  for ( j = 0 ; j < 2*bw ; j ++ )
    {
      tmpsum = 0.0 ;
      for ( k = 0 ; k < bw ; k ++ )
        tmpsum += 1./((double)(2*k+1)) *
          sin((double)((2*j+1)*(2*k+1))*fudge);
      tmpsum *= sin((double)(2*j+1)*fudge);
      tmpsum *= 2./((double) bw) ;
      
      weights[j] = tmpsum ;
      weights[j + 2*bw] = tmpsum * sin((double)(2*j+1)*fudge);
    }

}

/********************************************************/
/********************************************************/
/********************************************************/
/********************************************************/

/* just a hack to test the above function */
/*

#include <stdio.h>
#include <stdlib.h>

int main( int argc,
          char **argv )
{
  int i, bw ;
  double *weights ;


  bw = atoi(argv[1]);

  weights = (double *)malloc(sizeof(double)*4*bw);

  makeweights(bw, weights);

  for ( i = 0 ; i < 2*bw ; i ++ )
    printf("%d\t%f\t%f\n",
           i, weights[i], weights[i+2*bw]);

  free( weights );

  return 0 ;

}

*/
Revision:	1287
Committed:	Wed Jun 23 20:18:48 2004 UTC (21 years, 4 months ago) by chrisfen
Content type:	text/plain
File size:	3690 byte(s)
Log Message:	Major progress towards inclusion of spherical harmonic transform capability - still having some build issues...
#	User	Rev	Content
1	chrisfen	1287	/***************************************************************************
2			**************************************************************************
3
4			S2kit 1.0
5
6			A lite version of Spherical Harmonic Transform Kit
7
8			Peter Kostelec, Dan Rockmore
9			{geelong,rockmore}@cs.dartmouth.edu
10
11			Contact: Peter Kostelec
12			geelong@cs.dartmouth.edu
13
14			Copyright 2004 Peter Kostelec, Dan Rockmore
15
16			This file is part of S2kit.
17
18			S2kit is free software; you can redistribute it and/or modify
19			it under the terms of the GNU General Public License as published by
20			the Free Software Foundation; either version 2 of the License, or
21			(at your option) any later version.
22
23			S2kit is distributed in the hope that it will be useful,
24			but WITHOUT ANY WARRANTY; without even the implied warranty of
25			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26			GNU General Public License for more details.
27
28			You should have received a copy of the GNU General Public License
29			along with S2kit; if not, write to the Free Software
30			Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
31
32			See the accompanying LICENSE file for details.
33
34			************************************************************************
35			************************************************************************/
36
37
38			/*
39			source file which contains the function that generates the
40			weights for a bandwidth bw legendre transform. Basically,
41			it contains the implementation of the formula as defined in
42			the tensor paper, and also given in the so(3) paper. It's
43			just mentioned in the s^2 paper!
44
45			This formula is slightly different from the one given in the
46			original DH paper because they were sampling at the poles in
47			that paper, and now we're not.
48
49			In pseudo-TeX, the formula for the bandwidth B weights is
50
51			w_B(j) = 2/B sin((pi(2j+1))/(4B))
52			sum_{k=0}^{B-1} 1/(2k+1)*sin((2j+1)(2k+1)pi/(4B))
53
54			where j = 0, 1, ..., 2 B - 1
55
56
57			Note that if you want to use these weights for an odd
58			order transform, given the way the code is set up, you
59			have to MULTIPLY the j-th weight by sin(pi*(2j+1)/(4B))
60
61			*/
62
63			#include <math.h>
64
65			/*
66			makeweights: given a bandwidth bw, make weights for
67			both even and odd order Legendre transforms.
68
69			bw -> bandwidth of transform
70			weights -> pointer to double array of length 4*bw; this
71			array will contain the even and odd weights;
72			even weights start at weight[0], and odd weights
73			start at weights[2*bw]
74
75			*/
76
77			void makeweights( int bw,
78			double *weights )
79			{
80			int j, k ;
81			double fudge ;
82			double tmpsum ;
83
84			fudge = M_PI/((double)(4*bw)) ;
85
86
87			for ( j = 0 ; j < 2*bw ; j ++ )
88			{
89			tmpsum = 0.0 ;
90			for ( k = 0 ; k < bw ; k ++ )
91			tmpsum += 1./((double)(2k+1))
92			sin((double)((2j+1)(2k+1))fudge);
93			tmpsum = sin((double)(2j+1)*fudge);
94			tmpsum *= 2./((double) bw) ;
95
96			weights[j] = tmpsum ;
97			weights[j + 2bw] = tmpsum sin((double)(2j+1)fudge);
98			}
99
100			}
101
102			/********************************************************/
103			/********************************************************/
104			/********************************************************/
105			/********************************************************/
106
107			/* just a hack to test the above function */
108			/*
109
110			#include <stdio.h>
111			#include <stdlib.h>
112
113			int main( int argc,
114			char **argv )
115			{
116			int i, bw ;
117			double *weights ;
118
119
120			bw = atoi(argv[1]);
121
122			weights = (double )malloc(sizeof(double)4*bw);
123
124			makeweights(bw, weights);
125
126			for ( i = 0 ; i < 2*bw ; i ++ )
127			printf("%d\t%f\t%f\n",
128			i, weights[i], weights[i+2*bw]);
129
130			free( weights );
131
132			return 0 ;
133
134			}
135
136			*/