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...
#	Content
1	/***************************************************************************
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	*/