trunk/SHAPES/seminaive.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 for computing the Legendre transform where
   projections are carried out in cosine space, i.e., the
   "seminaive" algorithm.

   For a description, see the related paper or Sean's thesis.

*/

#include <math.h>
#include <stdio.h>
#include <string.h>   /** for memcpy **/
#include "fftw3.h"

#include "cospmls.h"


/************************************************************************/
/* InvSemiNaiveReduced computes the inverse Legendre transform
   using the transposed seminaive algorithm.  Note that because
   the Legendre transform is orthogonal, the inverse can be
   computed by transposing the matrix formulation of the
   problem.

   The forward transform looks like

   l = PCWf

   where f is the data vector, W is a quadrature matrix,
   C is a cosine transform matrix, P is a matrix
   full of coefficients of the cosine series representation
   of each Pml function P(m,m) P(m,m+1) ... P(m,bw-1),
   and l is the (associated) Legendre series representation
   of f.

   So to do the inverse, you do

   f = trans(C) trans(P) l

   so you need to transpose the matrix P from the forward transform
   and then do a cosine series evaluation.  No quadrature matrix
   is necessary.  If order m is odd, then there is also a sin
   factor that needs to be accounted for.

   Note that this function was written to be part of a full
   spherical harmonic transform, so a lot of precomputation
   has been assumed.

   Input argument description

   coeffs - a double pointer to an array of length
            (bw-m) containing associated
            Legendre series coefficients.  Assumed
            that first entry contains the P(m,m)
            coefficient.

   bw - problem bandwidth

   m - order of the associated Legendre functions

   result - a double pointer to an array of (2*bw) samples
            representing the evaluation of the Legendre
            series at (2*bw) Chebyshev nodes.

   trans_cos_pml_table - double pointer to array representing
                         the linearized form of trans(P) above.
                         See cospmls.{h,c} for a description
                         of the function Transpose_CosPmlTableGen()
                         which generates this array.

   sin_values - when m is odd, need to factor in the sin(x) that
                is factored out of the generation of the values
                in trans(P).

   workspace - a double array of size 2*bw -> temp space involving
          intermediate array

   fplan - pointer to fftw plan with input array being fcos
           and output being result; I'll probably be using the
           guru interface to execute - that way I can apply the
           same plan to different arrays; the plan should be
   
           fftw_plan_r2r_1d( 2*bw, fcos, result,
                             FFTW_REDFT01, FFTW_ESTIMATE );

*/

void InvSemiNaiveReduced(double *coeffs,
                         int bw, 
                         int m, 
                         double *result, 
                         double *trans_cos_pml_table, 
                         double *sin_values,
                         double *workspace,
                         fftw_plan *fplan )
{
  double *trans_tableptr;
  double *assoc_offset;
  int i, j, rowsize;
  double *p;
  double *fcos, fcos0, fcos1, fcos2, fcos3;
  double fudge ;

  fcos = workspace ;

  /* for paranoia, zero out arrays */
  memset( fcos, 0, sizeof(double) * 2 * bw );
  memset( result, 0, sizeof(double) * 2 * bw );

  trans_tableptr = trans_cos_pml_table;
  p = trans_cos_pml_table;

  /* main loop - compute each value of fcos

  Note that all zeroes have been stripped out of the
  trans_cos_pml_table, so indexing is somewhat complicated.
  */

  for (i=0; i<bw; i++)
    {
      if (i == (bw-1))
        {
          if ( m % 2 )
            {
              fcos[bw-1] = 0.0;
              break;
            }
        }

      rowsize = Transpose_RowSize(i, m, bw);
      if (i > m)
        assoc_offset = coeffs + (i - m) + (m % 2);
      else
        assoc_offset = coeffs + (i % 2);

      fcos0 = 0.0 ; fcos1 = 0.0; fcos2 = 0.0; fcos3 = 0.0;
          
      for (j = 0; j < rowsize % 4; ++j)
        fcos0 += assoc_offset[2*j] * trans_tableptr[j];
          
      for ( ; j < rowsize; j += 4){
        fcos0 += assoc_offset[2*j] * trans_tableptr[j];
        fcos1 += assoc_offset[2*(j+1)] * trans_tableptr[j+1];
        fcos2 += assoc_offset[2*(j+2)] * trans_tableptr[j+2];
        fcos3 += assoc_offset[2*(j+3)] * trans_tableptr[j+3];
      }
      fcos[i] = fcos0 + fcos1 + fcos2 + fcos3 ;

      trans_tableptr += rowsize;
    }
    

  /*
    now we have the cosine series for the result,
    so now evaluate the cosine series at 2*bw Chebyshev nodes 
  */

  /* scale coefficients prior to taking inverse DCT */
  fudge = 0.5 / sqrt((double) bw) ;
  for ( j = 1 ; j < 2*bw ; j ++ )
    fcos[j] *= fudge ;
  fcos[0] /= sqrt(2. * ((double) bw));

  /* now take the inverse dct */
  /* NOTE that I am using the guru interface */
  fftw_execute_r2r( *fplan,
                    fcos, result );

  /* if m is odd, then need to multiply by sin(x) at Chebyshev nodes */
  if ( m % 2 )
    {
      for (j=0; j<(2*bw); j++)
        result[j] *= sin_values[j];
    }

  trans_tableptr = p;

  /* amscray */

}

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

