00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058 #ifdef HAVE_CONFIG_H
00059 # include <config.h>
00060 #endif
00061 #define POSIX_SOURCE 1
00062 #include "sinfo_vltPort.h"
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 #include "sinfo_new_bezier.h"
00073 #include "sinfo_msg.h"
00081
00082
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
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
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
00261
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
00267
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
00274
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)]))
00302 {
00303
00304
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
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 newValue,old;
00337 double med, stdev;
00338
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
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
00407 rx_loop,
00408 ry_loop,
00409 rz_loop,
00410 &med,
00411 &stdev,
00412 factor );
00413
00414
00415 while ( newValue == ZERO && rx_loop < rx &&
00416 ry_loop < ry && rz_loop < rz )
00417 {
00418 rx_loop++ ; ry_loop++; rz_loop++;
00419
00420
00421 newValue = sinfo_new_c_bezier_correct_pixel_2D( i, j, im,
00422 mask,
00423 sc_im,
00424 drs_sc_mask,
00425
00426 rx_loop,
00427 ry_loop,
00428 rz_loop,
00429 &med,
00430 &stdev,
00431 factor );
00432 }
00433 if ( isnan(newValue))
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
00441
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,is,js,ks;
00470 short i, j, k, indexJ, indexI, lx, ly, lz, szx, szy, szz;
00471
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
00509
00510
00511
00512
00513
00514
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
00524
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
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
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
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
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
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
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,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
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
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
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
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
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
00791
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
00829
00830
00831
00832
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
00920
00921
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
00939
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
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
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
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
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
01122
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
01134
01135
01136 pZdata[sinfo_im_xy(l->Z ,indexI,indexJ)] = z;
01137
01138
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 newValue,old;
01245 double med, stdev;
01246
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
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
01328
01329 while ( newValue == cubeNONEIGHBOR && rx_loop < rx &&
01330 ry_loop < ry && rz_loop < rz )
01331 {
01332 rx_loop++ ; ry_loop++; rz_loop++;
01333
01334
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
01348 if ( isnan(newValue) || newValue == cubeNONEIGHBOR )
01349
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
01357
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,is,js,ks;
01407 short i, j, k, indexJ, indexI, lx, ly, lz, szx, szy, szz;
01408
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
01444
01445
01446
01447
01448
01449 if ( phXdata[sinfo_im_xy( hX, ipos, jpos)] < 1 )
01450 {
01451
01452
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
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
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
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
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
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
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