merge.c

00001 /*******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 * 
00005 *
00006 * who       when      what
00007 * --------  --------  ----------------------------------------------
00008 * schreib  04/07/00  created
00009 */
00010 
00011 /************************************************************************
00012 *   NAME              
00013 *        merge.c - merges the rows of two image data frames into
00014 *                  one frame with doubled column length     
00015 *        
00016 *   SYNOPSIS
00017 *   #include "merge.h"
00018 *
00019 *   a) OneImage * mergeImages ( OneImage * im1, 
00020 *                               OneImage * im2, 
00021 *                               OneImage * res_image )
00022 *
00023 *   1) OneImage * removeGeneralOffset( OneImage * im1, 
00024 *                                      OneImage * im2, 
00025 *                                      OneImage * res_image, 
00026 *                                      int n )
00027 *
00028 *   2) OneImage * removeRegionalTilt ( OneImage * im1,
00029 *                                      OneImage * im2, 
00030 *                                      OneImage * res_image )
00031 *
00032 *   3) OneImage * removeColumnOffset ( OneImage * im1, 
00033 *                                      OneImage * im2, 
00034 *                                      OneImage * res_image )
00035 *
00036 *   4) OneImage * removeResidualTilt ( OneImage * im2, OneImage * res_image )
00037 *
00038 *   5) OneImage * removeResidualOffset( OneImage * im2, OneImage * res_image )
00039 *
00040 *   DESCRIPTION
00041 *   a) merges the rows of two image frames in a way that the resulting
00042 *      image has double length in y-direction
00043 *
00044 *        The procedures are used in the SPIFFI data reduction to merge two 
00045 *        data frames. In order to fully match the two input frames there
00046 *        are five steps (procedures) foreseen:
00047 *        1) remove general offset between the frames, created by e.g. different
00048 *           air masses. 
00049 *        2) remove regional tilt between frames, created by e.g. different 
00050 *           emissivities.
00051 *        3) remove individual column offsets, created e.g. by imperfect 
00052 *           guiding, offset is divided out.
00053 *        4) remove residual individual column tilts, created by previous 
00054 *           operations.
00055 *        5) remove residual column offsets by subtracting the median.
00056 *
00057 *   FILES
00058 *
00059 *   ENVIRONMENT
00060 *
00061 *   RETURN VALUES 
00062 *   always the pointer to the image data structure OneImage
00063 *
00064 *   CAUTIONS 
00065 *
00066 *   EXAMPLES
00067 *
00068 *   SEE ALSO
00069 *   Python script merging.py
00070 *
00071 *   BUGS   
00072 *
00073 *------------------------------------------------------------------------
00074 */
00075 
00076 #include "vltPort.h"
00077 
00078 /* 
00079  * System Headers
00080  */
00081 
00082 /* 
00083  * Local Headers
00084  */
00085 
00086 #include "merge.h"
00087 
00088 /*----------------------------------------------------------------------------
00089    Function     :       mergeImages()
00090    In           :       two images to merge, dummy for residual image
00091                         note: first must have smaller wavelength than the 
00092                               second for the same pixel row.
00093    Out          :       resulting merged image, final residual image
00094    Job          :       merges the rows of two image frames in a way that the resulting
00095                         image has double length in y-direction
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     /* allocate memory */
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     /* now compute the final residual image */
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     /* now merge the two images */
00146     for ( i = 0 ; i < im1 -> ly ; i ++ )
00147     {
00148         for ( j = 0 ; j < im1 -> lx ; j ++ )
00149         {
00150             /* transfer rows to output */
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    Function     :       removeGeneralOffset()
00164    In           :       two images, number of rows from which the offset is
00165                         determined, residual image.
00166    Out          :       changed second image, residual image if wanted or needed
00167    Job          :       adds general offset between two frames to the second image
00168                         and delivers the residual image, assuming that the background 
00169                         cancellation did not work perfectly well.
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     /* allocate memory */
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      * first we determine the "good" pixels and subtract the two images 
00212      * then we determine the mean and 3 sigma
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         /* stdev is 3 sigma */
00243         stdev = 3 * sqrt(( sqr_sum - sum*mean ) / (pixelvalue)(npix - 1)) ;
00244     }
00245     
00246     /* exclude everything > 3 sigma */
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     /* now subtract the general offset which is determined as mean of the first n rows */
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     /* now apply this to the second input image */
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     /* now determine the residual image if available */
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    Function     :       removeRegionalTilt()
00311    In           :       both images to merge, residual image (obligatory no NullImage).
00312    Out          :       changed second image, residual image 
00313    Job          :       removes a general tilt from the spectra , created e.g. by
00314                         different emissivities of the telescope itself
00315                         and delivers the residual image 
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     /* allocate memory */
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      * Now work in the given difference image res_image on each column separately
00351      * This image is first smoothed columnwise by a running box
00352      */
00353     
00354     nrunning = 31 ; /* # of points in the running box, odd number required */
00355 
00356     for ( j = 0 ; j < res_image -> ly ; j ++ ) /* select a row */
00357     {
00358         for ( i = 0 ; i < res_image -> lx ; i ++ ) /* go through one row */
00359         {
00360             npix = 0 ;
00361             sum = 0. ;
00362             for (k = i - (nrunning-1)/2 ; k < i + (nrunning+1)/2; k ++ )
00363             {
00364                 /* marginal pixels are not considered */
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      * now determine the tilt in each column and remove it in such a way
00394      * that the first rows are used as references that are not changed
00395      * a free regression fit is (index i means the sum over i):
00396      * ax + b: a = [<xiyi>-<xi><yi>]/[<xi^2>-<xi>^2]
00397      * =>    : a = [xiyi - xi<yi>]/[xi^2 - xi<xi>]
00398      *         b = <yi> - a<xi>
00399      */
00400     
00401     for ( i = 0 ; i < filtered -> lx ; i ++ ) /* one column selected */
00402     {
00403         sumy = 0. ;                   /* yi   */
00404         sumc = 0. ;                   /* xiyi */
00405         sumx = 0. ;                   /* xi   */
00406         sum2 = 0. ;                   /* xi^2 */
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          * now correct the second input image im2 and the res_image.
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 ++ ) /* the same column */
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     /* now compute the final residual image */
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    Function     :       removeColumnOffset()
00469    In           :       first image and already corrected second image to merge, 
00470                         residual image (obligatory no NullImage).
00471    Out          :       changed second image, residual image 
00472    Job          :       removes individual column offset, created e.g. by imperfect
00473                         guiding. The offset is divided out. The ratio is derived
00474                         from the medians of the contributions.
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      *  now we deal with a constant offset in every column. We assume that it is
00503      *  due to redistribution of the flux. So we should divide the offset out.
00504      *  The ratio is derived from the medians of the contributions rather than the
00505      *  means.
00506      */
00507 
00508     for ( i = 0 ; i < im2 -> lx ; i ++ ) /* select a column */
00509     {
00510         /* statistics on columns */
00511         sum  = 0.  ;
00512         sum2 = 0. ;
00513         npix = 0  ;
00514         for ( j = 0 ; j < im2 -> ly ; j ++ )
00515         {
00516             /* first select only the good pixels */
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                 /* 2 sigma */
00541                 stdev = 2 * sqrt ( (sum2 - sum*mean)/(pixelvalue)(npix - 1) ) ;
00542             }
00543         }
00544 
00545         /* do it only if the S/N is high enough */
00546         if ( fabs(mean)/stdev < 0.5 )
00547         {
00548             continue ;
00549         }
00550 
00551         /* exclude everything > 2 sigma */
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         /* now deal with the offset */
00562         median1 = 0. ;
00563         median2 = 0. ;
00564         nrunning = 0 ;
00565         /* allocate memory for the column buffers */
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++ ) /* go through one column */
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         /* change the second input image only if there are more then 10 % good pixels in a column */
00585         if ( nrunning > 0.1*res_image -> ly )
00586         {
00587             /* ----------------------------------------------------------------------------
00588              * determine the medians of the columns of both images and compute the ratio, 
00589              * the columns of the second input image are multiplied by this ratio to
00590              * adjust the column offsets. 
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++ ) /* go through one column */
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     /* now compute the final residual image */
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    Function     :       removeResidualTilt()
00636    In           :       second image to merge, residual image (obligatory no NullImage).
00637    Out          :       changed second image, residual image 
00638    Job          :       removes a residual column tilt (determined from the residual
00639                         image) created by previous operations
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++ ) /* select one column */
00665     {
00666         sum  = 0. ;
00667         sum2 = 0. ;
00668         npix = 0  ;
00669         for ( j = 0 ; j < im2 -> ly ; j++ ) 
00670         {
00671             /* first select good pixels and derive the mean and sigma of each column */
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         /* exclude everything > 1.5 sigma */
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         /* now determine the tilt, see function removeRegionalTilt for explanation */
00703         sumy = 0. ;                   /* yi   */
00704         sumc = 0. ;                   /* xiyi */
00705         sumx = 0. ;                   /* xi   */
00706         sum2 = 0. ;                   /* xi^2 */
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          * now correct the second input image im2 and the res_image.
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 ++ ) /* the same column */
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    Function     :       removeResidualOffset()
00758    In           :       second image that will be changed,
00759                         residual image must be given.
00760    Out          :       changed second image, residual image
00761    Job          :       removes the residual offset by adding the median 
00762                         of the residual image to each column 
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++ ) /* select one column */
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++ ) /* go through one column */
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         /* determine the median of a column of the residual image */
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++ ) /* go through one column */
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 /*___oOo___*/

Generated on Wed Oct 26 13:08:53 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001