sinfo_resampling.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 /*----------------------------------------------------------------------------
00020    
00021    File name    :   resampling.c
00022    Author       :   Nicolas Devillard
00023    Created on   :   Jan 04, 1996
00024    Description  :   resampling routines
00025 
00026  ---------------------------------------------------------------------------*/
00027 
00028 /*
00029     $Id: sinfo_resampling.c,v 1.3 2006/12/01 12:45:36 amodigli Exp $
00030     $Author: amodigli $
00031     $Date: 2006/12/01 12:45:36 $
00032     $Revision: 1.3 $
00033  */
00034 
00035 #ifdef HAVE_CONFIG_H
00036 #  include <config.h>
00037 #endif
00038 /*---------------------------------------------------------------------------
00039                                 Includes
00040  ---------------------------------------------------------------------------*/
00041 #include "sinfo_resampling.h"
00042 #include <math.h>
00043 #include "sinfo_globals.h"
00044 /*---------------------------------------------------------------------------
00045                             Private functions
00046  ---------------------------------------------------------------------------*/
00047 static void reverse_tanh_kernel(double * data, int nn) ;
00048 static double * sinfo_generate_tanh_kernel(double steep);
00049 
00059 /*---------------------------------------------------------------------------
00060                             Function codes
00061  ---------------------------------------------------------------------------*/
00062 
00087 double   *
00088 sinfo_generate_interpolation_kernel(const char * kernel_type)
00089 {
00090     double  *   tab ;
00091     int         i ;
00092     double      x ;
00093     double      alpha ;
00094     double      inv_norm ;
00095     int         samples = KERNEL_SAMPLES ;
00096 
00097     if (kernel_type==NULL) {
00098         tab = sinfo_generate_interpolation_kernel("tanh") ;
00099     } else if (!strcmp(kernel_type, "default")) {
00100         tab = sinfo_generate_interpolation_kernel("tanh") ;
00101     } else if (!strcmp(kernel_type, "sinc")) {
00102         tab = cpl_malloc(samples * sizeof(double)) ;
00103         tab[0] = 1.0 ;
00104         tab[samples-1] = 0.0 ;
00105         for (i=1 ; i<samples ; i++) {
00106             x = (double)KERNEL_WIDTH * (double)i/(double)(samples-1) ;
00107             tab[i] = sinfo_sinc(x) ;
00108         }
00109     } else if (!strcmp(kernel_type, "sinc2")) {
00110         tab = cpl_malloc(samples * sizeof(double)) ;
00111         tab[0] = 1.0 ;
00112         tab[samples-1] = 0.0 ;
00113         for (i=1 ; i<samples ; i++) {
00114             x = 2.0 * (double)i/(double)(samples-1) ;
00115             tab[i] = sinfo_sinc(x) ;
00116             tab[i] *= tab[i] ;
00117         }
00118     } else if (!strcmp(kernel_type, "lanczos")) {
00119         tab = cpl_malloc(samples * sizeof(double)) ;
00120         for (i=0 ; i<samples ; i++) {
00121             x = (double)KERNEL_WIDTH * (double)i/(double)(samples-1) ;
00122             if (fabs(x)<2) {
00123                 tab[i] = sinfo_sinc(x) * sinfo_sinc(x/2) ;
00124             } else {
00125                 tab[i] = 0.00 ;
00126             }
00127         }
00128     } else if (!strcmp(kernel_type, "hamming")) {
00129         tab = cpl_malloc(samples * sizeof(double)) ;
00130         alpha = 0.54 ;
00131         inv_norm  = 1.00 / (double)(samples - 1) ;
00132         for (i=0 ; i<samples ; i++) {
00133             x = (double)i ;
00134             if (i<(samples-1)/2) {
00135                 tab[i] = alpha + (1-alpha) * cos(2.0*PI_NUMB*x*inv_norm) ;
00136             } else {
00137                 tab[i] = 0.0 ;
00138             }
00139         }
00140     } else if (!strcmp(kernel_type, "hann")) {
00141         tab = cpl_malloc(samples * sizeof(double)) ;
00142         alpha = 0.50 ;
00143         inv_norm  = 1.00 / (double)(samples - 1) ;
00144         for (i=0 ; i<samples ; i++) {
00145             x = (double)i ;
00146             if (i<(samples-1)/2) {
00147                 tab[i] = alpha + (1-alpha) * cos(2.0*PI_NUMB*x*inv_norm) ;
00148             } else {
00149                 tab[i] = 0.0 ;
00150             }
00151         }
00152     } else if (!strcmp(kernel_type, "tanh")) {
00153         tab = sinfo_generate_tanh_kernel(TANH_STEEPNESS) ;
00154     } else {
00155         sinfo_msg_error("unrecognized kernel type [%s]: aborting generation",
00156                 kernel_type) ;
00157         return NULL ;
00158     }
00159 
00160     return tab ;
00161 }
00162 
00174 double
00175 sinfo_sinc(double x)
00176 {
00177     if (fabs(x)<1e-4)
00178         return (double)1.00 ;
00179     else
00180         return ((sin(x * (double)PI_NUMB)) / (x * (double)PI_NUMB)) ;
00181 }
00182 
00204 #define hk_gen(x,s) (((tanh(s*(x+0.5))+1)/2)*((tanh(s*(-x+0.5))+1)/2))
00205 
00206 static double * 
00207 sinfo_generate_tanh_kernel(double steep)
00208 {
00209     double  *   kernel ;
00210     double  *   x ;
00211     double      width ;
00212     double      inv_np ;
00213     double      ind ;
00214     int         i ;
00215     int         np ;
00216     int         samples ;
00217 
00218     width   = (double)TABSPERPIX / 2.0 ; 
00219     samples = KERNEL_SAMPLES ;
00220     np      = 32768 ; /* Hardcoded: should never be changed */
00221     inv_np  = 1.00 / (double)np ;
00222 
00223     /*
00224      * Generate the kernel expression in Fourier space
00225      * with a correct frequency ordering to allow standard FT
00226      */
00227     x = cpl_malloc((2*np+1)*sizeof(double)) ;
00228     for (i=0 ; i<np/2 ; i++) {
00229         ind      = (double)i * 2.0 * width * inv_np ;
00230         x[2*i]   = hk_gen(ind, steep) ;
00231         x[2*i+1] = 0.00 ;
00232     }
00233     for (i=np/2 ; i<np ; i++) {
00234         ind      = (double)(i-np) * 2.0 * width * inv_np ;
00235         x[2*i]   = hk_gen(ind, steep) ;
00236         x[2*i+1] = 0.00 ;
00237     }
00238 
00239     /* 
00240      * Reverse Fourier to come back to image space
00241      */
00242     reverse_tanh_kernel(x, np) ;
00243 
00244     /*
00245      * Allocate and fill in returned array
00246      */
00247     kernel = cpl_malloc(samples * sizeof(double)) ;
00248     for (i=0 ; i<samples ; i++) {
00249         kernel[i] = 2.0 * width * x[2*i] * inv_np ;
00250     }
00251     cpl_free(x) ;
00252     return kernel ;
00253 }
00254 
00267 #define KERNEL_SW(a,b) tempr=(a);(a)=(b);(b)=tempr
00268 static void 
00269 reverse_tanh_kernel(double * data, int nn)
00270 {
00271     unsigned long   n,
00272                     mmax,
00273                     m,
00274                     i, j,
00275                     istep ;
00276     double  wtemp,
00277             wr,
00278             wpr,
00279             wpi,
00280             wi,
00281             theta;
00282     double  tempr,
00283             tempi;
00284 
00285     n = (unsigned long)nn << 1;
00286     j = 1;
00287     for (i=1 ; i<n ; i+=2) {
00288         if (j > i) {
00289             KERNEL_SW(data[j-1],data[i-1]);
00290             KERNEL_SW(data[j],data[i]);
00291         }
00292         m = n >> 1;
00293         while (m>=2 && j>m) {
00294             j -= m;
00295             m >>= 1;
00296         }
00297         j += m;
00298     }
00299     mmax = 2;
00300     while (n > mmax) {
00301         istep = mmax << 1;
00302         theta = 2 * PI_NUMB / mmax;
00303         wtemp = sin(0.5 * theta);
00304         wpr = -2.0 * wtemp * wtemp;
00305         wpi = sin(theta);
00306         wr  = 1.0;
00307         wi  = 0.0;
00308         for (m=1 ; m<mmax ; m+=2) {
00309             for (i=m ; i<=n ; i+=istep) {
00310                 j = i + mmax;
00311                 tempr = wr * data[j-1] - wi * data[j];
00312                 tempi = wr * data[j]   + wi * data[j-1];
00313                 data[j-1] = data[i-1] - tempr;
00314                 data[j]   = data[i]   - tempi;
00315                 data[i-1] += tempr;
00316                 data[i]   += tempi;
00317             }
00318             wr = (wtemp = wr) * wpr - wi * wpi + wr;
00319             wi = wi * wpr + wtemp * wpi + wi;
00320         }
00321         mmax = istep;
00322     }
00323 }
00324 #undef KERNEL_SW
00325 
00363 double *
00364 sinfo_invert_linear_transform(double *trans)
00365 {
00366     double   *  i_trans ;
00367     double      det ;
00368 
00369     if (trans==NULL) return NULL ;
00370     det = (trans[0] * trans[4]) - (trans[1] * trans[3]) ;
00371     if (fabs(det) < 1e-6) {
00372         sinfo_msg_error("NULL determinant: cannot sinfo_invert transform") ;
00373         return NULL ;
00374     }
00375 
00376     i_trans = cpl_calloc(6, sizeof(double)) ;
00377 
00378     i_trans[0] =  trans[4] / det ; 
00379     i_trans[1] = -trans[1] / det ;
00380     i_trans[2] = ((trans[1] * trans[5]) - (trans[2] * trans[4])) / det ;
00381     i_trans[3] = -trans[3] / det ;
00382     i_trans[4] =  trans[0] / det ;
00383     i_trans[5] = ((trans[2] * trans[3]) - (trans[0] * trans[5])) / det ;
00384 
00385     return i_trans ;
00386 }

Generated on Wed Jan 17 08:33:44 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4