sinfo_function_1d.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    :   function_1d.c
00022    Author       :   Nicolas Devillard
00023    Created on   :   Tue, Sept 23 1997   
00024    Description  :   1d signal processing related routines   
00025 
00026  ---------------------------------------------------------------------------*/
00027 /*
00028  $Id: sinfo_function_1d.c,v 1.3 2006/11/06 08:11:36 amodigli Exp $
00029  $Author: amodigli $
00030  $Date: 2006/11/06 08:11:36 $
00031  $Revision: 1.3 $
00032  */
00033 
00034 #ifdef HAVE_CONFIG_H
00035 #  include <config.h>
00036 #endif
00037 /*----------------------------------------------------------------------------
00038                                 Includes
00039  ---------------------------------------------------------------------------*/
00040 #include <math.h>
00041 #include "sinfo_function_1d.h"
00042 #include "sinfo_fit_curve.h"
00043 #include "sinfo_median.h"
00044 
00053 /*----------------------------------------------------------------------------
00054                                 Defines
00055  ---------------------------------------------------------------------------*/
00056 /*
00057  * This parameter sets up the half size of the domain around which a
00058  * centroid position will be computed.
00059  */
00060 #define HALF_CENTROID_DOMAIN    5
00061 /*----------------------------------------------------------------------------
00062                         Private function prototypes
00063  ---------------------------------------------------------------------------*/
00064 static double * function1d_generate_smooth_kernel(int filt_type, int hw);
00065 
00066 static int
00067 function1d_search_value(
00068     pixelvalue  *   x,
00069     int             len,
00070     pixelvalue      key,
00071     int         *   foundPtr
00072 ) ;
00073 
00074 /*----------------------------------------------------------------------------
00075                             Function codes
00076  ---------------------------------------------------------------------------*/
00093 pixelvalue *
00094 sinfo_function1d_new(int nsamples)
00095 {
00096     if (nsamples<1) return NULL ;
00097     return cpl_calloc(nsamples, sizeof(pixelvalue)) ;
00098 }
00099 
00100 
00110 void 
00111 sinfo_function1d_del(pixelvalue * s)
00112 {
00113     if (s)
00114         cpl_free(s);
00115     return ;
00116 }
00117 
00134 pixelvalue * 
00135 sinfo_function1d_dup(pixelvalue * arr, int ns)
00136 {
00137     pixelvalue  *   n_arr ;
00138 
00139     n_arr = sinfo_function1d_new(ns);
00140     memcpy(n_arr, arr, ns * sizeof(pixelvalue));
00141     return n_arr ;
00142 }
00143 
00159 double
00160 sinfo_function1d_find_centroid(
00161     pixelvalue  *   line,   /*  the input line  */
00162     int npix    /*  number of pixels in this line   */
00163 )
00164 {
00165     pixelvalue  max ;
00166     double      centroid ;
00167     double      weights ;
00168     int         i, maxpos ;
00169 
00170     /*
00171      * Search for the maximum pixel value on the line
00172      */
00173 
00174     max = line[0] ;
00175     maxpos = 0 ;
00176     for (i=1 ; i<npix ; i++) {
00177         if (line[i]>max) {
00178             max = line[i] ;
00179             maxpos = i ;
00180         }
00181     }
00182 
00183     /*
00184      * The centroid position is the weighted average over the maximum 
00185      * pixel neighborhood.
00186      */
00187 
00188     centroid = 0.0 ;
00189     weights  = 0.0 ;
00190     for (i=maxpos-HALF_CENTROID_DOMAIN; 
00191             i<=maxpos+HALF_CENTROID_DOMAIN; i++) {
00192         centroid += (double)line[i] * (double)i ;
00193         weights  += (double)line[i] ;
00194     }
00195 
00196     centroid /= weights ;
00197     return centroid ;   
00198 }
00199 
00222 double
00223 sinfo_function1d_find_locmax(
00224     pixelvalue  *   line,
00225     int             npix,
00226     int             where,
00227     int             hs
00228 )
00229 {
00230     pixelvalue  max ;
00231     double      centroid ;
00232     double      weights ;
00233     int         i, maxpos ;
00234 
00235 
00236     if ((where<hs) || (where>(npix-hs-1))) {
00237         return (double)-1.0 ;
00238     }
00239     
00240     /*
00241      * Search for the closest local maximal around the requested range.
00242      */
00243     max = line[where] ;
00244     maxpos = where ;
00245     for (i=-hs ; i<=hs ; i++) {
00246         if (line[where+i]>max) {
00247             max = line[where+i] ;
00248             maxpos = where+i ;
00249         }
00250     }
00251 
00252     /*
00253      * The centroid position is the weighted average over the maximum 
00254      * pixel neighborhood.
00255      */
00256 
00257     centroid = 0.0 ;
00258     weights  = 0.0 ;
00259     for (i=maxpos-hs; i<=maxpos+hs; i++) {
00260         centroid += (double)line[i] * (double)i ;
00261         weights  += (double)line[i] ;
00262     }
00263     if (fabs(weights)>1e-6) {
00264         centroid /= weights ;
00265     } else {
00266         centroid = -1.0 ;
00267     }
00268     return centroid ;   
00269 }
00270 
00295 pixelvalue *
00296 sinfo_function1d_filter_lowpass(
00297     pixelvalue  *   input_sig,
00298     int             samples,
00299     int             filter_type,
00300     int             hw
00301 )
00302 {
00303     pixelvalue  *   out_sig ;
00304     int             i, j ;
00305     double          replace ;
00306     double      *   kernel ;
00307 
00308     /* allocate output signal */
00309     out_sig = sinfo_function1d_new(samples);
00310 
00311     /* generate low-pass filter kernel */
00312     kernel = function1d_generate_smooth_kernel(filter_type, hw) ;
00313 
00314     /* compute sinfo_edge effects for the first hw elements */
00315     for (i=0 ; i<hw ; i++) {
00316         replace = 0.0 ;
00317         for (j=-hw ; j<=hw ; j++) {
00318             if (i+j<0) {
00319                 replace += kernel[hw+j] * (double)input_sig[0] ;
00320             } else {
00321                 replace += kernel[hw+j] * (double)input_sig[i+j] ;
00322             }
00323         }
00324         out_sig[i] = (pixelvalue)replace ;
00325     }
00326 
00327     /* compute sinfo_edge effects for the last hw elements */
00328     for (i=samples-hw ; i<samples ; i++) {
00329         replace = 0.0 ;
00330         for (j=-hw ; j<=hw ; j++) {
00331             if (i+j>samples-1) {
00332                 replace += kernel[hw+j] * (double)input_sig[samples-1] ;
00333             } else {
00334                 replace += kernel[hw+j] * (double)input_sig[i+j] ;
00335             }
00336         }
00337         out_sig[i] = (pixelvalue)replace ;
00338     }
00339 
00340     /* compute all other elements */
00341     for (i=hw ; i<samples-hw ; i++) {
00342         replace = 0.0 ;
00343         for (j=-hw ; j<=hw ; j++) {
00344             replace += kernel[hw+j] * (double)input_sig[i+j] ;
00345         }
00346         out_sig[i] = (pixelvalue)replace ;
00347     }
00348 
00349     cpl_free(kernel) ;
00350     return out_sig ;
00351 
00352 }
00353 
00370 static double * 
00371 function1d_generate_smooth_kernel(int filt_type, int hw)
00372 {
00373     double  *   kernel ;
00374     double      norm ;
00375     int         i ;
00376 
00377     kernel = (double*)cpl_calloc(2*hw+1, sizeof(double)) ;
00378 
00379     switch(filt_type) {
00380 
00381         case LOW_PASS_LINEAR:
00382         for (i=-hw ; i<=hw ; i++) {
00383             /* flat kernel */
00384             kernel[hw+i] = 1.0 / (double)(2*hw+1) ;
00385         }
00386         break ;
00387 
00388         case LOW_PASS_GAUSSIAN:
00389         norm = 0.00 ;
00390         for (i=-hw ; i<=hw ; i++) {
00391             /* sinfo_gaussian kernel */
00392             kernel[hw+i] = exp(-(double)(i*i)) ;
00393             norm += kernel[hw+i] ;
00394         }
00395         for (i=0 ; i<2*hw+1 ; i++) {
00396             kernel[i] /= norm ;
00397         }
00398         break ;
00399 
00400         default:
00401         sinfo_msg_error("unrecognized low pass filter: "
00402                                 "cannot generate kernel") ;
00403         return (double*)NULL ;
00404         break ;
00405     }
00406 
00407     return kernel ;
00408 }
00409 
00428 pixelvalue * 
00429 sinfo_function1d_median_smooth(
00430     pixelvalue * list,
00431     int          np,
00432     int          hw)
00433 {
00434     int             i,j ;
00435     pixelvalue  *   row ;
00436     pixelvalue  *   smoothed ;
00437 
00438     /* simply copy first 3 and last 3 items */
00439     smoothed = sinfo_function1d_new(np);
00440     for (i=0 ; i<hw ; i++) {
00441         smoothed[i] = list[i] ;
00442     }
00443     for (i=np-hw ; i<np ; i++) {
00444         smoothed[i] = list[i] ;
00445     }
00446 
00447     /* sinfo_median filter on all central items */
00448     row = sinfo_function1d_new(2*hw+1);
00449     for (i=hw ; i<np-hw ; i++) {
00450         for (j=-hw ; j<=hw ; j++) {
00451             row[j+hw] = list[i+j] ;
00452         }
00453         smoothed[i] = sinfo_median_pixelvalue(row, 2*hw+1) ; 
00454     }
00455     sinfo_function1d_del(row) ;
00456     return smoothed ;
00457 }
00458 
00475 #define LOWFREQ_PASSES      5
00476 
00477 pixelvalue * 
00478 sinfo_function1d_remove_lowfreq(
00479     pixelvalue * signal,
00480     int          ns)
00481 {
00482     pixelvalue  *   sig_in ;
00483     pixelvalue  *   smooth ;
00484     int             i ;
00485     
00486     
00487     /* Apply severe low-pass filter several times */
00488     sig_in = sinfo_function1d_dup(signal, ns);
00489     for (i=0 ; i<LOWFREQ_PASSES ; i++) {
00490         smooth = sinfo_function1d_filter_lowpass( sig_in, ns, 
00491                          LOW_PASS_LINEAR, 5);
00492         cpl_free(sig_in);
00493         sig_in = smooth ;
00494     }
00495 
00496     /* Subtract smoothed signal from input signal */
00497     for (i=0 ; i<ns ; i++) {
00498         smooth[i] = signal[i] - smooth[i];
00499     }
00500     return smooth ;
00501 }
00502 
00503 #undef LOWFREQ_PASSES
00504 
00526 #define SAMPLE_BORDER   10
00527 
00528 pixelvalue * 
00529 sinfo_function1d_remove_thermalbg(
00530     pixelvalue * signal,
00531     int          ns)
00532 {
00533     pixelvalue  *   smooth ;
00534     int             i ;
00535     
00536     int             nmin ;
00537     pixelvalue      lef[2], rig[2];
00538 
00539     pixelvalue  *   x,
00540                 *   y,
00541                 *   spl_x,
00542                 *   spl_y ;
00543     double          med_y ;
00544     double          avg2med ;
00545     double          dist ;
00546 
00547     
00548     /* Detect all local minima */
00549     nmin = 0 ;
00550     x = sinfo_function1d_new(ns);
00551     y = sinfo_function1d_new(ns);
00552 
00553     for (i=SAMPLE_BORDER ; i<(ns-SAMPLE_BORDER) ; i++) {
00554         lef[0] = signal[i-2];
00555         lef[1] = signal[i-1];
00556         rig[0] = signal[i+1];
00557         rig[1] = signal[i+2];
00558 
00559         if ( (signal[i] < lef[0]) &&
00560              (signal[i] < lef[1]) &&
00561              (signal[i] < rig[0]) &&
00562              (signal[i] < rig[1])) {
00563             x[nmin] = (pixelvalue)i ;
00564             y[nmin] = signal[i];
00565             nmin ++ ;
00566         }
00567     }
00568 
00569 
00570     /* Interpolate linearly missing values */
00571     spl_x = sinfo_function1d_new(ns);
00572     spl_y = sinfo_function1d_new(ns);
00573     for (i=0 ; i<ns ; i++) {
00574         spl_x[i] = (pixelvalue)i ;
00575     }
00576     sinfo_function1d_interpolate_linear(x, y, nmin, spl_x, spl_y, ns);
00577 
00578     sinfo_function1d_del(x) ;
00579     sinfo_function1d_del(y) ;
00580     sinfo_function1d_del(spl_x);
00581 
00582     /* Compute sinfo_median and average distance to the sinfo_median */
00583     med_y = (double)sinfo_median_pixelvalue(signal, ns);
00584     avg2med = 0.0 ;
00585     for (i=0 ; i<ns ; i++) {
00586         avg2med += fabs((double)signal[i] - med_y) ;
00587     }
00588     avg2med /= (double)ns ;
00589 
00590     /* Reset all pixels out of sinfo_median + 2 * avg2med to zero. */
00591     for (i=0 ; i<ns ; i++) {
00592         dist = fabs((double)signal[i] - med_y);
00593         if (dist > (2.0*avg2med)) {
00594             spl_y[i] = (pixelvalue)0 ;
00595         }
00596     }
00597 
00598 
00599     smooth = sinfo_function1d_new(ns);
00600     for (i=0 ; i<ns ; i++) {
00601         if (spl_y[i]>1e-4) {
00602             smooth[i] = signal[i] - spl_y[i];
00603         } else {
00604             smooth[i] = 0.0 ;
00605         }
00606     }
00607     sinfo_function1d_del(spl_y);
00608     return smooth ;
00609 }
00610 
00611 #undef LOWFREQ_PASSES
00612 
00634 void 
00635 sinfo_function1d_interpolate_linear(
00636     pixelvalue  *   x,
00637     pixelvalue  *   y,
00638     int             len,
00639     pixelvalue  *   spl_x,
00640     pixelvalue  *   spl_y,
00641     int             spl_len
00642 )
00643 {
00644     double      a, b ;
00645     int         i, j ;
00646     int         found ;
00647 
00648     for (i=0 ; i<spl_len ; i++) {
00649         /* Find (x1,y1) on the left of the current point */
00650         found = 0 ;
00651         for (j=0 ; j<(len-1) ; j++) {
00652             if ((spl_x[i]>=x[j]) && (spl_x[i]<=x[j+1])) {
00653                 found++ ;
00654                 break ;
00655             }
00656         }
00657         if (!found) {
00658             spl_y[i] = 0.0;
00659         } else {
00660             a = ((double)y[j+1]-(double)y[j]) /
00661                 ((double)x[j+1]-(double)x[j]);
00662             b = (double)y[j] - a * (double)x[j] ;
00663             spl_y[i] = (pixelvalue)(a * (double)spl_x[i] + b) ;
00664         }
00665     }
00666     return ;
00667 }
00668 
00685 static int
00686 function1d_search_value(
00687     pixelvalue  *   x,
00688     int             len,
00689     pixelvalue      key,
00690     int         *   foundPtr
00691 )
00692 {
00693     int high,
00694         low,
00695         middle;
00696 
00697     low  = 0;
00698     high = len - 1;
00699 
00700     while (high >= low) {
00701         middle = (high + low) / 2;
00702         if (key > x[middle]) {
00703             low = middle + 1;
00704         } else if (key < x[middle]) {
00705             high = middle - 1;
00706         } else {
00707             *foundPtr = 1;
00708             return (middle);
00709         }
00710     }
00711     *foundPtr = 0;
00712     return (low);
00713 }
00714 
00740 int
00741 sinfo_function1d_natural_spline(
00742     pixelvalue  *   x,
00743     pixelvalue  *   y,
00744     int             len,
00745     pixelvalue  *   splX,
00746     pixelvalue  *   splY,
00747     int             splLen
00748 )
00749 {
00750     int             end;
00751     int             loc,
00752                     found;
00753     register int    i,
00754                     j,
00755                     n;
00756     double      *   h; /* sinfo_vector of deltas in x */
00757     double      *   alpha;
00758     double      *   l,
00759                 *   mu,
00760                 *   z,
00761                 *   a,
00762                 *   b,
00763                 *   c,
00764                 *   d,
00765                     v;
00766 
00767     end = len - 1;
00768 
00769     a = cpl_malloc(sizeof(double) * splLen * 9) ;
00770     b = a + len;
00771     c = b + len;
00772     d = c + len;
00773     h = d + len;
00774     l = h + len;
00775     z = l + len;
00776     mu = z + len;
00777     alpha = mu + len;
00778 
00779     for (i = 0; i < len; i++) {
00780         a[i] = (double)y[i];
00781     }
00782 
00783     /* Calculate sinfo_vector of differences */
00784     for (i = 0; i < end; i++) {
00785         h[i] = (double)x[i + 1] - (double)x[i];
00786         if (h[i] < 0.0) {
00787             cpl_free(a) ;
00788             return -1;
00789         }
00790     }
00791 
00792     /* Calculate alpha sinfo_vector */
00793     for (n = 0, i = 1; i < end; i++, n++) {
00794         /* n = i - 1 */
00795         alpha[i] = 3.0 * ((a[i+1] / h[i]) - (a[i] / h[n]) - (a[i] / h[i]) +
00796                   (a[n] / h[n]));
00797     }
00798 
00799     /* Vectors to solve the tridiagonal sinfo_matrix */
00800     l[0] = l[end] = 1.0;
00801     mu[0] = mu[end] = 0.0;
00802     z[0] = z[end] = 0.0;
00803     c[0] = c[end] = 0.0;
00804 
00805     /* Calculate the intermediate results */
00806     for (n = 0, i = 1; i < end; i++, n++) {
00807         /* n = i-1 */
00808         l[i] = 2 * (h[i] + h[n]) - h[n] * mu[n];
00809         mu[i] = h[i] / l[i];
00810         z[i] = (alpha[i] - h[n] * z[n]) / l[i];
00811     }
00812     for (n = end, j = end - 1; j >= 0; j--, n--) {
00813         /* n = j + 1 */
00814         c[j] = z[j] - mu[j] * c[n];
00815         b[j] = (a[n] - a[j]) / h[j] - h[j] * (c[n] + 2.0 * c[j]) / 3.0;
00816         d[j] = (c[n] - c[j]) / (3.0 * h[j]);
00817     }
00818 
00819     /* Now calculate the new values */
00820     for (j = 0; j < splLen; j++) {
00821          v = (double)splX[j];
00822      splY[j] = (pixelvalue)0;
00823 
00824      /* Is it outside the interval? */
00825      if ((v < (double)x[0]) || (v > (double)x[end])) {
00826         continue;
00827      }
00828      /* Search for the interval containing v in the x sinfo_vector */
00829      loc = function1d_search_value(x, len, (pixelvalue)v, &found);
00830      if (found) {
00831         splY[j] = y[loc];
00832      } else {
00833         loc--;
00834         v -= (double)x[loc];
00835         splY[j] = (pixelvalue)( a[loc] +
00836                        v * (b[loc] +
00837                    v * (c[loc] +
00838                     v * d[loc])));
00839         }
00840     }
00841     cpl_free(a) ;
00842     return 0;
00843 }
00844 
00863 pixelvalue
00864 sinfo_function1d_average_reject(
00865     pixelvalue  *   line,
00866     int             npix,
00867     int             pix_low,
00868     int             pix_high)
00869 {
00870     pixelvalue  *   sorted ;
00871     int             i ;
00872     double          avg ;
00873 
00874     /* Sanity tests */
00875     if ((line==NULL) || (npix<1)) return (pixelvalue)0 ;
00876     if ((pix_low+pix_high)>=npix) return (pixelvalue)0 ;
00877 
00878     /* Copy input line and sort it */
00879     sorted = cpl_malloc(npix * sizeof(pixelvalue)) ;
00880     memcpy(sorted, line, npix * sizeof(pixelvalue)) ;
00881     sinfo_pixel_qsort(sorted, npix);
00882 
00883     /* Find out average of remaining values */
00884     avg = 0.00 ;
00885     for (i=pix_low+1 ; i<(npix-pix_high) ; i++) {
00886         avg += (double)sorted[i] ;
00887     }
00888     cpl_free(sorted);
00889     avg /= (double)(npix - pix_high - pix_low) ;
00890 
00891     return (pixelvalue)avg ;
00892 }
00893 
00919 #define STEP_MIN    (-half_search)
00920 #define STEP_MAX    (half_search)
00921 
00922 double 
00923 sinfo_function1d_xcorrelate(
00924     pixelvalue *    line_i,
00925     int             width_i,
00926     pixelvalue *    line_t,
00927     int             width_t,
00928     int             half_search,
00929     double     *    delta
00930 )
00931 {
00932     double * xcorr ;
00933     double   xcorr_max ;
00934     double   mean_i, mean_t ;
00935     double   rms_i, rms_t ;
00936     double   sum, sqsum ;
00937     double   norm ;
00938     int      maxpos ;
00939     int      nsteps ;
00940     int      i ;
00941     int      step ;
00942     int      nval ;
00943 
00944 
00945     /* Compute normalization factors */
00946     sum = sqsum = 0.00 ;
00947     for (i=0 ; i<width_i ; i++) {
00948         sum += (double)line_i[i] ;
00949         sqsum += (double)line_i[i] * (double)line_i[i];
00950     }
00951     mean_i = sum / (double)width_i ;
00952     sqsum /= (double)width_i ;
00953     rms_i = sqsum - mean_i*mean_i ;
00954 
00955     sum = sqsum = 0.00 ;
00956     for (i=0 ; i<width_t ; i++) {
00957         sum += (double)line_t[i] ;
00958         sqsum += (double)line_t[i] * (double)line_t[i];
00959     }
00960     mean_t = sum / (double)width_t ;
00961     sqsum /= (double)width_t ;
00962     rms_t = sqsum - mean_t*mean_t ;
00963 
00964     norm = 1.00 / sqrt(rms_i * rms_t);
00965 
00966     nsteps = (STEP_MAX - STEP_MIN) +1 ;
00967     xcorr = cpl_malloc(nsteps * sizeof(double));
00968     for (step=STEP_MIN ; step<=STEP_MAX ; step++) {
00969         xcorr[step-STEP_MIN] = 0.00 ;
00970         nval = 0 ;
00971         for (i=0 ; i<width_t ; i++) {
00972             if ((i+step > 0) &&
00973                 (i+step < width_i)) {
00974             xcorr[step-STEP_MIN] += ((double)line_t[i] - mean_t) *
00975                                     ((double)line_i[i+step] - mean_i) *
00976                                     norm ;
00977                 nval++ ;
00978             }
00979         }
00980         xcorr[step-STEP_MIN] /= (double)nval ;
00981     }
00982     xcorr_max = xcorr[0] ;
00983     maxpos    = 0 ;
00984     for (i=0 ; i<nsteps ; i++) {
00985         if (xcorr[i]>xcorr_max) {
00986             maxpos = i ;
00987             xcorr_max = xcorr[i];
00988         }
00989     }
00990     cpl_free(xcorr);
00991     (*delta) = + ((double)STEP_MIN + (double)maxpos);
00992     return xcorr_max ;
00993 }

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