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
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076 #include "vltPort.h"
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 #include "merge.h"
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098 OneImage * mergeImages ( OneImage * im1,
00099 OneImage * im2,
00100 OneImage * res_image )
00101 {
00102 OneImage * out_image ;
00103 OneImage * residual ;
00104 int i, j ;
00105
00106 if ( im1 == NullImage || im2 == NullImage )
00107 {
00108 cpl_msg_error ("mergeImages:"," null image as input\n") ;
00109 return NullImage ;
00110 }
00111
00112 if ( im1 -> lx != im2 -> lx || im1 -> ly != im2 -> ly )
00113 {
00114 cpl_msg_error ("mergeImages:"," input images are not compatible in size\n") ;
00115 return NullImage ;
00116 }
00117
00118
00119 if ( NULL == (out_image = new_image (im1 -> lx, 2 * im1 ->ly)) )
00120 {
00121 cpl_msg_error ("mergeImages:"," cannot allocate new image \n") ;
00122 return NullImage ;
00123 }
00124
00125 if ( NULL == (residual = new_image (im1 -> lx, im1 ->ly)) )
00126 {
00127 cpl_msg_error ("mergeImages:"," cannot allocate new image \n") ;
00128 return NullImage ;
00129 }
00130
00131
00132 for ( i = 0 ; i < (int) im1 -> nbpix ; i ++ )
00133 {
00134 if ( qfits_isnan(im1 -> data[i]) || qfits_isnan(im2 -> data[i]) )
00135 {
00136 residual -> data[i] = ZERO ;
00137 }
00138 else
00139 {
00140 residual -> data[i] = im1 -> data[i] - im2 -> data[i] ;
00141 }
00142 res_image -> data[i] = residual -> data[i] ;
00143 }
00144
00145
00146 for ( i = 0 ; i < im1 -> ly ; i ++ )
00147 {
00148 for ( j = 0 ; j < im1 -> lx ; j ++ )
00149 {
00150
00151 out_image -> data[2*i*im1->lx + j] = im1 -> data[i*im1->lx + j] ;
00152 out_image -> data[(2*i+1) * im1->lx + j] = im2 -> data[i*im1->lx + j] ;
00153 }
00154 }
00155
00156 destroy_image (residual) ;
00157
00158 return out_image ;
00159 }
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172 OneImage * removeGeneralOffset( OneImage * im1,
00173 OneImage * im2,
00174 OneImage * res_image,
00175 int n )
00176 {
00177 OneImage * out_image ;
00178 OneImage * residual ;
00179 pixelvalue sum, sqr_sum ;
00180 pixelvalue mean, stdev ;
00181 int i, npix ;
00182
00183 if ( im1 == NullImage || im2 == NullImage )
00184 {
00185 cpl_msg_error ("removeGeneralOffset:"," null image as input\n") ;
00186 return NullImage ;
00187 }
00188
00189 if ( im1 -> lx != im2 -> lx || im1 -> ly != im2 -> ly )
00190 {
00191 cpl_msg_error ("removeGeneralOffset:"," input images are not compatible in size\n") ;
00192 return NullImage ;
00193 }
00194
00195 if ( n <= 0 )
00196 {
00197 cpl_msg_error ("removeGeneralOffset:"," number of rows for offset determination is 0 or smaller \n") ;
00198 return NullImage ;
00199 }
00200
00201
00202 if ( NULL == (residual = new_image (im1 -> lx, im1 ->ly)) )
00203 {
00204 cpl_msg_error ("removeGeneralOffset:"," cannot allocate new image \n") ;
00205 return NullImage ;
00206 }
00207
00208 out_image = copy_image( im2 ) ;
00209
00210
00211
00212
00213
00214
00215 sum = 0. ;
00216 sqr_sum = 0. ;
00217 npix = 0 ;
00218 for ( i = 0 ; i < (int) im1 -> nbpix ; i ++ )
00219 {
00220 if ( qfits_isnan(im1 -> data[i]) || qfits_isnan(im2 -> data[i]) )
00221 {
00222 residual -> data[i] = ZERO ;
00223 continue ;
00224 }
00225 else
00226 {
00227 residual -> data[i] = im1 -> data[i] - im2 -> data[i] ;
00228 }
00229
00230 sum += residual -> data[i] ;
00231 sqr_sum += (residual -> data[i]) * (residual -> data[i]) ;
00232 npix ++ ;
00233 }
00234 if ( npix <= 1 )
00235 {
00236 mean = 0. ;
00237 stdev = 0. ;
00238 }
00239 else
00240 {
00241 mean = sum / (pixelvalue) npix ;
00242
00243 stdev = 3 * sqrt(( sqr_sum - sum*mean ) / (pixelvalue)(npix - 1)) ;
00244 }
00245
00246
00247
00248 for ( i = 0 ; i < (int) residual -> nbpix ; i++ )
00249 {
00250 if ( fabs( residual -> data[i] - mean ) > stdev )
00251 {
00252 residual -> data[i] = ZERO ;
00253 }
00254 }
00255
00256
00257
00258 sum = 0. ;
00259 npix = 0 ;
00260 for ( i = 0 ; i < n * residual -> lx ; i++ )
00261 {
00262 if ( qfits_isnan(residual -> data[i]) )
00263 {
00264 continue ;
00265 }
00266
00267 sum += residual -> data[i] ;
00268 npix ++ ;
00269 }
00270 if ( npix == 0 )
00271 {
00272 mean = 0. ;
00273 }
00274 else
00275 {
00276 mean = sum / (pixelvalue) npix ;
00277 }
00278
00279
00280 for ( i = 0 ; i < (int) im2 -> nbpix ; i++ )
00281 {
00282 if ( qfits_isnan(im2 -> data[i]) )
00283 {
00284 out_image -> data[i] = ZERO ;
00285 continue ;
00286 }
00287 out_image -> data[i] = im2 -> data[i] + mean ;
00288 }
00289
00290
00291 if ( res_image != NullImage )
00292 {
00293 for ( i = 0 ; i < (int) residual -> nbpix ; i++ )
00294 {
00295 if ( qfits_isnan(residual -> data[i]) )
00296 {
00297 res_image -> data[i] = ZERO ;
00298 continue ;
00299 }
00300 res_image -> data[i] = residual -> data[i] - mean ;
00301 }
00302 }
00303
00304 destroy_image (residual) ;
00305
00306 return out_image ;
00307 }
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318 OneImage * removeRegionalTilt ( OneImage * im1,
00319 OneImage * im2,
00320 OneImage * res_image )
00321 {
00322 OneImage * out_image ;
00323 OneImage * filtered ;
00324 int i, j, k, npix, nrunning ;
00325 pixelvalue a, b, sum, sumx, sumy, sumc, sum2 ;
00326
00327 if ( im1 == NullImage || im2 == NullImage || res_image == NullImage )
00328 {
00329 cpl_msg_error ("removeRegionalTilt:"," null image as input\n") ;
00330 return NullImage ;
00331 }
00332
00333 if ( im1 -> lx != im2 -> lx || im1 -> ly != im2 -> ly ||
00334 im2 -> lx != res_image -> lx || im2 -> ly != res_image -> ly )
00335 {
00336 cpl_msg_error ("removeRegionalTilt:"," input images are not compatible in size\n") ;
00337 return NullImage ;
00338 }
00339
00340
00341 if ( NULL == ( filtered = new_image (im2 -> lx, im2 ->ly)) )
00342 {
00343 cpl_msg_error ("removeRegionalTilt:"," cannot allocate new image \n") ;
00344 return NullImage ;
00345 }
00346
00347 out_image = copy_image( im2 ) ;
00348
00349
00350
00351
00352
00353
00354 nrunning = 31 ;
00355
00356 for ( j = 0 ; j < res_image -> ly ; j ++ )
00357 {
00358 for ( i = 0 ; i < res_image -> lx ; i ++ )
00359 {
00360 npix = 0 ;
00361 sum = 0. ;
00362 for (k = i - (nrunning-1)/2 ; k < i + (nrunning+1)/2; k ++ )
00363 {
00364
00365 if ( k < 2 )
00366 {
00367 continue ;
00368 }
00369 if ( k > (res_image -> lx) - 2 )
00370 {
00371 break ;
00372
00373 }
00374 if ( qfits_isnan(res_image -> data[j*res_image->lx + k]) )
00375 {
00376 continue ;
00377 }
00378 npix ++ ;
00379 sum += res_image -> data[j*res_image->lx + k] ;
00380 }
00381 if ( npix != 0 )
00382 {
00383 filtered -> data[j*res_image->lx + i] = sum/npix ;
00384 }
00385 else
00386 {
00387 filtered -> data[j*res_image->lx + i] = ZERO ;
00388 }
00389 }
00390 }
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 for ( i = 0 ; i < filtered -> lx ; i ++ )
00402 {
00403 sumy = 0. ;
00404 sumc = 0. ;
00405 sumx = 0. ;
00406 sum2 = 0. ;
00407 npix = 0 ;
00408
00409 for ( j = 0 ; j < filtered -> ly ; j ++ )
00410 {
00411 if ( qfits_isnan(filtered -> data[i + j*filtered -> lx]) )
00412 {
00413 continue ;
00414 }
00415 sumy += filtered -> data[i + j*filtered -> lx] ;
00416 sumc += (filtered -> data[i + j*filtered -> lx]) * j ;
00417 sum2 += j*j ;
00418 sumx += j ;
00419 npix ++ ;
00420 }
00421 if ( npix > 2 && fabs(sum2 - sumx*sumx/npix) >= 1e-6 )
00422 {
00423 a = ( sumc - sumx*sumy/npix ) / ( sum2 - sumx*sumx/npix ) ;
00424 b = ( sumy - a*sumx ) / npix ;
00425 }
00426 else
00427 {
00428 a = ZERO ;
00429 b = ZERO ;
00430 }
00431
00432
00433
00434
00435
00436 if ( !qfits_isnan(a) && !qfits_isnan(b) && fabs(a) < 1e8 && fabs(b) < 1e8 )
00437 {
00438 for ( j = 0 ; j < filtered -> ly ; j ++ )
00439 {
00440 if ( !qfits_isnan(out_image->data[i + j*filtered->lx]) )
00441 {
00442 out_image -> data[i + j*filtered -> lx] += a*j+b ;
00443 }
00444 }
00445 }
00446 }
00447
00448
00449 for ( i = 0 ; i < (int) im1 -> nbpix ; i ++ )
00450 {
00451 if ( qfits_isnan(im1 -> data[i]) || qfits_isnan(out_image -> data[i]) )
00452 {
00453 res_image -> data[i] = ZERO ;
00454 }
00455 else
00456 {
00457 res_image -> data[i] = im1 -> data[i] - out_image -> data[i] ;
00458 }
00459 }
00460
00461 destroy_image (filtered) ;
00462
00463 return out_image ;
00464 }
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 OneImage * removeColumnOffset ( OneImage * im1,
00478 OneImage * im2,
00479 OneImage * res_image )
00480 {
00481 OneImage * out_image ;
00482 int i, j, npix, nrunning ;
00483 pixelvalue sum, sum2, mean, stdev, median1, median2, ratio ;
00484 pixelvalue * column1, * column2 ;
00485
00486 if ( im1 == NullImage || im2 == NullImage || res_image == NullImage )
00487 {
00488 cpl_msg_error ("removeColumnOffset:"," null image as input\n") ;
00489 return NullImage ;
00490 }
00491
00492 if ( im1 -> lx != im2 -> lx || im1 -> ly != im2 -> ly ||
00493 im2 -> lx != res_image -> lx || im2 -> ly != res_image -> ly )
00494 {
00495 cpl_msg_error ("removeColumnOffset:"," input images are not compatible in size\n") ;
00496 return NullImage ;
00497 }
00498
00499 out_image = copy_image( im2 ) ;
00500
00501
00502
00503
00504
00505
00506
00507
00508 for ( i = 0 ; i < im2 -> lx ; i ++ )
00509 {
00510
00511 sum = 0. ;
00512 sum2 = 0. ;
00513 npix = 0 ;
00514 for ( j = 0 ; j < im2 -> ly ; j ++ )
00515 {
00516
00517 if ( qfits_isnan(res_image -> data[i + j*res_image->lx]) )
00518 {
00519 continue ;
00520 }
00521 sum += res_image -> data[i + j*res_image ->lx] ;
00522 sum2 += res_image -> data[i + j*res_image ->lx] *
00523 res_image -> data[i + j*res_image ->lx] ;
00524 npix ++ ;
00525 }
00526 if ( npix <= 1 )
00527 {
00528 continue ;
00529 }
00530 else
00531 {
00532 mean = sum/(pixelvalue) npix ;
00533 if ( (sum2 - sum * mean) < 0 )
00534 {
00535 cpl_msg_error ("removeColumnOffset:"," variance is negative\n") ;
00536 continue ;
00537 }
00538 else
00539 {
00540
00541 stdev = 2 * sqrt ( (sum2 - sum*mean)/(pixelvalue)(npix - 1) ) ;
00542 }
00543 }
00544
00545
00546 if ( fabs(mean)/stdev < 0.5 )
00547 {
00548 continue ;
00549 }
00550
00551
00552 for ( j = 0 ; j < res_image -> ly ; j ++ )
00553 {
00554 if ( res_image -> data[i + j*res_image -> lx] < mean - stdev ||
00555 res_image -> data[i + j*res_image -> lx] > mean + stdev )
00556 {
00557 res_image -> data[i + j*res_image -> lx] = ZERO ;
00558 }
00559 }
00560
00561
00562 median1 = 0. ;
00563 median2 = 0. ;
00564 nrunning = 0 ;
00565
00566 column1 = (pixelvalue *) cpl_calloc ( im1 -> ly , sizeof (pixelvalue *) ) ;
00567 column2 = (pixelvalue *) cpl_calloc ( im2 -> ly , sizeof (pixelvalue *) ) ;
00568
00569 for ( j = 0 ; j < res_image -> ly ; j++ )
00570 {
00571 if ( qfits_isnan(res_image -> data[i + j*res_image -> lx]) )
00572 {
00573 continue ;
00574 }
00575 if ( qfits_isnan(im1->data[i+j*im1->lx]) || qfits_isnan(im2->data[i+j*im2->lx]) )
00576 {
00577 continue ;
00578 }
00579 column1[nrunning] = im1 -> data[i + j*im1 -> lx] ;
00580 column2[nrunning] = im2 -> data[i + j*im2 -> lx] ;
00581 nrunning ++ ;
00582 }
00583
00584
00585 if ( nrunning > 0.1*res_image -> ly )
00586 {
00587
00588
00589
00590
00591
00592 median2 = median( column2, nrunning ) ;
00593 if ( median2 != 0. )
00594 {
00595 median1 = median( column1, nrunning ) ;
00596 ratio = median1 / median2 ;
00597 if ( ratio > 0 )
00598 {
00599 for ( j = 0 ; j < im2 -> ly ; j++ )
00600 {
00601 if ( !qfits_isnan(im2->data[i + j*im2->lx]) )
00602 {
00603 out_image -> data[i + j*im2 -> lx] = im2 -> data[i + j*im2 -> lx] * ratio ;
00604 }
00605 else
00606 {
00607 out_image -> data[i + j*im2 -> lx] = ZERO ;
00608 }
00609 }
00610 }
00611 }
00612 }
00613 cpl_free ( column1 ) ;
00614 cpl_free ( column2 ) ;
00615 }
00616
00617
00618 for ( i = 0 ; i < (int) im1 -> nbpix ; i ++ )
00619 {
00620 if ( qfits_isnan(im1 -> data[i]) || qfits_isnan(out_image -> data[i]) )
00621 {
00622 res_image -> data[i] = ZERO ;
00623 }
00624 else
00625 {
00626 res_image -> data[i] = im1 -> data[i] - out_image -> data[i] ;
00627 }
00628 }
00629
00630 return out_image ;
00631 }
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642 OneImage * removeResidualTilt ( OneImage * im2, OneImage * res_image )
00643 {
00644 OneImage * out_image ;
00645 OneImage * residual ;
00646 int i, j, npix ;
00647 pixelvalue a, b, sum, sumx, sumy, sumc, sum2, mean, stdev ;
00648
00649 if ( im2 == NullImage || res_image == NullImage )
00650 {
00651 cpl_msg_error ("removeResidualTilt:"," null image as input\n") ;
00652 return NullImage ;
00653 }
00654
00655 if ( im2 -> lx != res_image -> lx || im2 -> ly != res_image -> ly )
00656 {
00657 cpl_msg_error ("removeResidualTilt:"," input images are not compatible in size\n") ;
00658 return NullImage ;
00659 }
00660
00661 out_image = copy_image( im2 ) ;
00662 residual = copy_image( res_image ) ;
00663
00664 for ( i = 0 ; i < im2 -> lx; i++ )
00665 {
00666 sum = 0. ;
00667 sum2 = 0. ;
00668 npix = 0 ;
00669 for ( j = 0 ; j < im2 -> ly ; j++ )
00670 {
00671
00672 if ( qfits_isnan(res_image -> data[i + j*res_image -> lx]) )
00673 {
00674 continue ;
00675 }
00676 sum += res_image -> data[i + j*res_image -> lx] ;
00677 sum2 += res_image -> data[i + j*res_image -> lx] *
00678 res_image -> data[i + j*res_image -> lx] ;
00679 npix ++ ;
00680 }
00681
00682 if ( npix <= 1 )
00683 {
00684 continue ;
00685 }
00686 else
00687 {
00688 mean = sum / (pixelvalue) npix ;
00689 stdev = 1.5 * sqrt( (sum2 - sum*mean) / (pixelvalue)(npix - 1) ) ;
00690 }
00691
00692
00693 for ( j = 0 ; j < im2 -> ly ; j++ )
00694 {
00695 if ( res_image -> data[i + j*res_image -> lx] < mean - stdev ||
00696 res_image -> data[i + j*res_image -> lx] > mean + stdev )
00697 {
00698 res_image -> data[i + j*res_image -> lx] = ZERO ;
00699 }
00700 }
00701
00702
00703 sumy = 0. ;
00704 sumc = 0. ;
00705 sumx = 0. ;
00706 sum2 = 0. ;
00707 npix = 0 ;
00708
00709 for ( j = 0 ; j < res_image -> ly ; j ++ )
00710 {
00711 if ( qfits_isnan(res_image -> data[i + j*res_image -> lx]) )
00712 {
00713 continue ;
00714 }
00715 sumy += res_image -> data[i + j*res_image -> lx] ;
00716 sumc += (res_image -> data[i + j*res_image -> lx]) * j ;
00717 sum2 += j*j ;
00718 sumx += j ;
00719 npix ++ ;
00720 }
00721 if ( npix > 2 && fabs(sum2 - sumx*sumx/npix) >= 1e-6 )
00722 {
00723 a = ( sumc - sumx*sumy/npix ) / ( sum2 - sumx*sumx/npix ) ;
00724 b = ( sumy - a*sumx ) / npix ;
00725 }
00726 else
00727 {
00728 a = ZERO ;
00729 b = ZERO ;
00730 }
00731
00732
00733
00734
00735
00736 if ( !qfits_isnan(a) && !qfits_isnan(b) && fabs(a) < 1e8 && fabs(b) < 1e8 )
00737 {
00738 for ( j = 0 ; j < im2 -> ly ; j ++ )
00739 {
00740 if ( !qfits_isnan(out_image->data[i+j*im2->lx]) )
00741 {
00742 out_image -> data[i + j*im2 -> lx] += a*j+b ;
00743 res_image -> data[i + j*im2 -> lx] = residual->data[i + j*im2->lx] -
00744 (a*j+b) ;
00745 }
00746 }
00747 }
00748 }
00749
00750 destroy_image (residual) ;
00751
00752 return out_image ;
00753 }
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765 OneImage * removeResidualOffset( OneImage * im2, OneImage * res_image )
00766 {
00767 OneImage * out_image ;
00768 int i, j, npix ;
00769 pixelvalue res_median ;
00770 pixelvalue * column ;
00771
00772 if ( im2 == NullImage || res_image == NullImage )
00773 {
00774 cpl_msg_error ("removeResidualOffset:"," null image as input\n") ;
00775 return NullImage ;
00776 }
00777
00778 if ( im2 -> lx != res_image -> lx || im2 -> ly != res_image -> ly )
00779 {
00780 cpl_msg_error ("removeResidualOffset:"," input images are not compatible in size\n") ;
00781 return NullImage ;
00782 }
00783
00784 out_image = copy_image( im2 ) ;
00785
00786 column = (pixelvalue *) cpl_calloc ( im2 -> ly , sizeof (pixelvalue *) ) ;
00787
00788 for ( i = 0 ; i < im2 -> lx ; i++ )
00789 {
00790 npix = 0 ;
00791 for (j=0;j<im2->ly;j++)
00792 column[j]=0;
00793
00794 for ( j = 0 ; j < res_image -> ly ; j++ )
00795 {
00796 if ( qfits_isnan(res_image -> data[i + j*res_image -> lx]) )
00797 {
00798 continue ;
00799 }
00800
00801 column[npix] = res_image -> data[i + j*res_image -> lx] ;
00802 npix ++ ;
00803 }
00804
00805
00806 if ( npix > 0.1 * res_image -> ly )
00807 {
00808 res_median = median( column, npix ) ;
00809 }
00810 else
00811 {
00812 continue ;
00813 }
00814
00815 for ( j = 0 ; j < im2 -> ly ; j++ )
00816 {
00817 if ( !qfits_isnan(im2->data[i+j*im2->lx]))
00818 {
00819 out_image -> data[i + j*im2 -> lx] = im2 -> data[i + j*im2 -> lx] + res_median ;
00820 }
00821 else
00822 {
00823 out_image -> data[i + j*im2 -> lx] = ZERO ;
00824 }
00825 if ( !qfits_isnan(res_image->data[i + j*res_image->lx]) )
00826 {
00827 res_image -> data[i + j*res_image -> lx] -= res_median ;
00828 }
00829 }
00830 }
00831 cpl_free ( column ) ;
00832 return out_image ;
00833 }
00834
00835