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...
#	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
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			}