sinfo_new_bezier.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 * M.P.E. - SPIFFI project
00021 *
00022 *
00023 *
00024 * who       when      what
00025 
00026 * --------  --------  ----------------------------------------------
00027 * rabuter 2004-11-04 Fixed error on Find Cosmic when no valid neighboors 
00028                      where found in the subcube
00029 * rabuter 10/07/03  created
00030 */
00031 
00032 /************************************************************************
00033 *   NAME
00034 *        sinfo_new_bezier.c -
00035 *        procedures to correct for bad pixels using bezier splines 
00036 *
00037 *   SYNOPSIS
00038 *
00039 *   DESCRIPTION
00040 *
00041 *   FILES
00042 *
00043 *   ENVIRONMENT
00044 *
00045 *   RETURN VALUES
00046 *
00047 *   CAUTIONS
00048 *
00049 *   EXAMPLES
00050 *
00051 *   SEE ALSO
00052 *
00053 *   BUGS
00054 *
00055 *------------------------------------------------------------------------
00056 */
00057 
00058 #ifdef HAVE_CONFIG_H
00059 #  include <config.h>
00060 #endif
00061 #define POSIX_SOURCE 1
00062 #include "sinfo_vltPort.h"
00063 
00064 /*
00065  * System Headers
00066  */
00067 
00068 /*
00069  * Local Headers
00070  */
00071 
00072 #include "sinfo_new_bezier.h"
00073 #include "sinfo_msg.h"
00081 /*----------------------------------------------------------------------------
00082  *                          Function codes
00083  *--------------------------------------------------------------------------*/
00084 
00085 void sinfo_new_destroy_Lookup( new_Lookup *l ); 
00094 int sinfo_im_xy(cpl_image* im, int X, int Y)
00095 {
00096   int res=0;
00097 
00098    res=X+Y*cpl_image_get_size_x(im);
00099   return res;
00100 }
00108 int sinfo_im_xyz(cpl_image* im, int X, int Y, int Z)
00109 {
00110   int res=0;
00111    res = X+
00112          Y*cpl_image_get_size_x(im)+
00113          Z*cpl_image_get_size_x(im)*
00114            cpl_image_get_size_y(im);
00115   return res;
00116 }
00117 
00118 
00119 
00128 int sinfo_cu_xy(cpl_imagelist* cu, int X, int Y)
00129 {
00130   int res=0;
00131 
00132    res=X+Y*cpl_image_get_size_x(cpl_imagelist_get(cu,0));
00133   return res;
00134 }
00145 int sinfo_cu_xyz(cpl_imagelist* cu, int X, int Y, int Z)
00146 {
00147   int res=0;
00148    res = X+
00149          Y*cpl_image_get_size_x(cpl_imagelist_get(cu,0))+
00150          Z*cpl_image_get_size_x(cpl_imagelist_get(cu,0))*
00151            cpl_image_get_size_y(cpl_imagelist_get(cu,0));
00152   return res;
00153 }
00154 
00169 cpl_image * 
00170 sinfo_new_c_bezier_interpolate_image(cpl_image *im, 
00171                                      cpl_image *mask, 
00172                                      new_Lookup *look, 
00173                                            short rx, 
00174                                            short ry, 
00175                                            short rz, 
00176                            int max_rad , 
00177                                            float   ** slit_edges )
00178 {
00179     
00180     int i,j,count;
00181     cpl_imagelist * sc_im,* drs_sc_mask;
00182     cpl_image *auxImage;
00183     cpl_image *tempMask;
00184     short szx,szy,szz;
00185     short rx_loop, ry_loop, rz_loop; 
00186     /*float ant,new,dif;*/
00187 
00188     int ilx=0;
00189     int ily=0;
00190     int mlx=0;
00191     int mly=0;
00192 
00193     float* pidata=NULL;
00194     float* pmdata=NULL;
00195     float* ptdata=NULL;
00196     float* padata=NULL;
00197  
00198     cpl_image* sc_img=NULL;
00199     cpl_image* drs_img=NULL;
00200 
00201     mlx=cpl_image_get_size_x(mask);
00202     mly=cpl_image_get_size_y(mask);
00203     ilx=cpl_image_get_size_x(im);
00204     ily=cpl_image_get_size_y(im);
00205 
00206     pmdata=cpl_image_get_data_float(mask);
00207     pidata=cpl_image_get_data_float(im);
00208    
00209     if ( mlx != ilx || mly != ily )
00210     {
00211         sinfo_msg_error(" data & mask images not compatible in size\n") ;
00212         return NULL ;
00213     }
00214 
00215     /* allocate memory for sub cubes*/
00216     szx = (rx * 2 ) + 1;
00217     szy = (ry * 2 ) + 1;
00218     szz = (rz * 2 ) + 1;
00219 
00220     if ( NULL == ( sc_im = cpl_imagelist_new() ) )
00221     {
00222         sinfo_msg_error(" could not allocate memory for data subcube\n") ;
00223         return NULL ;
00224     }
00225 
00226     for(i=0;i<szz;i++) {
00227       sc_img=cpl_image_new(szx,szy,CPL_TYPE_FLOAT);
00228       cpl_imagelist_set(sc_im,sc_img,i);
00229     }
00230 
00231     if ( NULL == ( drs_sc_mask = cpl_imagelist_new() ) ) 
00232     {
00233         sinfo_msg_error(" could not allocate memory for mask subcube\n") ;
00234         return NULL ;
00235     }
00236     for(i=0;i<szz;i++) {
00237       drs_img=cpl_image_new(szx,szy,CPL_TYPE_FLOAT);
00238       cpl_imagelist_set(drs_sc_mask,drs_img,i);
00239     }
00240 
00241     if ( NULL == ( tempMask = cpl_image_new(mlx, mly, CPL_TYPE_FLOAT) ) ) 
00242     {
00243         sinfo_msg_error("could not allocate memory for temporary "
00244                         "dead pixel mask\n") ;
00245         return NULL ;
00246     }
00247     ptdata=cpl_image_get_data_float(tempMask);
00248    
00249     count=0;
00250     for ( i = 0 ; i < mlx; i++ )
00251     {
00252         for ( j = 0 ; j < mly ; j++ )
00253         {
00254         if ( pmdata[sinfo_im_xy(im,i,j)] ==  cubePT_BADPIXEL )
00255         {
00256           rx_loop = 1 ; ry_loop = 1 ; rz_loop = 1 ;
00257           pidata[sinfo_im_xy(im,i,j)] = 
00258                    sinfo_new_c_bezier_correct_pixel( i, j, im, mask, sc_im, 
00259                                drs_sc_mask, look, rx_loop, ry_loop, rz_loop );
00260           /* if not enough neighbors found, increase size of sub 
00261                      cube until max radius is reached */ 
00262           while ( pidata[sinfo_im_xy(im,i,j)] == cubeNONEIGHBOR && 
00263                           rx_loop < rx && ry_loop < ry && rz_loop < rz )
00264             {
00265               rx_loop++ ; ry_loop++; rz_loop++;  
00266               /* sinfo_msg_warning("Increasing radius to %d, in %d %d",
00267                                             rx_loop, i, j) ; */  
00268               
00269               pidata[sinfo_im_xy(im,i,j)] = 
00270                         sinfo_new_c_bezier_correct_pixel( i, j, im, mask, sc_im,
00271                                 drs_sc_mask, look, rx_loop, ry_loop, rz_loop );
00272             }
00273           /* If still not enough neighbors, make result NaN = ZERO 
00274                      in spred convention */     
00275           if ( pidata[sinfo_im_xy(im,i,j)] == cubeNONEIGHBOR ) 
00276             { 
00277               pidata[sinfo_im_xy(im,i,j)] = ZERO ;
00278             }
00279         count++;
00280         }
00281         if ( pidata[sinfo_im_xy(im,i,j)] == ZERO )
00282           {
00283         ptdata[sinfo_im_xy(tempMask,i,j)] = 0 ;
00284           }
00285          else
00286           {
00287         ptdata[sinfo_im_xy(tempMask,i,j)] = 1 ;
00288           }
00289         }
00290     }
00291 
00292     
00293     sinfo_msg("Replacing NaN\n"); 
00294     auxImage = sinfo_interpol_source_image ( im, tempMask, max_rad, slit_edges );
00295     padata=cpl_image_get_data_float(auxImage);
00296     for ( i = 0 ; i < mlx; i++ )
00297     {
00298         for ( j = 0 ; j < mly ; j++ )
00299         {
00300 
00301         if ( isnan(pidata[sinfo_im_xy(im,i,j)])) /*<= -2e10ZERO )*/
00302         {
00303           /* sinfo_msg_warning("Replacing NaN -> %d %d %f\n", 
00304                                        i,j, padata[sinfo_im_xy(im,i,j)] ); */
00305         pidata[sinfo_im_xy(im,i,j)] = padata[sinfo_im_xy(im,i,j)];
00306         }
00307         }
00308     }
00309     cpl_image_delete(auxImage);
00310     cpl_imagelist_delete(sc_im);
00311     cpl_imagelist_delete(drs_sc_mask);
00312     
00313     sinfo_msg("bad pixels count: %d\n",count);
00314 
00315     
00316     return im;
00317 }
00318 
00319 cpl_image * 
00320 sinfo_new_c_bezier_find_bad( cpl_image *im, 
00321                                    cpl_image *mask,
00322                                    /* Lookup *look,*/ 
00323                                    short rx, 
00324                                    short ry, 
00325                                    short rz,
00326                        short lowerI, 
00327                                    short highI, 
00328                                    short lowerJ, 
00329                                    short highJ, 
00330                                    float factor )
00331 {
00332     
00333     int i,j,count;
00334     cpl_imagelist * sc_im,* drs_sc_mask;
00335     short szx,szy,szz;
00336     float /*ant,*/newValue,old/*,dif,porcentage,distance*/;
00337     double med, stdev;
00338     /*cpl_image *out;*/
00339     short rx_loop, ry_loop, rz_loop; 
00340 
00341     int ilx=0;
00342     int ily=0;
00343     int mlx=0;
00344     int mly=0;
00345 
00346     float* pidata=NULL;
00347     float* pmdata=NULL;
00348  
00349     cpl_image* sc_img=NULL;
00350     cpl_image* drs_img=NULL;
00351 
00352 
00353     mlx=cpl_image_get_size_x(mask);
00354     mly=cpl_image_get_size_y(mask);
00355     ilx=cpl_image_get_size_x(im);
00356     ily=cpl_image_get_size_y(im);
00357 
00358     pmdata=cpl_image_get_data_float(mask);
00359     pidata=cpl_image_get_data_float(im);
00360 
00361 
00362     if ( mlx != ilx || mly != ily )
00363     {
00364         sinfo_msg_error(" data & mask images not compatible in size\n") ;
00365         return NULL ;
00366     }
00367 
00368     /* allocate memory for sub cubes*/
00369     szx = (rx * 2 ) + 1;
00370     szy = (ry * 2 ) + 1;
00371     szz = (rz * 2 ) + 1;
00372 
00373     if ( NULL == ( sc_im = cpl_imagelist_new() ) ) 
00374     {
00375         sinfo_msg_error(" could not allocate memory for data subcube\n") ;
00376         return NULL ;
00377     }
00378     for(i=0;i<szz;i++) {
00379       sc_img=cpl_image_new(szx,szy,CPL_TYPE_FLOAT);
00380       cpl_imagelist_set(sc_im,sc_img,i);
00381     }
00382 
00383     if ( NULL == ( drs_sc_mask = cpl_imagelist_new() ) ) 
00384     {
00385         sinfo_msg_error(" could not allocate memory for mask subcube\n") ;
00386         return NULL ;
00387     }
00388     for(i=0;i<szz;i++) {
00389       drs_img=cpl_image_new(szx,szy,CPL_TYPE_FLOAT);
00390       cpl_imagelist_set(drs_sc_mask,drs_img,i);
00391     }
00392     
00393     count=0;
00394     for ( i = 0 ; i < mlx; i++ )
00395     {
00396         for ( j = 0 ; j < mly ; j++ )
00397         {
00398         if ( i >= lowerI && i < highI && j >= lowerJ && j < highJ )
00399         {
00400 
00401         rx_loop = 1 ; ry_loop = 1 ; rz_loop = 1 ;
00402         newValue = sinfo_new_c_bezier_correct_pixel_2D( i, j, im, 
00403                                                                 mask, 
00404                                                                 sc_im, 
00405                                                                 drs_sc_mask,
00406                                                                 /* look,*/
00407                                                                 rx_loop, 
00408                                                                 ry_loop, 
00409                                                                 rz_loop, 
00410                                                                 &med, 
00411                                                                 &stdev,
00412                                                                 factor );
00413             /* if NaN returned, increase size of sub cube 
00414                    until max radius is reached */ 
00415         while ( newValue == ZERO && rx_loop < rx && 
00416                         ry_loop < ry && rz_loop < rz )
00417             {
00418               rx_loop++ ; ry_loop++; rz_loop++;  
00419               /*sinfo_msg_warning("Increasing radius to %d, 
00420                         in %d %d", rx_loop, i, j) ;  */
00421                   newValue = sinfo_new_c_bezier_correct_pixel_2D( i, j, im,
00422                                                                       mask, 
00423                                                                       sc_im, 
00424                                                                     drs_sc_mask,
00425                                                                     /*, look*/
00426                                                                     rx_loop, 
00427                                                                     ry_loop, 
00428                                                                     rz_loop, 
00429                                                                     &med, 
00430                                                                     &stdev,
00431                                                                     factor );
00432             }
00433         if ( isnan(newValue)) /*<= -3.e10 ZERO )*/
00434             continue;
00435         
00436         old = pidata[sinfo_im_xy(im,i,j)];
00437         if ( newValue != old )
00438             {
00439             pidata[sinfo_im_xy(im,i,j)] = newValue;
00440             /*sinfo_msg_warning("[%d,%d]=%f -> %f, med= %f, 
00441                       stdev=%f\n",i,j, old, newValue, med, stdev );*/
00442             count++;
00443             }
00444         }
00445         }
00446     }
00447     
00448     
00449     sinfo_msg("bad pixels count: %d\n",count);
00450 
00451     
00452     cpl_imagelist_delete(sc_im);
00453     cpl_imagelist_delete(drs_sc_mask);
00454     return im;
00455 }
00456 
00457 float 
00458 sinfo_new_c_bezier_correct_pixel(int ipos, 
00459                                  int jpos, 
00460                                  cpl_image * im, 
00461                                  cpl_image * mask, 
00462                      cpl_imagelist * sc_im, 
00463                                  cpl_imagelist * drs_sc_mask, 
00464                                  new_Lookup * look, 
00465                                  short rx, 
00466                                  short ry, 
00467                                  short rz )
00468 {
00469     short ic, jc, kc, ii, jj, kk/*, sjj, skk*/,is,js,ks;
00470     short i, j, k, indexJ, indexI, lx, ly, lz, szx, szy, szz;
00471     /*float indexIf,indexJf,sp;*/
00472     cpl_image * X, * Y, * Z, * hX;
00473     cpl_imagelist  * id, * jd;
00474 
00475     int idlx=0;
00476     int idly=0;
00477     int idnp=0;
00478 
00479     int drslx=0;
00480     int drsly=0;
00481     int drsnp=0;
00482 
00483 
00484     float* pXdata=NULL;
00485     float* pYdata=NULL;
00486     float* pZdata=NULL;
00487     float* phXdata=NULL;
00488     float* pidata=NULL;
00489     float* pmdata=NULL;
00490     float* piddata=NULL;
00491     float* pjddata=NULL;
00492     float* pscdata=NULL;
00493     float* pdrsdata=NULL;
00494 
00495     cpl_image* id_img=NULL;
00496     cpl_image* jd_img=NULL;
00497     cpl_image* sc_img=NULL;
00498     cpl_image* drs_img=NULL;
00499 
00500     X  = look -> X;
00501     Y  = look -> Y;
00502     Z  = look -> Z;
00503     hX = look -> hX;
00504     id = look -> id;
00505     jd = look -> jd;
00506 
00507     /*
00508       phXdata=cpl_image_get_data_float(hX);
00509       if ( phXdata[sinfo_im_xy( hX, ipos, jpos)] > 1 )
00510     {
00511     sinfo_msg_error(" double hit in position [%d,%d]=%f, "
00512                         "can not correct\n", 
00513         ipos,jpos,phXdata[sinfo_im_xy(hX,ipos,jpos)]) ;
00514     return ( -2e10 );
00515     }*/
00516 
00517     pidata=cpl_image_get_data_float(im);
00518     pmdata=cpl_image_get_data_float(mask);
00519 
00520     phXdata=cpl_image_get_data_float(hX);
00521     if ( phXdata[sinfo_im_xy( hX, ipos, jpos)] < 1 )
00522     {
00523     /*sinfo_msg_error("no lookup in position [%d,%d]=%f, can not correct", 
00524       ipos,jpos,phXdata[sinfo_im_xy(hX,ipos,jpos)]) ;*/
00525     return ( ZERO );
00526     }
00527     pXdata=cpl_image_get_data_float(X);
00528     pYdata=cpl_image_get_data_float(Y);
00529     pZdata=cpl_image_get_data_float(Z);
00530     
00531 
00532     ic = pXdata[sinfo_im_xy( X, ipos, jpos)];
00533     jc = pYdata[sinfo_im_xy( Y, ipos, jpos)];
00534     kc = pZdata[sinfo_im_xy( Z, ipos, jpos)];
00535     /*if ( !(ipos % 16 )  )*/
00536 #ifdef DEBUG    
00537     sinfo_msg_debug("Correcting bad pixel : ipos=%d,jpos=%d, "
00538                     "in Cube -> ic=%d, jc=%d, kc=%d\n", 
00539                     ipos,jpos, ic, jc, kc );
00540 #endif
00541     /*limit to start not before the beginning of the cube*/
00542     ii = ic - rx; if ( ii < 0 ) ii = 0;
00543     jj = jc - ry; if ( jj < 0 ) jj = 0;
00544     kk = kc - rz; if ( kk < 0 ) kk = 0;
00545 
00546 #ifdef DEBUG
00547     sinfo_msg_debug("Start Point in Cube -> ii=%d,jj=%d,kk=%d\n", ii, jj, kk );
00548 #endif
00549     
00550     /*limit to end not outside of the cube */
00551     szx = (rx * 2 ) + 1;
00552     szy = (ry * 2 ) + 1;
00553     szz = (rz * 2 ) + 1;
00554 
00555     idlx=cpl_image_get_size_x(cpl_imagelist_get(id,0));
00556     idly=cpl_image_get_size_y(cpl_imagelist_get(id,0));
00557     idnp=cpl_imagelist_get_size(id);
00558 
00559     lx = idlx;
00560     ly = idly;
00561     lz = idnp;
00562 
00563     if ( ( ic + rx ) >= idlx )
00564     szx = szx - ( (ic+rx)-(lx-1) );
00565 
00566     if ( ( jc + ry ) >= idly )
00567     szy = szy - ( (jc+ry)-(ly-1) );
00568 
00569     if ( ( kc + rz ) >= idnp )
00570     szz = szz - ( (kc+rz)-(lz-1) );
00571 
00572     drslx=cpl_image_get_size_x(cpl_imagelist_get(drs_sc_mask,0));
00573     drsly=cpl_image_get_size_y(cpl_imagelist_get(drs_sc_mask,0));
00574     drsnp=cpl_imagelist_get_size(drs_sc_mask);
00575 #ifdef DEBUG
00576 
00577 
00578     sinfo_msg_error("Size of subcube: szx=%d,szy=%d,szz=%d\n", szx, szy, szz );
00579     /*fill whole mask with not available*/
00580     sinfo_msg_error("Fill Mask subcube of size: %d,%d,%d, with NOINFO\n", 
00581                      drslx, drsly,  drsnp);
00582 #endif
00583     for( i = 0; i < drslx; i++) {
00584       for( j = 0; j < drsly; j++) {
00585     for( k = 0; k < drsnp; k++) {
00586       drs_img=cpl_imagelist_get(drs_sc_mask,k);
00587           pdrsdata=cpl_image_get_data_float(drs_img);
00588       pdrsdata[sinfo_cu_xy(drs_sc_mask,i,j)] = cubePT_NOINFO;
00589     }
00590       }
00591     }
00592 
00593      for( i = ii,is=0;  i < ii+szx; i++,is++)
00594      {
00595      for( j = jj,js=0;  j < jj+szy; j++,js++)
00596          {
00597          for( k = kk,ks=0;  k < kk+szz; k++,ks++)
00598          {
00599 #ifdef DEBUG
00600          sinfo_msg_debug("i=%d j=%d k=%d is=%d ij=%d ik=%d",
00601                                   i,j,k,is,js,ks);
00602 #endif
00603                  id_img=cpl_imagelist_get(id,k);
00604                  jd_img=cpl_imagelist_get(jd,k);
00605                  piddata=cpl_image_get_data_float(id_img);
00606                  pjddata=cpl_image_get_data_float(jd_img);
00607 
00608                  drs_img=cpl_imagelist_get(drs_sc_mask,ks);
00609                  pdrsdata=cpl_image_get_data_float(drs_img);
00610                  sc_img=cpl_imagelist_get(sc_im,ks);
00611                  pscdata=cpl_image_get_data_float(sc_img);
00612 
00613          indexI = sinfo_new_nint( piddata[sinfo_cu_xy(id,i,j)] ); 
00614          indexJ = sinfo_new_nint( pjddata[sinfo_cu_xy(jd,i,j)] );
00615          if ( indexJ <= -1 || indexJ>=2048 || indexI == -1 )
00616            {
00617              pdrsdata[sinfo_cu_xy(drs_sc_mask,is,js)] = cubePT_NOINFO;
00618              continue;
00619            }
00620          pscdata[sinfo_cu_xy(sc_im,is,js)]  = 
00621                          pidata[sinfo_im_xy(im,indexI,indexJ)];
00622          pdrsdata[sinfo_cu_xy(drs_sc_mask,is,js)]  = 
00623                          pmdata[sinfo_im_xy(mask,indexI,indexJ)];
00624 #ifdef DEBUG
00625         sinfo_msg_debug("Cube i=%d, j=%d, k=%d  ;"
00626                                 "  Sub is=%d, js=%d, ks=%d  ;"
00627                                 "  Plane I=%d,J=%d ; mask %f ; im %f", 
00628                         i, j, k, is, js, ks, indexI, indexJ,  
00629                                 mask -> data[sinfo_im_xy(mask,indexI,indexJ)], 
00630                                 im   -> data[sinfo_im_xy(im,indexI,indexJ)]);
00631 #endif
00632 
00633          }
00634          }
00635      }
00636 
00637 
00638     /*signal to correct this pixel*/
00639      drs_img=cpl_imagelist_get(drs_sc_mask,rz);
00640      pdrsdata=cpl_image_get_data_float(drs_img);
00641     pdrsdata[sinfo_cu_xy(drs_sc_mask,rx,ry)] = cubePT_FIND;
00642     return ( sinfo_new_c_bezier_interpol( sc_im, drs_sc_mask ) );
00643 }
00644 
00645 float 
00646 sinfo_new_c_bezier_correct_pixel_2D(int ipos, 
00647                                     int jpos, 
00648                                     cpl_image * im, 
00649                                     cpl_image * mask, 
00650                         cpl_imagelist * sc_im, 
00651                                     cpl_imagelist * drs_sc_mask,
00652                                     /* Lookup * look,*/ 
00653                                     short rx, short ry, 
00654                                     short rz , 
00655                                     double *med , 
00656                         double *stdev, 
00657                                     float factor )
00658 {
00659     short ic, jc, kc, ii, jj, kk/*, sjj, skk*/,is,js,ks;
00660     short i, j, k, indexJ, indexI, lx, ly, lz, szx, szy, szz;
00661     double sum,aux;
00662     int counter;
00663     float sumarr[100];
00664 
00665 
00666     int ilx=0;
00667     int ily=0;
00668 
00669     int drslx=0;
00670     int drsly=0;
00671     int drsnp=0;
00672 
00673 
00674     float* pidata=0;
00675     float* pmdata=0;
00676     float* pscdata=0;
00677     float* pdrsdata=0;
00678 
00679     cpl_image* drs_img=NULL;
00680     cpl_image* sc_img=NULL;
00681 
00682     jc = 0;
00683     ic = ipos;
00684     kc = jpos;
00685     sinfo_msg_debug("Correcting bad pixel : ipos=%d,jpos=%d, "
00686                     "in Cube -> ic=%d, jc=%d, kc=%d", ipos,jpos, ic, jc, kc );
00687     /*limit to start not before the beginning of the cube*/
00688     ii = ic - rx; if ( ii < 0 ) ii = 0;
00689     jj = jc - ry; if ( jj < 0 ) jj = 0;
00690     kk = kc - rz; if ( kk < 0 ) kk = 0;
00691 
00692     sinfo_msg_debug("Start Point in Cube -> ii=%d,jj=%d,kk=%d", ii, jj, kk );
00693     
00694     ilx=cpl_image_get_size_x(im);
00695     ily=cpl_image_get_size_y(im);
00696 
00697     /*limit to end not outside of the cube */
00698     szx = (rx * 2 ) + 1;
00699     szy = (ry * 2 ) + 1;
00700     szz = (rz * 2 ) + 1;
00701     lx = ilx;
00702     ly = ily;
00703     lz = ily;
00704     if ( ( ic + rx ) >= ilx )
00705     szx = szx - ( (ic+rx)-(lx-1) );
00706 
00707     if ( ( jc + ry ) >= ily )
00708     szy = szy - ( (jc+ry)-(ly-1) );
00709 
00710     if ( ( kc + rz ) >= ily )
00711     szz = szz - ( (kc+rz)-(lz-1) );
00712 
00713 #ifdef DEBUG
00714     drslx=cpl_image_get_size_x(cpl_imagelist_get(drs_sc_mask,0));
00715     drsly=cpl_image_get_size_y(cpl_imagelist_get(drs_sc_mask,0));
00716     drsnp=cpl_imagelist_get_size(drs_sc_mask);
00717     sinfo_msg_debug("Size of subcube : szx=%d,szy=%d,szz=%d", szx, szy, szz );
00718     /*fill whole mask with not available*/
00719     sinfo_msg_debug("Fill Mask subcube of size:%d,%d,%d, with NOINFO", 
00720                     drslx, drsly,  drsnp);
00721 #endif
00722     for( i = 0; i < drslx; i++) {
00723       for( j = 0; j < drsly; j++) {
00724     for( k = 0; k < drsnp; k++) {
00725       drs_img=cpl_imagelist_get(drs_sc_mask,k);
00726       pdrsdata=cpl_image_get_data_float(drs_img);
00727       pdrsdata[sinfo_cu_xy(drs_sc_mask,i,j)] = cubePT_NOINFO;
00728     }
00729       }
00730     }
00731     counter = 0;
00732     sum=0;
00733     memset(sumarr,0x00,sizeof(sumarr));
00734     pidata=cpl_image_get_data(im);
00735     pmdata=cpl_image_get_data(mask);
00736 
00737      for( i = ii,is=0;  i < ii+szx; i++,is++)
00738      {
00739      for( j = jj,js=0;  j < jj+szy; j++,js++)
00740          {
00741          for( k = kk,ks=0;  k < kk+szz; k++,ks++)
00742          {
00743 #ifdef DEBUG
00744          sinfo_msg_debug("i=%d j=%d k=%d is=%d ij=%d ik=%d",
00745                                   i,j,k,is,js,ks);
00746 #endif
00747          indexI = i; 
00748          indexJ = k;
00749          if ( isnan(pidata[sinfo_im_xy(mask,indexI,indexJ)]) )
00750            pmdata[sinfo_im_xy(mask,indexI,indexJ)] = 0;
00751 
00752          if ( pmdata[sinfo_im_xy(mask,indexI,indexJ)] == 1 && 
00753                             ( indexI != ipos || indexJ != jpos) )
00754            {
00755              /*sumarr[counter] = pidata[sinfo_im_xy(im,indexI,indexJ)];*/
00756              sum = sum + pidata[sinfo_im_xy(im,indexI,indexJ)];
00757              counter++;
00758            }
00759          sc_img=cpl_imagelist_get(sc_im,ks);
00760          pscdata[sinfo_cu_xy(sc_im,is,js)]=
00761                         pidata[sinfo_im_xy(im,indexI,indexJ)];
00762                  drs_img=cpl_imagelist_get(drs_sc_mask,ks);
00763                  pdrsdata=cpl_image_get_data_float(drs_img);
00764          pdrsdata[sinfo_cu_xy(drs_sc_mask,is,js)]=
00765                          pmdata[sinfo_im_xy(mask,indexI,indexJ)];
00766 #ifdef DEBUG
00767          sinfo_msg_debug("Cube i=%d, j=%d, k=%d  ;  "
00768                                  "Sub is=%d, js=%d, ks=%d  ; "
00769                                 " Plane I=%d,J=%d ; mask %f ; im %f", 
00770                         i, j, k, is, js, ks, indexI, indexJ,
00771                                 pmdata[sinfo_im_xy(mask,indexI,indexJ)], 
00772                                 pidata[sinfo_im_xy(im,indexI,indexJ)]);
00773 #endif
00774 
00775          }
00776          }
00777      }
00778 
00779 
00780     /*signal to correct this pixel*/
00781     drs_img=cpl_imagelist_get(drs_sc_mask,rz);
00782     pdrsdata=cpl_image_get_data_float(drs_img);
00783     pdrsdata[sinfo_cu_xy(drs_sc_mask,rx,ry)] = cubePT_FIND;
00784     if ( counter )
00786     *med = sum/counter;
00787     else
00788     return(pidata[sinfo_im_xy(im,ipos,jpos)]);
00789 
00790     /*sinfo_msg_debug("%f %f %d\n",
00791                       sum ,pidata[sinfo_im_xy(im,ipos,jpos)], counter);*/
00792 
00793 
00794 
00795     sum =0;
00796     counter=0;
00797     for( i = ii,is=0;  i < ii+szx; i++,is++)
00798       {
00799     for( j = jj,js=0;  j < jj+szy; j++,js++)
00800       {
00801         for( k = kk,ks=0;  k < kk+szz; k++,ks++)
00802           {
00803         drs_img=cpl_imagelist_get(drs_sc_mask,ks);
00804         pdrsdata=cpl_image_get_data_float(drs_img);
00805         indexI = i;
00806         indexJ = k;
00807         if ( pdrsdata[sinfo_cu_xy(drs_sc_mask,is,js)] == 1 && 
00808                          ( indexI != ipos || indexJ != jpos) )
00809           {
00810               sc_img=cpl_imagelist_get(sc_im,ks);
00811               pscdata=cpl_image_get_data_float(sc_img);
00812 
00813             sum=sum+((pscdata[sinfo_cu_xy(drs_sc_mask,is,js)]- *med) * 
00814              (pscdata[sinfo_cu_xy(drs_sc_mask,is,js)] - *med ) );
00815             counter++;
00816           }
00817           }
00818       }
00819       }
00820     
00821     aux = sum;
00822     sum   = sum / (counter - 1);
00823     *stdev = sqrt( sum );
00824     
00825     if ( (fabs( pidata[sinfo_im_xy(im,ipos,jpos)] - *med ) > factor * *stdev) || 
00826     isnan(pidata[sinfo_im_xy(im,ipos,jpos)]) )
00827     {
00828     /*sinfo_msg_debug("[%d,%d]: distance to mean = %f,"
00829                           " thres =%f sum=%f, stdev=%f, counter=%d, aux= %f", 
00830                 ipos,jpos, fabs( pidata[sinfo_im_xy(im,ipos,jpos)] - *med),
00831                     factor * *stdev, sum,*stdev, counter,aux );
00832     pmdata[sinfo_im_xy(mask,ipos,jpos)] = 0;*/
00833     return ( sinfo_new_c_bezier_interpol( sc_im, drs_sc_mask ) );
00834     }
00835     return(pidata[sinfo_im_xy(im,ipos,jpos)]); 
00836 }
00837 
00838 
00839 
00840 float 
00841 sinfo_new_c_bezier_interpol( cpl_imagelist * im, cpl_imagelist * action )
00842 {
00843     short pos;
00844     unsigned short i,j,k;
00845     new_XYZW  indata[1000],res,selected;
00846     float step,cumstep,distance,selected_distance;
00847     new_Dim point;
00848     double munk;
00849     int ilx=0;
00850     int ily=0;
00851     int inp=0;
00852    
00853     float* padata=NULL;
00854     float* pidata=NULL;
00855     cpl_image* i_img=NULL;
00856     cpl_image* a_img=NULL;
00857 
00858     memset(indata,0x00,1000*sizeof(new_XYZW));
00859     ilx=cpl_image_get_size_x(cpl_imagelist_get(im,0));
00860     ily=cpl_image_get_size_y(cpl_imagelist_get(im,0));
00861     inp=cpl_imagelist_get_size(im);
00862 
00863     pos=0;
00864     for( i=0; i < ilx; i++)
00865     {
00866     for( j=0; j < ily; j++)
00867         {
00868         for( k=0; k < inp; k++)
00869         {
00870           a_img=cpl_imagelist_get(action,k);
00871                   padata=cpl_image_get_data_float(a_img);
00872           i_img=cpl_imagelist_get(action,k);
00873                   pidata=cpl_image_get_data_float(i_img);
00874         if ( padata[sinfo_cu_xy(action,i,j)] == cubePT_USE )
00875             {
00876 #ifdef DEBUG
00877             sinfo_msg_debug("Used im[%d,%d,%d]=%lf\n",
00878                                       i,j,k,pidata[sinfo_im_xy(im,i,j)]);
00879 #endif
00880             indata[pos].x = i;
00881             indata[pos].y = j;
00882             indata[pos].z = k;
00883             indata[pos].w = pidata[sinfo_cu_xy(im,i,j)];
00884             pos++;
00885             }
00886         else
00887             {
00888             if ( padata[sinfo_cu_xy(action,i,j)] == cubePT_FIND )
00889             {
00890             point.x = i;
00891             point.y = j;
00892             point.z = k;
00893 #ifdef DEBUG
00894             sinfo_msg_debug("Find for im[%d,%d,%d]=%lf reason:%f",
00895                         i,j,k,pidata[sinfo_im_xy(im,i,j)],
00896                         padata[sinfo_cu_xy(action,i,j)]);
00897 #endif
00898             }
00899             else
00900             {
00901 #ifdef DEBUG
00902             sinfo_msg_debug("Ignored im[%d,%d,%d]=%lf reason:%f",
00903                         i,j,k,pidata[sinfo_im_xy(im,i,j)],
00904                         padata[sinfo_im_xy(action,i,j)]);
00905 #endif
00906             }
00907             }
00908         }   
00909         }   
00910     }
00911     
00912     
00913     if ( pos < 2 )
00914     {
00915 #ifdef DEBUG
00916        sinfo_msg_debug("subcube contains no valid pixels "
00917                        "to use in iterpolation");
00918 #endif
00919     /*i_img=cpl_imagelist_get(point.z);
00920           pidata=cpl_image_get_data_float(i_img);
00921           return( pidata[sinfo_im_xy(im,point.x,point.y)] );*/
00922     return( cubeNONEIGHBOR );
00923 
00924     }
00925     
00926 
00927     step    = 0.01;
00928     cumstep = 0.0;
00929     selected_distance=1000;
00930     munk = pow( 1.0-cumstep, (double)pos - 1 );
00931     for ( i = 0 ; ( i < 100 ) && ( munk != 0.0 ); i++ )
00932     {
00933     memset( &res, 0x00, sizeof(new_XYZW) );
00934     sinfo_new_bezier( indata, pos-1, cumstep, munk, &res);
00935     distance = sqrt( pow( (point.x-res.x), 2)+ 
00936                          pow( (point.y-res.y), 2)+ 
00937                          pow( (point.z-res.z), 2) );
00938     /*sinfo_msg_debug("%lf %lf %lf %lf %lf\n",
00939                           res.x,res.y,res.z,res.w,distance);*/
00940     if ( distance < selected_distance )
00941         {
00942         selected_distance = distance;
00943         selected.x = res.x;
00944         selected.y = res.y;
00945         selected.z = res.z;
00946         selected.w = res.w;
00947         }
00948     cumstep = cumstep + step;
00949     munk = pow( 1.0 - cumstep, (double)pos - 1 );
00950     
00951     }
00952     
00953 #ifdef DEBUG
00954    sinfo_msg_debug("Selected %lf %lf %lf %lf, distance=%lf",
00955                     selected.x,selected.y,selected.z,
00956                     selected.w,selected_distance);
00957 #endif
00958     i_img=cpl_imagelist_get(im,point.z);
00959     pidata=cpl_image_get_data_float(i_img);
00960     pidata[sinfo_cu_xy(im,point.x,point.y)] = selected.w;
00961     
00962     return selected.w;
00963 }
00964 
00965 
00966 
00967 int 
00968 sinfo_new_bezier(new_XYZW *p,int n,double mu,double munk,new_XYZW *res )
00969 {
00970    int k, kn, nn, nkn;
00971    double blend, muk;
00972 
00973    muk = 1;
00974    for ( k = 0; k <= n; k++ ) {
00975       nn = n;
00976       kn = k;
00977       nkn = n - k;
00978       blend = muk * munk;
00979       muk *= mu;
00980       munk /= ( 1.0 - mu );
00981       while ( nn >= 1 ) {
00982          blend *= (double)nn;
00983          nn--;
00984          if ( kn > 1 ) {
00985             blend /= (double)kn;
00986             kn--;
00987          }
00988          if ( nkn > 1 ) {
00989             blend /= (double)nkn;
00990             nkn--;
00991          }
00992       }
00993       res -> x += p[k].x * blend;
00994       res -> y += p[k].y * blend;
00995       res -> z += p[k].z * blend;
00996       res -> w += p[k].w * blend;
00997    }
00998    return( 0 );
00999 }
01000 
01001 int 
01002 sinfo_new_c_create_XYZ( new_Lookup *l )
01003 {
01004     cpl_image *imX,*imY,*imZ,*imcX;
01005     short i,j,k,indexI,indexJ,x,y,z;
01006     int size;
01007     int idlx=0;
01008     int idly=0;
01009     int idnp=0;
01010     float* piddata=NULL;
01011     float* pjddata=NULL;
01012     float* pXdata=NULL;
01013     float* pYdata=NULL;
01014     float* pZdata=NULL;
01015     float* phXdata=NULL;
01016 
01017     cpl_image* i_img=NULL;
01018     cpl_image* j_img=NULL;
01019 
01020     idlx=cpl_image_get_size_x(cpl_imagelist_get(l->id,0));
01021     idly=cpl_image_get_size_y(cpl_imagelist_get(l->id,0));
01022     idnp=cpl_imagelist_get_size(l->id);
01023 
01024     size = idlx*idly;
01025     /* allocate memory */
01026     if ( NULL == (imX = cpl_image_new(size, size, CPL_TYPE_FLOAT)) ) 
01027         {
01028     sinfo_msg_error(" could not allocate memory for X !\n") ;
01029     return -1 ;
01030         }
01031     if ( NULL == (imY = cpl_image_new(size, size, CPL_TYPE_FLOAT)) ) 
01032         {
01033     sinfo_msg_error(" could not allocate memory for Y !\n") ;
01034     return -1 ;
01035         }
01036     if ( NULL == (imZ = cpl_image_new(size, size, CPL_TYPE_FLOAT)) ) 
01037         {
01038     sinfo_msg_error(" could not allocate memory for Z !\n") ;
01039     return -1 ;
01040         }
01041     if ( NULL == (imcX = cpl_image_new(size, size, CPL_TYPE_FLOAT)) ) 
01042         {
01043     sinfo_msg_error(" could not allocate memory for cX !\n") ;
01044     return -1 ;
01045         }
01046 
01047     l -> X  = imX; 
01048     l -> Y  = imY; 
01049     l -> Z  = imZ; 
01050     l -> hX = imcX; 
01051 
01052     /*Round id*/
01053     for( i = 0; i < idlx; i++)
01054     {
01055     for( j = 0; j < idly; j++)
01056         {
01057         for( k = 0; k < idnp; k++)
01058         {
01059           i_img=cpl_imagelist_get(l->id,k);
01060           piddata=cpl_image_get_data_float(i_img);
01061           piddata[sinfo_cu_xy(l->id,i,j)] = 
01062                   (float)sinfo_new_nint(piddata[sinfo_cu_xy(l->id,i,j)]);
01063         }
01064         }
01065     }
01066 
01067     /*Round jd*/
01068     for( i = 0; i < idlx; i++)
01069     {
01070     for( j = 0; j < idly; j++)
01071         {
01072         for( k = 0; k < idnp; k++)
01073         {
01074           j_img=cpl_imagelist_get(l->jd,k);
01075           pjddata=cpl_image_get_data_float(j_img);
01076           pjddata[sinfo_cu_xy(l->jd,i,j)] = 
01077                  (float)sinfo_new_nint(pjddata[sinfo_cu_xy(l->jd,i,j)]);
01078         }
01079         }
01080     }
01081 
01082     /*Fill with -1 X Y Z*/
01083     for( i = 0; i < cpl_image_get_size_x(l -> X); i++)
01084     {
01085     for( j = 0; j < cpl_image_get_size_y(l -> X); j++)
01086         {
01087           pXdata=cpl_image_get_data_float(l->X);
01088           pYdata=cpl_image_get_data_float(l->Y);
01089           pZdata=cpl_image_get_data_float(l->Z);
01090 
01091           pXdata[sinfo_im_xy(l->X,i,j)] = ZERO;
01092           pYdata[sinfo_im_xy(l->Y,i,j)] = ZERO;
01093           pZdata[sinfo_im_xy(l->Z,i,j)] = ZERO;
01094         }
01095     }
01096 #define  FORW   
01097 #ifdef BACK
01098     for( x = idlx - 1;x>=0;x--)
01099     {
01100     for( y = idly - 1;y>=0;y--)
01101         {
01102         for( z = idnp - 1;z>=0;z--)
01103 #endif
01104 #ifdef FORW
01105     for( x = 0; x <  idlx; x++)
01106     {
01107     for( y = 0; y < idly; y++)
01108         {
01109         for( z = 0; z < idnp; z++)
01110 #endif
01111           {
01112         i_img=cpl_imagelist_get(l->id,z);
01113         piddata=cpl_image_get_data_float(i_img);
01114         j_img=cpl_imagelist_get(l->jd,z);
01115         pjddata=cpl_image_get_data_float(j_img);
01116         indexI = piddata [sinfo_cu_xy(l->id,x,y)];
01117         indexJ = pjddata [sinfo_cu_xy(l->jd,x,y)];
01118         if ( indexI > 0.0  && indexI < size  &&
01119              indexJ > 0.0  && indexJ < size )
01120             {
01121             /*sinfo_msg_debug("%d %d %d = %f, %f\n",
01122                       x,y,z,(float)ICube(x,y,z),(float)JCube(x,y,z));*/
01123               pXdata=cpl_image_get_data_float(l->X);
01124               pYdata=cpl_image_get_data_float(l->Y);
01125               pZdata=cpl_image_get_data_float(l->Z);
01126               phXdata=cpl_image_get_data_float(l->hX);
01127 
01128               pXdata[sinfo_im_xy(l->X ,indexI,indexJ)] = x;
01129               phXdata[sinfo_im_xy(l->hX,indexI,indexJ)] = 
01130                          phXdata[sinfo_im_xy(l->hX,indexI,indexJ)] + 1;
01131 
01132             pYdata[sinfo_im_xy(l->Y ,indexI,indexJ)] = y;
01133             /*phXdata[sinfo_im_xy(l->hX,indexI,indexJ)] = 
01134                       phXdata[sinfo_im_xy(l->hX,indexI,indexJ)] + 1;*/
01135 
01136             pZdata[sinfo_im_xy(l->Z ,indexI,indexJ)] = z;
01137             /*phXdata[sinfo_im_xy(l->hX,indexI,indexJ)] = 
01138                       phXdata[sinfo_im_xy(l->hX,indexI,indexJ)] + 1;*/
01139             }
01140         }
01141         }
01142     }
01143 
01144 
01145     sinfo_msg("Filled X Y Z , cX cY cZ 2D frames\n");
01146     return(0);
01147 }
01148 
01154 new_Lookup * 
01155 sinfo_new_lookup( void )
01156 {
01157     new_Lookup *l;
01158     l = (new_Lookup*)cpl_calloc(1, sizeof(new_Lookup));
01159     return (l);
01160 }
01167 void 
01168 sinfo_new_destroy_Lookup( new_Lookup *l ) 
01169 {
01170     if ( l )
01171     cpl_free(l);
01172 }
01180 int 
01181 sinfo_new_change_mask (cpl_image * mask, cpl_image * im)
01182 {
01183     int i ;
01184     int mlx=0;
01185     int mly=0;
01186     int ilx=0;
01187     int ily=0;
01188     float* pidata=NULL;
01189     float* pmdata=NULL;
01190 
01191     if (mask == NULL || im == NULL) return -1 ;
01192     ilx=cpl_image_get_size_x(im);
01193     ily=cpl_image_get_size_y(im);
01194     pidata=cpl_image_get_data_float(im);
01195 
01196     mlx=cpl_image_get_size_x(mask);
01197     mly=cpl_image_get_size_y(mask);
01198     pmdata=cpl_image_get_data_float(mask);
01199 
01200     for ( i = 0 ; i < (int) ilx*ily ; i++ )
01201     {
01202         if (isnan(pidata[i]))
01203         {
01204             pmdata[i] = 0. ;
01205         }
01206     }
01207     return 0 ;
01208 }
01209 
01210 
01227 cpl_image * 
01228 sinfo_new_c_bezier_find_cosmic( cpl_image *im, 
01229                                       cpl_image *mask, 
01230                                       new_Lookup *look, 
01231                                       short rx, 
01232                                       short ry, 
01233                                       short rz,
01234                           short lowerI, 
01235                                       short highI, 
01236                                       short lowerJ, 
01237                                       short highJ, 
01238                                       float factor )
01239 {
01240     
01241     int i,j,count;
01242     cpl_imagelist * sc_im,* drs_sc_mask;
01243     short szx,szy,szz;
01244     float /*ant,*/newValue,old/*,dif,porcentage,distance*/;
01245     double med, stdev;
01246     /*cpl_image *out;*/
01247     short rx_loop, ry_loop, rz_loop; 
01248 
01249 
01250     cpl_image* o_img=NULL;
01251     float* pmdata=NULL;
01252     float* pidata=NULL;
01253 
01254 
01255 
01256     int ilx=0;
01257     int ily=0;
01258 
01259     int mlx=0;
01260     int mly=0;
01261 
01262 
01263 
01264 
01265     mlx=cpl_image_get_size_x(mask);
01266     mly=cpl_image_get_size_y(mask);
01267     pmdata=cpl_image_get_data_float(mask);
01268 
01269     ilx=cpl_image_get_size_x(im);
01270     ily=cpl_image_get_size_y(im);
01271     pidata=cpl_image_get_data_float(im);
01272 
01273     if ( mlx != ilx || mly != ily )
01274     {
01275         sinfo_msg_error(" data & mask images not compatible in size\n") ;
01276         return NULL ;
01277     }
01278 
01279     /* allocate memory for sub cubes*/
01280     szx = (rx * 2 ) + 1;
01281     szy = (ry * 2 ) + 1;
01282     szz = (rz * 2 ) + 1;
01283 
01284     if ( NULL == ( sc_im = cpl_imagelist_new() ) ) 
01285     {
01286         sinfo_msg_error(" could not allocate memory for data subcube\n") ;
01287         return NULL ;
01288     }
01289     for(i=0;i<szz;i++){
01290       o_img=cpl_image_new(szx,szy,CPL_TYPE_FLOAT);
01291       cpl_imagelist_set(sc_im,o_img,i);
01292     }
01293 
01294     if ( NULL == ( drs_sc_mask = cpl_imagelist_new() ) ) 
01295     {
01296         sinfo_msg_error(" could not allocate memory for mask subcube\n") ;
01297         return NULL ;
01298     }
01299     
01300     for(i=0;i<szz;i++){
01301       o_img=cpl_image_new(szx,szy,CPL_TYPE_FLOAT);
01302       cpl_imagelist_set(drs_sc_mask,o_img,i);
01303     }
01304 
01305 
01306 
01307 
01308     count=0;
01309     for ( i = 0 ; i < mlx; i++ )
01310     {
01311         for ( j = 0 ; j < mly ; j++ )
01312         {
01313         if ( i >= lowerI && i < highI && j >= lowerJ && j < highJ )
01314         {
01315 
01316         rx_loop = 1 ; ry_loop = 1 ; rz_loop = 1 ;
01317         newValue = sinfo_new_c_bezier_correct_cosmic( i, j, im, 
01318                                                               mask, sc_im, 
01319                                                               drs_sc_mask, 
01320                                                               look, 
01321                                                               rx_loop, 
01322                                                               ry_loop, 
01323                                                               rz_loop, 
01324                                                               &med, 
01325                                                               &stdev,
01326                                                               factor );
01327             /* if no valid neighboors are found, increase size of 
01328                    sub cube until max radius is reached */ 
01329         while ( newValue == cubeNONEIGHBOR && rx_loop < rx && 
01330                         ry_loop < ry && rz_loop < rz )
01331             {
01332               rx_loop++ ; ry_loop++; rz_loop++;  
01333               /*sinfo_msg_debug("Increasing radius to %d, in %d %d", 
01334                           rx_loop, i, j) ;  */
01335                   newValue = sinfo_new_c_bezier_correct_cosmic( i, j, im, 
01336                                                                     mask, 
01337                                                                     sc_im, 
01338                                                                     drs_sc_mask,
01339                                                                     look, 
01340                                                                     rx_loop, 
01341                                                                     ry_loop, 
01342                                                                     rz_loop, 
01343                                                                     &med, 
01344                                                                     &stdev,
01345                                                                      factor );
01346             }
01347         /*give up on increasing the size*/
01348         if ( isnan(newValue) || newValue == cubeNONEIGHBOR ) 
01349                 /*<= -3.e10 ZERO )*/
01350             continue;
01351         
01352         old = pidata[sinfo_im_xy(im,i,j)];
01353         if ( newValue != old )
01354             {
01355             pidata[sinfo_im_xy(im,i,j)] = newValue;
01356             /*sinfo_msg_debug("[%d,%d]=%f -> %f, med= %f, stdev=%f",
01357                       i,j, old, newValue, med, stdev ); */
01358             count++;
01359             }
01360         }
01361         }
01362     }
01363     
01364     
01365     sinfo_msg_debug("bad pixels count: %d",count);
01366 
01367     
01368     cpl_imagelist_delete(sc_im);
01369     cpl_imagelist_delete(drs_sc_mask);
01370     return im;
01371 }
01372 
01373 
01391 float 
01392 sinfo_new_c_bezier_correct_cosmic(int ipos, 
01393                                   int jpos, 
01394                                   cpl_image * im, 
01395                                   cpl_image * mask, 
01396                       cpl_imagelist * sc_im, 
01397                                   cpl_imagelist * drs_sc_mask, 
01398                                   new_Lookup * look, 
01399                                   short rx, 
01400                                   short ry, 
01401                                   short rz , 
01402                                   double *med , 
01403                       double *stdev, 
01404                                   float factor )
01405 {
01406     short ic, jc, kc, ii, jj, kk/*, sjj, skk*/,is,js,ks;
01407     short i, j, k, indexJ, indexI, lx, ly, lz, szx, szy, szz;
01408     /*float indexIf,indexJf,sp;*/
01409     cpl_image * X, * Y, * Z, * hX;
01410     cpl_imagelist  * id, * jd;
01411     short counter;
01412     double sum;
01413     float* phXdata=NULL;
01414     float* pXdata=NULL;
01415     float* pYdata=NULL;
01416     float* pZdata=NULL;
01417 
01418     float* pimdata=NULL;
01419     float* pscdata=NULL;
01420     float* pdrsdata=NULL;
01421     float* piddata=NULL;
01422     float* pjddata=NULL;
01423     float* pmaskdata=NULL;
01424 
01425 
01426     int idlx=0;
01427     int idly=0;
01428     int idnp=0;
01429 
01430     int drslx=0;
01431     int drsly=0;
01432     int drsnp=0;
01433 
01434 
01435     X  = look -> X;
01436     Y  = look -> Y;
01437     Z  = look -> Z;
01438     hX = look -> hX;
01439     id = look -> id;
01440     jd = look -> jd;
01441 
01442     phXdata=cpl_image_get_data_float(hX);
01443     /*if ( phXdata[sinfo_im_xy( hX, ipos, jpos)] > 1 )
01444     {
01445     sinfo_msg_error("double hit in position [%d,%d]=%f, can not correct", 
01446         ipos,jpos,phXdata[sinfo_im_xy(hX,ipos,jpos)]) ;
01447     return ( -2e10 );
01448     }*/
01449     if ( phXdata[sinfo_im_xy( hX, ipos, jpos)] < 1 )
01450     {
01451     /*sinfo_msg_error("no lookup  in position [%d,%d]=%f, can not correct", 
01452       ipos,jpos,phXdata[sinfo_im_xy(hX,ipos,jpos)]) ;*/
01453     return ( ZERO );
01454     }
01455     
01456     pXdata=cpl_image_get_data_float(X);
01457     pYdata=cpl_image_get_data_float(Y);
01458     pZdata=cpl_image_get_data_float(Z);
01459 
01460     ic = pXdata[sinfo_im_xy( X, ipos, jpos)];
01461     jc = pYdata[sinfo_im_xy( Y, ipos, jpos)];
01462     kc = pZdata[sinfo_im_xy( Z, ipos, jpos)];
01463     /*if ( !(ipos % 16 )  )*/
01464 #ifdef DEBUG    
01465     sinfo_msg_debug("Correcting bad pixel : ipos=%d,jpos=%d, "
01466                     "in Cube -> ic=%d, jc=%d, kc=%d", ipos,jpos, ic, jc, kc );
01467 #endif
01468     /*limit to start not before the beginning of the cube*/
01469     ii = ic - rx; if ( ii < 0 ) ii = 0;
01470     jj = jc - ry; if ( jj < 0 ) jj = 0;
01471     kk = kc - rz; if ( kk < 0 ) kk = 0;
01472 
01473 #ifdef DEBUG
01474     sinfo_msg_debug("Start Point in Cube -> ii=%d,jj=%d,kk=%d", ii, jj, kk );
01475 #endif
01476     
01477     /*limit to end not outside of the cube */
01478     szx = (rx * 2 ) + 1;
01479     szy = (ry * 2 ) + 1;
01480     szz = (rz * 2 ) + 1;
01481 
01482     idlx = cpl_image_get_size_x(cpl_imagelist_get(id,0));
01483     idly = cpl_image_get_size_y(cpl_imagelist_get(id,0));
01484     idnp = cpl_imagelist_get_size(id);
01485 
01486     lx = idlx;
01487     ly = idly;
01488     lz = idnp;
01489     if ( ( ic + rx ) >= idlx )
01490     szx = szx - ( (ic+rx)-(lx-1) );
01491 
01492     if ( ( jc + ry ) >= idly )
01493     szy = szy - ( (jc+ry)-(ly-1) );
01494 
01495     if ( ( kc + rz ) >= idnp )
01496     szz = szz - ( (kc+rz)-(lz-1) );
01497 
01498 #ifdef DEBUG
01499     sinfo_msg_error("Size of subcube : szx=%d,szy=%d,szz=%d\n", szx, szy, szz );
01500     /*fill whole mask with not available*/
01501     drsnp=cpl_imagelist_get_size(drs_sc_mask);
01502     drslx=cpl_image_get_size_x(cpl_imagelist_get(drs_sc_mask,0));
01503     drsly=cpl_image_get_size_y(cpl_imagelist_get(drs_sc_mask,0));
01504 
01505     sinfo_msg_error("Fill Mask subcube of size: %d,%d,%d, with NOINFO", 
01506                      drslx, drsly,  drsnp);
01507 #endif
01508     for( i = 0; i < drslx; i++) {
01509       for( j = 0; j < drsly; j++) {
01510     for( k = 0; k < drsnp; k++) {
01511       pdrsdata=cpl_image_get_data_float(cpl_imagelist_get(drs_sc_mask,k));
01512       pdrsdata[sinfo_cu_xy(drs_sc_mask,i,j)] = cubePT_NOINFO;
01513     }
01514       }
01515     }
01516     pimdata=cpl_image_get_data_float(im);
01517     pmaskdata=cpl_image_get_data_float(mask);
01518     for( i = ii,is=0;  i < ii+szx; i++,is++)
01519       {
01520     for( j = jj,js=0;  j < jj+szy; j++,js++)
01521       {
01522         for( k = kk,ks=0;  k < kk+szz; k++,ks++)
01523           {
01524 #ifdef DEBUG
01525         sinfo_msg_debug("i=%d j=%d k=%d is=%d ij=%d ik=%d",
01526                 i,j,k,is,js,ks);
01527 #endif
01528         piddata=cpl_image_get_data_float(cpl_imagelist_get(id,k));
01529         pjddata=cpl_image_get_data_float(cpl_imagelist_get(id,k));
01530           pdrsdata=cpl_image_get_data_float(cpl_imagelist_get(drs_sc_mask,ks));
01531                 pscdata=cpl_image_get_data_float(cpl_imagelist_get(sc_im,ks));
01532 
01533 
01534         indexI = sinfo_new_nint( piddata[sinfo_cu_xy(id,i,j)] ); 
01535         indexJ = sinfo_new_nint( pjddata[sinfo_cu_xy(jd,i,j)] );
01536         if ( indexJ <= -1 || indexJ>=2048 || indexI == -1 )
01537           {
01538             pdrsdata[sinfo_cu_xy(drs_sc_mask,is,js)] = cubePT_NOINFO;
01539             continue;
01540           }
01541         pscdata[sinfo_cu_xy(sc_im,is,js)]=
01542                          pimdata[sinfo_im_xy(im,indexI,indexJ)];
01543         pdrsdata[sinfo_cu_xy(drs_sc_mask,is,js)]=
01544                          pmaskdata[sinfo_im_xy(mask,indexI,indexJ)];
01545 #ifdef DEBUG
01546         sinfo_msg_debug("Cube i=%d, j=%d, k=%d  ; "
01547                                 " Sub is=%d, js=%d, ks=%d  ; "
01548                                 " Plane I=%d,J=%d ; mask %f ; im %f\n", 
01549                         i, j, k, is, js, ks, indexI, indexJ,  
01550                                 pmaskdata[sinfo_im_xy(mask,indexI,indexJ)], 
01551                                 pimdata[sinfo_im_xy(im,indexI,indexJ)]);
01552 #endif
01553 
01554           }
01555       }
01556       }
01557 
01558     /* ignoring the elements in the slitlet of the tested pixel */
01559 
01560     for( i = 0; i < szx; i++) {
01561       for( k = 0;  k < szz; k++) {
01562     pdrsdata=cpl_image_get_data_float(cpl_imagelist_get(drs_sc_mask,k));
01563     pdrsdata[sinfo_cu_xy(drs_sc_mask,i,ry)] = cubePT_NOINFO;
01564       }
01565     }
01566 /* now calculate mean and stdev in subcube */
01567 
01568     counter = 0;
01569     sum=0;
01570     for( i = 0; i < szx; i++)
01571       {
01572     for( j = 0;  j < szy; j++)
01573       {
01574         for( k = 0;  k < szz; k++)
01575           {
01576         pdrsdata=cpl_image_get_data_float(cpl_imagelist_get(drs_sc_mask,k));
01577         pscdata=cpl_image_get_data_float(cpl_imagelist_get(sc_im,k));
01578         if (pscdata[sinfo_cu_xy(sc_im ,i ,j)] != ZERO && 
01579                     pdrsdata[sinfo_cu_xy(drs_sc_mask,i ,j)] != cubePT_NOINFO)
01580           {
01581             sum = sum + pscdata[sinfo_cu_xy(sc_im ,i ,j)];
01582             counter++;
01583           }
01584           }
01585       }
01586       }
01587 
01588     *med = sum / counter ; 
01589 
01590     counter = 0;
01591     sum=0;
01592     for( i = 0; i < szx; i++)
01593       {
01594     for( j = 0;  j < szy; j++)
01595       {
01596         for( k = 0;  k < szz; k++)
01597           {
01598         pscdata=cpl_image_get_data_float(cpl_imagelist_get(sc_im,k));
01599         pdrsdata=cpl_image_get_data_float(cpl_imagelist_get(drs_sc_mask,k));
01600         if (pscdata[sinfo_cu_xy(sc_im ,i ,j)] != ZERO && 
01601                     pdrsdata[sinfo_cu_xy(drs_sc_mask,i ,j)] != cubePT_NOINFO)
01602           {
01603                sum = sum + (pscdata[sinfo_cu_xy(sc_im ,i ,j)] - *med) * 
01604                                    (pscdata[sinfo_cu_xy(sc_im ,i ,j)] - *med) ;
01605                counter++;
01606           }
01607           }
01608       }
01609       }
01610 
01611     *stdev = sqrt( sum / ( counter - 1 ) );
01612     
01613 
01614     if ( (fabs( pimdata[sinfo_im_xy(im,ipos,jpos)] - *med ) > 
01615           factor * *stdev) || 
01616           isnan(pimdata[sinfo_im_xy(im,ipos,jpos)]) )
01617     {
01618     pdrsdata=cpl_image_get_data_float(cpl_imagelist_get(drs_sc_mask,rz));
01619     pdrsdata[sinfo_cu_xy(drs_sc_mask,rx,ry)] = cubePT_FIND;
01620     return ( sinfo_new_c_bezier_interpol( sc_im, drs_sc_mask ) );
01621     } 
01622     else 
01623         return(pimdata[sinfo_im_xy(im,ipos,jpos)]); 
01624     
01625 
01626 }
01627 

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