trunk/SHAPES/naive_synthesis.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 code to synthesize functions using a naive method
   based on recurrence.  This is slow but does not require any
   precomputed functions, and is also stable. 
*/


#include <math.h>
#include <string.h>


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

   Naive_AnalysisX: computing the discrete Legendre transform of
                    a function via summing naively. I.e. This is
                    the FORWARD discrete Legendre transform.

   bw - bandwidth
   m - order

   data - a pointer to double array of size (2*bw) containing
          the sample points

   result - a pointer to double array of size (bw-m) which, at the
            conclusion of the routine, will contains the coefficients

   plmtable - a pointer to a double array of size (2*bw*(bw-m));
              contains the PRECOMPUTED plms, i.e. associated Legendre
              functions.  E.g. Should be generated by a call to

              PmlTableGen()

              (see pmls.c)

              NOTE that these Legendres are normalized with norm
              equal to 1 !!!

   workspace - array of size 2 * bw;


*/


void Naive_AnalysisX(double *data,
                     int bw,
                     int m,
                     double *weights,
                     double *result,
                     double *plmtable,
                     double *workspace)
{
  int i, j;
  double result0, result1, result2, result3;
  register double *wdata;

  wdata = workspace;

  /* make sure result is zeroed out */
  memset( result, 0, sizeof(double) * (bw - m) );

  /* apply quadrature weights */
  /*
    I only have to differentiate between even and odd
    weights when doing something like seminaive, something
    which involves the dct. In this naive case, the parity of
    the order of the transform doesn't matter because I'm not
    dividing by sin(x) when precomputing the Legendres (because
    I'm not taking their dct). The plain ol' weights are just
    fine. 
  */
  
  for(i = 0; i < 2 * bw; i++)
    wdata[i] = data[i] * weights[i];

  /* unrolling seems to work */
  if ( 1 )
    {
      for (i = 0; i < bw - m; i++)
        {
          result0 = 0.0; result1 = 0.0;
          result2 = 0.0; result3 = 0.0;
          
          for(j = 0; j < (2 * bw) % 4; ++j)
            result0 += wdata[j] * plmtable[j];

          for( ; j < (2 * bw); j += 4)
            {
              result0 += wdata[j] * plmtable[j];
              result1 += wdata[j + 1] * plmtable[j + 1];
              result2 += wdata[j + 2] * plmtable[j + 2];
              result3 += wdata[j + 3] * plmtable[j + 3];
            }
          result[i] = result0 + result1 + result2 + result3;

          plmtable += (2 * bw);
        }
    }
  else
    {
      for (i = 0; i < bw - m; i++)
        {
          result0 = 0.0 ;         
          for(j = 0; j < (2 * bw) ; j++)
            result0 += wdata[j] * plmtable[j];
          result[i] = result0 ;

          plmtable += (2 * bw);
        }
    }
}


/************************************************************************/
/* This is the procedure that synthesizes a function from a list
   of coefficients of a Legendre series. I.e. this is the INVERSE
   discrete Legendre transform.

   Function is synthesized at the (2*bw) Chebyshev nodes. Associated
   Legendre functions are assumed to be precomputed.
   
   bw - bandwidth

   m - order

   coeffs - a pointer to double array of size (bw-m).  First coefficient is
            coefficient for Pmm
   result - a pointer to double array of size (2*bw); at the conclusion
            of the routine, this array will contain the
            synthesized function

   plmtable - a pointer to a double array of size (2*bw*(bw-m));
              contains the PRECOMPUTED plms, i.e. associated Legendre
              functions. E.g. Should be generated by a call to

              PmlTableGen(),

              (see pmls.c)


              NOTE that these Legendres are normalized with norm
              equal to 1 !!!
*/