/* SemiNaiveReduced computes the Legendre transform of data.
   This function has been designed to be a component in
   a full spherical harmonic transform.  

   data - pointer to double array of size (2*bw) containing
          function to be transformed.  Assumes sampling at Chebyshev nodes

   bw   - bandwidth of the problem
   m   - order of the problem.  0 <= m < bw

   result - pointer to double array of length bw for returning computed
            Legendre coefficients.  Contains 
            bw-m coeffs, with the <f,P(m,m)> coefficient
            located in result[0]

   cos_pml_table - a pointer to an array containing the cosine
                   series coefficients of the Pmls (or Gmls)
                   for this problem.  This table can be computed
                   using the CosPmlTableGen() function, and
                   the offset for a particular Pml can be had
                   by calling the function NewTableOffset().
                   The size of the table is computed using
                   the TableSize() function.  Note that
                   since the cosine series are always zero-striped,
                   the zeroes have been removed.

   weights -> ptr to double array of size 4*bw - this array holds
           the weights for both even (starting at weights[0])
           and odd (weights[2*bw]) transforms


   workspace -> tmp space: ptr to double array of size 4*bw

   fplan -> pointer to fftw plan with input array being weighted_data
           and output being cos_data; I'll probably be using the
           guru interface to execute; the plan should be

           fftw_plan_r2r_1d( 2*bw, weighted_data, cos_data,
                        FFTW_REDFT10, FFTW_ESTIMATE ) ;


*/

void SemiNaiveReduced(double *data, 
                      int bw, 
                      int m, 
                      double *result,
                      double *workspace,
                      double *cos_pml_table, 
                      double *weights,
                      fftw_plan *fplan )
{
  int i, j, n;
  double result0, result1, result2, result3;
  double fudge ;
  double d_bw;
  int toggle ;
  double *pml_ptr, *weighted_data, *cos_data ;

  n = 2*bw;
  d_bw = (double) bw;

  weighted_data = workspace ;
  cos_data = weighted_data + (2*bw) ;

  /* for paranoia, zero out the result array */
  memset( result, 0, sizeof(double)*(bw-m));
 
  /*
    need to apply quadrature weights to the data and compute
    the cosine transform
  */
  if ( m % 2 )
    for ( i = 0; i < n    ; ++i )
      weighted_data[i] = data[ i ] * weights[ 2*bw + i ];
  else
    for ( i = 0; i < n    ; ++i )
      weighted_data[i] = data[ i ] * weights[ i ];

  /*
    smooth the weighted signal
  */

  fftw_execute_r2r( *fplan,
                    weighted_data,
                    cos_data );

  /* need to normalize */
  cos_data[0] *= 0.707106781186547 ;
  fudge = 1./sqrt(2. * ((double) n ) );
  for ( j = 0 ; j < n ; j ++ )
    cos_data[j] *= fudge ;

  /*
    do the projections; Note that the cos_pml_table has
    had all the zeroes stripped out so the indexing is
    complicated somewhat
  */
  

  /******** this is the original loop

  toggle = 0 ;
  for (i=m; i<bw; i++)
  {
  pml_ptr = cos_pml_table + NewTableOffset(m,i);

  if ((m % 2) == 0)
  {
  for (j=0; j<(i/2)+1; j++)
  result[i-m] += cos_data[(2*j)+toggle] * pml_ptr[j];
  }
  else
  {
  if (((i-m) % 2) == 0)
  {
  for (j=0; j<(i/2)+1; j++)
  result[i-m] += cos_data[(2*j)+toggle] * pml_ptr[j];
  }
  else
  {
  for (j=0; j<(i/2); j++)
  result[i-m] += cos_data[(2*j)+toggle] * pml_ptr[j];
  }
  } 
      
  toggle = (toggle+1) % 2;
  }

  *****/
 
  /******** this is the new loop *********/
  toggle = 0 ;
  for ( i=m; i<bw; i++ )
    {
      pml_ptr = cos_pml_table + NewTableOffset(m,i);

      result0 = 0.0 ; result1 = 0.0 ;
      result2 = 0.0 ; result3 = 0.0 ; 

      for ( j = 0 ; j < ( (i/2) % 4 ) ; ++j )
        result0 += cos_data[(2*j)+toggle] * pml_ptr[j];

      for ( ; j < (i/2) ; j += 4 )
        {
          result0 += cos_data[(2*j)+toggle] * pml_ptr[j];
          result1 += cos_data[(2*(j+1))+toggle] * pml_ptr[j+1];
          result2 += cos_data[(2*(j+2))+toggle] * pml_ptr[j+2];
          result3 += cos_data[(2*(j+3))+toggle] * pml_ptr[j+3];
        }

      if ((((i-m) % 2) == 0 ) || ( (m % 2) == 0 ))
        result0 += cos_data[(2*(i/2))+toggle] * pml_ptr[(i/2)];

      result[i-m] = result0 + result1 + result2 + result3 ;
          
      toggle = (toggle + 1)%2 ;
          
    }
}


