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 #ifdef HAVE_CONFIG_H
00030 # include <config.h>
00031 #endif
00032 #define POSIX_SOURCE 1
00033 #include "sinfo_vltPort.h"
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #include "sinfo_detlin.h"
00044 #include "sinfo_recipes.h"
00045 #include "sinfo_fit_curve.h"
00052
00053
00054
00055
00077 cpl_imagelist *
00078 sinfo_new_fit_intensity_course(cpl_imagelist * flatStack,
00079 int order,
00080 float loReject,
00081 float hiReject )
00082 {
00083 cpl_imagelist * ret_iml ;
00084 dpoint * points ;
00085 int i, z ;
00086 double * coeffs ;
00087 Stats ** stats=NULL ;
00088 int sx;
00089 int sy;
00090 int sz;
00091 float* psrcdata;
00092 float* presdata;
00093 cpl_image* img_tmp=NULL;
00094 sx=cpl_image_get_size_x(cpl_imagelist_get(flatStack,0));
00095 sy=cpl_image_get_size_y(cpl_imagelist_get(flatStack,0));
00096 sz=cpl_imagelist_get_size(flatStack);
00097
00098 stats=(Stats**) cpl_calloc(sz,sizeof(Stats*)) ;
00099
00100 if ( NULL == flatStack )
00101 {
00102 sinfo_msg_error("no input cube given!") ;
00103 return NULL ;
00104 }
00105 if ( order <= 0 )
00106 {
00107 sinfo_msg_error("wrong order of polynomial given!") ;
00108 return NULL ;
00109 }
00110
00111 ret_iml = cpl_imagelist_new();
00112 for ( z = 0 ; z < order+1 ; z++ )
00113 {
00114 img_tmp=cpl_image_new(sx,sy,CPL_TYPE_FLOAT);
00115 cpl_imagelist_set(ret_iml,img_tmp,z);
00116 }
00117
00118 for ( z = 0 ; z < sz ; z++ )
00119 {
00120 stats[z]=
00121 sinfo_new_image_stats_on_rectangle(cpl_imagelist_get(flatStack,z),
00122 loReject,
00123 hiReject,
00124 0,
00125 0,
00126 sx-1,
00127 sy-1) ;
00128 if ( stats[z] == NULL )
00129 {
00130 sinfo_msg_error("could not compute image statistics "
00131 "in plane: %d", z) ;
00132 cpl_imagelist_delete(ret_iml) ;
00133 return NULL ;
00134 }
00135 }
00136
00137
00138
00139
00140 for ( i = 0 ; i < sx*sy ; i++ )
00141 {
00142
00143 if ( NULL == ( points = (dpoint*) cpl_calloc(sz, sizeof(dpoint)) ) )
00144 {
00145 sinfo_msg_error("could not allocate memory!\n") ;
00146 cpl_imagelist_delete(ret_iml) ;
00147 return NULL ;
00148 }
00149
00150 for ( z = 0 ; z < sz ; z++ )
00151 {
00152 if(NULL==(img_tmp = cpl_imagelist_get(flatStack,z))) {
00153 sinfo_msg_error("could not get image!");
00154 cpl_imagelist_delete(ret_iml) ;
00155 return NULL;
00156 } else {
00157 psrcdata=cpl_image_get_data_float(img_tmp);
00158 points[z].x = (double)stats[z]->cleanmean ;
00159 points[z].y = (double)psrcdata[i] ;
00160 }
00161 }
00162
00163
00164 if ( NULL == ( coeffs = sinfo_fit_1d_poly(order, points, sz, NULL) ) )
00165 {
00166 sinfo_msg_warning("could not fit spectrum of pixel: %d\n", i) ;
00167 for ( z = 0 ; z < order+1 ; z++ )
00168 {
00169 presdata=cpl_image_get_data_float(cpl_imagelist_get(ret_iml,z));
00170 presdata[i] = ZERO ;
00171 }
00172 }
00173 else
00174 {
00175 for ( z = 0 ; z < order+1 ; z++ )
00176 {
00177 if(NULL==(img_tmp = cpl_imagelist_get(ret_iml,z))) {
00178 sinfo_msg_error("could not get image!");
00179 cpl_imagelist_delete(ret_iml) ;
00180 return NULL;
00181 } else {
00182 presdata=cpl_image_get_data_float(img_tmp);
00183 presdata[i] = coeffs[z] ;
00184 }
00185 }
00186 }
00187 cpl_free(points) ;
00188 cpl_free(coeffs) ;
00189 }
00190
00191 for ( z = 0 ; z < sz ; z++ )
00192 {
00193 cpl_free (stats[z]) ;
00194 }
00195 cpl_free(stats);
00196 return ret_iml ;
00197 }
00198
00199
00223 cpl_image * sinfo_new_search_bad_pixels( cpl_imagelist * coeffs,
00224 double threshSigmaFactor,
00225 double nonlinearThresh,
00226 float loReject,
00227 float hiReject )
00228 {
00229 int i, z ;
00230 Stats * stats ;
00231 int sx=0;
00232 int sy=0;
00233 int sz=0;
00234
00235 cpl_image * img_res ;
00236 cpl_image* img_src=NULL;
00237
00238 float* psrcdata=NULL;
00239 float* presdata=NULL;
00240
00241 if ( NULL == coeffs )
00242 {
00243 sinfo_msg_error("no input cube given!\n") ;
00244 return NULL ;
00245 }
00246 if ( threshSigmaFactor <= 0. )
00247 {
00248 sinfo_msg_error("wrong sigma factor given, 0 or negativ!\n") ;
00249 return NULL ;
00250 }
00251 if ( nonlinearThresh <= 0. )
00252 {
00253 sinfo_msg_error("wrong nonlinear threshold value given, "
00254 "0 or negative!") ;
00255 return NULL ;
00256 }
00257
00258 sz=cpl_imagelist_get_size(coeffs);
00259
00260 if ( sz <= 1 )
00261 {
00262 sinfo_msg_error("no cube given, only one plane!\n") ;
00263 return NULL ;
00264 }
00265
00266
00267 img_src=cpl_imagelist_get(coeffs,1);
00268 sx=cpl_image_get_size_x(img_src);
00269 sy=cpl_image_get_size_y(img_src);
00270
00271
00272 if ( NULL == (img_res = cpl_image_new(sx, sy,CPL_TYPE_FLOAT)) )
00273 {
00274 sinfo_msg_error("could not allocate memory!\n") ;
00275 return NULL ;
00276 }
00277
00278
00279
00280
00281
00282
00283 stats = sinfo_new_image_stats_on_rectangle(img_src,
00284 loReject,
00285 hiReject, 0, 0,
00286 sx-1, sy-1) ;
00287 if ( NULL == stats )
00288 {
00289 sinfo_msg_error("could not determine image statistics!\n") ;
00290 cpl_image_delete(img_res) ;
00291 return NULL ;
00292 }
00293
00294
00295 psrcdata=cpl_image_get_data_float(img_src);
00296 presdata=cpl_image_get_data_float(img_res);
00297 for ( i = 0 ; i < (int) sx*sy ; i++ )
00298 {
00299
00300 if ( isnan(psrcdata[i]) )
00301 {
00302 presdata[i] = 0. ;
00303 }
00304 else if ( stats->cleanmean - psrcdata[i] >
00305 threshSigmaFactor*stats->cleanstdev )
00306 {
00307 presdata[i] = 0. ;
00308 }
00309 else
00310 {
00311 presdata[i] = 1. ;
00312 }
00313 }
00314 cpl_free(stats) ;
00315
00316
00317
00318
00319
00320
00321
00322
00323 if (sz > 1)
00324 {
00325 for ( z = 2 ; z < sz ; z++ )
00326 {
00327 img_src=cpl_imagelist_get(coeffs,z);
00328 sx=cpl_image_get_size_x(img_src);
00329 sy=cpl_image_get_size_y(img_src);
00330
00331 psrcdata=cpl_image_get_data_float(img_src);
00332 stats = sinfo_new_image_stats_on_rectangle(img_src, loReject,
00333 hiReject, 0, 0, sx-1, sy-1) ;
00334 if ( NULL == stats )
00335 {
00336 sinfo_msg_error("could not determine image statistics!\n") ;
00337 cpl_image_delete(img_res) ;
00338 return NULL ;
00339 }
00340 presdata=cpl_image_get_data_float(img_res);
00341 for ( i = 0 ; i < (int) sx*sy ; i++ )
00342 {
00343 if ( presdata[i] == 1. &&
00344 (fabs(psrcdata[i] - stats->cleanmean) >
00345 threshSigmaFactor*stats->cleanstdev ||
00346 fabs(psrcdata[i]) > nonlinearThresh ) )
00347 {
00348 presdata[i] = 0. ;
00349 }
00350 }
00351 cpl_free(stats) ;
00352 }
00353 }
00354
00355 return img_res ;
00356 }
00357
00358
00359
00360
00381 cpl_image * sinfo_new_search_bad_pixels_via_noise(cpl_imagelist * darks,
00382 float threshSigmaFactor,
00383 float loReject,
00384 float hiReject )
00385 {
00386 cpl_image * bp_map ;
00387 int z, n, i ;
00388 int lx, ly ;
00389 int row, col ;
00390 int low_n, high_n ;
00391 float * spectrum ;
00392 double pix_sum ;
00393 double sqr_sum ;
00394 Stats * stats ;
00395 cpl_image* img_src=NULL;
00396
00397 float* psrcdata=NULL;
00398 float* pbpdata=NULL;
00399
00400 int lz=0;
00401
00402 if ( NULL == darks )
00403 {
00404 sinfo_msg_error("no input cube given!\n") ;
00405 return NULL ;
00406 }
00407
00408 if ( threshSigmaFactor <= 0. )
00409 {
00410 sinfo_msg_error("factor is smaller or equal zero!\n") ;
00411 return NULL ;
00412 }
00413 if ( loReject < 0. || hiReject < 0. || (loReject + hiReject) >= 100. )
00414 {
00415 sinfo_msg_error("wrong reject percentage values!\n") ;
00416 return NULL ;
00417 }
00418
00419 lz=cpl_imagelist_get_size(darks);
00420 if ( lz < 1 )
00421 {
00422 sinfo_msg_error("not enough dark frames given for good statistics!") ;
00423 return NULL ;
00424 }
00425 img_src=cpl_imagelist_get(darks,0);
00426
00427 lx = cpl_image_get_size_x(img_src) ;
00428 ly = cpl_image_get_size_y(img_src) ;
00429
00430 low_n = (int)(loReject/100. *(float)lz) ;
00431 high_n = (int)(hiReject/100. *(float)lz) ;
00432 if (NULL == (bp_map = cpl_image_new (lx, ly,CPL_TYPE_FLOAT) ) )
00433 {
00434 sinfo_msg_error("could not allocate new memory!\n") ;
00435 return NULL ;
00436 }
00437 pbpdata=cpl_image_get_data(bp_map);
00438 if (NULL == (spectrum = (float*) cpl_calloc(lz, sizeof(float)) ) )
00439 {
00440 sinfo_msg_error("could not allocate new memory!\n") ;
00441 return NULL ;
00442 }
00443 for ( row = 0 ; row < ly ; row++ ) {
00444
00445 for ( col = 0 ; col < lx ; col++ ) {
00446
00447 for ( z = 0 ; z < lz ; z++ ) {
00448 img_src=cpl_imagelist_get(darks,z);
00449 psrcdata=cpl_image_get_data(img_src);
00450 spectrum[z] = psrcdata[col+lx*row] ;
00451 }
00452 sinfo_pixel_qsort(spectrum, lz) ;
00453 n = 0 ;
00454 pix_sum = 0.;
00455 sqr_sum = 0.;
00456 for ( i = low_n ; i < lz - high_n ; i++ ) {
00457 pix_sum += (double)spectrum[i] ;
00458 sqr_sum += ((double)spectrum[i]*(double)spectrum[i]) ;
00459 n++ ;
00460 }
00461
00462 pix_sum /= (double)n ;
00463 sqr_sum /= (double)n ;
00464
00465 pbpdata[col+lx*row] = (float)sqrt(sqr_sum - pix_sum*pix_sum) ;
00466 }
00467 }
00468 cpl_free(spectrum) ;
00469 if ( NULL == (stats = sinfo_new_image_stats_on_rectangle (bp_map, loReject,
00470 hiReject, 200, 200, 800, 800) ) )
00471 {
00472 sinfo_msg_error("could not get image statistics!\n") ;
00473 cpl_image_delete (bp_map) ;
00474 return NULL ;
00475 }
00476
00477
00478
00479 for ( row = 0 ; row < ly ; row++ ) {
00480 for ( col = 0 ; col < lx ; col++ ) {
00481 if (pbpdata[col+lx*row] >
00482 stats->cleanmean+threshSigmaFactor*stats->cleanstdev ||
00483 pbpdata[col+lx*row] <
00484 stats->cleanmean-threshSigmaFactor*stats->cleanstdev)
00485 {
00486 pbpdata[col+lx*row] = 0. ;
00487 }
00488 else
00489 {
00490 pbpdata[col+lx*row] = 1. ;
00491 }
00492 }
00493 }
00494 cpl_free (stats) ;
00495 return bp_map ;
00496 }
00497
00498
00499
00508 int sinfo_new_count_bad_pixels (cpl_image * bad )
00509 {
00510 int i, n ;
00511 int sx=cpl_image_get_size_x(bad);
00512 int sy=cpl_image_get_size_y(bad);
00513 float* pbpdata=cpl_image_get_data(bad);
00514
00515 n = 0 ;
00516 for ( i = 0 ; i < (int) sx*sy ; i++ )
00517 {
00518 if ( pbpdata[i] == 0 || isnan(pbpdata[i]) )
00519 {
00520 n++ ;
00521 }
00522 }
00523 return n ;
00524 }
00525
00526
00554 cpl_image * sinfo_new_abs_dist_image(cpl_image * im, float fmedian )
00555 {
00556
00557 cpl_image * image ;
00558 pixelvalue * value ;
00559 pixelvalue dist ;
00560 pixelvalue median_dist ;
00561 pixelvalue* pix_dist=NULL ;
00562 int * position ;
00563 int nposition ;
00564 int n, m, i, j ;
00565 double sum, sum2 ;
00566 double stdev ;
00567 float* pdata=NULL;
00568 int lx=0;
00569 int ly=0;
00570
00571 if ( im == NULL )
00572 {
00573 sinfo_msg_error ("no image input\n") ;
00574 return NULL ;
00575 }
00576
00577 image = cpl_image_duplicate ( im ) ;
00578
00579
00580
00581
00582
00583 sum = 0. ;
00584 sum2 = 0. ;
00585 m = 0 ;
00586
00587 pdata = cpl_image_get_data(im);
00588 lx=cpl_image_get_size_x(im);
00589 ly=cpl_image_get_size_y(im);
00590 pix_dist=(pixelvalue*)cpl_calloc(lx*ly,sizeof(pixelvalue)) ;
00591
00592 for ( i = 0 ; i < (int) lx*ly ; i++ )
00593 {
00594
00595 if ( isnan(pdata[i]) )
00596 {
00597 continue ;
00598 }
00599
00600
00601 value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ;
00602 position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ;
00603
00604
00605
00606
00607
00608 position[0] = i + lx - 1 ;
00609 position[1] = i + lx ;
00610 position[2] = i + lx + 1 ;
00611 position[3] = i + 1 ;
00612 position[4] = i - lx + 1 ;
00613 position[5] = i - lx ;
00614 position[6] = i - lx - 1 ;
00615 position[7] = i - 1 ;
00616
00617
00618
00619
00620
00621
00622
00623 if ( i >= 0 && i < lx )
00624 {
00625 position[4] += 2 * lx ;
00626 position[5] += 2 * lx ;
00627 position[6] += 2 * lx ;
00628 }
00629 else if ( i >= ((int) lx*ly - lx ) && i < (int) lx*ly )
00630 {
00631 position[0] -= 2 * lx ;
00632 position[1] -= 2 * lx ;
00633 position[2] -= 2 * lx ;
00634 }
00635 else if ( i % lx == 0 )
00636 {
00637 position[0] += 2 ;
00638 position[6] += 2 ;
00639 position[7] += 2 ;
00640 }
00641 else if ( i % lx == lx - 1 )
00642 {
00643 position[2] -= 2 ;
00644 position[3] -= 2 ;
00645 position[4] -= 2 ;
00646 }
00647
00648
00649
00650
00651
00652
00653 nposition = 8 ;
00654 n = 0 ;
00655 for ( j = 0 ; j < nposition ; j ++ )
00656 {
00657 if ( !isnan(pdata[position[j]]) )
00658 {
00659 value[n] = pdata[position[j]] ;
00660 n ++ ;
00661 }
00662 }
00663 nposition = n ;
00664
00665 if ( nposition <= 1 )
00666 {
00667 pdata[i] = ZERO ;
00668 cpl_free(value) ;
00669 cpl_free(position) ;
00670 continue ;
00671 }
00672
00673
00674 dist = 0. ;
00675 for ( n = 0 ; n < nposition ; n++ )
00676 {
00677 dist += (pdata[i] - value[n])*(pdata[i] - value[n]) ;
00678 }
00679 dist = sqrt(dist)/(float) nposition ;
00680 pix_dist[m] = dist ;
00681 m++ ;
00682 sum += (double)dist ;
00683 sum2 += (double)dist * (double)dist ;
00684 cpl_free(value) ;
00685 cpl_free(position) ;
00686 }
00687 sum /= (double)m ;
00688 sum2 /= (double)m ;
00689 stdev = sqrt(sum2 - sum*sum) ;
00690
00691 median_dist = sinfo_new_median(pix_dist, m) ;
00692
00693 for ( i = 0 ; i < (int) lx*ly ; i++ )
00694 {
00695
00696 if ( isnan(pdata[i]) )
00697 {
00698 continue ;
00699 }
00700
00701
00702 value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ;
00703 position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ;
00704
00705
00706
00707
00708
00709 position[0] = i + lx - 1 ;
00710 position[1] = i + lx ;
00711 position[2] = i + lx + 1 ;
00712 position[3] = i + 1 ;
00713 position[4] = i - lx + 1 ;
00714 position[5] = i - lx ;
00715 position[6] = i - lx - 1 ;
00716 position[7] = i - 1 ;
00717
00718
00719
00720
00721
00722
00723
00724 if ( i >= 0 && i < lx )
00725 {
00726 position[4] += 2 * lx ;
00727 position[5] += 2 * lx ;
00728 position[6] += 2 * lx ;
00729 }
00730 else if ( i >= ((int) lx*ly - lx ) && i < (int) lx*ly )
00731 {
00732 position[0] -= 2 * lx ;
00733 position[1] -= 2 * lx ;
00734 position[2] -= 2 * lx ;
00735 }
00736 else if ( i % lx == 0 )
00737 {
00738 position[0] += 2 ;
00739 position[6] += 2 ;
00740 position[7] += 2 ;
00741 }
00742 else if ( i % lx == lx - 1 )
00743 {
00744 position[2] -= 2 ;
00745 position[3] -= 2 ;
00746 position[4] -= 2 ;
00747 }
00748
00749
00750
00751
00752
00753
00754 nposition = 8 ;
00755 n = 0 ;
00756 for ( j = 0 ; j < nposition ; j ++ )
00757 {
00758 if ( !isnan(pdata[position[j]]) )
00759 {
00760 value[n] = pdata[position[j]] ;
00761 n ++ ;
00762 }
00763 }
00764 nposition = n ;
00765
00766 if ( nposition <= 1 )
00767 {
00768 pdata[i] = ZERO ;
00769 cpl_free(value) ;
00770 cpl_free(position) ;
00771 continue ;
00772 }
00773
00774
00775 dist = 0. ;
00776 for ( n = 0 ; n < nposition ; n++ )
00777 {
00778 dist += (pdata[i] - value[n])*(pdata[i] - value[n]) ;
00779 }
00780 dist = sqrt(dist)/(float) nposition ;
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796 if ( fmedian == 0 )
00797 {
00798 pdata[i] = dist ;
00799 }
00800 else if ( fmedian < 0 &&
00801 fabs ( median_dist - dist ) >= -fmedian*stdev )
00802 {
00803 pdata[i] = dist ;
00804 }
00805 else if ( fmedian > 0 &&
00806 fabs ( median_dist - dist ) >=
00807 fmedian*stdev * sqrt(fabs(dist)) )
00808 {
00809 pdata[i] = dist ;
00810 }
00811 else
00812 {
00813 cpl_free (value) ;
00814 cpl_free (position) ;
00815 continue ;
00816 }
00817
00818 cpl_free (value) ;
00819 cpl_free (position) ;
00820 }
00821 cpl_free(pix_dist);
00822 return image ;
00823 }