bezier.c

00001 /*******************************************************************************
00002 * M.P.E. - SPIFFI project
00003 *
00004 *
00005 *
00006 * who       when      what
00007 * --------  --------  ----------------------------------------------
00008 * rabuter 2004-11-04 Fixed error on Find Cosmic when no valid neighboors where found
00009 *                    in the subcube
00010 * rabuter 10/07/03  created
00011 */
00012 
00013 /************************************************************************
00014 *   NAME
00015 *        bezier.c -
00016 *        procedures to correct for bad pixels using bezier splines 
00017 *
00018 *   SYNOPSIS
00019 *
00020 *   DESCRIPTION
00021 *
00022 *   FILES
00023 *
00024 *   ENVIRONMENT
00025 *
00026 *   RETURN VALUES
00027 *
00028 *   CAUTIONS
00029 *
00030 *   EXAMPLES
00031 *
00032 *   SEE ALSO
00033 *
00034 *   BUGS
00035 *
00036 *------------------------------------------------------------------------
00037 */
00038 
00039 #define POSIX_SOURCE 1
00040 #include "vltPort.h"
00041 
00042 /*
00043  * System Headers
00044  */
00045 
00046 /*
00047  * Local Headers
00048  */
00049 
00050 #include "bezier.h"
00051 
00052 /*----------------------------------------------------------------------------
00053  *                          Function codes
00054  *--------------------------------------------------------------------------*/
00055 
00056 
00057 OneImage * cbezierInterpolateImage( OneImage *im, OneImage *mask, Lookup *look, short rx, short ry, short rz, 
00058                     int max_rad , float   ** slit_edges )
00059 {
00060     
00061     int i,j,count;
00062     OneCube * sc_im,* drs_sc_mask;
00063     OneImage *auxImage;
00064     OneImage *tempMask;
00065     short szx,szy,szz;
00066     short rx_loop, ry_loop, rz_loop; 
00067     /*float ant,new,dif;*/
00068 
00069     if ( mask -> lx != im -> lx || mask -> ly != im -> ly )
00070     {
00071         cpl_msg_error("cbezierInterpolate:"," data & mask images not compatible in size\n") ;
00072         return NullImage ;
00073     }
00074 
00075     /* allocate memory for sub cubes*/
00076     szx = (rx * 2 ) + 1;
00077     szy = (ry * 2 ) + 1;
00078     szz = (rz * 2 ) + 1;
00079 
00080     if ( NullCube == ( sc_im = newCube( szx, szy, szz) ) ) 
00081     {
00082         cpl_msg_error("cbezierInterpolate:"," could not allocate memory for data subcube\n") ;
00083         return NullImage ;
00084     }
00085     if ( NullCube == ( drs_sc_mask = newCube( szx, szy, szz) ) ) 
00086     {
00087         cpl_msg_error("cbezierInterpolate:"," could not allocate memory for mask subcube\n") ;
00088         return NullImage ;
00089     }
00090     if ( NullImage == ( tempMask = new_image(mask -> lx, mask -> ly) ) ) 
00091     {
00092         cpl_msg_error("cbezierInterpolate:"," could not allocate memory for temporary dead pixel mask\n") ;
00093         return NullImage ;
00094     }
00095    
00096     count=0;
00097     for ( i = 0 ; i < mask -> lx; i++ )
00098     {
00099         for ( j = 0 ; j < mask -> ly ; j++ )
00100         {
00101         if ( mask -> data[i2(im,i,j)] ==  cubePT_BADPIXEL )
00102         {
00103           rx_loop = 1 ; ry_loop = 1 ; rz_loop = 1 ;
00104           im -> data[i2(im,i,j)] = cbezierCorrectPixel( i, j, im, mask, sc_im, drs_sc_mask, look, rx_loop, ry_loop, rz_loop );
00105           /* if not enough neighbors found, increase size of sub cube until max radius is reached */ 
00106           while ( im -> data[i2(im,i,j)] == cubeNONEIGHBOR && rx_loop < rx && ry_loop < ry && rz_loop < rz )
00107             {
00108               rx_loop++ ; ry_loop++; rz_loop++;  
00109               /*              printf (" cbezierInterpolateImage: Increasing radius to %d, in %d %d \n", rx_loop, i, j) ; */  
00110               
00111               im -> data[i2(im,i,j)] = cbezierCorrectPixel( i, j, im, mask, sc_im, drs_sc_mask, look, rx_loop, ry_loop, rz_loop );
00112             }
00113           /* If still not enough neighbors, make result NaN = ZERO in spred convention */     
00114           if ( im -> data[i2(im,i,j)] == cubeNONEIGHBOR ) 
00115             { 
00116               im -> data[i2(im,i,j)] = ZERO ;
00117             }
00118         count++;
00119         }
00120         if ( im -> data[i2(im,i,j)] == ZERO )
00121           {
00122         tempMask -> data[i2(tempMask,i,j)] = 0 ;
00123           }
00124          else
00125           {
00126         tempMask -> data[i2(tempMask,i,j)] = 1 ;
00127           }
00128         }
00129     }
00130 
00131     
00132     printf("Replacing NaN\n"); 
00133     auxImage = interpolSourceImage ( im, tempMask, max_rad, slit_edges );
00134     for ( i = 0 ; i < mask -> lx; i++ )
00135     {
00136         for ( j = 0 ; j < mask -> ly ; j++ )
00137         {
00138 
00139         if ( qfits_isnan(im -> data[i2(im,i,j)])) /*<= -2e10ZERO )*/
00140         {
00141           /*        printf("Replacing NaN -> %d %d %f\n", i,j, auxImage -> data[i2(im,i,j)] ); */
00142         im -> data[i2(im,i,j)] = auxImage -> data[i2(im,i,j)];
00143         }
00144         }
00145     }
00146     destroy_image(auxImage);
00147     destroy_cube(sc_im);
00148     destroy_cube(drs_sc_mask);
00149     
00150     printf("bad pixels count: %d\n",count);
00151 
00152     
00153     return im;
00154 }
00155 
00156 OneImage * cbezierFindBad( OneImage *im, OneImage *mask/*, Lookup *look*/, short rx, short ry, short rz,
00157                short lowerI, short highI, short lowerJ, short highJ, float factor )
00158 {
00159     
00160     int i,j,count;
00161     OneCube * sc_im,* drs_sc_mask;
00162     short szx,szy,szz;
00163     float /*ant,*/newValue,old/*,dif,porcentage,distance*/;
00164     double med, stdev;
00165     /*OneImage *out;*/
00166     short rx_loop, ry_loop, rz_loop; 
00167 
00168     if ( mask -> lx != im -> lx || mask -> ly != im -> ly )
00169     {
00170         cpl_msg_error("cbezierInterpolate:"," data & mask images not compatible in size\n") ;
00171         return NullImage ;
00172     }
00173 
00174     /* allocate memory for sub cubes*/
00175     szx = (rx * 2 ) + 1;
00176     szy = (ry * 2 ) + 1;
00177     szz = (rz * 2 ) + 1;
00178 
00179     if ( NullCube == ( sc_im = newCube( szx, szy, szz) ) ) 
00180     {
00181         cpl_msg_error("cbezierInterpolate:"," could not allocate memory for data subcube\n") ;
00182         return NullImage ;
00183     }
00184     if ( NullCube == ( drs_sc_mask = newCube( szx, szy, szz) ) ) 
00185     {
00186         cpl_msg_error("cbezierInterpolate:"," could not allocate memory for mask subcube\n") ;
00187         return NullImage ;
00188     }
00189     
00190     count=0;
00191     for ( i = 0 ; i < mask -> lx; i++ )
00192     {
00193         for ( j = 0 ; j < mask -> ly ; j++ )
00194         {
00195         if ( i >= lowerI && i < highI && j >= lowerJ && j < highJ )
00196         {
00197 
00198         rx_loop = 1 ; ry_loop = 1 ; rz_loop = 1 ;
00199         newValue = cbezierCorrectPixel2D( i, j, im, mask, sc_im, drs_sc_mask/*, look*/, rx_loop, ry_loop, rz_loop, &med, &stdev,factor );
00200             /* if NaN returned, increase size of sub cube until max radius is reached */ 
00201         while ( newValue == ZERO && rx_loop < rx && ry_loop < ry && rz_loop < rz )
00202             {
00203               rx_loop++ ; ry_loop++; rz_loop++;  
00204               /*              printf (" cbezierFindBad: Increasing radius to %d, in %d %d \n", rx_loop, i, j) ;  */
00205                   newValue = cbezierCorrectPixel2D( i, j, im, mask, sc_im, drs_sc_mask/*, look*/, rx_loop, ry_loop, rz_loop, &med, &stdev,factor );
00206             }
00207         if ( qfits_isnan(newValue)) /*<= -3.e10 ZERO )*/
00208             continue;
00209         
00210         old = im -> data[i2(im,i,j)];
00211         if ( newValue != old )
00212             {
00213             im -> data[i2(im,i,j)] = newValue;
00214             /*printf("[%d,%d]=%f -> %f, med= %f, stdev=%f\n",i,j, old, newValue, med, stdev );*/
00215             count++;
00216             }
00217         }
00218         }
00219     }
00220     
00221     
00222     printf("bad pixels count: %d\n",count);
00223 
00224     
00225     destroy_cube(sc_im);
00226     destroy_cube(drs_sc_mask);
00227     return im;
00228 }
00229 
00230 float cbezierCorrectPixel(int ipos, int jpos, OneImage * im, OneImage * mask, 
00231               OneCube * sc_im, OneCube * drs_sc_mask, Lookup * look, short rx, short ry, short rz )
00232 {
00233     short ic, jc, kc, ii, jj, kk/*, sjj, skk*/,is,js,ks;
00234     short i, j, k, indexJ, indexI, lx, ly, lz, szx, szy, szz;
00235     /*float indexIf,indexJf,sp;*/
00236     OneImage * X, * Y, * Z, * hX;
00237     OneCube  * id, * jd;
00238 
00239     X  = look -> X;
00240     Y  = look -> Y;
00241     Z  = look -> Z;
00242     hX = look -> hX;
00243     id = look -> id;
00244     jd = look -> jd;
00245 
00246     /*if ( hX -> data[i2( hX, ipos, jpos)] > 1 )
00247     {
00248     cpl_msg_error("cbezierCorrect:"," double hit in position [%d,%d]=%f, can not correct\n", 
00249         ipos,jpos,hX->data[i2(hX,ipos,jpos)]) ;
00250     return ( -2e10 );
00251     }*/
00252     if ( hX -> data[i2( hX, ipos, jpos)] < 1 )
00253     {
00254     /*cpl_msg_error("cbezierCorrect:"," no lookup  in position [%d,%d]=%f, can not correct\n", 
00255       ipos,jpos,hX->data[i2(hX,ipos,jpos)]) ;*/
00256     return ( ZERO );
00257     }
00258     
00259 
00260     ic = X -> data[i2( X, ipos, jpos)];
00261     jc = Y -> data[i2( Y, ipos, jpos)];
00262     kc = Z -> data[i2( Z, ipos, jpos)];
00263     /*if ( !(ipos % 16 )  )*/
00264 #ifdef DEBUG    
00265     printf("Correcting bad pixel : ipos=%d,jpos=%d, in Cube -> ic=%d, jc=%d, kc=%d\n", ipos,jpos, ic, jc, kc );
00266 #endif
00267     /*limit to start not before the beginning of the cube*/
00268     ii = ic - rx; if ( ii < 0 ) ii = 0;
00269     jj = jc - ry; if ( jj < 0 ) jj = 0;
00270     kk = kc - rz; if ( kk < 0 ) kk = 0;
00271 
00272 #ifdef DEBUG
00273     printf("Start Point in Cube -> ii=%d,jj=%d,kk=%d\n", ii, jj, kk );
00274 #endif
00275     
00276     /*limit to end not outside of the cube */
00277     szx = (rx * 2 ) + 1;
00278     szy = (ry * 2 ) + 1;
00279     szz = (rz * 2 ) + 1;
00280     lx = id -> lx;
00281     ly = id -> ly;
00282     lz = id -> np;
00283     if ( ( ic + rx ) >= id -> lx )
00284     szx = szx - ( (ic+rx)-(lx-1) );
00285 
00286     if ( ( jc + ry ) >= id -> ly )
00287     szy = szy - ( (jc+ry)-(ly-1) );
00288 
00289     if ( ( kc + rz ) >= id -> np )
00290     szz = szz - ( (kc+rz)-(lz-1) );
00291 
00292 #ifdef DEBUG
00293     cpl_msg_error("Size of subcube :"," szx=%d,szy=%d,szz=%d\n", szx, szy, szz );
00294     /*fill whole mask with not available*/
00295     cpl_msg_error("Fill Mask subcube of size:","%d,%d,%d, with NOINFO\n", drs_sc_mask -> lx, drs_sc_mask -> ly,  drs_sc_mask -> np);
00296 #endif
00297     for( i = 0; i < drs_sc_mask -> lx; i++)
00298     for( j = 0; j < drs_sc_mask -> ly; j++)
00299         for( k = 0; k < drs_sc_mask -> np; k++)
00300         drs_sc_mask -> plane[k] -> data[i2(drs_sc_mask,i,j)] = cubePT_NOINFO;
00301 
00302 
00303      for( i = ii,is=0;  i < ii+szx; i++,is++)
00304      {
00305      for( j = jj,js=0;  j < jj+szy; j++,js++)
00306          {
00307          for( k = kk,ks=0;  k < kk+szz; k++,ks++)
00308          {
00309 #ifdef DEBUG
00310          printf("i=%d j=%d k=%d is=%d ij=%d ik=%d\n",i,j,k,is,js,ks);
00311 #endif
00312          indexI = nint( id -> plane[k] -> data[i2(id,i,j)] ); 
00313          indexJ = nint( jd -> plane[k] -> data[i2(jd,i,j)] );
00314          if ( indexJ <= -1 || indexJ>=2048 || indexI == -1 )
00315              {
00316              drs_sc_mask -> plane[ks] -> data[i2(drs_sc_mask,is,js)] = cubePT_NOINFO;
00317              continue;
00318              }
00319          sc_im   -> plane[ks] -> data[i2(sc_im  ,is,js)]  = im   -> data[i2(im,indexI,indexJ)];
00320          drs_sc_mask -> plane[ks] -> data[i2(drs_sc_mask,is,js)]  = mask -> data[i2(mask,indexI,indexJ)];
00321 #ifdef DEBUG
00322          printf("Cube i=%d, j=%d, k=%d  ;  Sub is=%d, js=%d, ks=%d  ;  Plane I=%d,J=%d ; mask %f ; im %f\n", 
00323            i, j, k, is, js, ks, indexI, indexJ,  mask -> data[i2(mask,indexI,indexJ)], im   -> data[i2(im,indexI,indexJ)]);
00324 #endif
00325 
00326          }
00327          }
00328      }
00329 
00330 
00331     /*signal to correct this pixel*/
00332     drs_sc_mask -> plane[rz] -> data[i2(drs_sc_mask,rx,ry)] = cubePT_FIND;
00333     return ( cbezierInterpol( sc_im, drs_sc_mask ) );
00334 }
00335 
00336 float cbezierCorrectPixel2D(int ipos, int jpos, OneImage * im, OneImage * mask, 
00337               OneCube * sc_im, OneCube * drs_sc_mask/*, Lookup * look*/, short rx, short ry, short rz , double *med , 
00338                 double *stdev, float factor )
00339 {
00340     short ic, jc, kc, ii, jj, kk/*, sjj, skk*/,is,js,ks;
00341     short i, j, k, indexJ, indexI, lx, ly, lz, szx, szy, szz;
00342     double sum,aux;
00343     int counter;
00344     float sumarr[100];
00345 
00346     ic = ipos;
00347     jc = 0;
00348     kc = jpos;
00349 
00350 #ifdef DEBUG    
00351     printf("Correcting bad pixel : ipos=%d,jpos=%d, in Cube -> ic=%d, jc=%d, kc=%d\n", ipos,jpos, ic, jc, kc );
00352 #endif
00353     /*limit to start not before the beginning of the cube*/
00354     ii = ic - rx; if ( ii < 0 ) ii = 0;
00355     jj = jc - ry; if ( jj < 0 ) jj = 0;
00356     kk = kc - rz; if ( kk < 0 ) kk = 0;
00357 
00358 #ifdef DEBUG
00359     printf("Start Point in Cube -> ii=%d,jj=%d,kk=%d\n", ii, jj, kk );
00360 #endif
00361     
00362 
00363     /*limit to end not outside of the cube */
00364     szx = (rx * 2 ) + 1;
00365     szy = (ry * 2 ) + 1;
00366     szz = (rz * 2 ) + 1;
00367     lx = im -> lx;
00368     ly = im -> ly;
00369     lz = im -> ly;
00370     if ( ( ic + rx ) >= im -> lx )
00371     szx = szx - ( (ic+rx)-(lx-1) );
00372 
00373     if ( ( jc + ry ) >= im -> ly )
00374     szy = szy - ( (jc+ry)-(ly-1) );
00375 
00376     if ( ( kc + rz ) >= im -> ly )
00377     szz = szz - ( (kc+rz)-(lz-1) );
00378 
00379 #ifdef DEBUG
00380     printf("Size of subcube : szx=%d,szy=%d,szz=%d\n", szx, szy, szz );
00381     /*fill whole mask with not available*/
00382     printf("Fill Mask subcube of size:%d,%d,%d, with NOINFO\n", drs_sc_mask -> lx, drs_sc_mask -> ly,  drs_sc_mask -> np);
00383 #endif
00384     for( i = 0; i < drs_sc_mask -> lx; i++)
00385     for( j = 0; j < drs_sc_mask -> ly; j++)
00386         for( k = 0; k < drs_sc_mask -> np; k++)
00387         drs_sc_mask -> plane[k] -> data[i2(drs_sc_mask,i,j)] = cubePT_NOINFO;
00388 
00389     counter = 0;
00390     sum=0;
00391     memset(sumarr,0x00,sizeof(sumarr));
00392      for( i = ii,is=0;  i < ii+szx; i++,is++)
00393      {
00394      for( j = jj,js=0;  j < jj+szy; j++,js++)
00395          {
00396          for( k = kk,ks=0;  k < kk+szz; k++,ks++)
00397          {
00398 #ifdef DEBUG
00399          printf("i=%d j=%d k=%d is=%d ij=%d ik=%d\n",i,j,k,is,js,ks);
00400 #endif
00401          indexI = i; 
00402          indexJ = k;
00403          if ( qfits_isnan(im -> data[i2(mask,indexI,indexJ)]) )
00404              mask -> data[i2(mask,indexI,indexJ)] = 0;
00405 
00406          if ( mask -> data[i2(mask,indexI,indexJ)] == 1 && ( indexI != ipos || indexJ != jpos) )
00407              {
00408              /*sumarr[counter] = im   -> data[i2(im,indexI,indexJ)];*/
00409              sum = sum + im   -> data[i2(im,indexI,indexJ)];
00410              counter++;
00411              }
00412          sc_im   -> plane[ks] -> data[i2(sc_im  ,is,js)]  = im   -> data[i2(im,indexI,indexJ)];
00413          drs_sc_mask -> plane[ks] -> data[i2(drs_sc_mask,is,js)]  = mask -> data[i2(mask,indexI,indexJ)];
00414 #ifdef DEBUG
00415          printf("Cube i=%d, j=%d, k=%d  ;  Sub is=%d, js=%d, ks=%d  ;  Plane I=%d,J=%d ; mask %f ; im %f\n", 
00416            i, j, k, is, js, ks, indexI, indexJ,  mask -> data[i2(mask,indexI,indexJ)], im   -> data[i2(im,indexI,indexJ)]);
00417 #endif
00418 
00419          }
00420          }
00421      }
00422 
00423 
00424     /*signal to correct this pixel*/
00425     drs_sc_mask -> plane[rz] -> data[i2(drs_sc_mask,rx,ry)] = cubePT_FIND;
00426     if ( counter )
00428     *med = sum/counter;
00429     else
00430     return(im -> data[i2(im,ipos,jpos)]);
00431 
00432     /*printf("%f %f %d\n",sum ,im -> data[i2(im,ipos,jpos)], counter);*/
00433 
00434 
00435 
00436     sum =0;
00437     counter=0;
00438     for( i = ii,is=0;  i < ii+szx; i++,is++)
00439     {
00440     for( j = jj,js=0;  j < jj+szy; j++,js++)
00441         {
00442         for( k = kk,ks=0;  k < kk+szz; k++,ks++)
00443         {
00444 
00445         indexI = i; 
00446         indexJ = k;
00447         if ( drs_sc_mask -> plane[ks] -> data[i2(drs_sc_mask,is,js)] == 1 && ( indexI != ipos || indexJ != jpos) )
00448             {
00449             sum= sum + ( (sc_im -> plane[ks] -> data[i2(drs_sc_mask,is,js)] - *med ) * 
00450                  (sc_im -> plane[ks] -> data[i2(drs_sc_mask,is,js)] - *med ) );
00451             counter++;
00452             }
00453         }
00454         }
00455     }
00456     
00457     aux = sum;
00458     sum   = sum / (counter - 1);
00459     *stdev = sqrt( sum );
00460     
00461     if ( (fabs( im -> data[i2(im,ipos,jpos)] - *med ) > factor * *stdev) || qfits_isnan(im -> data[i2(im,ipos,jpos)]) )
00462     {
00463     /*printf("[%d,%d]: distance to mean = %f, thres =%f sum=%f, stdev=%f, counter=%d, aux= %f\n", 
00464            ipos,jpos, fabs( im -> data[i2(im,ipos,jpos)] - *med), factor * *stdev, sum,*stdev, counter,aux );
00465     mask -> data[i2(mask,ipos,jpos)] = 0;*/
00466     return ( cbezierInterpol( sc_im, drs_sc_mask ) );
00467     }
00468     return(im -> data[i2(im,ipos,jpos)]); 
00469 }
00470 
00471 
00472 
00473 float cbezierInterpol( OneCube * im, OneCube * action )
00474 {
00475     short pos;
00476     unsigned short i,j,k;
00477     XYZW  indata[1000],res,selected;
00478     float step,cumstep,distance,selected_distance;
00479     Dim point;
00480     double munk;
00481 
00482 
00483     memset(indata,0x00,1000*sizeof(XYZW));
00484 
00485     pos=0;
00486     for( i=0; i < im -> lx; i++)
00487     {
00488     for( j=0; j < im -> ly; j++)
00489         {
00490         for( k=0; k < im ->np; k++)
00491         {
00492         if ( action -> plane[k] -> data[i2(action,i,j)] == cubePT_USE )
00493             {
00494 #ifdef DEBUG
00495             printf("Used im[%d,%d,%d]=%lf\n",i,j,k,im -> plane[k] -> data[i2(im,i,j)]);
00496 #endif
00497             indata[pos].x = i;
00498             indata[pos].y = j;
00499             indata[pos].z = k;
00500             indata[pos].w = im -> plane[k] -> data[i2(im,i,j)];
00501             pos++;
00502             }
00503         else
00504             {
00505             if ( action -> plane[k] -> data[i2(action,i,j)] == cubePT_FIND )
00506             {
00507             point.x = i;
00508             point.y = j;
00509             point.z = k;
00510 #ifdef DEBUG
00511             printf("Find for im[%d,%d,%d]=%lf reason:%f\n",i,j,k,im -> plane[k] -> data[i2(im,i,j)],
00512                    action -> plane[k] -> data[i2(action,i,j)]);
00513 #endif
00514             }
00515             else
00516             {
00517 #ifdef DEBUG
00518             printf("Ignored im[%d,%d,%d]=%lf reason:%f\n",i,j,k,im -> plane[k] -> data[i2(im,i,j)],
00519                    action -> plane[k] -> data[i2(action,i,j)]);
00520 #endif
00521             }
00522             }
00523         }   
00524         }   
00525     }
00526     
00527     
00528     if ( pos < 2 )
00529     {
00530 #ifdef DEBUG
00531     printf("subcube contains no valid pixels to use in iterpolation\n");
00532 #endif
00533     /*return( im -> plane[point.z] -> data[i2(im,point.x,point.y)] );*/
00534     return( cubeNONEIGHBOR );
00535 
00536     }
00537     
00538     /*printf("\n");*/
00539     /*printf("\"bezier interpolated data\n");*/
00540 
00541     step    = 0.01;
00542     cumstep = 0.0;
00543     selected_distance=1000;
00544     munk = pow( 1.0-cumstep, (double)pos - 1 );
00545     for ( i = 0 ; ( i < 100 ) && ( munk != 0.0 ); i++ )
00546     {
00547     memset( &res, 0x00, sizeof(XYZW) );
00548     Bezier( indata, pos-1, cumstep, munk, &res);
00549     distance = sqrt( pow( (point.x-res.x), 2)+ pow( (point.y-res.y), 2)+ pow( (point.z-res.z), 2) );
00550     /*printf("%lf %lf %lf %lf %lf\n",res.x,res.y,res.z,res.w,distance);*/
00551     if ( distance < selected_distance )
00552         {
00553         selected_distance = distance;
00554         selected.x = res.x;
00555         selected.y = res.y;
00556         selected.z = res.z;
00557         selected.w = res.w;
00558         }
00559     cumstep = cumstep + step;
00560     munk = pow( 1.0 - cumstep, (double)pos - 1 );
00561     
00562     }
00563     
00564 #ifdef DEBUG
00565     printf("Selected %lf %lf %lf %lf, distance=%lf\n",selected.x,selected.y,selected.z,selected.w,selected_distance);
00566 #endif
00567     im -> plane[point.z] -> data[i2(im,point.x,point.y)] = selected.w;
00568     
00569     return selected.w;
00570 }
00571 
00572 
00573 
00574 int Bezier( XYZW *p, int n, double mu, double munk, XYZW *res )
00575 {
00576    int k, kn, nn, nkn;
00577    double blend, muk;
00578 
00579    muk = 1;
00580    for ( k = 0; k <= n; k++ ) {
00581       nn = n;
00582       kn = k;
00583       nkn = n - k;
00584       blend = muk * munk;
00585       muk *= mu;
00586       munk /= ( 1.0 - mu );
00587       while ( nn >= 1 ) {
00588          blend *= (double)nn;
00589          nn--;
00590          if ( kn > 1 ) {
00591             blend /= (double)kn;
00592             kn--;
00593          }
00594          if ( nkn > 1 ) {
00595             blend /= (double)nkn;
00596             nkn--;
00597          }
00598       }
00599       res -> x += p[k].x * blend;
00600       res -> y += p[k].y * blend;
00601       res -> z += p[k].z * blend;
00602       res -> w += p[k].w * blend;
00603    }
00604    return( 0 );
00605 }
00606 int cCreateXYZ( Lookup *l )
00607 {
00608     OneImage *imX,*imY,*imZ,*imcX;
00609     short i,j,k,indexI,indexJ,x,y,z;
00610     int size;
00611     
00612     size = l -> id -> lx*l -> id -> ly;
00613     /* allocate memory */
00614     if ( NULL == (imX = new_image(size, size)) ) 
00615         {
00616     cpl_msg_error("cCreateXYZ:"," could not allocate memory for X !\n") ;
00617     return -1 ;
00618         }
00619     if ( NULL == (imY = new_image(size, size)) ) 
00620         {
00621     cpl_msg_error("cCreateXYZ:"," could not allocate memory for Y !\n") ;
00622     return -1 ;
00623         }
00624     if ( NULL == (imZ = new_image(size, size)) ) 
00625         {
00626     cpl_msg_error("cCreateXYZ:"," could not allocate memory for Z !\n") ;
00627     return -1 ;
00628         }
00629     if ( NULL == (imcX = new_image(size, size)) ) 
00630         {
00631     cpl_msg_error("cCreateXYZ:"," could not allocate memory for cX !\n") ;
00632     return -1 ;
00633         }
00634 
00635     l -> X  = imX; 
00636     l -> Y  = imY; 
00637     l -> Z  = imZ; 
00638     l -> hX = imcX; 
00639 
00640     /*Round id*/
00641     for( i = 0; i < l -> id -> lx; i++)
00642     {
00643     for( j = 0; j < l -> id -> ly; j++)
00644         {
00645         for( k = 0; k < l -> id -> np; k++)
00646         {
00647         l -> id -> plane[k] -> data[i2(l->id,i,j)] = (float)nint(l -> id -> plane[k] -> data[i2(l->id,i,j)]);
00648         }
00649         }
00650     }
00651 
00652     /*Round jd*/
00653     for( i = 0; i < l -> id -> lx; i++)
00654     {
00655     for( j = 0; j < l -> id -> ly; j++)
00656         {
00657         for( k = 0; k < l -> id -> np; k++)
00658         {
00659         l -> jd -> plane[k] -> data[i2(l->jd,i,j)] = (float)nint(l -> jd -> plane[k] -> data[i2(l->jd,i,j)]);
00660         }
00661         }
00662     }
00663 
00664     /*Fill with -1 X Y Z*/
00665     for( i = 0; i < l -> X -> lx; i++)
00666     {
00667     for( j = 0; j < l -> X -> ly; j++)
00668         {
00669         l -> X -> data[i2(l->X,i,j)] = ZERO;
00670         l -> Y -> data[i2(l->Y,i,j)] = ZERO;
00671         l -> Z -> data[i2(l->Z,i,j)] = ZERO;
00672         }
00673     }
00674 #define  FORW   
00675 #ifdef BACK
00676     for( x = l -> id -> lx - 1;x>=0;x--)
00677     {
00678     for( y = l -> id -> ly - 1;y>=0;y--)
00679         {
00680         for( z = l -> id -> np - 1;z>=0;z--)
00681 #endif
00682 #ifdef FORW
00683     for( x = 0; x <  l -> id -> lx; x++)
00684     {
00685     for( y = 0; y < l -> id -> ly; y++)
00686         {
00687         for( z = 0; z < l -> id -> np; z++)
00688 #endif
00689         {
00690         indexI = l -> id -> plane[z] -> data [i2(l->id,x,y)];
00691         indexJ = l -> jd -> plane[z] -> data [i2(l->jd,x,y)];
00692         if ( indexI > 0.0  && indexI < size  &&
00693              indexJ > 0.0  && indexJ < size )
00694             {
00695             /*printf("%d %d %d = %f, %f\n",x,y,z,(float)ICube(x,y,z),(float)JCube(x,y,z));*/
00696             
00697             l -> X  -> data[i2(l->X ,indexI,indexJ)] = x;
00698             l -> hX -> data[i2(l->hX,indexI,indexJ)] = l -> hX -> data[i2(l->hX,indexI,indexJ)] + 1;
00699 
00700             l -> Y  -> data[i2(l->Y ,indexI,indexJ)] = y;
00701             /*l -> hX -> data[i2(l->hX,indexI,indexJ)] = l -> hX -> data[i2(l->hX,indexI,indexJ)] + 1;*/
00702 
00703             l -> Z  -> data[i2(l->Z ,indexI,indexJ)] = z;
00704             /*l -> hX -> data[i2(l->hX,indexI,indexJ)] = l -> hX -> data[i2(l->hX,indexI,indexJ)] + 1;*/
00705             }
00706         }
00707         }
00708     }
00709 
00710 
00711     printf("Filled X Y Z , cX cY cZ 2D frames\n");
00712     printf("cCreateXYZ: Ok..\n");
00713     return(0);
00714 }
00715 
00716 Lookup * newLookup( void )
00717 {
00718     Lookup *l;
00719     l = (Lookup*)cpl_calloc(1, sizeof(Lookup));
00720     return (l);
00721 }
00722 
00723 void destroyLookup( Lookup *l ) 
00724 {
00725     if ( l )
00726     cpl_free(l);
00727 }
00728 
00729 int changeMask (OneImage * mask, OneImage * im)
00730 {
00731     int i ;
00732     if (mask == NullImage || im == NullImage) return -1 ;
00733 
00734     for ( i = 0 ; i < (int) im -> nbpix ; i++ )
00735     {
00736         if (qfits_isnan(im -> data[i]))
00737         {
00738             mask->data[i] = 0. ;
00739         }
00740     }
00741     return 0 ;
00742 }
00743 
00744 
00745 
00746 OneImage * cbezierFindCosmic( OneImage *im, OneImage *mask, Lookup *look, short rx, short ry, short rz,
00747                short lowerI, short highI, short lowerJ, short highJ, float factor )
00748 {
00749     
00750     int i,j,count;
00751     OneCube * sc_im,* drs_sc_mask;
00752     short szx,szy,szz;
00753     float /*ant,*/newValue,old/*,dif,porcentage,distance*/;
00754     double med, stdev;
00755     /*OneImage *out;*/
00756     short rx_loop, ry_loop, rz_loop; 
00757 
00758     if ( mask -> lx != im -> lx || mask -> ly != im -> ly )
00759     {
00760         cpl_msg_error("cbezierInterpolate:"," data & mask images not compatible in size\n") ;
00761         return NullImage ;
00762     }
00763 
00764     /* allocate memory for sub cubes*/
00765     szx = (rx * 2 ) + 1;
00766     szy = (ry * 2 ) + 1;
00767     szz = (rz * 2 ) + 1;
00768 
00769     if ( NullCube == ( sc_im = newCube( szx, szy, szz) ) ) 
00770     {
00771         cpl_msg_error("cbezierInterpolate:"," could not allocate memory for data subcube\n") ;
00772         return NullImage ;
00773     }
00774     if ( NullCube == ( drs_sc_mask = newCube( szx, szy, szz) ) ) 
00775     {
00776         cpl_msg_error("cbezierInterpolate:"," could not allocate memory for mask subcube\n") ;
00777         return NullImage ;
00778     }
00779     
00780     count=0;
00781     for ( i = 0 ; i < mask -> lx; i++ )
00782     {
00783         for ( j = 0 ; j < mask -> ly ; j++ )
00784         {
00785         if ( i >= lowerI && i < highI && j >= lowerJ && j < highJ )
00786         {
00787 
00788         rx_loop = 1 ; ry_loop = 1 ; rz_loop = 1 ;
00789         newValue = cbezierCorrectCosmic( i, j, im, mask, sc_im, drs_sc_mask, look, rx_loop, ry_loop, rz_loop, &med, &stdev,factor );
00790             /* if no valid neighboors are found, increase size of sub cube until max radius is reached */ 
00791         while ( newValue == cubeNONEIGHBOR && rx_loop < rx && ry_loop < ry && rz_loop < rz )
00792             {
00793               rx_loop++ ; ry_loop++; rz_loop++;  
00794               /*              printf (" cbezierFindCosmic: Increasing radius to %d, in %d %d \n", rx_loop, i, j) ;  */
00795                   newValue = cbezierCorrectCosmic( i, j, im, mask, sc_im, drs_sc_mask, look, rx_loop, ry_loop, rz_loop, &med, &stdev,factor );
00796             }
00797         /*give up on increasing the size*/
00798         if ( qfits_isnan(newValue) || newValue == cubeNONEIGHBOR ) /*<= -3.e10 ZERO )*/
00799             continue;
00800         
00801         old = im -> data[i2(im,i,j)];
00802         if ( newValue != old )
00803             {
00804             im -> data[i2(im,i,j)] = newValue;
00805             /*          printf("[%d,%d]=%f -> %f, med= %f, stdev=%f\n",i,j, old, newValue, med, stdev ); */
00806             count++;
00807             }
00808         }
00809         }
00810     }
00811     
00812     
00813     printf("bad pixels count: %d\n",count);
00814 
00815     
00816     destroy_cube(sc_im);
00817     destroy_cube(drs_sc_mask);
00818     return im;
00819 }
00820 
00821 
00822 
00823 float cbezierCorrectCosmic(int ipos, int jpos, OneImage * im, OneImage * mask, 
00824               OneCube * sc_im, OneCube * drs_sc_mask, Lookup * look, short rx, short ry, short rz , double *med , 
00825                 double *stdev, float factor )
00826 {
00827     short ic, jc, kc, ii, jj, kk/*, sjj, skk*/,is,js,ks;
00828     short i, j, k, indexJ, indexI, lx, ly, lz, szx, szy, szz;
00829     /*float indexIf,indexJf,sp;*/
00830     OneImage * X, * Y, * Z, * hX;
00831     OneCube  * id, * jd;
00832     short counter;
00833     double sum;
00834 
00835     X  = look -> X;
00836     Y  = look -> Y;
00837     Z  = look -> Z;
00838     hX = look -> hX;
00839     id = look -> id;
00840     jd = look -> jd;
00841 
00842     /*if ( hX -> data[i2( hX, ipos, jpos)] > 1 )
00843     {
00844     cpl_msg_error("cbezierCorrect:"," double hit in position [%d,%d]=%f, can not correct\n", 
00845         ipos,jpos,hX->data[i2(hX,ipos,jpos)]) ;
00846     return ( -2e10 );
00847     }*/
00848     if ( hX -> data[i2( hX, ipos, jpos)] < 1 )
00849     {
00850     /*cpl_msg_error("cbezierCorrect:"," no lookup  in position [%d,%d]=%f, can not correct\n", 
00851       ipos,jpos,hX->data[i2(hX,ipos,jpos)]) ;*/
00852     return ( ZERO );
00853     }
00854     
00855     ic = X -> data[i2( X, ipos, jpos)];
00856     jc = Y -> data[i2( Y, ipos, jpos)];
00857     kc = Z -> data[i2( Z, ipos, jpos)];
00858     /*if ( !(ipos % 16 )  )*/
00859 #ifdef DEBUG    
00860     printf("Correcting bad pixel : ipos=%d,jpos=%d, in Cube -> ic=%d, jc=%d, kc=%d\n", ipos,jpos, ic, jc, kc );
00861 #endif
00862     /*limit to start not before the beginning of the cube*/
00863     ii = ic - rx; if ( ii < 0 ) ii = 0;
00864     jj = jc - ry; if ( jj < 0 ) jj = 0;
00865     kk = kc - rz; if ( kk < 0 ) kk = 0;
00866 
00867 #ifdef DEBUG
00868     printf("Start Point in Cube -> ii=%d,jj=%d,kk=%d\n", ii, jj, kk );
00869 #endif
00870     
00871     /*limit to end not outside of the cube */
00872     szx = (rx * 2 ) + 1;
00873     szy = (ry * 2 ) + 1;
00874     szz = (rz * 2 ) + 1;
00875     lx = id -> lx;
00876     ly = id -> ly;
00877     lz = id -> np;
00878     if ( ( ic + rx ) >= id -> lx )
00879     szx = szx - ( (ic+rx)-(lx-1) );
00880 
00881     if ( ( jc + ry ) >= id -> ly )
00882     szy = szy - ( (jc+ry)-(ly-1) );
00883 
00884     if ( ( kc + rz ) >= id -> np )
00885     szz = szz - ( (kc+rz)-(lz-1) );
00886 
00887 #ifdef DEBUG
00888     cpl_msg_error("Size of subcube :"," szx=%d,szy=%d,szz=%d\n", szx, szy, szz );
00889     /*fill whole mask with not available*/
00890     cpl_msg_error("Fill Mask subcube of size:","%d,%d,%d, with NOINFO\n", drs_sc_mask -> lx, drs_sc_mask -> ly,  drs_sc_mask -> np);
00891 #endif
00892     for( i = 0; i < drs_sc_mask -> lx; i++)
00893     for( j = 0; j < drs_sc_mask -> ly; j++)
00894         for( k = 0; k < drs_sc_mask -> np; k++)
00895         drs_sc_mask -> plane[k] -> data[i2(drs_sc_mask,i,j)] = cubePT_NOINFO;
00896 
00897 
00898      for( i = ii,is=0;  i < ii+szx; i++,is++)
00899      {
00900      for( j = jj,js=0;  j < jj+szy; j++,js++)
00901          {
00902          for( k = kk,ks=0;  k < kk+szz; k++,ks++)
00903          {
00904 #ifdef DEBUG
00905          printf("i=%d j=%d k=%d is=%d ij=%d ik=%d\n",i,j,k,is,js,ks);
00906 #endif
00907          indexI = nint( id -> plane[k] -> data[i2(id,i,j)] ); 
00908          indexJ = nint( jd -> plane[k] -> data[i2(jd,i,j)] );
00909          if ( indexJ <= -1 || indexJ>=2048 || indexI == -1 )
00910              {
00911              drs_sc_mask -> plane[ks] -> data[i2(drs_sc_mask,is,js)] = cubePT_NOINFO;
00912              continue;
00913              }
00914          sc_im   -> plane[ks] -> data[i2(sc_im  ,is,js)]  = im   -> data[i2(im,indexI,indexJ)];
00915          drs_sc_mask -> plane[ks] -> data[i2(drs_sc_mask,is,js)]  = mask -> data[i2(mask,indexI,indexJ)];
00916 #ifdef DEBUG
00917          printf("Cube i=%d, j=%d, k=%d  ;  Sub is=%d, js=%d, ks=%d  ;  Plane I=%d,J=%d ; mask %f ; im %f\n", 
00918            i, j, k, is, js, ks, indexI, indexJ,  mask -> data[i2(mask,indexI,indexJ)], im   -> data[i2(im,indexI,indexJ)]);
00919 #endif
00920 
00921          }
00922          }
00923      }
00924 
00925      /* ignoring the elements in the slitlet of the tested pixel */
00926 
00927     for( i = 0; i < szx; i++)
00928      for( k = 0;  k < szz; k++)
00929            drs_sc_mask -> plane[k] -> data[i2(drs_sc_mask,i,ry)] = cubePT_NOINFO;
00930 
00931 /* now calculate mean and stdev in subcube */
00932 
00933     counter = 0;
00934     sum=0;
00935     for( i = 0; i < szx; i++)
00936      {
00937      for( j = 0;  j < szy; j++)
00938          {
00939          for( k = 0;  k < szz; k++)
00940          {
00941            if (sc_im -> plane[k] -> data[i2(sc_im ,i ,j)] != ZERO && drs_sc_mask -> plane[k] -> data[i2(drs_sc_mask,i ,j)] != cubePT_NOINFO)
00942              {
00943                sum = sum + sc_im   -> plane[k] -> data[i2(sc_im ,i ,j)];
00944                counter++;
00945              }
00946          }
00947          }
00948      }
00949 
00950     *med = sum / counter ;  
00951 
00952     counter = 0;
00953     sum=0;
00954     for( i = 0; i < szx; i++)
00955      {
00956      for( j = 0;  j < szy; j++)
00957          {
00958          for( k = 0;  k < szz; k++)
00959          {
00960            if (sc_im -> plane[k] -> data[i2(sc_im ,i ,j)] != ZERO && drs_sc_mask -> plane[k] -> data[i2(drs_sc_mask,i ,j)] != cubePT_NOINFO)
00961              {
00962                sum = sum + (sc_im   -> plane[k] -> data[i2(sc_im ,i ,j)] - *med) * (sc_im   -> plane[k] -> data[i2(sc_im ,i ,j)] - *med) ;
00963                counter++;
00964              }
00965          }
00966          }
00967      }
00968 
00969     *stdev = sqrt( sum / ( counter - 1 ) );
00970     
00971 
00972     if ( (fabs( im -> data[i2(im,ipos,jpos)] - *med ) > factor * *stdev) || qfits_isnan(im -> data[i2(im,ipos,jpos)]) )
00973     {
00974     drs_sc_mask -> plane[rz] -> data[i2(drs_sc_mask,rx,ry)] = cubePT_FIND;
00975     return ( cbezierInterpol( sc_im, drs_sc_mask ) );
00976     } 
00977     else 
00978         return(im -> data[i2(im,ipos,jpos)]); 
00979     
00980 
00981 }
00982 

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