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
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187 #include "vltPort.h"
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197 #include "cube_ops.h"
00198 #include <sys/types.h>
00199 #include <sys/times.h>
00200 #include <math.h>
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215 OneCube *
00216 newCube(
00217 int lx,
00218 int ly,
00219 int n_im
00220 )
00221 {
00222 OneCube * new_cube ;
00223 int i ;
00224
00225
00226
00227 if ((lx<=0)||(lx>MAX_COLUMN_NUMBER)||(ly<=0)||(ly>MAX_LINE_NUMBER) ||
00228 (n_im<=0) || (n_im>MAX_IMAGE_NUMBER))
00229 {
00230 cpl_msg_error("newCube:","error in requested cube size: [%d x %d x %d]", lx, ly, n_im) ;
00231 return NullCube ;
00232 }
00233
00234
00235 new_cube = (OneCube*) malloc(sizeof(OneCube)) ;
00236 new_cube->plane = (OneImage**) calloc(n_im , sizeof(OneImage*)) ;
00237 new_cube->data = (pixelvalue *) calloc(n_im * lx * ly, sizeof(pixelvalue));
00238 new_cube->planes= (OneImage *) calloc(n_im, sizeof(OneImage));
00239
00240 new_cube->lx = lx ;
00241 new_cube->ly = ly ;
00242 new_cube->np = n_im ;
00243 new_cube->nbpix = (ulong32)lx * (ulong32)ly * (ulong32)n_im ;
00244 new_cube->history = (char*)NULL ;
00245 new_cube->n_comments = 0 ;
00246 new_cube->orig_ptype = BPP_DEFAULT ;
00247 new_cube->filename = NULL ;
00248 new_cube->type = 1;
00249
00250 for (i=0 ; i<n_im ; i++)
00251 {
00252 fill_new_image(lx, ly, new_cube->planes+i, new_cube->data+(i*lx*ly));
00253 new_cube->plane[i] = new_cube->planes+i;
00254 }
00255
00256 return new_cube ;
00257 }
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271 OneCube *
00272 cubeOps(
00273 OneCube * cube1,
00274 OneCube * cube2,
00275 int operation)
00276 {
00277
00278 if (cube1==NullCube || cube2==NullCube)
00279 {
00280 cpl_msg_error("cubeOps:","null cubes");
00281 return NullCube ;
00282 }
00283
00284 switch(operation)
00285 {
00286 case '+':
00287 return cubeAdd(cube1, cube2) ;
00288 break ;
00289
00290 case '-':
00291 return cubeSub(cube1, cube2) ;
00292 break ;
00293
00294 case '*':
00295 return cubeMul(cube1, cube2) ;
00296 break ;
00297
00298 case '/':
00299 return cubeDiv(cube1, cube2) ;
00300 break ;
00301
00302 default:
00303 cpl_msg_error("cubeOps:","illegal requested operation: aborting cube arithmetic") ;
00304 return NullCube ;
00305 }
00306 }
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325 OneCube *
00326 cubeConstOps(
00327 OneCube * cube1,
00328 double constant,
00329 int operation)
00330 {
00331 OneCube * tmp_cube ;
00332 OneImage * processed ;
00333 int p ;
00334
00335 if (cube1 == NULL)
00336 {
00337 cpl_msg_error("cubeConstOps:","null cube") ;
00338 return NullCube ;
00339 }
00340
00341 if ((constant == 0.0) && (operation == '/'))
00342 {
00343 cpl_msg_error("cubeConstOps:","division by zero requested in cube/constant operation") ;
00344 return NullCube ;
00345 }
00346
00347 if ( NULL == (tmp_cube = newCube(cube1->lx, cube1->ly, cube1->np)) )
00348 {
00349 cpl_msg_error ("cubeConstOps:","cannot allocate new cube" ) ;
00350 return NullCube ;
00351 }
00352
00353 for (p=0 ; p<tmp_cube->np ; p++)
00354 {
00355 processed = cst_op_image( cube1->plane[p],
00356 constant,
00357 operation) ;
00358 tmp_cube->plane[p] = processed ;
00359 }
00360 return tmp_cube ;
00361 }
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371 OneCube *
00372 cubeSub(
00373 OneCube * c1,
00374 OneCube * c2
00375 )
00376 {
00377 OneCube * c3 ;
00378 ulong32 i ;
00379 int np ;
00380
00381 if ((c1->lx != c2->lx) ||
00382 (c1->ly != c2->ly))
00383 {
00384 cpl_msg_error("cubeSub:","incompatible size: cannot subtract") ;
00385 return NullCube ;
00386 }
00387
00388 if ((c2->np != c1->np) &&
00389 (c2->np != 1))
00390 {
00391 cpl_msg_error("cubeSub:","cannot compute with these number of planes") ;
00392 return NullCube ;
00393 }
00394
00395 if ( NULL == (c3 = newCube( c1->lx, c1->ly, c1->np)) )
00396 {
00397 cpl_msg_error ("cubeSub:","cannot allocate new cube" ) ;
00398 return NullCube ;
00399 }
00400
00401 for (np=0 ; np < c1->np ; np++)
00402 {
00403 for (i=0 ; i< c1->plane[0]->nbpix ; i++)
00404 {
00405 c3->plane[np]->data[i] = c1 -> plane[np] -> data[i] - c2->plane[np]->data[i] ;
00406 }
00407 }
00408
00409 return c3 ;
00410 }
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420 OneCube *
00421 cubeAdd(
00422 OneCube * c1,
00423 OneCube * c2
00424 )
00425 {
00426 OneCube * c3 ;
00427 ulong32 i ;
00428 int np ;
00429
00430 if ((c1->lx != c2->lx) ||
00431 (c1->ly != c2->ly))
00432 {
00433 cpl_msg_error("cubeAdd:","incompatible size: cannot add") ;
00434 return NullCube ;
00435 }
00436 if ((c2->np != c1->np) &&
00437 (c2->np != 1))
00438 {
00439 cpl_msg_error("cubeAdd:","cannot compute with these number of planes") ;
00440 return NullCube ;
00441 }
00442
00443 if (NULL == (c3 = newCube( c1->lx, c1->ly, c1->np)) )
00444 {
00445 cpl_msg_error ("cubeAdd:","cannot allocate new cube") ;
00446 return NullCube ;
00447 }
00448
00449 for (np=0 ; np < c1->np ; np++)
00450 {
00451 for (i=0 ; i< c1->plane[0]->nbpix ; i++)
00452 {
00453 c3 -> plane[np]->data[i] = c1->plane[np]->data[i] + c2->plane[np]->data[i] ;
00454 }
00455 }
00456
00457 return c3 ;
00458 }
00459
00460
00461
00462
00463
00464
00465
00466
00467 OneCube *
00468 cubeMul(
00469 OneCube * c1,
00470 OneCube * c2
00471 )
00472 {
00473 OneCube *c3 ;
00474 ulong32 i ;
00475 int np ;
00476
00477 if ((c1->lx != c2->lx) ||
00478 (c1->ly != c2->ly))
00479 {
00480 cpl_msg_error("cubeMul:","incompatible size: cannot multiply") ;
00481 return NullCube ;
00482 }
00483
00484 if ((c2->np != c1->np) &&
00485 (c2->np != 1))
00486 {
00487 cpl_msg_error("cubeMul:","cannot compute with these number of planes") ;
00488 return NullCube ;
00489 }
00490
00491 if ( NULL == (c3 = newCube( c1->lx, c1->ly, c1->np)) )
00492 {
00493 cpl_msg_error ( "cubeMul:","cannot allocate new cube" ) ;
00494 return NullCube ;
00495 }
00496
00497 for (np=0 ; np < c1->np ; np++)
00498 {
00499 for (i=0 ; i< c1->plane[0]->nbpix ; i++)
00500 {
00501 c3->plane[np]->data[i] = c1->plane[np]->data[i] * c2->plane[np]->data[i] ;
00502 }
00503 }
00504
00505 return c3 ;
00506 }
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516 OneCube *
00517 cubeDiv(
00518 OneCube * c1,
00519 OneCube * c2
00520 )
00521 {
00522 OneCube * c3 ;
00523 ulong32 i ;
00524 int np ;
00525
00526 if ((c1->lx != c2->lx) ||
00527 (c1->ly != c2->ly))
00528 {
00529 cpl_msg_error("cubeDiv:","incompatible size: cannot divide") ;
00530 return NullCube ;
00531 }
00532
00533 if ((c2->np != c1->np) &&
00534 (c2->np != 1))
00535 {
00536 cpl_msg_error("cubeDiv:","cannot compute with these number of planes") ;
00537 return NullCube ;
00538 }
00539
00540 if (NULL == (c3 = newCube( c1->lx, c1->ly, c1->np)) )
00541 {
00542 cpl_msg_error ("cubeDiv:","cannot allocate a new cube") ;
00543 return NullCube ;
00544 }
00545
00546 for (np=0 ; np < c1->np ; np++)
00547 {
00548 for (i=0 ; i< c1->plane[0]->nbpix ; i++)
00549 {
00550 if (fabs((double)c2->plane[np]->data[i]) < 1e-10)
00551 {
00552 c3->plane[np]->data[i] = 0.0 ;
00553 }
00554 else
00555 {
00556 c3->plane[np]->data[i] = c1->plane[np]->data[i] / c2->plane[np]->data[i] ;
00557 }
00558 }
00559 }
00560
00561 return c3 ;
00562 }
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573 OneCube * addImageToCube(OneCube * cu, OneImage * im)
00574 {
00575 OneCube * cube ;
00576 int i ;
00577
00578 if (cu==NULL || im==NULL)
00579 {
00580 cpl_msg_error ("addImageToCube:","null cube or null image") ;
00581 return NullCube ;
00582 }
00583
00584 if ((cu->lx != im->lx) || (cu->ly != im->ly))
00585 {
00586 cpl_msg_error("addImageToCube:","incompatible size: cannot add image to cube") ;
00587 return NullCube ;
00588 }
00589
00590 cube = copy_cube (cu) ;
00591
00592 for (i=0 ; i<cu->np ; i++)
00593 {
00594 add_img_local(cube->plane[i], im) ;
00595 }
00596
00597 return cube ;
00598 }
00599
00600
00601
00602
00603
00604
00605
00606
00607 OneCube * subImageFromCube (OneCube * cu, OneImage * im)
00608 {
00609 OneCube * cube ;
00610 int i ;
00611
00612 if (cu==NULL || im==NULL)
00613 {
00614 cpl_msg_error ("subImageFromCube:","null cube or null image") ;
00615 return NullCube ;
00616 }
00617
00618 if ((cu->lx != im->lx) || (cu->ly != im->ly))
00619 {
00620
00621 cpl_msg_error("subImageFromCube:","incompatible size: cannot subtract image from cube") ;
00622 return NullCube ;
00623 }
00624
00625 cube = copy_cube (cu) ;
00626
00627 for (i=0 ; i<cu->np ; i++)
00628 {
00629 sub_img_local(cube->plane[i], im) ;
00630 }
00631 return cube ;
00632 }
00633
00634
00635
00636
00637
00638
00639
00640
00641 OneCube * mulImageToCube(OneCube * cu, OneImage * im)
00642 {
00643 OneCube * cube ;
00644 int i ;
00645
00646 if (cu==NULL || im==NULL)
00647 {
00648 cpl_msg_error("mulImageToCube:","null cube or null image") ;
00649 return NullCube ;
00650 }
00651
00652 if ((cu->lx != im->lx) || (cu->ly != im->ly))
00653 {
00654 cpl_msg_error("mulImageToCube:","incompatible size: cannot multiply image by cube") ;
00655 return NullCube ;
00656 }
00657
00658 cube = copy_cube (cu) ;
00659
00660 for (i=0 ; i<cu->np ; i++)
00661 {
00662 mul_img_local(cube->plane[i], im) ;
00663 }
00664
00665 return cube ;
00666 }
00667
00668
00669
00670
00671
00672
00673
00674
00675 OneCube * divCubeByImage(OneCube * cu, OneImage * im)
00676 {
00677 OneCube * cube ;
00678 int i ;
00679
00680 if (cu==NULL || im==NULL)
00681 {
00682 cpl_msg_error ("divCubeByImage:","null cube or null image") ;
00683 return NullCube ;
00684 }
00685 if ((cu->lx != im->lx) || (cu->ly != im->ly))
00686 {
00687 cpl_msg_error("divCubeByImage:","incompatible size: cannot divide cube by image") ;
00688 return NullCube ;
00689 }
00690
00691 cube = copy_cube (cu) ;
00692
00693 for (i=0 ; i<cu->np ; i++)
00694 {
00695 div_img_local(cube->plane[i], im) ;
00696 }
00697
00698 return cube ;
00699 }
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710 OneCube * addSpectrumToCube(OneCube *cu, Vector *spec)
00711 {
00712 OneCube * cube ;
00713 int i ,j ;
00714
00715 if (cu == NULL || spec == NULL)
00716 {
00717 cpl_msg_error ("addSpectrumToCube:","null cube or null spectrum") ;
00718 return NullCube ;
00719 }
00720
00721 if ( cu -> np != spec -> n_elements )
00722 {
00723 cpl_msg_error ("addSpectrumToCube:","cube length and spectrum length are not compatible") ;
00724 return NullCube ;
00725 }
00726
00727 if ( NULL == (cube = newCube (cu->lx, cu->ly, cu->np)) )
00728 {
00729 cpl_msg_error ("addSpectrumToCube:","cannot allocate new cube" ) ;
00730 return NullCube ;
00731 }
00732
00733 for ( i = 0; i < cu -> np; i++)
00734 {
00735 for ( j = 0; j < (int) cu -> plane[i] -> nbpix; j++)
00736 {
00737 cube -> plane[i] -> data[j] = cu -> plane[i] -> data[j] + spec -> data[i] ;
00738 }
00739 }
00740
00741 return cube ;
00742 }
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752 OneCube * subSpectrumFromCube(OneCube *cu, Vector *spec)
00753 {
00754 OneCube * cube ;
00755 int i ,j ;
00756
00757 if (cu == NULL || spec == NULL)
00758 {
00759 cpl_msg_error ("subSpectrumFromCube:","null cube or null spectrum") ;
00760 return NullCube ;
00761 }
00762
00763 if ( cu -> np != spec -> n_elements )
00764 {
00765 cpl_msg_error ("subSpectrumFromCube:","cube length and spectrum length are not compatible") ;
00766 return NullCube ;
00767 }
00768
00769 if ( NULL == (cube = newCube (cu->lx, cu->ly, cu->np)) )
00770 {
00771 cpl_msg_error ("subSpectrumFromCube:","cannot allocate new cube" ) ;
00772 return NullCube ;
00773 }
00774
00775 for ( i = 0; i < cu -> np; i++)
00776 {
00777 for ( j = 0; j < (int) cu -> plane[i] -> nbpix; j++)
00778 {
00779 if ( qfits_isnan(cu -> plane[i] -> data[j]) || qfits_isnan(spec -> data[i]) )
00780 {
00781 cube->plane[i]->data[j] = ZERO ;
00782 }
00783 else
00784 {
00785 cube -> plane[i] -> data[j] = cu -> plane[i] -> data[j] - spec -> data[i] ;
00786 }
00787 }
00788 }
00789
00790 return cube ;
00791 }
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802 OneCube * mulSpectrumToCube(OneCube *cu, Vector *spec)
00803 {
00804 OneCube * cube ;
00805 int i ,j ;
00806
00807 if (cu == NULL || spec == NULL)
00808 {
00809 cpl_msg_error ("mulSpectrumToCube:","null cube or null spectrum") ;
00810 return NullCube ;
00811 }
00812
00813 if ( cu -> np != spec -> n_elements )
00814 {
00815 cpl_msg_error ("mulSpectrumToCube:","cube length and spectrum length are not compatible") ;
00816 return NullCube ;
00817 }
00818
00819 if ( NULL == (cube = newCube (cu->lx, cu->ly, cu->np)) )
00820 {
00821 cpl_msg_error ("mulSpectrumToCube:","cannot allocate new cube" ) ;
00822 return NullCube ;
00823 }
00824
00825 for ( i = 0; i < cu -> np; i++)
00826 {
00827 for ( j = 0; j < (int) cu -> plane[i] -> nbpix; j++)
00828 {
00829 if ( qfits_isnan(cu->plane[i]->data[j]) || qfits_isnan(spec->data[i]) )
00830 {
00831 cube -> plane[i] -> data[j] = ZERO ;
00832 }
00833 else
00834 {
00835 cube -> plane[i] -> data[j] = cu -> plane[i] -> data[j] * spec -> data[i] ;
00836 }
00837 }
00838 }
00839
00840 return cube ;
00841 }
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852 OneCube * divCubeBySpectrum(OneCube *cu, Vector *spec)
00853 {
00854 OneCube * cube ;
00855 float help ;
00856 int i ,j ;
00857
00858 if (cu == NULL || spec == NULL)
00859 {
00860 cpl_msg_error ("divCubeBySpectrum:","null cube or null spectrum") ;
00861 return NullCube ;
00862 }
00863
00864 if ( cu -> np != spec -> n_elements )
00865 {
00866 cpl_msg_error ("divCubeBySpectrum:","cube length and spectrum length are not compatible") ;
00867 return NullCube ;
00868 }
00869
00870 if (NULL == (cube = newCube (cu->lx, cu->ly, cu->np)) )
00871 {
00872 cpl_msg_error ("divCubeBySpectrum:","cannot allocate new cube") ;
00873 return NullCube ;
00874 }
00875
00876 for ( i = 0; i < cu -> np; i++)
00877 {
00878 for ( j = 0; j < (int) cu -> plane[i] -> nbpix; j++)
00879 {
00880 if (!qfits_isnan(spec->data[i]) && spec->data[i] != 0.)
00881 {
00882 help = 1/spec->data[i] ;
00883 if ( help > THRESH )
00884 {
00885 help = 1. ;
00886 }
00887 }
00888 else
00889 {
00890 help = ZERO ;
00891 }
00892
00893 if ( qfits_isnan(help) || qfits_isnan(cu -> plane[i] -> data[j]) )
00894 {
00895 cube -> plane[i] -> data[j] = ZERO ;
00896 }
00897 else
00898 {
00899 cube -> plane[i] -> data[j] = cu->plane[i]->data[j] * help ;
00900 }
00901 }
00902 }
00903 return cube ;
00904 }
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920 Vector * cleanMeanOfSpectra(OneCube * cube,
00921 int llx,
00922 int lly,
00923 int urx,
00924 int ury,
00925 double lo_reject,
00926 double hi_reject)
00927 {
00928 Vector * mean ;
00929 pixelvalue *rectangle ;
00930 int i, j, k, l, m ;
00931 int recsize, lo_n, hi_n, nv ;
00932
00933 if ( cube == NULL || cube -> np < 1 )
00934 {
00935 cpl_msg_error ("cleanMeanOfSpectra:","no cube to take the mean of his spectra") ;
00936 return NullVector ;
00937 }
00938
00939 if ((llx<1) || (llx>cube->lx) ||
00940 (urx<1) || (urx>cube->lx) ||
00941 (lly<1) || (lly>cube->ly) ||
00942 (ury<1) || (ury>cube->ly) ||
00943 (llx>=urx) || (lly>=ury))
00944 {
00945 cpl_msg_error("cleanMeanOfSpectra:","invalid rectangle coordinates:") ;
00946 cpl_msg_error("cleanMeanOfSpectra:","lower left is [%d %d] upper right is [%d %d]", llx, lly, urx, ury) ;
00947 return NullVector ;
00948 }
00949
00950 if ((lo_reject + hi_reject) > 0.9)
00951 {
00952 cpl_msg_error("cleanMeanOfSpectra:","illegal rejection thresholds: [%f] and [%f]", lo_reject, hi_reject) ;
00953 cpl_msg_error("cleanMeanOfSpectra:","threshold sum should not be over 0.9 aborting average") ;
00954 return NullVector ;
00955 }
00956
00957
00958 llx -- ;
00959 lly -- ;
00960 urx -- ;
00961 ury -- ;
00962
00963 recsize = (urx - llx + 1) * (ury - lly + 1) ;
00964
00965 lo_n = (int) (recsize * lo_reject + 0.5) ;
00966 hi_n = (int) (recsize * hi_reject + 0.5) ;
00967
00968 if (lo_n + hi_n >= recsize)
00969 {
00970 cpl_msg_error ("cleanMeanOfSpectra:","everything would be rejected") ;
00971 return NullVector;
00972 }
00973
00974
00975 if (NULL == (mean = newVector (cube -> np)) )
00976 {
00977 cpl_msg_error ("cleanMeanOfSpectra:","cannot allocate a new vector") ;
00978 return NullVector ;
00979 }
00980
00981
00982
00983
00984
00985 for ( i = 0 ; i < cube->np ; i++ )
00986 {
00987 m = 0 ;
00988 rectangle = (pixelvalue *) cpl_calloc (recsize, sizeof (pixelvalue*));
00989
00990 for ( j = lly ; j <= ury ; j++ )
00991 {
00992 for ( k = llx ; k <= urx ; k++ )
00993 {
00994 rectangle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
00995 m ++ ;
00996 }
00997 }
00998 pixel_qsort (rectangle, recsize) ;
00999
01000 nv = 0 ;
01001 for ( l = lo_n ; l < (recsize - hi_n) ; l++ )
01002 {
01003 mean -> data[i] += rectangle[l] ;
01004 nv ++;
01005 }
01006 mean -> data[i] /= nv ;
01007
01008 cpl_free ( rectangle ) ;
01009 }
01010 return mean ;
01011 }
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024 OneImage * medianCube(OneCube * cube)
01025 {
01026 OneImage * im ;
01027 pixelvalue * buffer ;
01028 int i, j, k, nz ;
01029
01030 if ( cube == NULL )
01031 {
01032 cpl_msg_error ("medianCube:","null cube") ;
01033 return NullImage ;
01034 }
01035
01036
01037 if (NULL == (im = new_image (cube -> lx, cube -> ly )) )
01038 {
01039 cpl_msg_error ("medianCube:","cannot allocate new image") ;
01040 return NullImage ;
01041 }
01042
01043
01044
01045
01046
01047
01048 for ( i = 0 ; i < (int) cube -> plane[0] -> nbpix ; i++ )
01049 {
01050 buffer = (pixelvalue *) cpl_calloc (cube -> np, sizeof (pixelvalue *));
01051 k = 0 ;
01052 for ( j = 0 ; j < cube -> np ; j ++ )
01053 {
01054 if ( !qfits_isnan(cube -> plane[j] -> data[i]) )
01055 {
01056 buffer[k] = cube -> plane[j] -> data[i] ;
01057 k ++ ;
01058 }
01059 }
01060 nz = k ;
01061
01062
01063 if ( nz > 2 )
01064 {
01065 im-> data[i] = median ( buffer, nz ) ;
01066 }
01067 else if (nz == 2)
01068 {
01069 im -> data[i] = (buffer[0] + buffer[1]) / 2. ;
01070 }
01071 else if (nz == 1)
01072 {
01073 im -> data[i] = buffer[0] ;
01074 }
01075 else if (nz == 0)
01076 {
01077 im -> data[i] = ZERO ;
01078 }
01079
01080 cpl_free ( buffer ) ;
01081 }
01082
01083 return im ;
01084 }
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097 OneImage * averageCubeToImage(OneCube * cube)
01098 {
01099 OneImage * im ;
01100 int i, j, nz ;
01101
01102 if ( cube == NULL )
01103 {
01104 cpl_msg_error ("averageCubeToImage:","null cube") ;
01105 return NullImage ;
01106 }
01107
01108
01109 if (NULL == (im = new_image (cube -> lx, cube -> ly )) )
01110 {
01111 cpl_msg_error ("averageCubeToImage:","cannot allocate new image") ;
01112 return NullImage ;
01113 }
01114
01115
01116
01117
01118
01119
01120 for ( i = 0 ; i < (int) cube -> plane[0] -> nbpix ; i++ )
01121 {
01122 nz = 0 ;
01123 for ( j = 0 ; j < cube -> np ; j ++ )
01124 {
01125 if ( !qfits_isnan(cube -> plane[j] -> data[i]) )
01126 {
01127 nz ++ ;
01128 im -> data[i] += cube -> plane[j] -> data[i] ;
01129 }
01130 }
01131
01132
01133 if ( nz >= 1 )
01134 {
01135 im -> data[i] /= nz ;
01136 }
01137 else if (nz == 0)
01138 {
01139 im -> data[i] = ZERO ;
01140 }
01141 }
01142
01143 return im ;
01144 }
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156 OneImage * sumCubeToImage(OneCube * cube)
01157 {
01158 OneImage * im ;
01159 int i, j, nz ;
01160
01161 if ( cube == NULL )
01162 {
01163 cpl_msg_error ("sumCubeToImage:","null cube") ;
01164 return NullImage ;
01165 }
01166
01167
01168 if (NULL == (im = new_image (cube -> lx, cube -> ly )) )
01169 {
01170 cpl_msg_error ("sumCubeToImage:","cannot allocate new image") ;
01171 return NullImage ;
01172 }
01173
01174
01175
01176
01177
01178
01179 for ( i = 0 ; i < (int) cube -> plane[0] -> nbpix ; i++ )
01180 {
01181 nz = 0 ;
01182 for ( j = 0 ; j < cube -> np ; j ++ )
01183 {
01184 if ( !qfits_isnan(cube -> plane[j] -> data[i]) )
01185 {
01186 nz++ ;
01187 im -> data[i] += cube -> plane[j] -> data[i] ;
01188 }
01189 }
01190
01191
01192 if (nz == 0)
01193 {
01194 im -> data[i] = ZERO ;
01195 }
01196 }
01197
01198 return im ;
01199 }
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221 OneImage * averageCubeToImageBetweenWaves (OneCube * cube,
01222 float dispersion,
01223 float centralWave,
01224 float initialLambda,
01225 float finalLambda)
01226 {
01227 OneImage * im ;
01228 int firstPlane ;
01229 int lastPlane ;
01230 int i, j, nz ;
01231 float minWave ;
01232
01233 if ( cube == NULL )
01234 {
01235 cpl_msg_error ("averageCubeToImageBetweenWaves:","null cube") ;
01236 return NullImage ;
01237 }
01238
01239 minWave = centralWave - (float) (cube->np / 2)*dispersion ;
01240
01241 if ( dispersion <= 0. || minWave <= 0. )
01242 {
01243 cpl_msg_error ("averageCubeToImageBetweenWaves:","wrong dispersion or minimum wavelength given") ;
01244 return NullImage ;
01245 }
01246
01247 if ( initialLambda < minWave ||
01248 (initialLambda >= minWave + dispersion * cube->np) )
01249 {
01250 cpl_msg_error ("averageCubeToImageBetweenWaves:","wrong initial wavelength given") ;
01251 return NullImage ;
01252 }
01253
01254 if ( finalLambda <= minWave ||
01255 (finalLambda > minWave + dispersion * cube->np) )
01256 {
01257 cpl_msg_error ("averageCubeToImageBetweenWaves:","wrong final wavelength given") ;
01258 return NullImage ;
01259 }
01260
01261
01262 if (NULL == (im = new_image (cube -> lx, cube -> ly )) )
01263 {
01264 cpl_msg_error ("averageCubeToImageBetweenWaves:","cannot allocate new image") ;
01265 return NullImage ;
01266 }
01267
01268
01269 firstPlane = nint ((double) ((initialLambda - minWave) / dispersion)) ;
01270 lastPlane = nint ((double) ((finalLambda - minWave) / dispersion)) ;
01271
01272 if ( firstPlane < 0 || firstPlane >= cube -> np ||
01273 lastPlane < 0 || lastPlane > cube -> np )
01274 {
01275 cpl_msg_error ("averageCubeToImageBetweenWaves:","wrong values given!") ;
01276 return NullImage ;
01277 }
01278
01279
01280
01281
01282
01283
01284 for ( i = 0 ; i < (int) cube -> plane[0] -> nbpix ; i++ )
01285 {
01286 nz = 0 ;
01287 for ( j = firstPlane ; j <= lastPlane ; j ++ )
01288 {
01289 if ( !qfits_isnan(cube -> plane[j] -> data[i]) )
01290 {
01291 nz ++ ;
01292 im -> data[i] += cube -> plane[j] -> data[i] ;
01293 }
01294 }
01295
01296
01297 if ( nz >= 1 )
01298 {
01299 im -> data[i] /= nz ;
01300 }
01301 else if (nz == 0)
01302 {
01303 im -> data[i] = ZERO ;
01304 }
01305 }
01306
01307 return im ;
01308 }
01309
01310
01311
01312
01313
01314
01315
01316
01317 OneImage * extractImageFromCube(OneCube * cube, int plane_index)
01318 {
01319 if ( cube == NULL )
01320 {
01321 cpl_msg_error ("extractImageFromCube:","null cube") ;
01322 return NullImage ;
01323 }
01324
01325 if ( plane_index < 0 || plane_index >= cube -> np )
01326 {
01327 cpl_msg_error ("extractImageFromCube:","wrong plane index for image to be extracted") ;
01328 return NullImage ;
01329 }
01330
01331 return cube -> plane[plane_index] ;
01332 }
01333
01334
01335
01336
01337
01338
01339
01340
01341 Vector * extractSpectrumFromCube( OneCube * cube, int x_pos, int y_pos )
01342 {
01343 Vector * returnedSpectrum ;
01344 int i ;
01345
01346 if ( cube == NullCube )
01347 {
01348 cpl_msg_error ("extractSpectrumFromCube:","no cube given!") ;
01349 return NullVector ;
01350 }
01351
01352 if ( x_pos < 0 || x_pos >= cube -> lx )
01353 {
01354 cpl_msg_error ("extractSpectrumFromCube:","wrong x-positon of spectrum given!") ;
01355 return NullVector ;
01356 }
01357
01358 if ( y_pos < 0 || y_pos >= cube -> ly )
01359 {
01360 cpl_msg_error ("extractSpectrumFromCube:","wrong y-positon of spectrum given!") ;
01361 return NullVector ;
01362 }
01363
01364
01365 if ( NULL == (returnedSpectrum = newVector ( cube -> np )) )
01366 {
01367 cpl_msg_error ("extractSpectrumFromCube:","cannot allocate new spectrum!") ;
01368 return NullVector ;
01369 }
01370
01371 for ( i = 0 ; i < cube -> np ; i++ )
01372 {
01373 returnedSpectrum -> data[i] = cube -> plane[i] -> data[x_pos + cube->lx*y_pos] ;
01374 }
01375
01376 return returnedSpectrum ;
01377 }
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414 OneCube * combineJitteredCubes ( OneCube ** cubes,
01415 OneCube * mergedCube,
01416 int n_cubes,
01417 float * cumoffsetx,
01418 float * cumoffsety,
01419 float * exptimes,
01420 char * kernel_type )
01421 {
01422
01423 int i ;
01424 int x, y, z ;
01425 int llx[n_cubes], lly[n_cubes] ;
01426 int llx0, lly0 ;
01427 int posx, posy ;
01428 float weight ;
01429 OneCube * mask ;
01430 float sub_offsetx[n_cubes] ;
01431 float sub_offsety[n_cubes] ;
01432 double * kernel ;
01433
01434 OneCube * tmpcubes[n_cubes] ;
01435 pixelvalue * tmpspace;
01436
01437 if ( cubes == NULL )
01438 {
01439 cpl_msg_error ("combineJitteredCubes:","no cube list given!") ;
01440 return NullCube ;
01441 }
01442 if ( n_cubes <= 0 )
01443 {
01444 cpl_msg_error ("combineJitteredCubes:","wrong number of data cubes in list!") ;
01445 return NullCube ;
01446 }
01447 if ( cumoffsetx == NULL || cumoffsety == NULL )
01448 {
01449 cpl_msg_error ("combineJitteredCubes:","no cumoffsetx/y given!") ;
01450 return NullCube ;
01451 }
01452 if ( exptimes == NULL )
01453 {
01454 cpl_msg_error ("combineJitteredCubes:","no exposure time array given!") ;
01455 return NullCube ;
01456 }
01457
01458 if ( NullCube == (mask = newCube(mergedCube->lx, mergedCube->ly, cubes[0]->np)) )
01459 {
01460 cpl_msg_error ("combineJitteredCubes:","could not allocate cube!") ;
01461 return NullCube ;
01462 }
01463
01464 for ( i = 0 ; i < n_cubes ; i++ )
01465 {
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475 }
01476
01477
01478
01479
01480
01481
01482
01483
01484 llx0 = mergedCube->lx/2 - cubes[0]->lx/2 ;
01485 lly0 = mergedCube->ly/2 - cubes[0]->ly/2 ;
01486
01487
01488
01489
01490
01491
01492 for ( i = 0 ; i < n_cubes ; i++ )
01493 {
01494 llx[i] = llx0 - nint(cumoffsetx[i]) ;
01495 sub_offsetx[i] = (float)nint(cumoffsetx[i]) - cumoffsetx[i] ;
01496 lly[i] = lly0 - nint(cumoffsety[i]) ;
01497 sub_offsety[i] = (float)nint(cumoffsety[i]) - cumoffsety[i] ;
01498 }
01499
01500
01501
01502
01503
01504
01505
01506 if ( NULL == (kernel = generate_interpolation_kernel( kernel_type )) )
01507 {
01508 cpl_msg_warning ("combineJitteredCubes:","could not generate desired interpolation kernel or no kernel_typ was given, the default kernel is used now!") ;
01509 }
01510
01511
01512
01513
01514 for ( i = 0 ; i < n_cubes ; i++ )
01515 {
01516 tmpspace = cpl_calloc(cubes[i]->lx, cubes[i]->ly*sizeof(pixelvalue)) ;
01517 tmpcubes[i] = new_cube(cubes[i]->lx,cubes[i]->ly,cubes[i]->np);
01518
01519 for ( z = 0 ; z < cubes[0]->np ; z++ )
01520 {
01521
01522
01523
01524
01525
01526
01527 tmpcubes[i]->plane[z]=shift_image(cubes[i]->plane[z], sub_offsetx[i], sub_offsety[i], kernel);
01528
01529 if (tmpcubes[i]->plane[z]==NULL)
01530 {
01531 cpl_msg_error ("combineJitteredCubes:","could not shift image plane no %d in cube no %d!", z, i) ;
01532 destroy_cube(mergedCube) ;
01533 destroy_cube(mask) ;
01534 cpl_free(kernel) ;
01535 return NullCube ;
01536 }
01537 }
01538 cpl_free(tmpspace);
01539
01540 }
01541
01542
01543
01544
01545
01546
01547
01548
01549 for ( i = 0 ; i < n_cubes ; i++ )
01550 {
01551
01552 for ( y = 0 ; y < mergedCube->ly ; y++ )
01553 {
01554 for ( x = 0 ; x < mergedCube->lx ; x++ )
01555 {
01556
01557 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
01558 x >= llx[i] && x < llx[i]+cubes[i]->lx )
01559 {
01560 posx = x - llx[i] ;
01561 posy = y - lly[i] ;
01562 for ( z = 0 ; z < mergedCube->np ; z++ )
01563 {
01564 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
01565 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.)
01566 {
01567 mask->plane[z]->data[x+y*mask->lx] += exptimes[i] ;
01568 }
01569 }
01570 }
01571 }
01572 }
01573 }
01574
01575
01576
01577
01578
01579
01580
01581
01582 for ( i = 0 ; i < n_cubes ; i++ )
01583 {
01584
01585
01586 for ( y = 0 ; y < mergedCube->ly ; y++ )
01587 {
01588
01589 for ( x = 0 ; x < mergedCube->lx ; x++ )
01590 {
01591
01592
01593 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
01594 x >= llx[i] && x < llx[i]+cubes[i]->lx )
01595 {
01596
01597 posx = x - llx[i] ;
01598 posy = y - lly[i] ;
01599 for ( z = 0 ; z < mergedCube->np ; z++ )
01600 {
01601
01602 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]))
01603 {
01604 if (mask->plane[z]->data[x+y*mask->lx] != 0.)
01605 {
01606
01607 weight = exptimes[0] / mask->plane[z]->data[x+y*mask->lx] ;
01608 }
01609 else
01610 {
01611 weight = 0. ;
01612 }
01613 mergedCube->plane[z]->data[x+y*mergedCube->lx] +=
01614 weight*tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] ;
01615 }
01616 }
01617 }
01618 }
01619 }
01620 }
01621
01622
01623
01624
01625
01626
01627 free(kernel) ;
01628 for( i = 0 ; i < n_cubes ; i++ )
01629 {
01630 destroy_cube (tmpcubes[i]) ;
01631 }
01632 return mask ;
01633 }
01634
01635
01636
01637
01638
01639 int combineJitteredCubesRange ( OneCube ** cubes,
01640 OneCube * mergedCube,
01641 OneCube * mask,
01642 int n_cubes,
01643 float * cumoffsetx,
01644 float * cumoffsety,
01645 float * exptimes,
01646 char * kernel_type,
01647 const int z_min, const int z_max )
01648 {
01649
01650 int i,m ;
01651 int x, y, z ;
01652 int llx[n_cubes], lly[n_cubes] ;
01653 int llx0, lly0 ;
01654 int posx, posy ;
01655 float weight ;
01656 float sub_offsetx[n_cubes] ;
01657 float sub_offsety[n_cubes] ;
01658 double * kernel ;
01659 const int z_siz=z_max-z_min;
01660
01661 OneCube * tmpcubes[n_cubes] ;
01662 pixelvalue * tmpspace;
01663
01664 if ( cubes == NULL )
01665 {
01666 cpl_msg_error ("combineJitteredCubes:","no cube list given!") ;
01667 return -1 ;
01668 }
01669 if ( n_cubes <= 0 )
01670 {
01671 cpl_msg_error ("combineJitteredCubes:","wrong number of data cubes in list!") ;
01672 return -1 ;
01673 }
01674 if ( cumoffsetx == NULL || cumoffsety == NULL )
01675 {
01676 cpl_msg_error ("combineJitteredCubes:","no cumoffsetx/y given!") ;
01677 return -1;
01678 }
01679 if ( exptimes == NULL )
01680 {
01681 cpl_msg_error ("combineJitteredCubes:","no exposure time array given!") ;
01682 return -1 ;
01683 }
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698 llx0 = mergedCube->lx/2 - cubes[0]->lx/2 ;
01699 lly0 = mergedCube->ly/2 - cubes[0]->ly/2 ;
01700
01701
01702
01703
01704
01705
01706 for ( i = 0 ; i < n_cubes ; i++ )
01707 {
01708 llx[i] = llx0 - nint(cumoffsetx[i]) ;
01709 sub_offsetx[i] = (float)nint(cumoffsetx[i]) - cumoffsetx[i] ;
01710 lly[i] = lly0 - nint(cumoffsety[i]) ;
01711 sub_offsety[i] = (float)nint(cumoffsety[i]) - cumoffsety[i] ;
01712 }
01713
01714
01715
01716
01717
01718
01719 if ( NULL == (kernel = generate_interpolation_kernel( kernel_type )) )
01720 {
01721 cpl_msg_warning ("combineJitteredCubes:","could not generate desired interpolation kernel or no kernel_typ was given, the default kernel is used now!") ;
01722 }
01723
01724
01725
01726 for ( i = 0 ; i < n_cubes ; i++ )
01727 {
01728 tmpspace = cpl_calloc(cubes[i]->lx, cubes[i]->ly*sizeof(pixelvalue)) ;
01729 tmpcubes[i] = new_cube(cubes[i]->lx,cubes[i]->ly,z_siz);
01730 for ( z = z_min, m=0 ; z< z_max ; z++, m++ )
01731 {
01732
01733
01734
01735 tmpcubes[i]->plane[m]=shift_image(cubes[i]->plane[z], sub_offsetx[i], sub_offsety[i], kernel);
01736 if (tmpcubes[i]->plane[m]==NULL)
01737 {
01738 cpl_msg_error ("combineJitteredCubes:","could not shift image plane no %d in cube no %d!", z, i) ;
01739 destroy_cube(mergedCube) ;
01740 destroy_cube(mask) ;
01741 cpl_free(kernel) ;
01742 return -1 ;
01743 }
01744 }
01745 cpl_free(tmpspace);
01746 }
01747
01748
01749
01750
01751
01752
01753
01754
01755 for ( i = 0 ; i < n_cubes ; i++ )
01756 {
01757
01758 for ( y = 0 ; y < mergedCube->ly ; y++ )
01759 {
01760 for ( x = 0 ; x < mergedCube->lx ; x++ )
01761 {
01762
01763 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
01764 x >= llx[i] && x < llx[i]+cubes[i]->lx )
01765 {
01766 posx = x - llx[i] ;
01767 posy = y - lly[i] ;
01768 for ( z = z_min, m=0 ; z < z_max ; z++, m++ )
01769 {
01770 if (!qfits_isnan(tmpcubes[i]->plane[m]->data[posx+posy*cubes[i]->lx]) &&
01771 tmpcubes[i]->plane[m]->data[posx+posy*cubes[i]->lx] != 0.)
01772 {
01773 mask->plane[z]->data[x+y*mask->lx] += exptimes[i] ;
01774 }
01775 }
01776 }
01777 }
01778 }
01779 }
01780
01781
01782
01783 for ( i = 0 ; i < n_cubes ; i++ )
01784 {
01785
01786
01787 for ( y = 0 ; y < mergedCube->ly ; y++ )
01788 {
01789
01790 for ( x = 0 ; x < mergedCube->lx ; x++ )
01791 {
01792
01793
01794 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
01795 x >= llx[i] && x < llx[i]+cubes[i]->lx )
01796 {
01797
01798 posx = x - llx[i] ;
01799 posy = y - lly[i] ;
01800 for ( z = z_min, m = 0 ; z < z_max ; z++, m++ )
01801 {
01802
01803 if (!qfits_isnan(tmpcubes[i]->plane[m]->data[posx+posy*cubes[i]->lx]))
01804 {
01805 if (mask->plane[z]->data[x+y*mask->lx] != 0.)
01806 {
01807
01808 weight = exptimes[0] / mask->plane[z]->data[x+y*mask->lx] ;
01809 }
01810 else
01811 {
01812 weight = 0. ;
01813 }
01814 mergedCube->plane[z]->data[x+y*mergedCube->lx] +=
01815 weight*tmpcubes[i]->plane[m]->data[posx+posy*cubes[i]->lx] ;
01816 }
01817 }
01818 }
01819 }
01820 }
01821 }
01822
01823
01824
01825 free(kernel) ;
01826 for( i = 0 ; i < n_cubes ; i++ )
01827 {
01828 destroy_cube (tmpcubes[i]) ;
01829 }
01830 return 0 ;
01831 }
01832
01833 int combineJitteredCubesSkyRange ( OneCube ** cubes,
01834 OneCube * mergedCube,
01835 OneCube * mask,
01836 OneCube * mergedSky,
01837 OneCube * msk_cube,
01838 OneCube * med_cube,
01839 OneCube * avg_cube,
01840 OneCube * sig_cube,
01841 OneCube * nc_cube,
01842 int n_cubes,
01843 float * cumoffsetx,
01844 float * cumoffsety,
01845 float * exptimes,
01846 char * kernel_type,
01847 const int z_min,
01848 const int z_max )
01849 {
01850
01851 const char* _id="combineJitteredCubesSky";
01852 int i ;
01853 int x, y, z ;
01854 int k=2;
01855 int llx[n_cubes], lly[n_cubes] ;
01856 int llx0, lly0 ;
01857 int posx, posy ;
01858 float weight ;
01859
01860 double med=0;
01861 double sig=0;
01862 double avg=0;
01863 int nc=0;
01864 float sub_offsetx[n_cubes] ;
01865 float sub_offsety[n_cubes] ;
01866 double * kernel ;
01867 cpl_vector* val=NULL;
01868 cpl_vector* msk=NULL;
01869 int nclip=0;
01870 int ks=0;
01871 int ovr=0;
01872
01873 double msk_sum=0;
01874 double val_msk_sum=0;
01875
01876 OneCube * tmpcubes[n_cubes] ;
01877 pixelvalue * tmpspace;
01878 const int z_siz=z_max-z_min;
01879 int m=0;
01880
01881
01882 if ( cubes == NULL )
01883 {
01884 cpl_msg_error ("combineJitteredCubes:","no cube list given!") ;
01885 return -1 ;
01886 }
01887 if ( n_cubes <= 0 )
01888 {
01889 cpl_msg_error ("combineJitteredCubes:","wrong number of data cubes in list!") ;
01890 return -1 ;
01891 }
01892 if ( cumoffsetx == NULL || cumoffsety == NULL )
01893 {
01894 cpl_msg_error ("combineJitteredCubes:","no cumoffsetx/y given!") ;
01895 return -1 ;
01896 }
01897 if ( exptimes == NULL )
01898 {
01899 cpl_msg_error ("combineJitteredCubes:","no exposure time array given!") ;
01900 return -1 ;
01901 }
01902 if (z_siz <= 0 )
01903 {
01904 cpl_msg_error ("combineJitteredCubes:","z_max <= z_min given!") ;
01905 return -1 ;
01906 }
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916 llx0 = mergedCube->lx/2 - cubes[0]->lx/2 ;
01917 lly0 = mergedCube->ly/2 - cubes[0]->ly/2 ;
01918
01919
01920
01921
01922
01923
01924 for ( i = 0 ; i < n_cubes ; i++ )
01925 {
01926 llx[i] = llx0 - nint(cumoffsetx[i]) ;
01927 sub_offsetx[i] = (float)nint(cumoffsetx[i]) - cumoffsetx[i] ;
01928 lly[i] = lly0 - nint(cumoffsety[i]) ;
01929 sub_offsety[i] = (float)nint(cumoffsety[i]) - cumoffsety[i] ;
01930 }
01931
01932
01933
01934
01935
01936
01937
01938 if ( NULL == (kernel = generate_interpolation_kernel( kernel_type )) )
01939 {
01940 cpl_msg_warning ("combineJitteredCubes:","could not generate desired interpolation kernel or no kernel_typ was given, the default kernel is used now!") ;
01941 }
01942
01943
01944
01945 cpl_msg_info(_id,"drs5");
01946 for ( i = 0 ; i < n_cubes ; i++ )
01947 {
01948 tmpspace = cpl_calloc(cubes[i]->lx, cubes[i]->ly*sizeof(pixelvalue)) ;
01949 tmpcubes[i]=new_cube(cubes[i]->lx,cubes[i]->ly,cubes[i]->np);
01950 for ( z = z_min ; z < z_max ; z++ )
01951 {
01952
01953
01954
01955
01956
01957 tmpcubes[i]->plane[z]=shift_image(cubes[i]->plane[z], sub_offsetx[i],
01958 sub_offsety[i], kernel);
01959 if (tmpcubes[i]->plane[z]==NULL)
01960 {
01961 cpl_msg_error ("combineJitteredCubes:","could not shift image plane no %d in cube no %d!", z, i) ;
01962 destroy_cube(mergedCube) ;
01963 destroy_cube(mask) ;
01964 cpl_free(kernel) ;
01965 return -1 ;
01966 }
01967 }
01968 cpl_free(tmpspace);
01969
01970 }
01971 cpl_msg_info(_id,"drs6");
01972
01973
01974
01975
01976
01977
01978
01979 for ( i = 0 ; i < n_cubes ; i++ )
01980 {
01981
01982 for ( y = 0 ; y < mergedCube->ly ; y++ )
01983 {
01984 for ( x = 0 ; x < mergedCube->lx ; x++ )
01985 {
01986
01987 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
01988 x >= llx[i] && x < llx[i]+cubes[i]->lx )
01989 {
01990 posx = x - llx[i] ;
01991 posy = y - lly[i] ;
01992 for ( z = z_min ; z < z_max ; z++ )
01993 {
01994 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
01995 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.)
01996 {
01997 mask->plane[z]->data[x+y*mask->lx] += exptimes[i] ;
01998 }
01999 }
02000 }
02001 }
02002 }
02003 }
02004
02005 cpl_msg_info(_id,"drs6a");
02006
02007 for ( i = 0 ; i < n_cubes ; i++ )
02008 {
02009
02010 for ( y = 0 ; y < mergedCube->ly ; y++ )
02011 {
02012 for ( x = 0 ; x < mergedCube->lx ; x++ )
02013 {
02014
02015 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02016 x >= llx[i] && x < llx[i]+cubes[i]->lx )
02017 {
02018 posx = x - llx[i] ;
02019 posy = y - lly[i] ;
02020 m=0;
02021 for ( z = z_min ; z < z_max ; z++ )
02022 {
02023 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02024 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.)
02025 {
02026 msk_cube->plane[m]->data[x+y*msk_cube->lx] += 1 ;
02027 }
02028 m++;
02029 }
02030 }
02031 }
02032 }
02033 }
02034 cpl_msg_info(_id,"drs7");
02035 m=0;
02036 for ( z = z_min ; z < z_max ; z++ ) {
02037
02038 for ( y = 0 ; y < mergedCube->ly ; y++ ) {
02039 for ( x = 0 ; x < mergedCube->lx ; x++ ) {
02040 avg=0;
02041 nc=0;
02042
02043 for ( i = 0 ; i < n_cubes ; i++ ) {
02044 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02045 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02046 posx = x - llx[i] ;
02047 posy = y - lly[i] ;
02048 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02049 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.) {
02050 nc++;
02051 }
02052 }
02053 }
02054 if( nc > 0 ) {
02055 msk=cpl_vector_new(n_cubes);
02056 for (i=0;i<n_cubes;i++) {
02057 cpl_vector_set(msk,i,1);
02058 }
02059
02060
02061 nclip=0;
02062
02063
02064 nc_cube->plane[m]->data[x+y*nc_cube->lx]=nc;
02065 for (ks=0;ks<nc;ks++) {
02066 sig=0;
02067 med=0;
02068 ovr=0;
02069 if(nc-nclip >0) {
02070 val=cpl_vector_new(nc-nclip);
02071 }
02072
02073 for ( i = 0 ; i < n_cubes ; i++ ) {
02074 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02075 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02076 posx = x - llx[i] ;
02077 posy = y - lly[i] ;
02078 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02079 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0. &&
02080 (cpl_vector_get(msk,i) != 0)) {
02081 cpl_vector_set(val,ovr,
02082 (double)tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]);
02083 ovr++;
02084 }
02085 }
02086 }
02087
02088
02089 if(ovr>0) {
02090 avg=cpl_vector_get_mean(val);
02091 med=cpl_vector_get_median(val);
02092 if(ovr>1) {
02093 sig=cpl_vector_get_stdev(val);
02094 } else {
02095 sig=0;
02096 }
02097 med_cube->plane[m]->data[x+y*med_cube->lx]=med;
02098 avg_cube->plane[m]->data[x+y*avg_cube->lx]=avg;
02099 sig_cube->plane[m]->data[x+y*sig_cube->lx]=sig;
02100 cpl_vector_delete(val);
02101 }
02102
02103 for ( i = 0 ; i < n_cubes ; i++ ) {
02104
02105 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02106 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02107 posx = x - llx[i] ;
02108 posy = y - lly[i] ;
02109 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02110 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0. &&
02111 (cpl_vector_get(msk,i) != 0)) {
02112 if(abs((tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]-med))> k*sig) {
02113 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]=0;
02114 msk_cube->plane[m]->data[x+y*msk_cube->lx]-=1;
02115 mask->plane[z]->data[x+y*mask->lx] -= exptimes[i] ;
02116 cpl_vector_set(msk,i,0);
02117 nclip++;
02118 }
02119 }
02120 }
02121 }
02122 }
02123
02124 msk_sum=0;
02125 val_msk_sum=0;
02126 for ( i = 0 ; i < n_cubes ; i++ ) {
02127
02128 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02129 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02130 posx = x - llx[i] ;
02131 posy = y - lly[i] ;
02132 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02133 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0. &&
02134 (cpl_vector_get(msk,i) != 0)) {
02135
02136
02137
02138
02139
02140 msk_sum+=mask->plane[z]->data[x+y*mask->lx];
02141 val_msk_sum+=tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]*
02142 mask->plane[z]->data[x+y*mask->lx];
02143 }
02144 }
02145 }
02146 mergedSky->plane[m]->data[x+y*mergedSky->lx]=val_msk_sum/msk_sum;
02147 mergedCube->plane[z]->data[x+y*mergedCube->lx]=val_msk_sum/msk_sum;
02148 cpl_vector_delete(msk);
02149 }
02150 }
02151 }
02152 m++;
02153 }
02154
02155 cpl_msg_info(_id,"drs8");
02156
02157
02158
02159
02160
02161
02162
02163
02164
02165
02166
02167
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185
02186
02187
02188
02189
02190
02191
02192
02193
02194
02195
02196
02197
02198 cpl_msg_info(_id,"drs9");
02199
02200
02201
02202 cpl_msg_info(_id,"drs9");
02203
02204 cpl_msg_info(_id,"drs9");
02205 free(kernel) ;
02206 for( i = 0 ; i < n_cubes ; i++ )
02207 {
02208 destroy_cube (tmpcubes[i]) ;
02209 }
02210 cpl_msg_info(_id,"drs10");
02211 return 0 ;
02212 }
02213
02214
02215
02216 int combineJitteredCubesThomasRange(OneCube ** cubes,
02217 OneCube * mergedCube,
02218 OneCube * mask,
02219 int n_cubes,
02220 float * cumoffsetx,
02221 float * cumoffsety,
02222 float * exptimes,
02223 char * kernel_type,
02224 const int z_min,
02225 const int z_max,
02226 const double kappa )
02227 {
02228
02229 int i ;
02230 int x, y, z ;
02231 int llx[n_cubes], lly[n_cubes] ;
02232 int llx0, lly0 ;
02233 int posx, posy ;
02234 float weight ;
02235
02236 double med=0;
02237 double sig=0;
02238 double avg=0;
02239 int nc=0;
02240 float sub_offsetx[n_cubes] ;
02241 float sub_offsety[n_cubes] ;
02242 double * kernel ;
02243 cpl_vector* val=NULL;
02244 cpl_vector* msk=NULL;
02245 int nclip=0;
02246 int ks=0;
02247 int ovr=0;
02248
02249 double msk_sum=0;
02250 double val_msk_sum=0;
02251
02252 OneCube * tmpcubes[n_cubes] ;
02253 pixelvalue * tmpspace;
02254 const int z_siz=z_max-z_min;
02255 int m=0;
02256
02257 if ( cubes == NULL )
02258 {
02259 cpl_msg_error ("combineJitteredCubes:","no cube list given!") ;
02260 return -1 ;
02261 }
02262 if ( n_cubes <= 0 )
02263 {
02264 cpl_msg_error ("combineJitteredCubes:","wrong number of data cubes in list!") ;
02265 return -1 ;
02266 }
02267 if ( cumoffsetx == NULL || cumoffsety == NULL )
02268 {
02269 cpl_msg_error ("combineJitteredCubes:","no cumoffsetx/y given!") ;
02270 return -1 ;
02271 }
02272 if ( exptimes == NULL )
02273 {
02274 cpl_msg_error ("combineJitteredCubes:","no exposure time array given!") ;
02275 return -1 ;
02276 }
02277 if (z_siz <= 0 )
02278 {
02279 cpl_msg_error ("combineJitteredCubes:","z_max <= z_min given!") ;
02280 return -1 ;
02281 }
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291 llx0 = mergedCube->lx/2 - cubes[0]->lx/2 ;
02292 lly0 = mergedCube->ly/2 - cubes[0]->ly/2 ;
02293
02294
02295
02296
02297
02298
02299 for ( i = 0 ; i < n_cubes ; i++ )
02300 {
02301 llx[i] = llx0 - nint(cumoffsetx[i]) ;
02302 sub_offsetx[i] = (float)nint(cumoffsetx[i]) - cumoffsetx[i] ;
02303 lly[i] = lly0 - nint(cumoffsety[i]) ;
02304 sub_offsety[i] = (float)nint(cumoffsety[i]) - cumoffsety[i] ;
02305 }
02306
02307
02308
02309
02310
02311
02312
02313 if ( NULL == (kernel = generate_interpolation_kernel( kernel_type )) )
02314 {
02315 cpl_msg_warning ("combineJitteredCubes:","could not generate desired interpolation kernel or no kernel_typ was given, the default kernel is used now!") ;
02316 }
02317
02318
02319
02320 for ( i = 0 ; i < n_cubes ; i++ )
02321 {
02322 tmpspace = cpl_calloc(cubes[i]->lx, cubes[i]->ly*sizeof(pixelvalue)) ;
02323 tmpcubes[i]=new_cube(cubes[i]->lx,cubes[i]->ly,cubes[i]->np);
02324 for ( z = z_min ; z < z_max ; z++ )
02325 {
02326
02327
02328
02329
02330
02331 tmpcubes[i]->plane[z]=shift_image(cubes[i]->plane[z], sub_offsetx[i],
02332 sub_offsety[i], kernel);
02333 if (tmpcubes[i]->plane[z]==NULL)
02334 {
02335 cpl_msg_error ("combineJitteredCubes:","could not shift image plane no %d in cube no %d!", z, i) ;
02336 destroy_cube(mergedCube) ;
02337 destroy_cube(mask) ;
02338 cpl_free(kernel) ;
02339 return -1 ;
02340 }
02341 }
02342 cpl_free(tmpspace);
02343
02344 }
02345
02346
02347
02348
02349
02350
02351
02352 for ( i = 0 ; i < n_cubes ; i++ )
02353 {
02354
02355 for ( y = 0 ; y < mergedCube->ly ; y++ )
02356 {
02357 for ( x = 0 ; x < mergedCube->lx ; x++ )
02358 {
02359
02360 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02361 x >= llx[i] && x < llx[i]+cubes[i]->lx )
02362 {
02363 posx = x - llx[i] ;
02364 posy = y - lly[i] ;
02365 for ( z = z_min ; z < z_max ; z++ )
02366 {
02367 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02368 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.)
02369 {
02370 mask->plane[z]->data[x+y*mask->lx] += exptimes[i] ;
02371 }
02372 }
02373 }
02374 }
02375 }
02376 }
02377
02378 m=0;
02379 for ( z = z_min ; z < z_max ; z++ ) {
02380
02381 for ( y = 0 ; y < mergedCube->ly ; y++ ) {
02382 for ( x = 0 ; x < mergedCube->lx ; x++ ) {
02383 avg=0;
02384 nc=0;
02385
02386 for ( i = 0 ; i < n_cubes ; i++ ) {
02387 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02388 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02389 posx = x - llx[i] ;
02390 posy = y - lly[i] ;
02391 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02392 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.) {
02393 nc++;
02394 }
02395 }
02396 }
02397 if( nc > 0 ) {
02398 msk=cpl_vector_new(n_cubes);
02399 for (i=0;i<n_cubes;i++) {
02400 cpl_vector_set(msk,i,1);
02401 }
02402
02403
02404 nclip=0;
02405
02406
02407 for (ks=0;ks<nc;ks++) {
02408 sig=0;
02409 med=0;
02410 ovr=0;
02411 if(nc-nclip >0) {
02412 val=cpl_vector_new(nc-nclip);
02413 }
02414
02415 for ( i = 0 ; i < n_cubes ; i++ ) {
02416 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02417 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02418 posx = x - llx[i] ;
02419 posy = y - lly[i] ;
02420 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02421 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0. &&
02422 (cpl_vector_get(msk,i) != 0)) {
02423 cpl_vector_set(val,ovr,
02424 (double)tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]);
02425 ovr++;
02426 }
02427 }
02428 }
02429
02430
02431 if(ovr>0) {
02432 avg=cpl_vector_get_mean(val);
02433 med=cpl_vector_get_median(val);
02434 if(ovr>1) {
02435 sig=cpl_vector_get_stdev(val);
02436 } else {
02437 sig=0;
02438 }
02439 cpl_vector_delete(val);
02440 }
02441
02442 for ( i = 0 ; i < n_cubes ; i++ ) {
02443
02444 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02445 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02446 posx = x - llx[i] ;
02447 posy = y - lly[i] ;
02448 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02449 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0. &&
02450 (cpl_vector_get(msk,i) != 0)) {
02451 if(abs((tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]-med))> kappa*sig) {
02452 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]=0;
02453 mask->plane[z]->data[x+y*mask->lx] -= exptimes[i] ;
02454 cpl_vector_set(msk,i,0);
02455 nclip++;
02456 }
02457 }
02458 }
02459 }
02460 }
02461
02462 msk_sum=0;
02463 val_msk_sum=0;
02464 for ( i = 0 ; i < n_cubes ; i++ ) {
02465
02466 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02467 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02468 posx = x - llx[i] ;
02469 posy = y - lly[i] ;
02470 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02471 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0. &&
02472 (cpl_vector_get(msk,i) != 0)) {
02473
02474
02475
02476
02477
02478 msk_sum+=mask->plane[z]->data[x+y*mask->lx];
02479 val_msk_sum+=tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]*
02480 mask->plane[z]->data[x+y*mask->lx];
02481 }
02482 }
02483 }
02484 mergedCube->plane[z]->data[x+y*mergedCube->lx]=val_msk_sum/msk_sum;
02485 cpl_vector_delete(msk);
02486 }
02487 }
02488 }
02489 m++;
02490 }
02491
02492
02493
02494
02495
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512
02513
02514
02515
02516
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526
02527
02528
02529
02530
02531
02532
02533
02534
02535
02536
02537
02538 free(kernel) ;
02539 for( i = 0 ; i < n_cubes ; i++ )
02540 {
02541 destroy_cube (tmpcubes[i]) ;
02542 }
02543 return 0 ;
02544 }
02545
02546
02547
02548
02549
02550
02551 int combineJitteredCubesSkyRange2 ( OneCube ** cubes,
02552 OneCube * mergedCube,
02553 OneCube * mask,
02554 OneCube * mergedSky,
02555 OneCube * msk_cube,
02556 OneCube * med_cube,
02557 OneCube * avg_cube,
02558 OneCube * sig_cube,
02559 OneCube * nc_cube,
02560 int n_cubes,
02561 float * cumoffsetx,
02562 float * cumoffsety,
02563 float * exptimes,
02564 char * kernel_type,
02565 const int z_min,
02566 const int z_max )
02567 {
02568
02569 const char* _id="combineJitteredCubesSky";
02570 int i ;
02571 int x, y, z ;
02572 int k=1;
02573 int llx[n_cubes], lly[n_cubes] ;
02574 int llx0, lly0 ;
02575 int posx, posy ;
02576 float weight ;
02577
02578 double med=0;
02579 double sig=0;
02580 double avg=0;
02581 int nc=0;
02582 float sub_offsetx[n_cubes] ;
02583 float sub_offsety[n_cubes] ;
02584 double * kernel ;
02585 cpl_vector* val=NULL;
02586 double msk_sum=0;
02587 double val_msk_sum=0;
02588
02589 OneCube * tmpcubes[n_cubes] ;
02590 pixelvalue * tmpspace;
02591 const int z_siz=z_max-z_min;
02592 int m=0;
02593
02594 if ( cubes == NULL )
02595 {
02596 cpl_msg_error ("combineJitteredCubes:","no cube list given!") ;
02597 return -1 ;
02598 }
02599 if ( n_cubes <= 0 )
02600 {
02601 cpl_msg_error ("combineJitteredCubes:","wrong number of data cubes in list!") ;
02602 return -1 ;
02603 }
02604 if ( cumoffsetx == NULL || cumoffsety == NULL )
02605 {
02606 cpl_msg_error ("combineJitteredCubes:","no cumoffsetx/y given!") ;
02607 return -1 ;
02608 }
02609 if ( exptimes == NULL )
02610 {
02611 cpl_msg_error ("combineJitteredCubes:","no exposure time array given!") ;
02612 return -1 ;
02613 }
02614 if (z_siz <= 0 )
02615 {
02616 cpl_msg_error ("combineJitteredCubes:","z_max <= z_min given!") ;
02617 return -1 ;
02618 }
02619
02620
02621
02622
02623
02624
02625
02626
02627
02628 llx0 = mergedCube->lx/2 - cubes[0]->lx/2 ;
02629 lly0 = mergedCube->ly/2 - cubes[0]->ly/2 ;
02630
02631
02632
02633
02634
02635
02636 for ( i = 0 ; i < n_cubes ; i++ )
02637 {
02638 llx[i] = llx0 - nint(cumoffsetx[i]) ;
02639 sub_offsetx[i] = (float)nint(cumoffsetx[i]) - cumoffsetx[i] ;
02640 lly[i] = lly0 - nint(cumoffsety[i]) ;
02641 sub_offsety[i] = (float)nint(cumoffsety[i]) - cumoffsety[i] ;
02642 }
02643
02644
02645
02646
02647
02648
02649
02650 if ( NULL == (kernel = generate_interpolation_kernel( kernel_type )) )
02651 {
02652 cpl_msg_warning ("combineJitteredCubes:","could not generate desired interpolation kernel or no kernel_typ was given, the default kernel is used now!") ;
02653 }
02654
02655
02656
02657 cpl_msg_info(_id,"drs5");
02658 for ( i = 0 ; i < n_cubes ; i++ )
02659 {
02660 tmpspace = cpl_calloc(cubes[i]->lx, cubes[i]->ly*sizeof(pixelvalue)) ;
02661 tmpcubes[i]=new_cube(cubes[i]->lx,cubes[i]->ly,cubes[i]->np);
02662 for ( z = z_min ; z < z_max ; z++ )
02663 {
02664
02665
02666
02667
02668
02669 tmpcubes[i]->plane[z]=shift_image(cubes[i]->plane[z], sub_offsetx[i],
02670 sub_offsety[i], kernel);
02671 if (tmpcubes[i]->plane[z]==NULL)
02672 {
02673 cpl_msg_error ("combineJitteredCubes:","could not shift image plane no %d in cube no %d!", z, i) ;
02674 destroy_cube(mergedCube) ;
02675 destroy_cube(mask) ;
02676 cpl_free(kernel) ;
02677 return -1 ;
02678 }
02679 }
02680 cpl_free(tmpspace);
02681
02682 }
02683 cpl_msg_info(_id,"drs6");
02684
02685
02686
02687
02688
02689
02690
02691 for ( i = 0 ; i < n_cubes ; i++ )
02692 {
02693
02694 for ( y = 0 ; y < mergedCube->ly ; y++ )
02695 {
02696 for ( x = 0 ; x < mergedCube->lx ; x++ )
02697 {
02698
02699 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02700 x >= llx[i] && x < llx[i]+cubes[i]->lx )
02701 {
02702 posx = x - llx[i] ;
02703 posy = y - lly[i] ;
02704 for ( z = z_min ; z < z_max ; z++ )
02705 {
02706 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02707 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.)
02708 {
02709 mask->plane[z]->data[x+y*mask->lx] += exptimes[i] ;
02710 }
02711 }
02712 }
02713 }
02714 }
02715 }
02716
02717 cpl_msg_info(_id,"drs6a");
02718
02719 for ( i = 0 ; i < n_cubes ; i++ )
02720 {
02721
02722 for ( y = 0 ; y < mergedCube->ly ; y++ )
02723 {
02724 for ( x = 0 ; x < mergedCube->lx ; x++ )
02725 {
02726
02727 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02728 x >= llx[i] && x < llx[i]+cubes[i]->lx )
02729 {
02730 posx = x - llx[i] ;
02731 posy = y - lly[i] ;
02732 m=0;
02733 for ( z = z_min ; z < z_max ; z++ )
02734 {
02735 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02736 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.)
02737 {
02738 msk_cube->plane[m]->data[x+y*msk_cube->lx] += 1 ;
02739 }
02740 m++;
02741 }
02742 }
02743 }
02744 }
02745 }
02746 cpl_msg_info(_id,"drs7");
02747 m=0;
02748 for ( z = z_min ; z < z_max ; z++ ) {
02749
02750 for ( y = 0 ; y < mergedCube->ly ; y++ ) {
02751 for ( x = 0 ; x < mergedCube->lx ; x++ ) {
02752 avg=0;
02753 nc=0;
02754
02755 for ( i = 0 ; i < n_cubes ; i++ ) {
02756 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02757 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02758 posx = x - llx[i] ;
02759 posy = y - lly[i] ;
02760 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02761 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.) {
02762 nc++;
02763 }
02764 }
02765 }
02766 if(nc>0) {
02767 val=cpl_vector_new(nc);
02768 }
02769 nc_cube->plane[m]->data[x+y*nc_cube->lx]=nc;
02770 sig=0;
02771 med=0;
02772 nc=0;
02773
02774 for ( i = 0 ; i < n_cubes ; i++ ) {
02775 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02776 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02777 posx = x - llx[i] ;
02778 posy = y - lly[i] ;
02779 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02780 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.) {
02781 cpl_vector_set(val,nc,(double)tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]);
02782 nc++;
02783 }
02784 }
02785 }
02786
02787 if(nc>0) {
02788 avg=cpl_vector_get_mean(val);
02789 med=cpl_vector_get_median(val);
02790 if(nc>1) {
02791 sig=cpl_vector_get_stdev(val);
02792 } else {
02793 sig=0;
02794 }
02795 med_cube->plane[m]->data[x+y*med_cube->lx]=med;
02796 avg_cube->plane[m]->data[x+y*avg_cube->lx]=avg;
02797 sig_cube->plane[m]->data[x+y*sig_cube->lx]=sig;
02798 cpl_vector_delete(val);
02799 }
02800 for ( i = 0 ; i < n_cubes ; i++ ) {
02801
02802 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02803 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02804 posx = x - llx[i] ;
02805 posy = y - lly[i] ;
02806 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02807 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.) {
02808 if(abs((tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]-med))> k*sig) {
02809 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]=0;
02810 msk_cube->plane[m]->data[x+y*msk_cube->lx]-=1;
02811 }
02812 }
02813 }
02814 }
02815 msk_sum=0;
02816 val_msk_sum=0;
02817 for ( i = 0 ; i < n_cubes ; i++ ) {
02818
02819 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02820 x >= llx[i] && x < llx[i]+cubes[i]->lx ) {
02821 posx = x - llx[i] ;
02822 posy = y - lly[i] ;
02823 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) &&
02824 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] != 0.) {
02825 msk_sum+=msk_cube->plane[m]->data[x+y*msk_cube->lx];
02826 val_msk_sum+=tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]*
02827 msk_cube->plane[m]->data[x+y*msk_cube->lx];
02828 }
02829 }
02830 }
02831 mergedSky->plane[m]->data[x+y*mergedSky->lx]=val_msk_sum/msk_sum;
02832
02833 }
02834 }
02835 m++;
02836 }
02837
02838 cpl_msg_info(_id,"drs8");
02839
02840
02841
02842 for ( i = 0 ; i < n_cubes ; i++ )
02843 {
02844
02845
02846 for ( y = 0 ; y < mergedCube->ly ; y++ )
02847 {
02848 for ( x = 0 ; x < mergedCube->lx ; x++ )
02849 {
02850
02851 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
02852 x >= llx[i] && x < llx[i]+cubes[i]->lx )
02853 {
02854 posx = x - llx[i] ;
02855 posy = y - lly[i] ;
02856 for ( z = z_min ; z < z_max ; z++ )
02857 {
02858 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]))
02859 {
02860 if (mask->plane[z]->data[x+y*mask->lx] != 0.)
02861 {
02862
02863 weight = exptimes[0] / mask->plane[z]->data[x+y*mask->lx] ;
02864 }
02865 else
02866 {
02867 weight = 0. ;
02868 }
02869
02870 mergedCube->plane[z]->data[x+y*mergedCube->lx] +=
02871 weight*tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] ;
02872
02873 }
02874 }
02875 }
02876 }
02877 }
02878 }
02879 cpl_msg_info(_id,"drs9");
02880
02881
02882
02883 cpl_msg_info(_id,"drs9");
02884
02885 cpl_msg_info(_id,"drs9");
02886 free(kernel) ;
02887 for( i = 0 ; i < n_cubes ; i++ )
02888 {
02889 destroy_cube (tmpcubes[i]) ;
02890 }
02891 cpl_msg_info(_id,"drs10");
02892 return 0 ;
02893 }
02894
02895
02896
02897
02898
02899
02900
02901
02902
02903
02904
02905
02906
02907
02908
02909 OneCube * interpolCubeSimple( OneCube * cube,
02910 OneCube * badcube,
02911 int maxdist )
02912 {
02913 OneCube * intercube ;
02914 float goodNeighbors[(2*maxdist+1)*(2*maxdist+1)*(2*maxdist+1) -1] ;
02915 int z, row, col ;
02916 int nx, ny, nz ;
02917 int llx, lly, llz ;
02918 int zi, coli, rowi ;
02919 int n ;
02920
02921 if ( cube == NULL || badcube == NULL )
02922 {
02923 cpl_msg_error("interpolCubeSimple:","no cube given!") ;
02924 return NullCube ;
02925 }
02926 if ( maxdist < 1 )
02927 {
02928 cpl_msg_error("interpolCubeSimple:","wrong maxrad given!") ;
02929 return NullCube ;
02930 }
02931 intercube = copy_cube(cube) ;
02932
02933 for ( z = 0 ; z < cube->np ; z++ )
02934 {
02935 for ( row = 0 ; row < cube->ly ; row++ )
02936 {
02937 for ( col = 0 ; col < cube->lx ; col++ )
02938 {
02939 if ( badcube->plane[z]->data[col+row*cube->lx] == 0 )
02940 {
02941
02942 llx = col - maxdist ;
02943 nx = 2*maxdist +1 ;
02944 if (llx < 0)
02945 {
02946 nx += llx ;
02947 llx = 0 ;
02948 }
02949 if ( llx + nx > cube->lx )
02950 {
02951 nx -= (llx + nx - cube->lx) ;
02952 }
02953
02954 lly = row - maxdist ;
02955 ny = 2*maxdist +1 ;
02956 if (lly < 0)
02957 {
02958 ny += lly ;
02959 lly = 0 ;
02960 }
02961 if ( lly + ny > cube->ly )
02962 {
02963 ny -= (lly + ny - cube->ly) ;
02964 }
02965
02966 llz = z - maxdist ;
02967 nz = 2*maxdist +1 ;
02968 if (llz < 0)
02969 {
02970 nz += llz ;
02971 llz = 0 ;
02972 }
02973 if ( llz + nz > cube->np )
02974 {
02975 nz -= (llz + nz - cube->np) ;
02976 }
02977 n = 0 ;
02978 for ( zi = llz ; zi < llz+nz ; zi++ )
02979 {
02980 for ( rowi = lly ; rowi < lly+ny ; rowi++ )
02981 {
02982 for ( coli = llx ; coli < llx+nx ; coli++ )
02983 {
02984 if ( badcube->plane[zi]->data[coli+rowi*badcube->lx] == 1 )
02985 {
02986 goodNeighbors[n] = cube->plane[zi]->data[coli+rowi*cube->lx] ;
02987 n++ ;
02988 }
02989 }
02990 }
02991 }
02992 if ( n > 0 )
02993 {
02994 intercube->plane[z]->data[col+row*cube->lx] = median(goodNeighbors, n) ;
02995 badcube->plane[z]->data[col+row*cube->lx] = 1 ;
02996 }
02997 else
02998 {
02999 continue ;
03000 }
03001 }
03002 }
03003 }
03004 }
03005 return intercube ;
03006 }
03007
03008
03009
03010
03011
03012
03013
03014
03015
03016
03017
03018
03019
03020
03021
03022
03023
03024
03025
03026
03027
03028
03029
03030
03031
03032
03033
03034
03035
03036
03037
03038
03039
03040
03041
03042
03043
03044 OneCube * combineCubes ( OneCube ** cubes,
03045 OneCube * mergedCube,
03046 int n_cubes,
03047 float * cumoffsetx,
03048 float * cumoffsety,
03049 float factor,
03050 char * kernel_type )
03051 {
03052 int i ;
03053 int x, y, z ;
03054 int llx[n_cubes], lly[n_cubes] ;
03055 int llx0, lly0 ;
03056 int posx, posy ;
03057 OneCube * mask ;
03058 float sub_offsetx[n_cubes] ;
03059 float sub_offsety[n_cubes] ;
03060 double * kernel ;
03061 OneImage * shiftedImage ;
03062 int n, ns ;
03063 OneCube * tmpcubes[n_cubes] ;
03064 double sum , sum2 ;
03065 double mean, sigma ;
03066 float cubedata[n_cubes] ;
03067
03068 if ( cubes == NULL )
03069 {
03070 cpl_msg_error ("combineCubes:","no cube list given!") ;
03071 return NullCube ;
03072 }
03073 if ( n_cubes <= 0 )
03074 {
03075 cpl_msg_error ("combineCubes:","wrong number of data cubes in list!") ;
03076 return NullCube ;
03077 }
03078 if ( cumoffsetx == NULL || cumoffsety == NULL )
03079 {
03080 cpl_msg_error ("combineCubes:","no cumoffsetx/y given!") ;
03081 return NullCube ;
03082 }
03083 if ( factor <= 0. )
03084 {
03085 cpl_msg_error ("combineCubes:","wrong factor given!") ;
03086 return NullCube ;
03087 }
03088
03089 if ( NullCube == (mask = newCube(mergedCube->lx, mergedCube->ly, cubes[0]->np)) )
03090 {
03091 cpl_msg_error ("combineCubes:","could not allocate cube!") ;
03092 return NullCube ;
03093 }
03094 for ( i = 0 ; i < n_cubes ; i++ )
03095 {
03096
03097 tmpcubes[i] = (OneCube*)cpl_malloc(sizeof(OneCube)) ;
03098 tmpcubes[i]->plane = (OneImage**)cpl_calloc(cubes[0]->np , sizeof(OneImage*)) ;
03099
03100 tmpcubes[i]->lx = cubes[0]->lx ;
03101 tmpcubes[i]->ly = cubes[0]->ly ;
03102 tmpcubes[i]->np = cubes[0]->np ;
03103 tmpcubes[i]->nbpix = (ulong32)cubes[0]->lx * (ulong32)cubes[0]->ly * (ulong32)cubes[0]->np ;
03104 tmpcubes[i]->history = (char*)NULL ;
03105 tmpcubes[i]->n_comments = 0 ;
03106 tmpcubes[i]->orig_ptype = BPP_DEFAULT ;
03107 tmpcubes[i]->filename = NULL ;
03108 }
03109 for ( z = 0 ; z < cubes[0]->np ; z++ )
03110 {
03111 tmpcubes[0]->plane[z] = copy_image(cubes[0]->plane[z]) ;
03112 }
03113
03114
03115
03116
03117
03118
03119
03120 llx0 = mergedCube->lx/2 - cubes[0]->lx/2 ;
03121 lly0 = mergedCube->ly/2 - cubes[0]->ly/2 ;
03122
03123
03124
03125
03126
03127
03128 for ( i = 0 ; i < n_cubes ; i++ )
03129 {
03130 llx[i] = llx0 - nint(cumoffsetx[i]) ;
03131 sub_offsetx[i] = (float)nint(cumoffsetx[i]) - cumoffsetx[i] ;
03132 lly[i] = lly0 - nint(cumoffsety[i]) ;
03133 sub_offsety[i] = (float)nint(cumoffsety[i]) - cumoffsety[i] ;
03134 }
03135
03136
03137
03138
03139
03140
03141 if ( NULL == (kernel = generate_interpolation_kernel( kernel_type )) )
03142 {
03143 cpl_msg_warning ("combineJitteredCubes:","could not generate desired interpolation kernel or no kernel_typ was given, the default kernel is used now!") ;
03144 }
03145
03146 for ( i = 0 ; i < n_cubes ; i++ )
03147 {
03148
03149 for ( z = 0 ; z < cubes[0]->np ; z++ )
03150 {
03151 if ( NullImage == (shiftedImage = shiftImage( cubes[i]->plane[z], sub_offsetx[i],
03152 sub_offsety[i], kernel ) ) )
03153 {
03154 cpl_msg_error ("combineJitteredCubes:","could not shift image plane no %d in cube no %d!", z, i) ;
03155 destroy_cube(mergedCube) ;
03156 destroy_cube(mask) ;
03157 cpl_free(kernel) ;
03158 return NullCube ;
03159 }
03160 tmpcubes[i]->plane[z] = copy_image(shiftedImage) ;
03161 destroy_image (shiftedImage) ;
03162 }
03163 }
03164
03165 for ( y = 0 ; y < mergedCube->ly ; y++ )
03166 {
03167 for ( x = 0 ; x < mergedCube->lx ; x++ )
03168 {
03169 for ( z = 0 ; z < mergedCube->np ; z++ )
03170 {
03171 sum = 0. ;
03172 sum2 = 0. ;
03173 n = 0 ;
03174 for ( i = 0 ; i < n_cubes ; i++ )
03175 {
03176
03177 if ( y >= lly[i] && y < lly[i]+cubes[i]->ly &&
03178 x >= llx[i] && x < llx[i]+cubes[i]->lx )
03179 {
03180 posx = x - llx[i] ;
03181 posy = y - lly[i] ;
03182 if (!qfits_isnan(tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]))
03183 {
03184 sum += tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] ;
03185 sum2 += (tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] *
03186 tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx]) ;
03187 cubedata[n] = tmpcubes[i]->plane[z]->data[posx+posy*cubes[i]->lx] ;
03188 n++ ;
03189 }
03190 }
03191 }
03192
03193 if ( n == 0 )
03194 {
03195 mean = 0. ;
03196 sigma = 0. ;
03197 mergedCube->plane[z]->data[x+y*mergedCube->lx] = 0. ;
03198 mask->plane[z]->data[x+y*mergedCube->lx] = 0 ;
03199 }
03200 else if ( n == 1 )
03201 {
03202 mean = sum ;
03203 sigma = 0. ;
03204 mergedCube->plane[z]->data[x+y*mergedCube->lx] = mean ;
03205 mask->plane[z]->data[x+y*mergedCube->lx] = 1 ;
03206 }
03207 else
03208 {
03209 mean = sum/(double)n ;
03210 sigma = sqrt( (sum2 - sum*mean) / (double)(n - 1) ) ;
03211 ns = 0 ;
03212 for ( i = 0 ; i < n ; i++ )
03213 {
03214 if ( cubedata[i] > mean+factor*sigma || cubedata[i] < mean-factor*sigma )
03215 {
03216 continue ;
03217 }
03218 else
03219 {
03220 mergedCube->plane[z]->data[x+y*mergedCube->lx] += cubedata[i] ;
03221 ns++ ;
03222 }
03223 }
03224 if ( ns == 0 )
03225 {
03226 mergedCube->plane[z]->data[x+y*mergedCube->lx] = 0. ;
03227 }
03228 else
03229 {
03230 mergedCube->plane[z]->data[x+y*mergedCube->lx] /= (float)ns ;
03231 }
03232 mask->plane[z]->data[x+y*mergedCube->lx] = (float)ns ;
03233 }
03234 }
03235 }
03236 }
03237
03238 for( i = 0 ; i < n_cubes ; i++ )
03239 {
03240 destroy_cube (tmpcubes[i]) ;
03241 }
03242
03243
03244 convert0ToZEROForCubes(mergedCube) ;
03245 cpl_free(kernel) ;
03246 return mask ;
03247 }
03248
03249 OneCube * binCube(OneCube *cu, int xscale, int yscale, int xmin, int xmax, int ymin, int ymax)
03250 {
03251 int i,j,k;
03252 OneCube * cube;
03253
03254 if (NULL == (cube = newCube (xmax-xmin+1,ymax-ymin+1, cu->np)) )
03255 {
03256 cpl_msg_error ("binCube:","cannot allocate new cube") ;
03257 return NullCube ;
03258 }
03259
03260 for (i=0;i<cu->np;i++){
03261 for (j=0 ; j < cube->lx ; j++) {
03262 for (k=0 ; k< cube->ly ; k++) {
03263 cube->plane[i]->data[j+k*cube->lx]=cu->plane[i]->data[((int) (j+xmin)/xscale)+((int) (k+ymin)/yscale)*cu->lx]/(xscale*yscale);
03264 }
03265 }
03266 }
03267
03268
03269
03270
03271
03272
03273
03274
03275
03276
03277
03278
03279
03280
03281
03282 return cube;
03283 }
03284
03285
03286 OneCube * scaleCube(OneCube *cu, float xscale, float yscale, char * kernel_type)
03287 {
03288 OneCube * cube ;
03289 int i, j, k, l ;
03290 int lx_out, ly_out ;
03291 double cur ;
03292 double * invert_transform ;
03293 double neighbors[16] ;
03294 double rsc[8],
03295 sumrs ;
03296 double param[6];
03297 double x, y ;
03298 int px, py ;
03299 int pos ;
03300 int tabx, taby ;
03301 double * kernel ;
03302 int leaps[16] ;
03303
03304 if (cu == NULL)
03305 {
03306 cpl_msg_error ("scalerecCube:","null cube") ;
03307 return NullCube ;
03308 }
03309
03310 param[0]=xscale;
03311 param[1]=0;
03312 param[2]=0;
03313 param[3]=0;
03314 param[4]=yscale;
03315 param[5]=0;
03316
03317
03318 invert_transform = invert_linear_transform(param) ;
03319 if (invert_transform == NULL) {
03320 cpl_msg_error("scalerecCube:","cannot compute invert transform: aborting warping") ;
03321 return NULL ;
03322 }
03323
03324
03325 kernel = generate_interpolation_kernel(kernel_type) ;
03326 if (kernel == NULL) {
03327 cpl_msg_error("scalerecCube:","cannot generate kernel: aborting resampling") ;
03328 return NULL ;
03329 }
03330
03331
03332
03333 lx_out = (int) cu->lx*xscale ;
03334 ly_out = (int) cu->ly*yscale ;
03335
03336 if (NULL == (cube = newCube (lx_out, ly_out, cu->np)) )
03337 {
03338 cpl_msg_error ("scalerecCube:"," cannot allocate new cube") ;
03339 return NullCube ;
03340 }
03341
03342
03343 for (l=0;l<cu->np;l++){
03344
03345
03346 leaps[0] = -1 - cu->plane[l]->lx ;
03347 leaps[1] = - cu->plane[l]->lx ;
03348 leaps[2] = 1 - cu->plane[l]->lx ;
03349 leaps[3] = 2 - cu->plane[l]->lx ;
03350
03351 leaps[4] = -1 ;
03352 leaps[5] = 0 ;
03353 leaps[6] = 1 ;
03354 leaps[7] = 2 ;
03355
03356 leaps[8] = -1 + cu->plane[l]->lx ;
03357 leaps[9] = cu->plane[l]->lx ;
03358 leaps[10]= 1 + cu->plane[l]->lx ;
03359 leaps[11]= 2 + cu->plane[l]->lx ;
03360
03361 leaps[12]= -1 + 2*cu->plane[l]->lx ;
03362 leaps[13]= 2*cu->plane[l]->lx ;
03363 leaps[14]= 1 + 2*cu->plane[l]->lx ;
03364 leaps[15]= 2 + 2*cu->plane[l]->lx ;
03365
03366
03367 for (j=0 ; j < ly_out ; j++) {
03368 for (i=0 ; i< lx_out ; i++) {
03369
03370
03371 x = invert_transform[0] * (double)i +
03372 invert_transform[1] * (double)j +
03373 invert_transform[2] ;
03374
03375 y = invert_transform[3] * (double)i +
03376 invert_transform[4] * (double)j +
03377 invert_transform[5] ;
03378
03379
03380 px = (int)x ;
03381 py = (int)y ;
03382
03383 if ((px < 1) ||
03384 (px > (cu->plane[l]->lx-2)) ||
03385 (py < 1) ||
03386 (py > (cu->plane[l]->ly-2)))
03387 cube->plane[l]->data[i+j*lx_out] = (pixelvalue)0.0 ;
03388 else {
03389
03390 pos = px + py * cu->plane[l]->lx ;
03391 for (k=0 ; k<16 ; k++){
03392 if(!qfits_isnan(cu->plane[l]->data[(int)(pos+leaps[k])])) neighbors[k] = (double)(cu->plane[l]->data[(int)(pos+leaps[k])]) ;
03393 else neighbors[k]=0;
03394 }
03395
03396
03397 tabx = (x - (double)px) * (double)(TABSPERPIX) ;
03398 taby = (y - (double)py) * (double)(TABSPERPIX) ;
03399
03400
03401
03402
03403 rsc[0] = kernel[TABSPERPIX + tabx] ;
03404 rsc[1] = kernel[tabx] ;
03405 rsc[2] = kernel[TABSPERPIX - tabx] ;
03406 rsc[3] = kernel[2 * TABSPERPIX - tabx] ;
03407 rsc[4] = kernel[TABSPERPIX + taby] ;
03408 rsc[5] = kernel[taby] ;
03409 rsc[6] = kernel[TABSPERPIX - taby] ;
03410 rsc[7] = kernel[2 * TABSPERPIX - taby] ;
03411
03412 sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) *
03413 (rsc[4]+rsc[5]+rsc[6]+rsc[7]) ;
03414
03415
03416 cur = rsc[4] * ( rsc[0]*neighbors[0] +
03417 rsc[1]*neighbors[1] +
03418 rsc[2]*neighbors[2] +
03419 rsc[3]*neighbors[3] ) +
03420 rsc[5] * ( rsc[0]*neighbors[4] +
03421 rsc[1]*neighbors[5] +
03422 rsc[2]*neighbors[6] +
03423 rsc[3]*neighbors[7] ) +
03424 rsc[6] * ( rsc[0]*neighbors[8] +
03425 rsc[1]*neighbors[9] +
03426 rsc[2]*neighbors[10] +
03427 rsc[3]*neighbors[11] ) +
03428 rsc[7] * ( rsc[0]*neighbors[12] +
03429 rsc[1]*neighbors[13] +
03430 rsc[2]*neighbors[14] +
03431 rsc[3]*neighbors[15] ) ;
03432
03433
03434 cube->plane[l]->data[i+j*lx_out] = (pixelvalue)(cur/sumrs) ;
03435
03436 }
03437 }
03438 }
03439 }
03440 cpl_free(kernel) ;
03441 cpl_free(invert_transform) ;
03442 return cube ;
03443 }
03444
03445