00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #ifdef HAVE_CONFIG_H
00035 # include <config.h>
00036 #endif
00037
00038
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
00055
00056
00057
00058
00059
00060 #define HALF_CENTROID_DOMAIN 5
00061
00062
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
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,
00162 int npix
00163 )
00164 {
00165 pixelvalue max ;
00166 double centroid ;
00167 double weights ;
00168 int i, maxpos ;
00169
00170
00171
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
00185
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
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
00254
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
00309 out_sig = sinfo_function1d_new(samples);
00310
00311
00312 kernel = function1d_generate_smooth_kernel(filter_type, hw) ;
00313
00314
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
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
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
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
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
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
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
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
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
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
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
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
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
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;
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
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
00793 for (n = 0, i = 1; i < end; i++, n++) {
00794
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
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
00806 for (n = 0, i = 1; i < end; i++, n++) {
00807
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
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
00820 for (j = 0; j < splLen; j++) {
00821 v = (double)splX[j];
00822 splY[j] = (pixelvalue)0;
00823
00824
00825 if ((v < (double)x[0]) || (v > (double)x[end])) {
00826 continue;
00827 }
00828
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
00875 if ((line==NULL) || (npix<1)) return (pixelvalue)0 ;
00876 if ((pix_low+pix_high)>=npix) return (pixelvalue)0 ;
00877
00878
00879 sorted = cpl_malloc(npix * sizeof(pixelvalue)) ;
00880 memcpy(sorted, line, npix * sizeof(pixelvalue)) ;
00881 sinfo_pixel_qsort(sorted, npix);
00882
00883
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
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 }