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 #define POSIX_SOURCE 1
00040 #include "vltPort.h"
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 #include "bezier.h"
00051
00052
00053
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
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
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
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
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
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)]))
00140 {
00141
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, 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 newValue,old;
00164 double med, stdev;
00165
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
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, rx_loop, ry_loop, rz_loop, &med, &stdev,factor );
00200
00201 while ( newValue == ZERO && rx_loop < rx && ry_loop < ry && rz_loop < rz )
00202 {
00203 rx_loop++ ; ry_loop++; rz_loop++;
00204
00205 newValue = cbezierCorrectPixel2D( i, j, im, mask, sc_im, drs_sc_mask, rx_loop, ry_loop, rz_loop, &med, &stdev,factor );
00206 }
00207 if ( qfits_isnan(newValue))
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
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,is,js,ks;
00234 short i, j, k, indexJ, indexI, lx, ly, lz, szx, szy, szz;
00235
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
00247
00248
00249
00250
00251
00252 if ( hX -> data[i2( hX, ipos, jpos)] < 1 )
00253 {
00254
00255
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
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
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
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
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
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, short rx, short ry, short rz , double *med ,
00338 double *stdev, float factor )
00339 {
00340 short ic, jc, kc, ii, jj, kk,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
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
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
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
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
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
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
00464
00465
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
00534 return( cubeNONEIGHBOR );
00535
00536 }
00537
00538
00539
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
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
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
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
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
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
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
00702
00703 l -> Z -> data[i2(l->Z ,indexI,indexJ)] = z;
00704
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 newValue,old;
00754 double med, stdev;
00755
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
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
00791 while ( newValue == cubeNONEIGHBOR && rx_loop < rx && ry_loop < ry && rz_loop < rz )
00792 {
00793 rx_loop++ ; ry_loop++; rz_loop++;
00794
00795 newValue = cbezierCorrectCosmic( i, j, im, mask, sc_im, drs_sc_mask, look, rx_loop, ry_loop, rz_loop, &med, &stdev,factor );
00796 }
00797
00798 if ( qfits_isnan(newValue) || newValue == cubeNONEIGHBOR )
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
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,is,js,ks;
00828 short i, j, k, indexJ, indexI, lx, ly, lz, szx, szy, szz;
00829
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
00843
00844
00845
00846
00847
00848 if ( hX -> data[i2( hX, ipos, jpos)] < 1 )
00849 {
00850
00851
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
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
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
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
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
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
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