Revision:	1287
Committed:	Wed Jun 23 20:18:48 2004 UTC (20 years ago) by chrisfen
Content type:	text/plain
File size:	10526 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			/* Source code for computing the Legendre transform where
39			projections are carried out in cosine space, i.e., the
40			"seminaive" algorithm.
41
42			For a description, see the related paper or Sean's thesis.
43
44			*/
45
46			#include <math.h>
47			#include <stdio.h>
48			#include <string.h> / for memcpy /
49			#include "fftw3.h"
50
51			#include "cospmls.h"
52
53
54			/************************************************************************/
55			/* InvSemiNaiveReduced computes the inverse Legendre transform
56			using the transposed seminaive algorithm. Note that because
57			the Legendre transform is orthogonal, the inverse can be
58			computed by transposing the matrix formulation of the
59			problem.
60
61			The forward transform looks like
62
63			l = PCWf
64
65			where f is the data vector, W is a quadrature matrix,
66			C is a cosine transform matrix, P is a matrix
67			full of coefficients of the cosine series representation
68			of each Pml function P(m,m) P(m,m+1) ... P(m,bw-1),
69			and l is the (associated) Legendre series representation
70			of f.
71
72			So to do the inverse, you do
73
74			f = trans(C) trans(P) l
75
76			so you need to transpose the matrix P from the forward transform
77			and then do a cosine series evaluation. No quadrature matrix
78			is necessary. If order m is odd, then there is also a sin
79			factor that needs to be accounted for.
80
81			Note that this function was written to be part of a full
82			spherical harmonic transform, so a lot of precomputation
83			has been assumed.
84
85			Input argument description
86
87			coeffs - a double pointer to an array of length
88			(bw-m) containing associated
89			Legendre series coefficients. Assumed
90			that first entry contains the P(m,m)
91			coefficient.
92
93			bw - problem bandwidth
94
95			m - order of the associated Legendre functions
96
97			result - a double pointer to an array of (2*bw) samples
98			representing the evaluation of the Legendre
99			series at (2*bw) Chebyshev nodes.
100
101			trans_cos_pml_table - double pointer to array representing
102			the linearized form of trans(P) above.
103			See cospmls.{h,c} for a description
104			of the function Transpose_CosPmlTableGen()
105			which generates this array.
106
107			sin_values - when m is odd, need to factor in the sin(x) that
108			is factored out of the generation of the values
109			in trans(P).
110
111			workspace - a double array of size 2*bw -> temp space involving
112			intermediate array
113
114			fplan - pointer to fftw plan with input array being fcos
115			and output being result; I'll probably be using the
116			guru interface to execute - that way I can apply the
117			same plan to different arrays; the plan should be
118
119			fftw_plan_r2r_1d( 2*bw, fcos, result,
120			FFTW_REDFT01, FFTW_ESTIMATE );
121
122			*/
123
124			void InvSemiNaiveReduced(double *coeffs,
125			int bw,
126			int m,
127			double *result,
128			double *trans_cos_pml_table,
129			double *sin_values,
130			double *workspace,
131			fftw_plan *fplan )
132			{
133			double *trans_tableptr;
134			double *assoc_offset;
135			int i, j, rowsize;
136			double *p;
137			double *fcos, fcos0, fcos1, fcos2, fcos3;
138			double fudge ;
139
140			fcos = workspace ;
141
142			/* for paranoia, zero out arrays */
143			memset( fcos, 0, sizeof(double) * 2 * bw );
144			memset( result, 0, sizeof(double) * 2 * bw );
145
146			trans_tableptr = trans_cos_pml_table;
147			p = trans_cos_pml_table;
148
149			/* main loop - compute each value of fcos
150
151			Note that all zeroes have been stripped out of the
152			trans_cos_pml_table, so indexing is somewhat complicated.
153			*/
154
155			for (i=0; i<bw; i++)
156			{
157			if (i == (bw-1))
158			{
159			if ( m % 2 )
160			{
161			fcos[bw-1] = 0.0;
162			break;
163			}
164			}
165
166			rowsize = Transpose_RowSize(i, m, bw);
167			if (i > m)
168			assoc_offset = coeffs + (i - m) + (m % 2);
169			else
170			assoc_offset = coeffs + (i % 2);
171
172			fcos0 = 0.0 ; fcos1 = 0.0; fcos2 = 0.0; fcos3 = 0.0;
173
174			for (j = 0; j < rowsize % 4; ++j)
175			fcos0 += assoc_offset[2j] trans_tableptr[j];
176
177			for ( ; j < rowsize; j += 4){
178			fcos0 += assoc_offset[2j] trans_tableptr[j];
179			fcos1 += assoc_offset[2(j+1)] trans_tableptr[j+1];
180			fcos2 += assoc_offset[2(j+2)] trans_tableptr[j+2];
181			fcos3 += assoc_offset[2(j+3)] trans_tableptr[j+3];
182			}
183			fcos[i] = fcos0 + fcos1 + fcos2 + fcos3 ;
184
185			trans_tableptr += rowsize;
186			}
187
188
189			/*
190			now we have the cosine series for the result,
191			so now evaluate the cosine series at 2*bw Chebyshev nodes
192			*/
193
194			/* scale coefficients prior to taking inverse DCT */
195			fudge = 0.5 / sqrt((double) bw) ;
196			for ( j = 1 ; j < 2*bw ; j ++ )
197			fcos[j] *= fudge ;
198			fcos[0] /= sqrt(2. * ((double) bw));
199
200			/* now take the inverse dct */
201			/* NOTE that I am using the guru interface */
202			fftw_execute_r2r( *fplan,
203			fcos, result );
204
205			/* if m is odd, then need to multiply by sin(x) at Chebyshev nodes */
206			if ( m % 2 )
207			{
208			for (j=0; j<(2*bw); j++)
209			result[j] *= sin_values[j];
210			}
211
212			trans_tableptr = p;
213
214			/* amscray */
215
216			}
217
218			/************************************************************************/
219
220			/* SemiNaiveReduced computes the Legendre transform of data.
221			This function has been designed to be a component in
222			a full spherical harmonic transform.
223
224			data - pointer to double array of size (2*bw) containing
225			function to be transformed. Assumes sampling at Chebyshev nodes
226
227			bw - bandwidth of the problem
228			m - order of the problem. 0 <= m < bw
229
230			result - pointer to double array of length bw for returning computed
231			Legendre coefficients. Contains
232			bw-m coeffs, with the <f,P(m,m)> coefficient
233			located in result[0]
234
235			cos_pml_table - a pointer to an array containing the cosine
236			series coefficients of the Pmls (or Gmls)
237			for this problem. This table can be computed
238			using the CosPmlTableGen() function, and
239			the offset for a particular Pml can be had
240			by calling the function NewTableOffset().
241			The size of the table is computed using
242			the TableSize() function. Note that
243			since the cosine series are always zero-striped,
244			the zeroes have been removed.
245
246			weights -> ptr to double array of size 4*bw - this array holds
247			the weights for both even (starting at weights[0])
248			and odd (weights[2*bw]) transforms
249
250
251			workspace -> tmp space: ptr to double array of size 4*bw
252
253			fplan -> pointer to fftw plan with input array being weighted_data
254			and output being cos_data; I'll probably be using the
255			guru interface to execute; the plan should be
256
257			fftw_plan_r2r_1d( 2*bw, weighted_data, cos_data,
258			FFTW_REDFT10, FFTW_ESTIMATE ) ;
259
260
261			*/
262
263			void SemiNaiveReduced(double *data,
264			int bw,
265			int m,
266			double *result,
267			double *workspace,
268			double *cos_pml_table,
269			double *weights,
270			fftw_plan *fplan )
271			{
272			int i, j, n;
273			double result0, result1, result2, result3;
274			double fudge ;
275			double d_bw;
276			int toggle ;
277			double pml_ptr, weighted_data, *cos_data ;
278
279			n = 2*bw;
280			d_bw = (double) bw;
281
282			weighted_data = workspace ;
283			cos_data = weighted_data + (2*bw) ;
284
285			/* for paranoia, zero out the result array */
286			memset( result, 0, sizeof(double)*(bw-m));
287
288			/*
289			need to apply quadrature weights to the data and compute
290			the cosine transform
291			*/
292			if ( m % 2 )
293			for ( i = 0; i < n ; ++i )
294			weighted_data[i] = data[ i ] * weights[ 2*bw + i ];
295			else
296			for ( i = 0; i < n ; ++i )
297			weighted_data[i] = data[ i ] * weights[ i ];
298
299			/*
300			smooth the weighted signal
301			*/
302
303			fftw_execute_r2r( *fplan,
304			weighted_data,
305			cos_data );
306
307			/* need to normalize */
308			cos_data[0] *= 0.707106781186547 ;
309			fudge = 1./sqrt(2. * ((double) n ) );
310			for ( j = 0 ; j < n ; j ++ )
311			cos_data[j] *= fudge ;
312
313			/*
314			do the projections; Note that the cos_pml_table has
315			had all the zeroes stripped out so the indexing is
316			complicated somewhat
317			*/
318
319
320			/******** this is the original loop
321
322			toggle = 0 ;
323			for (i=m; i<bw; i++)
324			{
325			pml_ptr = cos_pml_table + NewTableOffset(m,i);
326
327			if ((m % 2) == 0)
328			{
329			for (j=0; j<(i/2)+1; j++)
330			result[i-m] += cos_data[(2j)+toggle] pml_ptr[j];
331			}
332			else
333			{
334			if (((i-m) % 2) == 0)
335			{
336			for (j=0; j<(i/2)+1; j++)
337			result[i-m] += cos_data[(2j)+toggle] pml_ptr[j];
338			}
339			else
340			{
341			for (j=0; j<(i/2); j++)
342			result[i-m] += cos_data[(2j)+toggle] pml_ptr[j];
343			}
344			}
345
346			toggle = (toggle+1) % 2;
347			}
348
349			*****/
350
351			/****** this is the new loop *******/
352			toggle = 0 ;
353			for ( i=m; i<bw; i++ )
354			{
355			pml_ptr = cos_pml_table + NewTableOffset(m,i);
356
357			result0 = 0.0 ; result1 = 0.0 ;
358			result2 = 0.0 ; result3 = 0.0 ;
359
360			for ( j = 0 ; j < ( (i/2) % 4 ) ; ++j )
361			result0 += cos_data[(2j)+toggle] pml_ptr[j];
362
363			for ( ; j < (i/2) ; j += 4 )
364			{
365			result0 += cos_data[(2j)+toggle] pml_ptr[j];
366			result1 += cos_data[(2(j+1))+toggle] pml_ptr[j+1];
367			result2 += cos_data[(2(j+2))+toggle] pml_ptr[j+2];
368			result3 += cos_data[(2(j+3))+toggle] pml_ptr[j+3];
369			}
370
371			if ((((i-m) % 2) == 0 ) \|\| ( (m % 2) == 0 ))
372			result0 += cos_data[(2(i/2))+toggle] pml_ptr[(i/2)];
373
374			result[i-m] = result0 + result1 + result2 + result3 ;
375
376			toggle = (toggle + 1)%2 ;
377
378			}
379			}
380
381