void Naive_SynthesizeX(double *coeffs,
                       int bw,
                       int m,
                       double *result,
                       double *plmtable)
{
  int i, j;
  double tmpcoef;

  /* make sure result is zeroed out */
  memset( result, 0, sizeof(double) * 2 * bw );

  for ( i = 0 ; i < bw - m ; i ++ )
    {
      tmpcoef = coeffs[i] ;
      if ( tmpcoef != 0.0 )
        for (j=0; j<(2*bw); j++)
          result[j] += (tmpcoef * plmtable[j]);
      plmtable += (2 * bw ) ;
    }
}
Revision:	1287
Committed:	Wed Jun 23 20:18:48 2004 UTC (21 years, 4 months ago) by chrisfen
Content type:	text/plain
File size:	5589 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
40	/* Source code to synthesize functions using a naive method
41	based on recurrence. This is slow but does not require any
42	precomputed functions, and is also stable.
43	*/
44
45
46	#include <math.h>
47	#include <string.h>
48
49
50	/************************************************************************/
51	/*
52
53	Naive_AnalysisX: computing the discrete Legendre transform of
54	a function via summing naively. I.e. This is
55	the FORWARD discrete Legendre transform.
56
57	bw - bandwidth
58	m - order
59
60	data - a pointer to double array of size (2*bw) containing
61	the sample points
62
63	result - a pointer to double array of size (bw-m) which, at the
64	conclusion of the routine, will contains the coefficients
65
66	plmtable - a pointer to a double array of size (2bw(bw-m));
67	contains the PRECOMPUTED plms, i.e. associated Legendre
68	functions. E.g. Should be generated by a call to
69
70	PmlTableGen()
71
72	(see pmls.c)
73
74	NOTE that these Legendres are normalized with norm
75	equal to 1 !!!
76
77	workspace - array of size 2 * bw;
78
79
80	*/
81
82
83	void Naive_AnalysisX(double *data,
84	int bw,
85	int m,
86	double *weights,
87	double *result,
88	double *plmtable,
89	double *workspace)
90	{
91	int i, j;
92	double result0, result1, result2, result3;
93	register double *wdata;
94
95	wdata = workspace;
96
97	/* make sure result is zeroed out */
98	memset( result, 0, sizeof(double) * (bw - m) );
99
100	/* apply quadrature weights */
101	/*
102	I only have to differentiate between even and odd
103	weights when doing something like seminaive, something
104	which involves the dct. In this naive case, the parity of
105	the order of the transform doesn't matter because I'm not
106	dividing by sin(x) when precomputing the Legendres (because
107	I'm not taking their dct). The plain ol' weights are just
108	fine.
109	*/
110
111	for(i = 0; i < 2 * bw; i++)
112	wdata[i] = data[i] * weights[i];
113
114	/* unrolling seems to work */
115	if ( 1 )
116	{
117	for (i = 0; i < bw - m; i++)
118	{
119	result0 = 0.0; result1 = 0.0;
120	result2 = 0.0; result3 = 0.0;
121
122	for(j = 0; j < (2 * bw) % 4; ++j)
123	result0 += wdata[j] * plmtable[j];
124
125	for( ; j < (2 * bw); j += 4)
126	{
127	result0 += wdata[j] * plmtable[j];
128	result1 += wdata[j + 1] * plmtable[j + 1];
129	result2 += wdata[j + 2] * plmtable[j + 2];
130	result3 += wdata[j + 3] * plmtable[j + 3];
131	}
132	result[i] = result0 + result1 + result2 + result3;
133
134	plmtable += (2 * bw);
135	}
136	}
137	else
138	{
139	for (i = 0; i < bw - m; i++)
140	{
141	result0 = 0.0 ;
142	for(j = 0; j < (2 * bw) ; j++)
143	result0 += wdata[j] * plmtable[j];
144	result[i] = result0 ;
145
146	plmtable += (2 * bw);
147	}
148	}
149	}
150
151
152	/************************************************************************/
153	/* This is the procedure that synthesizes a function from a list
154	of coefficients of a Legendre series. I.e. this is the INVERSE
155	discrete Legendre transform.
156
157	Function is synthesized at the (2*bw) Chebyshev nodes. Associated
158	Legendre functions are assumed to be precomputed.
159
160	bw - bandwidth
161
162	m - order
163
164	coeffs - a pointer to double array of size (bw-m). First coefficient is
165	coefficient for Pmm
166	result - a pointer to double array of size (2*bw); at the conclusion
167	of the routine, this array will contain the
168	synthesized function
169
170	plmtable - a pointer to a double array of size (2bw(bw-m));
171	contains the PRECOMPUTED plms, i.e. associated Legendre
172	functions. E.g. Should be generated by a call to
173
174	PmlTableGen(),
175
176	(see pmls.c)
177
178
179	NOTE that these Legendres are normalized with norm
180	equal to 1 !!!
181	*/
182
183	void Naive_SynthesizeX(double *coeffs,
184	int bw,
185	int m,
186	double *result,
187	double *plmtable)
188	{
189	int i, j;
190	double tmpcoef;
191
192	/* make sure result is zeroed out */
193	memset( result, 0, sizeof(double) * 2 * bw );
194
195	for ( i = 0 ; i < bw - m ; i ++ )
196	{
197	tmpcoef = coeffs[i] ;
198	if ( tmpcoef != 0.0 )
199	for (j=0; j<(2*bw); j++)
200	result[j] += (tmpcoef * plmtable[j]);
201	plmtable += (2 * bw ) ;
202	}
203	}