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
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228 #ifdef HAVE_CONFIG_H
00229 # include <config.h>
00230 #endif
00231
00232 #include "sinfo_vltPort.h"
00233
00234
00235
00236
00237
00238
00239
00240
00241 #include "sinfo_dfs.h"
00242 #include "sinfo_new_cube_ops.h"
00243 #include <sys/types.h>
00244 #include <sys/times.h>
00245 #include <math.h>
00246 #include "sinfo_resampling.h"
00247 #include "sinfo_function_1d.h"
00248 #include "sinfo_error.h"
00249 #include "sinfo_globals.h"
00250 #include "sinfo_utils_wrappers.h"
00251
00252
00253
00277 cpl_imagelist *
00278 sinfo_new_cube_ops(
00279 cpl_imagelist * cube1,
00280 cpl_imagelist * cube2,
00281 int operation)
00282 {
00283
00284 if (cube1==NULL || cube2==NULL)
00285 {
00286 sinfo_msg_error("null cubes");
00287 return NULL ;
00288 }
00289
00290 switch(operation)
00291 {
00292 case '+':
00293 return sinfo_new_cube_add(cube1, cube2) ;
00294 break ;
00295 case '-':
00296 return sinfo_new_cube_sub(cube1, cube2) ;
00297 break ;
00298
00299 case '*':
00300 return sinfo_new_cube_mul(cube1, cube2) ;
00301 break ;
00302
00303 case '/':
00304 return sinfo_new_cube_div(cube1, cube2) ;
00305 break ;
00306
00307 default:
00308 sinfo_msg_error("illegal requested operation: aborting cube arithmetic") ;
00309 return NULL ;
00310 }
00311 }
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330 cpl_imagelist *
00331 sinfo_new_cube_const_ops(
00332 cpl_imagelist * c1,
00333 double constant,
00334 int operation)
00335 {
00336 int ilx1=0;
00337 int ily1=0;
00338 int inp1=0;
00339 cpl_imagelist* c2=NULL;
00340 cpl_image* img1=NULL;
00341
00342
00343
00344 if (c1 == NULL)
00345 {
00346 sinfo_msg_error("null cube") ;
00347 return NULL ;
00348 }
00349 inp1=cpl_imagelist_get_size(c1);
00350 img1=cpl_imagelist_get(c1,0);
00351 ilx1=cpl_image_get_size_x(img1);
00352 ily1=cpl_image_get_size_y(img1);
00353
00354
00355
00356
00357
00358 if ((constant == 0.0) && (operation == '/'))
00359 {
00360 sinfo_msg_error("division by zero requested "
00361 "in cube/constant operation") ;
00362 return NULL ;
00363 }
00364
00365 if ( NULL == (c2 = cpl_imagelist_new()) )
00366 {
00367 sinfo_msg_error ("cannot allocate new cube" ) ;
00368 return NULL ;
00369 }
00370
00371 c2=cpl_imagelist_duplicate(c1);
00372 if(operation == '+') {
00373 cpl_imagelist_add_scalar(c2,constant);
00374 } else if (operation == '-') {
00375 cpl_imagelist_subtract_scalar(c2,constant);
00376 } else if (operation == '*') {
00377 cpl_imagelist_multiply_scalar(c2,constant);
00378 } else if (operation == '/') {
00379 cpl_imagelist_divide_scalar(c2,constant);
00380
00381 } else {
00382 sinfo_msg_error("operation not supported");
00383 return NULL;
00384 }
00385 return c2 ;
00386 }
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396 cpl_imagelist *
00397 sinfo_new_cube_sub(
00398 cpl_imagelist * c1,
00399 cpl_imagelist * c2
00400 )
00401 {
00402 cpl_imagelist * c3 ;
00403 ulong32 i ;
00404 int np ;
00405 int ilx1=0;
00406 int ily1=0;
00407 int inp1=0;
00408 int ilx2=0;
00409 int ily2=0;
00410 int inp2=0;
00411
00412
00413 cpl_image* i_img=NULL;
00414 cpl_image* img1=NULL;
00415 cpl_image* img2=NULL;
00416 cpl_image* img3=NULL;
00417 float* p1data=NULL;
00418 float* p2data=NULL;
00419 float* p3data=NULL;
00420
00421
00422
00423 inp1=cpl_imagelist_get_size(c1);
00424 i_img=cpl_imagelist_get(c1,0);
00425 ilx1=cpl_image_get_size_x(i_img);
00426 ily1=cpl_image_get_size_y(i_img);
00427
00428
00429 inp2=cpl_imagelist_get_size(c2);
00430 i_img=cpl_imagelist_get(c2,0);
00431 ilx2=cpl_image_get_size_x(i_img);
00432 ily2=cpl_image_get_size_y(i_img);
00433
00434 if ((ilx1 != ilx2) ||
00435 (ily1 != ily2))
00436 {
00437 sinfo_msg_error("incompatible size: cannot subtract") ;
00438 return NULL ;
00439 }
00440
00441 if ((inp2 != inp1) &&
00442 (inp2 != 1))
00443 {
00444 sinfo_msg_error("cannot compute with these number of planes") ;
00445 return NULL ;
00446 }
00447
00448 if ( NULL == (c3 = cpl_imagelist_new()) )
00449 {
00450 sinfo_msg_error ("cannot allocate new cube" ) ;
00451 return NULL ;
00452 }
00453
00454 for (np=0 ; np < inp1 ; np++)
00455 {
00456 img3=cpl_image_new(ilx1,ily1,CPL_TYPE_FLOAT);
00457 cpl_imagelist_set(c3,img3,np);
00458 }
00459
00460
00461 for (np=0 ; np < inp1 ; np++)
00462 {
00463 img1=cpl_imagelist_get(c1,np);
00464 p1data=cpl_image_get_data_float(img1);
00465 img2=cpl_imagelist_get(c2,np);
00466 p2data=cpl_image_get_data_float(img2);
00467 img3=cpl_imagelist_get(c3,np);
00468 p3data=cpl_image_get_data_float(img3);
00469
00470 for (i=0 ; i< (ulong32)ilx1*ily1 ; i++)
00471 {
00472 p3data[i] = p1data[i] - p2data[i] ;
00473 }
00474 }
00475
00476 return c3 ;
00477 }
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487 cpl_imagelist *
00488 sinfo_new_cube_add(
00489 cpl_imagelist * c1,
00490 cpl_imagelist * c2
00491 )
00492 {
00493 cpl_imagelist * c3 ;
00494 ulong32 i ;
00495 int np ;
00496 int ilx1=0;
00497 int ily1=0;
00498 int inp1=0;
00499 int ilx2=0;
00500 int ily2=0;
00501 int inp2=0;
00502
00503
00504 cpl_image* i_img=NULL;
00505 cpl_image* img1=NULL;
00506 cpl_image* img2=NULL;
00507 cpl_image* img3=NULL;
00508 float* p1data=NULL;
00509 float* p2data=NULL;
00510 float* p3data=NULL;
00511
00512
00513
00514 inp1=cpl_imagelist_get_size(c1);
00515 i_img=cpl_imagelist_get(c1,0);
00516 ilx1=cpl_image_get_size_x(i_img);
00517 ily1=cpl_image_get_size_y(i_img);
00518
00519
00520 inp2=cpl_imagelist_get_size(c2);
00521 i_img=cpl_imagelist_get(c2,0);
00522 ilx2=cpl_image_get_size_x(i_img);
00523 ily2=cpl_image_get_size_y(i_img);
00524 if ((ilx1 != ilx2) || (ily1 != ily2))
00525 {
00526 sinfo_msg_error("incompatible size: cannot add") ;
00527 return NULL ;
00528 }
00529 if ((inp2 != inp1) && (inp2 != 1))
00530 {
00531 sinfo_msg_error("cannot compute with these number of planes") ;
00532 return NULL ;
00533 }
00534
00535 if (NULL == (c3 = cpl_imagelist_new()) )
00536 {
00537 sinfo_msg_error ("cannot allocate new cube") ;
00538 return NULL ;
00539 }
00540
00541 for (np=0 ; np < inp1 ; np++)
00542 {
00543 img3=cpl_image_new(ilx1,ily1,CPL_TYPE_FLOAT);
00544 cpl_imagelist_set(c3,img3,np);
00545 }
00546
00547 for (np=0 ; np < inp1 ; np++)
00548 {
00549 img1=cpl_imagelist_get(c1,np);
00550 p1data=cpl_image_get_data_float(img1);
00551 img2=cpl_imagelist_get(c2,np);
00552 p2data=cpl_image_get_data_float(img2);
00553 img3=cpl_imagelist_get(c3,np);
00554 p3data=cpl_image_get_data_float(img3);
00555 for (i=0 ; i< (ulong32)ilx1*ily1 ; i++)
00556 {
00557 p3data[i] = p1data[i] + p2data[i] ;
00558 }
00559 }
00560
00561 return c3 ;
00562 }
00563
00564
00565
00566
00567
00568
00569
00570
00571 cpl_imagelist *
00572 sinfo_new_cube_mul(
00573 cpl_imagelist * c1,
00574 cpl_imagelist * c2
00575 )
00576 {
00577 cpl_imagelist *c3 ;
00578 ulong32 i ;
00579 int np ;
00580 int ilx1=0;
00581 int ily1=0;
00582 int inp1=0;
00583 int ilx2=0;
00584 int ily2=0;
00585 int inp2=0;
00586
00587
00588 cpl_image* i_img=NULL;
00589 cpl_image* img1=NULL;
00590 cpl_image* img2=NULL;
00591 cpl_image* img3=NULL;
00592 float* p1data=NULL;
00593 float* p2data=NULL;
00594 float* p3data=NULL;
00595
00596
00597
00598
00599 inp1=cpl_imagelist_get_size(c1);
00600 i_img=cpl_imagelist_get(c1,0);
00601 ilx1=cpl_image_get_size_x(i_img);
00602 ily1=cpl_image_get_size_y(i_img);
00603
00604
00605 inp2=cpl_imagelist_get_size(c2);
00606 i_img=cpl_imagelist_get(c2,0);
00607 ilx2=cpl_image_get_size_x(i_img);
00608 ily2=cpl_image_get_size_y(i_img);
00609
00610 if ((ilx1 != ilx2) || (ily1 != ily2))
00611 {
00612 sinfo_msg_error("incompatible size: cannot multiply") ;
00613 return NULL ;
00614 }
00615
00616 if ((inp2 != inp1) && (inp2 != 1))
00617 {
00618 sinfo_msg_error("cannot compute with these number of planes") ;
00619 return NULL ;
00620 }
00621
00622 if ( NULL == (c3 = cpl_imagelist_new()) )
00623 {
00624 sinfo_msg_error ("cannot allocate new cube" ) ;
00625 return NULL ;
00626 }
00627
00628
00629 for (np=0 ; np < inp1 ; np++)
00630 {
00631 img3=cpl_image_new(ilx1,ily1,CPL_TYPE_FLOAT);
00632 cpl_imagelist_set(c3,img3,np);
00633 }
00634
00635 for (np=0 ; np < inp1 ; np++)
00636 {
00637 img1=cpl_imagelist_get(c1,np);
00638 p1data=cpl_image_get_data_float(img1);
00639 img2=cpl_imagelist_get(c2,np);
00640 p2data=cpl_image_get_data_float(img2);
00641 img3=cpl_imagelist_get(c3,np);
00642 p3data=cpl_image_get_data_float(img3);
00643 for (i=0 ; i< (ulong32)ilx1*ilx2 ; i++)
00644 {
00645 p3data[i] = p1data[i] * p2data[i] ;
00646 }
00647 }
00648
00649 return c3 ;
00650 }
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660 cpl_imagelist *
00661 sinfo_new_cube_div(
00662 cpl_imagelist * c1,
00663 cpl_imagelist * c2
00664 )
00665 {
00666 cpl_imagelist * c3 ;
00667 ulong32 i ;
00668 int np ;
00669 int ilx1=0;
00670 int ily1=0;
00671 int inp1=0;
00672 int ilx2=0;
00673 int ily2=0;
00674 int inp2=0;
00675
00676
00677 cpl_image* i_img=NULL;
00678 cpl_image* img1=NULL;
00679 cpl_image* img2=NULL;
00680 cpl_image* img3=NULL;
00681 float* p1data=NULL;
00682 float* p2data=NULL;
00683 float* p3data=NULL;
00684
00685
00686 inp1=cpl_imagelist_get_size(c1);
00687 i_img=cpl_imagelist_get(c1,0);
00688 ilx1=cpl_image_get_size_x(i_img);
00689 ily1=cpl_image_get_size_y(i_img);
00690
00691
00692 inp2=cpl_imagelist_get_size(c2);
00693 i_img=cpl_imagelist_get(c2,0);
00694 ilx2=cpl_image_get_size_x(i_img);
00695 ily2=cpl_image_get_size_y(i_img);
00696
00697 if ((ilx1 != ilx2) ||
00698 (ily1 != ily2))
00699 {
00700 sinfo_msg_error("incompatible size: cannot divide") ;
00701 return NULL ;
00702 }
00703
00704 if ((inp2 != inp1) && (inp2 != 1))
00705 {
00706 sinfo_msg_error("cannot compute with these number of planes") ;
00707 return NULL ;
00708 }
00709
00710 if (NULL == (c3 = cpl_imagelist_new()) )
00711 {
00712 sinfo_msg_error ("cannot allocate a new cube") ;
00713 return NULL ;
00714 }
00715
00716 for (np=0 ; np < inp1 ; np++)
00717 {
00718 img3=cpl_image_new(ilx1,ily1,CPL_TYPE_FLOAT);
00719 cpl_imagelist_set(c3,img3,np);
00720 }
00721
00722 for (np=0 ; np < inp1 ; np++)
00723 {
00724 img1=cpl_imagelist_get(c1,np);
00725 p1data=cpl_image_get_data_float(img1);
00726 img2=cpl_imagelist_get(c2,np);
00727 p2data=cpl_image_get_data_float(img2);
00728 img3=cpl_imagelist_get(c3,np);
00729 p3data=cpl_image_get_data_float(img3);
00730
00731
00732 for (i=0 ; i< (ulong32) ilx1*ily1 ; i++)
00733 {
00734 if (fabs((double)p2data[i]) < 1e-10)
00735 {
00736 p3data[i] = 0.0 ;
00737 }
00738 else
00739 {
00740 p3data[i] = p1data[i] / p2data[i] ;
00741 }
00742 }
00743 }
00744
00745 return c3 ;
00746 }
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757 cpl_imagelist *
00758 sinfo_new_add_image_to_cube(cpl_imagelist * cu, cpl_image * im)
00759 {
00760 cpl_imagelist * cube ;
00761 int i ;
00762 int clx=0;
00763 int cly=0;
00764 int cnp=0;
00765 int ilx=0;
00766 int ily=0;
00767
00768
00769 cpl_image* i_img=NULL;
00770
00771 if (cu==NULL || im==NULL)
00772 {
00773 sinfo_msg_error ("null cube or null image") ;
00774 return NULL ;
00775 }
00776 cnp=cpl_imagelist_get_size(cu);
00777 i_img=cpl_imagelist_get(cu,0);
00778 clx=cpl_image_get_size_x(i_img);
00779 cly=cpl_image_get_size_y(i_img);
00780
00781 ilx=cpl_image_get_size_x(im);
00782 ily=cpl_image_get_size_y(im);
00783
00784 if ((clx != ilx) || (cly != ily))
00785 {
00786 sinfo_msg_error("incompatible size: cannot add image to cube") ;
00787 return NULL ;
00788 }
00789
00790 cube = cpl_imagelist_duplicate (cu) ;
00791
00792 for (i=0 ; i<cnp ; i++)
00793 {
00794
00795
00796
00797 cpl_image_add(cpl_imagelist_get(cube,i), im) ;
00798 }
00799
00800 return cube ;
00801 }
00802
00803
00804
00805
00806
00807
00808
00809
00810 cpl_imagelist *
00811 sinfo_new_sub_image_from_cube (cpl_imagelist * cu, cpl_image * im)
00812 {
00813 cpl_imagelist * cube ;
00814 int i ;
00815 int clx=0;
00816 int cly=0;
00817 int cnp=0;
00818 int ilx=0;
00819 int ily=0;
00820
00821
00822 cpl_image* i_img=NULL;
00823
00824 if (cu==NULL || im==NULL)
00825 {
00826 sinfo_msg_error ("null cube or null image") ;
00827 return NULL ;
00828 }
00829 cnp=cpl_imagelist_get_size(cu);
00830 i_img=cpl_imagelist_get(cu,0);
00831 clx=cpl_image_get_size_x(i_img);
00832 cly=cpl_image_get_size_y(i_img);
00833
00834 ilx=cpl_image_get_size_x(im);
00835 ily=cpl_image_get_size_y(im);
00836
00837 if ((clx != ilx) || (cly != ily))
00838 {
00839
00840 sinfo_msg_error("incompatible size: cannot subtract image from cube") ;
00841 return NULL ;
00842 }
00843
00844 cube = cpl_imagelist_duplicate (cu) ;
00845
00846 for (i=0 ; i<cnp ; i++)
00847 {
00848
00849
00850
00851 cpl_image_subtract(cpl_imagelist_get(cube,i), im) ;
00852 }
00853 return cube ;
00854 }
00855
00856
00857
00858
00859
00860
00861
00862
00863 cpl_imagelist *
00864 sinfo_new_mul_image_to_cube(cpl_imagelist * cu, cpl_image * im)
00865 {
00866 cpl_imagelist * cube ;
00867 int i ;
00868 int clx=0;
00869 int cly=0;
00870 int cnp=0;
00871 int ilx=0;
00872 int ily=0;
00873
00874
00875 cpl_image* i_img=NULL;
00876
00877 if (cu==NULL || im==NULL)
00878 {
00879 sinfo_msg_error("null cube or null image") ;
00880 return NULL ;
00881 }
00882 cnp=cpl_imagelist_get_size(cu);
00883 i_img=cpl_imagelist_get(cu,0);
00884 clx=cpl_image_get_size_x(i_img);
00885 cly=cpl_image_get_size_y(i_img);
00886
00887 ilx=cpl_image_get_size_x(im);
00888 ily=cpl_image_get_size_y(im);
00889
00890 if ((clx != ilx) || (cly != ily))
00891 {
00892 sinfo_msg_error("incompatible size: cannot multiply image by cube") ;
00893 return NULL ;
00894 }
00895
00896 cube = cpl_imagelist_duplicate (cu) ;
00897
00898 for (i=0 ; i<cnp ; i++)
00899 {
00900
00901
00902
00903 cpl_image_multiply(cpl_imagelist_get(cube,i), im) ;
00904 }
00905
00906 return cube ;
00907 }
00908
00909
00910
00911
00912
00913
00914
00915
00916 cpl_imagelist *
00917 sinfo_new_div_cube_by_image(cpl_imagelist * cu, cpl_image * im)
00918 {
00919 cpl_imagelist * cube ;
00920 int i ;
00921 int clx=0;
00922 int cly=0;
00923 int cnp=0;
00924 int ilx=0;
00925 int ily=0;
00926
00927
00928 cpl_image* i_img=NULL;
00929
00930 if (cu==NULL || im==NULL)
00931 {
00932 sinfo_msg_error ("null cube or null image") ;
00933 return NULL ;
00934 }
00935 cnp=cpl_imagelist_get_size(cu);
00936 i_img=cpl_imagelist_get(cu,0);
00937 clx=cpl_image_get_size_x(i_img);
00938 cly=cpl_image_get_size_y(i_img);
00939
00940 ilx=cpl_image_get_size_x(im);
00941 ily=cpl_image_get_size_y(im);
00942
00943 if ((clx != ilx) || (cly != ily))
00944 {
00945 sinfo_msg_error("incompatible size: cannot divide cube by image") ;
00946 return NULL ;
00947 }
00948
00949 cube = cpl_imagelist_duplicate (cu) ;
00950
00951 for (i=0 ; i<cnp ; i++)
00952 {
00953
00954
00955
00956 cpl_image_divide(cpl_imagelist_get(cube,i), im) ;
00957 }
00958
00959 return cube ;
00960 }
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971 cpl_imagelist *
00972 sinfo_new_add_spectrum_to_cube(cpl_imagelist *cu, Vector *spec)
00973 {
00974 cpl_imagelist * cube ;
00975 int i ,j ;
00976 int ilx=0;
00977 int ily=0;
00978 int inp=0;
00979 float* pidata=NULL;
00980 float* podata=NULL;
00981 cpl_image* i_img=NULL;
00982 cpl_image* o_img=NULL;
00983
00984 if (cu == NULL || spec == NULL)
00985 {
00986 sinfo_msg_error ("null cube or null spectrum") ;
00987 return NULL ;
00988 }
00989 inp=cpl_imagelist_get_size(cu);
00990 i_img=cpl_imagelist_get(cu,0);
00991 ilx=cpl_image_get_size_x(i_img);
00992 ily=cpl_image_get_size_y(i_img);
00993
00994 if ( inp != spec -> n_elements )
00995 {
00996 sinfo_msg_error("cube length and spectrum length are not compatible") ;
00997 return NULL ;
00998 }
00999
01000 if ( NULL == (cube = cpl_imagelist_new ()) )
01001 {
01002 sinfo_msg_error ("cannot allocate new cube" ) ;
01003 return NULL ;
01004 }
01005 for ( i = 0; i < inp; i++)
01006 {
01007 o_img=cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
01008 cpl_imagelist_set(cube,o_img,i);
01009 }
01010
01011
01012 for ( i = 0; i < inp; i++)
01013 {
01014 i_img=cpl_imagelist_get(cu,i);
01015 pidata=cpl_image_get_data_float(i_img);
01016 o_img=cpl_imagelist_get(cube,i);
01017 podata=cpl_image_get_data_float(o_img);
01018 for ( j = 0; j < (int) ilx*ily; j++)
01019 {
01020 podata[j] = pidata[j] + spec -> data[i] ;
01021 }
01022 }
01023
01024 return cube ;
01025 }
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035 cpl_imagelist *
01036 sinfo_new_sub_spectrum_from_cube(cpl_imagelist *cu, Vector *spec)
01037 {
01038 cpl_imagelist * cube ;
01039 int i ,j ;
01040 int ilx=0;
01041 int ily=0;
01042 int inp=0;
01043 float* pidata=NULL;
01044 float* podata=NULL;
01045 cpl_image* i_img=NULL;
01046 cpl_image* o_img=NULL;
01047
01048 if (cu == NULL || spec == NULL)
01049 {
01050 sinfo_msg_error ("null cube or null spectrum") ;
01051 return NULL ;
01052 }
01053 inp=cpl_imagelist_get_size(cu);
01054 i_img=cpl_imagelist_get(cu,0);
01055 ilx=cpl_image_get_size_x(i_img);
01056 ily=cpl_image_get_size_y(i_img);
01057
01058 if ( inp != spec -> n_elements )
01059 {
01060 sinfo_msg_error("cube length and spectrum length are not compatible") ;
01061 return NULL ;
01062 }
01063
01064 if ( NULL == (cube = cpl_imagelist_new()) )
01065 {
01066 sinfo_msg_error ("cannot allocate new cube" ) ;
01067 return NULL ;
01068 }
01069 for ( i = 0; i < inp; i++)
01070 {
01071 o_img=cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
01072 cpl_imagelist_set(cube,o_img,i);
01073 }
01074
01075 for ( i = 0; i < inp; i++)
01076 {
01077 i_img=cpl_imagelist_get(cu,i);
01078 pidata=cpl_image_get_data_float(i_img);
01079 o_img=cpl_imagelist_get(cube,i);
01080 podata=cpl_image_get_data_float(o_img);
01081 for ( j = 0; j < (int) ilx*ily; j++)
01082 {
01083 if ( isnan(pidata[j]) || isnan(spec -> data[i]) )
01084 {
01085 podata[j] = ZERO ;
01086 }
01087 else
01088 {
01089 podata[j] = pidata[j] - spec -> data[i] ;
01090 }
01091 }
01092 }
01093
01094 return cube ;
01095 }
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106 cpl_imagelist *
01107 sinfo_new_mul_spectrum_to_cube(cpl_imagelist *cu, Vector *spec)
01108 {
01109 cpl_imagelist * cube ;
01110 int i ,j ;
01111 int ilx=0;
01112 int ily=0;
01113 int inp=0;
01114 float* pidata=NULL;
01115 float* podata=NULL;
01116 cpl_image* i_img=NULL;
01117 cpl_image* o_img=NULL;
01118
01119 if (cu == NULL || spec == NULL)
01120 {
01121 sinfo_msg_error ("null cube or null spectrum") ;
01122 return NULL ;
01123 }
01124 inp=cpl_imagelist_get_size(cu);
01125 i_img=cpl_imagelist_get(cu,0);
01126 ilx=cpl_image_get_size_x(i_img);
01127 ily=cpl_image_get_size_y(i_img);
01128
01129 if ( inp != spec -> n_elements )
01130 {
01131 sinfo_msg_error("cube length and spectrum length are not compatible") ;
01132 return NULL ;
01133 }
01134
01135 if ( NULL == (cube = cpl_imagelist_new ()) )
01136 {
01137 sinfo_msg_error ("cannot allocate new cube" ) ;
01138 return NULL ;
01139 }
01140
01141 for ( i = 0; i < inp; i++)
01142 {
01143 o_img=cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
01144 cpl_imagelist_set(cube,o_img,i);
01145 }
01146
01147 for ( i = 0; i < inp; i++)
01148 {
01149 i_img=cpl_imagelist_get(cu,i);
01150 pidata=cpl_image_get_data_float(i_img);
01151 o_img=cpl_imagelist_get(cube,i);
01152 podata=cpl_image_get_data_float(o_img);
01153 for ( j = 0; j < (int) ilx*ily; j++)
01154 {
01155 if ( isnan(pidata[j]) || isnan(spec->data[i]) )
01156 {
01157 podata[j] = ZERO ;
01158 }
01159 else
01160 {
01161 podata[j] = pidata[j] * spec -> data[i] ;
01162 }
01163 }
01164 }
01165
01166 return cube ;
01167 }
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178 cpl_imagelist *
01179 sinfo_new_div_cube_by_spectrum(cpl_imagelist *cu, Vector *spec)
01180 {
01181 cpl_imagelist * cube ;
01182 float help ;
01183 int i ,j ;
01184 int ilx=0;
01185 int ily=0;
01186 int inp=0;
01187 float* pidata=NULL;
01188 float* podata=NULL;
01189 cpl_image* i_img=NULL;
01190 cpl_image* o_img=NULL;
01191
01192 if (cu == NULL || spec == NULL)
01193 {
01194 sinfo_msg_error ("null cube or null spectrum") ;
01195 return NULL ;
01196 }
01197 inp=cpl_imagelist_get_size(cu);
01198 i_img=cpl_imagelist_get(cu,0);
01199 ilx=cpl_image_get_size_x(i_img);
01200 ily=cpl_image_get_size_y(i_img);
01201
01202 if ( inp != spec -> n_elements )
01203 {
01204 sinfo_msg_error("cube length and spectrum length are not compatible") ;
01205 return NULL ;
01206 }
01207
01208 if (NULL == (cube = cpl_imagelist_new ()) )
01209 {
01210 sinfo_msg_error ("cannot allocate new cube") ;
01211 return NULL ;
01212 }
01213
01214 for ( i = 0; i < inp; i++)
01215 {
01216 o_img=cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
01217 cpl_imagelist_set(cube,o_img,i);
01218 }
01219
01220 for ( i = 0; i < inp; i++)
01221 {
01222
01223 i_img=cpl_imagelist_get(cu,i);
01224 pidata=cpl_image_get_data_float(i_img);
01225 o_img=cpl_imagelist_get(cube,i);
01226 podata=cpl_image_get_data_float(o_img);
01227 for ( j = 0; j < (int) ilx*ily; j++)
01228 {
01229 if (!isnan(spec->data[i]) && spec->data[i] != 0.)
01230 {
01231 help = 1/spec->data[i] ;
01232 if ( help > THRESH )
01233 {
01234 help = 1. ;
01235 }
01236 }
01237 else
01238 {
01239 help = ZERO ;
01240 }
01241
01242 if ( isnan(help) || isnan(pidata[j]) )
01243 {
01244 podata[j] = ZERO ;
01245 }
01246 else
01247 {
01248 podata[j] = pidata[j] * help ;
01249 }
01250 }
01251 }
01252 return cube ;
01253 }
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269 Vector *
01270 sinfo_new_clean_mean_of_spectra(cpl_imagelist * cube,
01271 int llx,
01272 int lly,
01273 int urx,
01274 int ury,
01275 double lo_reject,
01276 double hi_reject)
01277 {
01278 Vector * mean ;
01279 pixelvalue *local_rectangle ;
01280 int i, j, k, l, m ;
01281 int recsize, lo_n, hi_n, nv ;
01282 int ilx=0;
01283 int ily=0;
01284 int inp=0;
01285 float* pidata=NULL;
01286 cpl_image* i_img=NULL;
01287
01288 if ( cube == NULL || cpl_imagelist_get_size(cube) < 1 )
01289 {
01290 sinfo_msg_error ("no cube to take the mean of his spectra") ;
01291 return NullVector ;
01292 }
01293 inp=cpl_imagelist_get_size(cube);
01294 i_img=cpl_imagelist_get(cube,0);
01295 ilx=cpl_image_get_size_x(i_img);
01296 ily=cpl_image_get_size_y(i_img);
01297
01298 if ((llx<1) || (llx>ilx) ||
01299 (urx<1) || (urx>ilx) ||
01300 (lly<1) || (lly>ily) ||
01301 (ury<1) || (ury>ily) ||
01302 (llx>=urx) || (lly>=ury))
01303 {
01304 sinfo_msg_error("invalid rectangle coordinates:") ;
01305 sinfo_msg_error("lower left is [%d %d] upper right is [%d %d]",
01306 llx, lly, urx, ury) ;
01307 return NullVector ;
01308 }
01309
01310 if ((lo_reject + hi_reject) > 0.9)
01311 {
01312 sinfo_msg_error("illegal rejection thresholds: [%f] and [%f]",
01313 lo_reject, hi_reject) ;
01314 sinfo_msg_error("threshold sum should not be over 0.9"
01315 " aborting average") ;
01316 return NullVector ;
01317 }
01318
01319
01320 llx -- ;
01321 lly -- ;
01322 urx -- ;
01323 ury -- ;
01324
01325 recsize = (urx - llx + 1) * (ury - lly + 1) ;
01326
01327 lo_n = (int) (recsize * lo_reject + 0.5) ;
01328 hi_n = (int) (recsize * hi_reject + 0.5) ;
01329
01330 if (lo_n + hi_n >= recsize)
01331 {
01332 sinfo_msg_error ("everything would be rejected") ;
01333 return NullVector;
01334 }
01335
01336
01337 if (NULL == (mean = sinfo_new_vector (inp)) )
01338 {
01339 sinfo_msg_error ("cannot allocate a new sinfo_vector") ;
01340 return NullVector ;
01341 }
01342
01343
01344
01345
01346
01347 for ( i = 0 ; i < inp ; i++ )
01348 {
01349 i_img=cpl_imagelist_get(cube,i);
01350 pidata=cpl_image_get_data_float(i_img);
01351 m = 0 ;
01352 local_rectangle=(pixelvalue *)cpl_calloc(recsize, sizeof (pixelvalue*));
01353
01354 for ( j = lly ; j <= ury ; j++ )
01355 {
01356 for ( k = llx ; k <= urx ; k++ )
01357 {
01358 local_rectangle[m] = pidata[k + j * ilx] ;
01359 m ++ ;
01360 }
01361 }
01362
01363 sinfo_pixel_qsort (local_rectangle, recsize) ;
01364
01365 nv = 0 ;
01366 for ( l = lo_n ; l < (recsize - hi_n) ; l++ )
01367 {
01368 mean -> data[i] += local_rectangle[l] ;
01369 nv ++;
01370 }
01371 mean -> data[i] /= nv ;
01372
01373 cpl_free ( local_rectangle ) ;
01374 }
01375 return mean ;
01376 }
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389 cpl_image *
01390 sinfo_new_median_cube(cpl_imagelist * cube)
01391 {
01392 cpl_image * im ;
01393 pixelvalue * buffer ;
01394 int i, j, k, nz ;
01395 int ilx=0;
01396 int ily=0;
01397 int inp=0;
01398 float* pidata=NULL;
01399 float* podata=NULL;
01400 cpl_image* i_img=NULL;
01401
01402 if ( cube == NULL )
01403 {
01404 sinfo_msg_error ("null cube") ;
01405 return NULL ;
01406 }
01407 inp=cpl_imagelist_get_size(cube);
01408 i_img=cpl_imagelist_get(cube,0);
01409 ilx=cpl_image_get_size_x(i_img);
01410 ily=cpl_image_get_size_y(i_img);
01411
01412
01413 if (NULL == (im = cpl_image_new (ilx, ily, CPL_TYPE_FLOAT )) )
01414 {
01415 sinfo_msg_error ("cannot allocate new image") ;
01416 return NULL ;
01417 }
01418
01419
01420
01421
01422
01423
01424 podata=cpl_image_get_data_float(im);
01425 for ( i = 0 ; i < (int) ilx*ily ; i++ )
01426 {
01427 buffer = (pixelvalue *) cpl_calloc (inp, sizeof (pixelvalue *));
01428 k = 0 ;
01429 for ( j = 0 ; j < inp ; j ++ )
01430 {
01431 i_img=cpl_imagelist_get(cube,j);
01432 pidata=cpl_image_get_data_float(i_img);
01433 if ( !isnan(pidata[i]) )
01434 {
01435 buffer[k] = pidata[i] ;
01436 k ++ ;
01437 }
01438 }
01439 nz = k ;
01440
01441
01442 if ( nz > 2 )
01443 {
01444 podata[i] = sinfo_new_median ( buffer, nz ) ;
01445 }
01446 else if (nz == 2)
01447 {
01448 podata[i] = (buffer[0] + buffer[1]) / 2. ;
01449 }
01450 else if (nz == 1)
01451 {
01452 podata[i] = buffer[0] ;
01453 }
01454 else if (nz == 0)
01455 {
01456 podata[i] = ZERO ;
01457 }
01458
01459 cpl_free ( buffer ) ;
01460 }
01461
01462 return im ;
01463 }
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476 cpl_image *
01477 sinfo_new_average_cube_to_image(cpl_imagelist * cube)
01478 {
01479 cpl_image * im ;
01480 int i, j, nz ;
01481 int ilx=0;
01482 int ily=0;
01483 int inp=0;
01484 float* pidata=NULL;
01485 float* podata=NULL;
01486 cpl_image* i_img=NULL;
01487
01488 if ( cube == NULL )
01489 {
01490 sinfo_msg_error ("null cube") ;
01491 return NULL ;
01492 }
01493 inp=cpl_imagelist_get_size(cube);
01494 i_img=cpl_imagelist_get(cube,0);
01495 ilx=cpl_image_get_size_x(i_img);
01496 ily=cpl_image_get_size_y(i_img);
01497
01498
01499 if (NULL == (im = cpl_image_new (ilx, ily,CPL_TYPE_FLOAT )) )
01500 {
01501 sinfo_msg_error ("cannot allocate new image") ;
01502 return NULL ;
01503 }
01504
01505
01506
01507
01508
01509
01510 podata=cpl_image_get_data_float(im);
01511 for ( i = 0 ; i < (int) ilx*ily ; i++ )
01512 {
01513 nz = 0 ;
01514 for ( j = 0 ; j < inp ; j ++ )
01515 {
01516 i_img=cpl_imagelist_get(cube,j);
01517 pidata=cpl_image_get_data_float(i_img);
01518 if ( !isnan(pidata[i]) )
01519 {
01520 nz ++ ;
01521 podata[i] += pidata[i] ;
01522 }
01523 }
01524
01525
01526 if ( nz >= 1 )
01527 {
01528 podata[i] /= nz ;
01529 }
01530 else if (nz == 0)
01531 {
01532 podata[i] = ZERO ;
01533 }
01534 }
01535
01536 return im ;
01537 }
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549 cpl_image *
01550 sinfo_new_sum_cube_to_image(cpl_imagelist * cube)
01551 {
01552 cpl_image * im ;
01553 int i, j, nz ;
01554 int ilx=0;
01555 int ily=0;
01556 int inp=0;
01557 float* pidata=NULL;
01558 float* podata=NULL;
01559 cpl_image* i_img=NULL;
01560
01561 if ( cube == NULL )
01562 {
01563 sinfo_msg_error ("null cube") ;
01564 return NULL ;
01565 }
01566 inp=cpl_imagelist_get_size(cube);
01567 i_img=cpl_imagelist_get(cube,0);
01568 ilx=cpl_image_get_size_x(i_img);
01569 ily=cpl_image_get_size_y(i_img);
01570
01571
01572 if (NULL == (im = cpl_image_new (ilx, ily, CPL_TYPE_FLOAT )) )
01573 {
01574 sinfo_msg_error ("cannot allocate new image") ;
01575 return NULL ;
01576 }
01577
01578
01579
01580
01581
01582
01583 podata=cpl_image_get_data_float(im);
01584 for ( i = 0 ; i < (int) ilx*ily ; i++ )
01585 {
01586 nz = 0 ;
01587 for ( j = 0 ; j < inp ; j ++ )
01588 {
01589 i_img=cpl_imagelist_get(cube,j);
01590 pidata=cpl_image_get_data_float(i_img);
01591 if ( !isnan(pidata[i]) )
01592 {
01593 nz++ ;
01594 podata[i] += pidata[i] ;
01595 }
01596 }
01597
01598
01599 if (nz == 0)
01600 {
01601 podata[i] = ZERO ;
01602 }
01603 }
01604
01605 return im ;
01606 }
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628 cpl_image *
01629 sinfo_new_average_cube_to_image_between_waves (cpl_imagelist * cube,
01630 float dispersion,
01631 float centralWave,
01632 float initialLambda,
01633 float finalLambda)
01634 {
01635 cpl_image * im ;
01636 int firstPlane ;
01637 int lastPlane ;
01638 int i, j, nz ;
01639 float minWave ;
01640 int ilx=0;
01641 int ily=0;
01642 int inp=0;
01643 float* pidata=NULL;
01644 float* podata=NULL;
01645 cpl_image* i_img=NULL;
01646
01647 if ( cube == NULL )
01648 {
01649 sinfo_msg_error ("null cube") ;
01650 return NULL ;
01651 }
01652 i_img=cpl_imagelist_get(cube,0);
01653 ilx=cpl_image_get_size_x(i_img);
01654 ily=cpl_image_get_size_y(i_img);
01655
01656 inp=cpl_imagelist_get_size(cube);
01657
01658 minWave = centralWave - (float) (inp / 2)*dispersion ;
01659
01660 if ( dispersion <= 0. || minWave <= 0. )
01661 {
01662 sinfo_msg_error ("wrong dispersion or minimum wavelength given") ;
01663 return NULL ;
01664 }
01665
01666 if ( initialLambda < minWave ||
01667 (initialLambda >= minWave + dispersion * inp) )
01668 {
01669 sinfo_msg_error ("wrong initial wavelength given") ;
01670 return NULL ;
01671 }
01672
01673 if ( finalLambda <= minWave ||
01674 (finalLambda > minWave + dispersion * inp) )
01675 {
01676 sinfo_msg_error ("wrong final wavelength given") ;
01677 return NULL ;
01678 }
01679
01680
01681 if (NULL == (im = cpl_image_new (ilx, ily, CPL_TYPE_FLOAT )) )
01682 {
01683 sinfo_msg_error ("cannot allocate new image") ;
01684 return NULL ;
01685 }
01686
01687
01688 firstPlane = sinfo_new_nint ((double) ((initialLambda - minWave) /
01689 dispersion)) ;
01690 lastPlane = sinfo_new_nint ((double) ((finalLambda - minWave) /
01691 dispersion)) ;
01692
01693 if ( firstPlane < 0 || firstPlane >= inp ||
01694 lastPlane < 0 || lastPlane > inp )
01695 {
01696 sinfo_msg_error ("wrong values given!") ;
01697 return NULL ;
01698 }
01699
01700
01701
01702
01703
01704
01705
01706
01707 podata=cpl_image_get_data_float(im);
01708 for ( i = 0 ; i < (int) ilx*ily ; i++ )
01709 {
01710 nz = 0 ;
01711
01712 for ( j = firstPlane ; j <= lastPlane ; j ++ )
01713 {
01714 i_img=cpl_imagelist_get(cube,j);
01715 pidata=cpl_image_get_data_float(i_img);
01716 if ( !isnan(pidata[i]) )
01717 {
01718 nz ++ ;
01719 podata[i] += pidata[i] ;
01720 }
01721 }
01722
01723
01724 if ( nz >= 1 )
01725 {
01726 podata[i] /= nz ;
01727 }
01728 else if (nz == 0)
01729 {
01730 podata[i] = ZERO ;
01731 }
01732 }
01733
01734 return im ;
01735 }
01736
01737
01738
01739
01740
01741
01742
01743
01744 cpl_image *
01745 sinfo_new_extract_image_from_cube(cpl_imagelist * cube, int plane_index)
01746 {
01747 if ( cube == NULL )
01748 {
01749 sinfo_msg_error ("null cube") ;
01750 return NULL ;
01751 }
01752
01753 if ( plane_index < 0 || plane_index >= cpl_imagelist_get_size(cube) )
01754 {
01755 sinfo_msg_error ("wrong plane index for image to be extracted") ;
01756 return NULL ;
01757 }
01758
01759 return cpl_imagelist_get(cube,plane_index) ;
01760 }
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770 Vector *
01771 sinfo_new_extract_spectrum_from_cube(cpl_imagelist * cube,
01772 int x_pos, int y_pos)
01773 {
01774 Vector * returnedSpectrum ;
01775 int i ;
01776 int ilx=0;
01777 int ily=0;
01778 int inp=0;
01779 float* pidata=NULL;
01780 cpl_image* i_img=NULL;
01781
01782 if ( cube == NULL )
01783 {
01784 sinfo_msg_error ("no cube given!") ;
01785 return NullVector ;
01786 }
01787 i_img=cpl_imagelist_get(cube,0);
01788 ilx=cpl_image_get_size_x(i_img);
01789 ily=cpl_image_get_size_y(i_img);
01790 inp=cpl_imagelist_get_size(cube);
01791
01792 if ( x_pos < 0 || x_pos >= ilx )
01793 {
01794 sinfo_msg_error ("wrong x-positon of spectrum given!") ;
01795 return NullVector ;
01796 }
01797
01798 if ( y_pos < 0 || y_pos >= ily )
01799 {
01800 sinfo_msg_error ("wrong y-positon of spectrum given!") ;
01801 return NullVector ;
01802 }
01803
01804
01805 if ( NULL == (returnedSpectrum = sinfo_new_vector ( inp )) )
01806 {
01807 sinfo_msg_error ("cannot allocate new spectrum!") ;
01808 return NullVector ;
01809 }
01810
01811 for ( i = 0 ; i < inp ; i++ )
01812 {
01813 i_img=cpl_imagelist_get(cube,i);
01814 pidata=cpl_image_get_data_float(i_img);
01815 returnedSpectrum -> data[i] = pidata[x_pos + ilx*y_pos] ;
01816 }
01817
01818 return returnedSpectrum ;
01819 }
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860 cpl_imagelist *
01861 sinfo_new_combine_jittered_cubes ( cpl_imagelist ** cubes,
01862 cpl_imagelist * mergedCube,
01863 int n_cubes,
01864 float * cumoffsetx,
01865 float * cumoffsety,
01866 float * exptimes,
01867 char * kernel_type )
01868 {
01869
01870 int i=0 ;
01871 int x=0;
01872 int y=0;
01873 int z=0;
01874 int llx0=0;
01875 int lly0=0;
01876 int posx=0;
01877 int posy=0;
01878 float weight=0;
01879 cpl_imagelist * mask=NULL;
01880 double * kernel=NULL;
01881
01882
01883 int* llx=NULL ;
01884 int* lly=NULL ;
01885
01886 float* sub_offsetx=NULL ;
01887 float* sub_offsety=NULL ;
01888
01889 cpl_imagelist ** tmpcubes=NULL ;
01890 pixelvalue * tmpspace=NULL;
01891
01892
01893 int ilx=0;
01894 int ily=0;
01895 int olx=0;
01896 int oly=0;
01897 int mlx=0;
01898 int onp=0;
01899 int inp=0;
01900
01901
01902
01903 float* podata=NULL;
01904 float* pmdata=NULL;
01905 float* ptdata=NULL;
01906
01907 cpl_image* i_img=NULL;
01908 cpl_image* o_img=NULL;
01909 cpl_image* m_img=NULL;
01910 cpl_image* t_img=NULL;
01911
01912
01913 if ( cubes == NULL )
01914 {
01915 sinfo_msg_error ("no cube list given!") ;
01916 return NULL ;
01917 }
01918 if ( n_cubes <= 0 )
01919 {
01920 sinfo_msg_error ("wrong number of data cubes in list!") ;
01921 return NULL ;
01922 }
01923 if ( cumoffsetx == NULL || cumoffsety == NULL )
01924 {
01925 sinfo_msg_error ("no cumoffsetx/y given!") ;
01926 return NULL ;
01927 }
01928 if ( exptimes == NULL )
01929 {
01930 sinfo_msg_error ("no exposure time array given!") ;
01931 return NULL ;
01932 }
01933
01934 o_img=cpl_imagelist_get(mergedCube,0);
01935 olx=cpl_image_get_size_x(o_img);
01936 oly=cpl_image_get_size_y(o_img);
01937 onp=cpl_imagelist_get_size(mergedCube);
01938 if ( NULL == (mask = cpl_imagelist_new()) )
01939 {
01940 sinfo_msg_error ("could not allocate cube!") ;
01941 return NULL ;
01942 }
01943 for(i=0;i<onp;i++){
01944 o_img=cpl_image_new(olx,oly,CPL_TYPE_FLOAT);
01945 cpl_imagelist_set(mergedCube,o_img,i);
01946 }
01947
01948 i_img=cpl_imagelist_get(cubes[0],0);
01949 ilx=cpl_image_get_size_x(i_img);
01950 ily=cpl_image_get_size_y(i_img);
01951
01952 inp=cpl_imagelist_get_size(cubes[0]);
01953
01954
01955
01956
01957
01958
01959
01960 llx0 = olx/2 - ilx/2 ;
01961 lly0 = oly/2 - ily/2 ;
01962
01963
01964
01965
01966
01967
01968
01969 llx=cpl_calloc(n_cubes,sizeof(int)); ;
01970 lly=cpl_calloc(n_cubes,sizeof(int)) ;
01971
01972 sub_offsetx=cpl_calloc(n_cubes,sizeof(float)) ;
01973 sub_offsety=cpl_calloc(n_cubes,sizeof(float)) ;
01974
01975 for ( i = 0 ; i < n_cubes ; i++ )
01976 {
01977 llx[i] = llx0 - sinfo_new_nint(cumoffsetx[i]) ;
01978 sub_offsetx[i] = (float)sinfo_new_nint(cumoffsetx[i]) - cumoffsetx[i] ;
01979 lly[i] = lly0 - sinfo_new_nint(cumoffsety[i]) ;
01980 sub_offsety[i] = (float)sinfo_new_nint(cumoffsety[i]) - cumoffsety[i] ;
01981 }
01982
01983
01984
01985
01986
01987
01988
01989 if ( NULL == (kernel = sinfo_generate_interpolation_kernel(kernel_type)))
01990 {
01991 sinfo_msg_warning ("could not generate desired interpolation kernel"
01992 " or no kernel_typ was given, the default kernel"
01993 " is used now!") ;
01994 }
01995
01996
01997
01998 tmpcubes=(cpl_imagelist**)cpl_calloc(n_cubes,sizeof(cpl_imagelist*)) ;
01999
02000 for ( i = 0 ; i < n_cubes ; i++ )
02001 {
02002 tmpspace = cpl_calloc(ilx, ily*sizeof(pixelvalue)) ;
02003 tmpcubes[i] = cpl_imagelist_new();
02004
02005 for ( z = 0 ; z < inp ; z++ )
02006 {
02007
02008
02009 t_img=sinfo_new_shift_image(cpl_imagelist_get(cubes[i],z),
02010 sub_offsetx[i], sub_offsety[i], kernel);
02011
02012 if (t_img==NULL)
02013 {
02014 sinfo_msg_error ("could not shift image plane no %d"
02015 " in cube no %d!", z, i) ;
02016 cpl_imagelist_delete(mergedCube) ;
02017 cpl_imagelist_delete(mask) ;
02018 cpl_free(kernel) ;
02019 return NULL ;
02020 }
02021 cpl_imagelist_set(tmpcubes[i],t_img,z);
02022 }
02023 cpl_free(tmpspace);
02024 }
02025
02026
02027
02028
02029
02030
02031
02032
02033 for ( i = 0 ; i < n_cubes ; i++ )
02034 {
02035
02036
02037 for ( y = 0 ; y < oly ; y++ )
02038 {
02039 for ( x = 0 ; x < olx ; x++ )
02040 {
02041
02042
02043 if ( y >= lly[i] && y < lly[i]+ily &&
02044 x >= llx[i] && x < llx[i]+ilx )
02045 {
02046 posx = x - llx[i] ;
02047 posy = y - lly[i] ;
02048 for ( z = 0 ; z < onp ; z++ )
02049 {
02050 t_img=cpl_imagelist_get(tmpcubes[i],z);
02051 ptdata=cpl_image_get_data_float(t_img);
02052 m_img=cpl_imagelist_get(mask,z);
02053 pmdata=cpl_image_get_data_float(m_img);
02054 if (!isnan(ptdata[posx+posy*ilx]) &&
02055 ptdata[posx+posy*ilx] != 0.)
02056 {
02057 pmdata[x+y*mlx] += exptimes[i] ;
02058 }
02059 }
02060 }
02061 }
02062 }
02063 }
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073 for ( i = 0 ; i < n_cubes ; i++ )
02074 {
02075
02076
02077 for ( y = 0 ; y < oly ; y++ )
02078 {
02079
02080 for ( x = 0 ; x < olx ; x++ )
02081 {
02082
02083
02084
02085 if ( y >= lly[i] && y < lly[i]+ily &&
02086 x >= llx[i] && x < llx[i]+ilx )
02087 {
02088
02089 posx = x - llx[i] ;
02090 posy = y - lly[i] ;
02091 for ( z = 0 ; z < onp ; z++ )
02092 {
02093
02094 t_img=cpl_imagelist_get(tmpcubes[i],z);
02095 ptdata=cpl_image_get_data_float(t_img);
02096 m_img=cpl_imagelist_get(mask,z);
02097 pmdata=cpl_image_get_data_float(m_img);
02098 mlx=cpl_image_get_size_x(m_img);
02099
02100 o_img=cpl_imagelist_get(mergedCube,z);
02101 podata=cpl_image_get_data_float(o_img);
02102 podata[x+y*olx]=0;
02103 if (!isnan(ptdata[posx+posy*ilx]))
02104 {
02105 if (pmdata[x+y*mlx] != 0.)
02106 {
02107
02108
02109 weight = exptimes[0] / pmdata[x+y*mlx] ;
02110 }
02111 else
02112 {
02113 weight = 0. ;
02114 }
02115 podata[x+y*olx] +=
02116 weight*ptdata[posx+posy*ilx] ;
02117 }
02118 }
02119 }
02120 }
02121 }
02122 }
02123
02124
02125
02126
02127
02128
02129 cpl_free(kernel) ;
02130 for( i = 0 ; i < n_cubes ; i++ )
02131 {
02132 cpl_imagelist_delete (tmpcubes[i]) ;
02133 }
02134
02135 cpl_free(tmpcubes); ;
02136 cpl_free(llx); ;
02137 cpl_free(lly) ;
02138
02139 cpl_free(sub_offsetx) ;
02140 cpl_free(sub_offsety) ;
02141
02142 return mask ;
02143 }
02144
02145
02146
02147
02148
02149 int
02150 sinfo_new_combine_jittered_cubes_range ( cpl_imagelist ** cubes,
02151 cpl_imagelist * mergedCube,
02152 cpl_imagelist * mask,
02153 int n_cubes,
02154 float * cumoffsetx,
02155 float * cumoffsety,
02156 double * exptimes,
02157 char * kernel_type,
02158 const int z_min, const int z_max )
02159 {
02160
02161 int i,m ;
02162 int x, y, z ;
02163 int llx0, lly0 ;
02164 int posx, posy ;
02165 float weight ;
02166 double * kernel ;
02167
02168
02169
02170 cpl_imagelist ** tmpcubes=NULL ;
02171 int* llx=NULL ;
02172 int* lly=NULL ;
02173 float* sub_offsetx=NULL ;
02174 float* sub_offsety=NULL ;
02175
02176 pixelvalue * tmpspace;
02177
02178
02179
02180 int ilx=0;
02181 int ily=0;
02182 int olx=0;
02183 int oly=0;
02184 int mlx=0;
02185 int mly=0;
02186
02187
02188
02189
02190 float* podata=NULL;
02191 float* pmdata=NULL;
02192 float* ptdata=NULL;
02193
02194 cpl_image* i_img=NULL;
02195 cpl_image* o_img=NULL;
02196 cpl_image* m_img=NULL;
02197 cpl_image* t_img=NULL;
02198
02199
02200 if ( cubes == NULL )
02201 {
02202 sinfo_msg_error ("no cube list given!") ;
02203 return -1 ;
02204 }
02205 if ( n_cubes <= 0 )
02206 {
02207 sinfo_msg_error ("wrong number of data cubes in list!") ;
02208 return -1 ;
02209 }
02210 if ( cumoffsetx == NULL || cumoffsety == NULL )
02211 {
02212 sinfo_msg_error ("no cumoffsetx/y given!") ;
02213 return -1;
02214 }
02215 if ( exptimes == NULL )
02216 {
02217 sinfo_msg_error ("no exposure time array given!") ;
02218 return -1 ;
02219 }
02220
02221
02222
02223
02224
02225
02226
02227
02228 o_img=cpl_imagelist_get(mergedCube,z_min);
02229 olx=cpl_image_get_size_x(o_img);
02230 oly=cpl_image_get_size_y(o_img);
02231 i_img=cpl_imagelist_get(cubes[0],0);
02232 ilx=cpl_image_get_size_x(i_img);
02233 ily=cpl_image_get_size_y(i_img);
02234 mlx=olx;
02235 mly=oly;
02236
02237
02238
02239
02240
02241
02242
02243
02244 llx0 = olx/2 - ilx/2 ;
02245 lly0 = oly/2 - ily/2 ;
02246
02247
02248
02249
02250
02251
02252
02253
02254 llx=cpl_calloc(n_cubes,sizeof(int)) ;
02255 lly=cpl_calloc(n_cubes,sizeof(int)) ;
02256 sub_offsetx=cpl_calloc(n_cubes,sizeof(float)) ;
02257 sub_offsety=cpl_calloc(n_cubes,sizeof(float)) ;
02258
02259 for ( i = 0 ; i < n_cubes ; i++ )
02260 {
02261 llx[i] = llx0 - sinfo_new_nint(cumoffsetx[i]) ;
02262 sub_offsetx[i] = (float)sinfo_new_nint(cumoffsetx[i]) - cumoffsetx[i] ;
02263 lly[i] = lly0 - sinfo_new_nint(cumoffsety[i]) ;
02264 sub_offsety[i] = (float)sinfo_new_nint(cumoffsety[i]) - cumoffsety[i] ;
02265 }
02266
02267
02268
02269
02270
02271
02272
02273 if ( NULL == (kernel = sinfo_generate_interpolation_kernel( kernel_type )) )
02274 {
02275 sinfo_msg_warning ("could not generate desired interpolation kernel"
02276 " or no kernel_typ was given, the default kernel"
02277 " is used now!") ;
02278 }
02279
02280
02281 tmpcubes=(cpl_imagelist**)cpl_calloc(n_cubes,sizeof(cpl_imagelist*)) ;
02282
02283 for ( i = 0 ; i < n_cubes ; i++ )
02284 {
02285 tmpspace = cpl_calloc(ilx, ily*sizeof(pixelvalue)) ;
02286
02287
02288
02289 tmpcubes[i] = cpl_imagelist_new();
02290 for ( z = z_min, m=0 ; z< z_max ; z++, m++ )
02291 {
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301 t_img=sinfo_new_shift_image(cpl_imagelist_get(cubes[i],z),
02302 sub_offsetx[i], sub_offsety[i], kernel);
02303
02304
02305 if (t_img==NULL)
02306 {
02307 sinfo_msg_error ("could not shift image plane no %d"
02308 "in cube no %d!", z, i) ;
02309 cpl_imagelist_delete(mergedCube) ;
02310 cpl_imagelist_delete(mask) ;
02311 cpl_free(kernel) ;
02312 return -1 ;
02313 }
02314 cpl_imagelist_set(tmpcubes[i],t_img,m);
02315 m_img=cpl_image_new(mlx,mly,CPL_TYPE_FLOAT);
02316 cpl_imagelist_set(mask,m_img,z);
02317 }
02318 cpl_free(tmpspace);
02319 }
02320
02321
02322
02323
02324
02325
02326
02327
02328
02329 for ( z = z_min, m=0 ; z < z_max ; z++, m++ )
02330 {
02331
02332
02333 for ( y = 0 ; y < oly ; y++ )
02334 {
02335 for ( x = 0 ; x < olx ; x++ )
02336 {
02337
02338
02339 for ( i = 0 ; i < n_cubes ; i++ )
02340 {
02341
02342 i_img=cpl_imagelist_get(cubes[i],0);
02343 ilx=cpl_image_get_size_x(i_img);
02344 ily=cpl_image_get_size_y(i_img);
02345
02346
02347
02348
02349 if ( y >= lly[i] && y < lly[i]+ily &&
02350 x >= llx[i] && x < llx[i]+ilx )
02351 {
02352 posx = x - llx[i] ;
02353 posy = y - lly[i] ;
02354
02355
02356 t_img=cpl_imagelist_get(tmpcubes[i],m);
02357 ptdata=cpl_image_get_data_float(t_img);
02358 m_img=cpl_imagelist_get(mask,z);
02359 pmdata=cpl_image_get_data_float(m_img);
02360 mlx=cpl_image_get_size_x(m_img);
02361
02362 if (!isnan(ptdata[posx+posy*ilx]) &&
02363 ptdata[posx+posy*ilx] != 0.)
02364 {
02365 pmdata[x+y*mlx] += (float)exptimes[i] ;
02366 } else if (isnan(ptdata[posx+posy*ilx])) {
02367 sinfo_msg_debug("ptdata %d, %d, %d is NAN\t",x,y,z);
02368 } else if (ptdata[posx+posy*ilx] == 0.) {
02369 sinfo_msg_debug("ptdata %d, %d, %d is 0\t",x,y,z);
02370 }
02371
02372 } else {
02373 sinfo_msg_debug("point %d, %d, %d outside range\n",x,y,z);
02374 }
02375 }
02376 }
02377 }
02378 }
02379
02380
02381
02382 for ( z = z_min, m = 0 ; z < z_max ; z++, m++ )
02383 {
02384 o_img=cpl_imagelist_get(mergedCube,z);
02385 podata=cpl_image_get_data_float(o_img);
02386 m_img=cpl_imagelist_get(mask,z);
02387 pmdata=cpl_image_get_data_float(m_img);
02388 mlx=cpl_image_get_size_x(m_img);
02389
02390
02391 for ( y = 0 ; y < oly ; y++ )
02392 {
02393 for ( x = 0 ; x < olx ; x++ )
02394 {
02395
02396
02397
02398
02399 for ( i = 0 ; i < n_cubes ; i++ )
02400 {
02401
02402 if ( y >= lly[i] && y < lly[i]+ily &&
02403 x >= llx[i] && x < llx[i]+ilx )
02404 {
02405 posx = x - llx[i] ;
02406 posy = y - lly[i] ;
02407
02408
02409
02410
02411 t_img=cpl_imagelist_get(tmpcubes[i],m);
02412 ptdata=cpl_image_get_data_float(t_img);
02413
02414
02415
02416 if (!isnan(ptdata[posx+posy*ilx]))
02417 {
02418 if (pmdata[x+y*mlx] != 0.)
02419 {
02420
02421
02422 weight = exptimes[0] / pmdata[x+y*mlx] ;
02423 }
02424 else
02425 {
02426 weight = 0. ;
02427 }
02428 podata[x+y*olx] += weight*ptdata[posx+posy*ilx] ;
02429
02430 }
02431 }
02432 }
02433 }
02434 }
02435 }
02436
02437
02438
02439
02440 cpl_free(kernel) ;
02441 for( i = 0 ; i < n_cubes ; i++ )
02442 {
02443 cpl_imagelist_delete (tmpcubes[i]) ;
02444 }
02445
02446
02447 cpl_free(tmpcubes) ;
02448 cpl_free(llx) ;
02449 cpl_free(lly) ;
02450 cpl_free(sub_offsetx) ;
02451 cpl_free(sub_offsety) ;
02452
02453 return 0 ;
02454 }
02455
02456
02457
02458 int
02459 sinfo_new_combine_jittered_cubes_thomas_range(cpl_imagelist ** cubes,
02460 cpl_imagelist * mergedCube,
02461 cpl_imagelist * mask,
02462 int n_cubes,
02463 float * cumoffsetx,
02464 float * cumoffsety,
02465 double * exptimes,
02466 char * kernel_type,
02467 const int z_min,
02468 const int z_max,
02469 const double kappa )
02470 {
02471
02472 int i ;
02473 int x, y, z ;
02474 int llx0, lly0 ;
02475 int posx, posy ;
02476
02477
02478 double med=0;
02479 double sig=0;
02480 double avg=0;
02481 int nc=0;
02482
02483
02484 double * kernel ;
02485
02486 cpl_vector* val=NULL;
02487 cpl_vector* msk=NULL;
02488 int nclip=0;
02489 int ks=0;
02490 int ovr=0;
02491
02492 double msk_sum=0;
02493 double val_msk_sum=0;
02494
02495
02496
02497 int* llx=NULL;
02498 int* lly=NULL ;
02499 float* sub_offsetx=NULL ;
02500 float* sub_offsety=NULL ;
02501 cpl_imagelist ** tmpcubes=NULL ;
02502
02503 pixelvalue * tmpspace;
02504 const int z_siz=z_max-z_min;
02505 int m=0;
02506
02507
02508
02509
02510
02511
02512 int ilx=0;
02513 int ily=0;
02514 int olx=0;
02515 int oly=0;
02516 int mlx=0;
02517 int mly=0;
02518
02519
02520 int inp=0;
02521 int onp=0;
02522
02523
02524 float* podata=NULL;
02525 float* pmdata=NULL;
02526 float* ptdata=NULL;
02527 float* pvdata=NULL;
02528
02529 cpl_image* i_img=NULL;
02530 cpl_image* o_img=NULL;
02531 cpl_image* m_img=NULL;
02532 cpl_image* t_img=NULL;
02533 cpl_image* v_img=NULL;
02534
02535
02536
02537 if ( cubes == NULL )
02538 {
02539 sinfo_msg_error ("no cube list given!") ;
02540 return -1 ;
02541 }
02542 if ( n_cubes <= 0 )
02543 {
02544 sinfo_msg_error ("wrong number of data cubes in list!") ;
02545 return -1 ;
02546 }
02547 if ( cumoffsetx == NULL || cumoffsety == NULL )
02548 {
02549 sinfo_msg_error ("no cumoffsetx/y given!") ;
02550 return -1 ;
02551 }
02552 if ( exptimes == NULL )
02553 {
02554 sinfo_msg_error ("no exposure time array given!") ;
02555 return -1 ;
02556 }
02557 if (z_siz <= 0 )
02558 {
02559 sinfo_msg_error ("z_max <= z_min given!") ;
02560 return -1 ;
02561 }
02562
02563 i_img=cpl_imagelist_get(cubes[0],0);
02564 o_img=cpl_imagelist_get(mergedCube,0);
02565 ilx=cpl_image_get_size_x(i_img);
02566 ily=cpl_image_get_size_y(i_img);
02567 olx=cpl_image_get_size_x(o_img);
02568 oly=cpl_image_get_size_y(o_img);
02569 mlx=olx;
02570 mly=oly;
02571
02572
02573
02574
02575
02576
02577
02578
02579 llx0 = olx/2 - ilx/2 ;
02580 lly0 = oly/2 - ily/2 ;
02581
02582
02583
02584
02585
02586
02587
02588 llx=cpl_calloc(n_cubes,sizeof(int));
02589 lly=cpl_calloc(n_cubes,sizeof(int)) ;
02590 sub_offsetx=cpl_calloc(n_cubes,sizeof(float)) ;
02591 sub_offsety=cpl_calloc(n_cubes,sizeof(float)) ;
02592
02593
02594
02595 for ( i = 0 ; i < n_cubes ; i++ )
02596 {
02597 llx[i] = llx0 - sinfo_new_nint(cumoffsetx[i]) ;
02598 sub_offsetx[i] = (float)sinfo_new_nint(cumoffsetx[i]) - cumoffsetx[i] ;
02599 lly[i] = lly0 - sinfo_new_nint(cumoffsety[i]) ;
02600 sub_offsety[i] = (float)sinfo_new_nint(cumoffsety[i]) - cumoffsety[i] ;
02601 }
02602
02603
02604
02605
02606
02607
02608
02609 if ( NULL == (kernel = sinfo_generate_interpolation_kernel(kernel_type)) )
02610 {
02611 sinfo_msg_warning ("could not generate desired interpolation kernel"
02612 "or no kernel_typ was given, the default kernel"
02613 "is used now!") ;
02614 }
02615
02616
02617 tmpcubes=(cpl_imagelist**)cpl_calloc(n_cubes,sizeof(cpl_imagelist*)) ;
02618
02619 for ( i = 0 ; i < n_cubes ; i++ )
02620 {
02621 i_img=cpl_imagelist_get(cubes[i],0);
02622 ilx=cpl_image_get_size_x(i_img);
02623 ily=cpl_image_get_size_y(i_img);
02624 inp=cpl_imagelist_get_size(cubes[i]);
02625 tmpspace = cpl_calloc(ilx, ily*sizeof(pixelvalue)) ;
02626 tmpcubes[i]=cpl_imagelist_new();
02627
02628
02629
02630 for ( z = z_min, m=0 ; z < z_max ; z++, m++ )
02631 {
02632 t_img=sinfo_new_shift_image(cpl_imagelist_get(cubes[i],z),
02633 sub_offsetx[i],
02634 sub_offsety[i],
02635 kernel);
02636
02637 if (t_img==NULL)
02638 {
02639 sinfo_msg_error("could not shift image plane no %d "
02640 "in cube no %d!", z, i) ;
02641 cpl_imagelist_delete(mergedCube) ;
02642 cpl_imagelist_delete(mask) ;
02643 cpl_free(kernel) ;
02644 return -1 ;
02645 }
02646
02647 cpl_imagelist_set(tmpcubes[i],t_img,m);
02648 m_img=cpl_image_new(mlx,mly,CPL_TYPE_FLOAT);
02649 cpl_imagelist_set(mask,m_img,z);
02650 }
02651 cpl_free(tmpspace);
02652
02653 }
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664 o_img=cpl_imagelist_get(mergedCube,0);
02665 olx=cpl_image_get_size_x(o_img);
02666 oly=cpl_image_get_size_y(o_img);
02667 onp=cpl_imagelist_get_size(mergedCube);
02668
02669 for ( i = 0 ; i < n_cubes ; i++ )
02670 {
02671
02672
02673 i_img=cpl_imagelist_get(cubes[i],0);
02674 ilx=cpl_image_get_size_x(i_img);
02675 ily=cpl_image_get_size_y(i_img);
02676 inp=cpl_imagelist_get_size(cubes[i]);
02677
02678
02679 for ( y = 0 ; y < oly ; y++ )
02680 {
02681
02682 for ( x = 0 ; x < olx ; x++ )
02683 {
02684
02685
02686
02687 if ( y >= lly[i] && y < lly[i]+ily &&
02688 x >= llx[i] && x < llx[i]+ilx )
02689 {
02690 posx = x - llx[i] ;
02691 posy = y - lly[i] ;
02692
02693 for ( z = z_min,m=0 ; z < z_max ; z++,m++ )
02694 {
02695 t_img=cpl_imagelist_get(tmpcubes[i],m);
02696 ptdata=cpl_image_get_data_float(t_img);
02697 m_img=cpl_imagelist_get(mask,z);
02698 pmdata=cpl_image_get_data_float(m_img);
02699 mlx=cpl_image_get_size_x(m_img);
02700
02701 if (!isnan(ptdata[posx+posy*ilx]) &&
02702 ptdata[posx+posy*ilx] != 0.)
02703 {
02704 pmdata[x+y*mlx] += (float)exptimes[i] ;
02705 }
02706 }
02707 }
02708 }
02709 }
02710 }
02711
02712 m=0;
02713 for ( z = z_min; z < z_max ; z++ ) {
02714 m_img=cpl_imagelist_get(mask,z);
02715 pmdata=cpl_image_get_data_float(m_img);
02716 o_img=cpl_imagelist_get(mergedCube,z);
02717 podata=cpl_image_get_data_float(o_img);
02718 mlx=cpl_image_get_size_x(m_img);
02719 mly=cpl_image_get_size_y(m_img);
02720
02721
02722
02723
02724 for ( y = 0 ; y < oly ; y++ ) {
02725 for ( x = 0 ; x < olx ; x++ ) {
02726 avg=0;
02727 nc=0;
02728
02729 for ( i = 0 ; i < n_cubes ; i++ ) {
02730 t_img=cpl_imagelist_get(tmpcubes[i],m);
02731 ptdata=cpl_image_get_data_float(t_img);
02732 if ( y >= lly[i] && y < lly[i]+ily &&
02733 x >= llx[i] && x < llx[i]+ilx ) {
02734 posx = x - llx[i] ;
02735 posy = y - lly[i] ;
02736 if (!isnan(ptdata[posx+posy*ilx]) &&
02737 ptdata[posx+posy*ilx] != 0.) {
02738 nc++;
02739 }
02740 }
02741 }
02742 if( nc > 0 ) {
02743 msk=cpl_vector_new(n_cubes);
02744 for (i=0;i<n_cubes;i++) {
02745 cpl_vector_set(msk,i,1);
02746 }
02747
02748
02749 nclip=0;
02750
02751
02752 for (ks=0;ks<nc;ks++) {
02753 sig=0;
02754 med=0;
02755 ovr=0;
02756 if(nc-nclip >0) {
02757 val=cpl_vector_new(nc-nclip);
02758 }
02759
02760 for ( i = 0 ; i < n_cubes ; i++ ) {
02761 t_img=cpl_imagelist_get(tmpcubes[i],m);
02762 ptdata=cpl_image_get_data_float(t_img);
02763 if ( y >= lly[i] && y < lly[i]+ily &&
02764 x >= llx[i] && x < llx[i]+ilx ) {
02765 posx = x - llx[i] ;
02766 posy = y - lly[i] ;
02767 if (!isnan(ptdata[posx+posy*ilx]) &&
02768 ptdata[posx+posy*ilx] != 0. &&
02769 (cpl_vector_get(msk,i) != 0)) {
02770 cpl_vector_set(val,ovr,
02771 (double)ptdata[posx+posy*ilx]);
02772 ovr++;
02773 }
02774 }
02775 }
02776
02777
02778 if(ovr>0) {
02779 avg=cpl_vector_get_mean(val);
02780 med=cpl_vector_get_median(val);
02781 if(ovr>1) {
02782 sig=cpl_vector_get_stdev(val);
02783 } else {
02784 sig=0;
02785 }
02786 cpl_vector_delete(val);
02787 }
02788
02789 for ( i = 0 ; i < n_cubes ; i++ ) {
02790 t_img=cpl_imagelist_get(tmpcubes[i],m);
02791 ptdata=cpl_image_get_data_float(t_img);
02792
02793 if ( y >= lly[i] && y < lly[i]+ily &&
02794 x >= llx[i] && x < llx[i]+ilx ) {
02795 posx = x - llx[i] ;
02796 posy = y - lly[i] ;
02797 if (!isnan(ptdata[posx+posy*ilx]) &&
02798 ptdata[posx+posy*ilx] != 0. &&
02799 (cpl_vector_get(msk,i) != 0)) {
02800 if(abs((ptdata[posx+posy*ilx]-med))> kappa*sig) {
02801 ptdata[posx+posy*ilx]=0;
02802 pmdata[x+y*mlx] -= exptimes[i] ;
02803 cpl_vector_set(msk,i,0);
02804 nclip++;
02805 }
02806 }
02807 }
02808 }
02809 }
02810
02811 msk_sum=0;
02812 val_msk_sum=0;
02813 for ( i = 0 ; i < n_cubes ; i++ ) {
02814 v_img=cpl_imagelist_get(tmpcubes[i],m);
02815 pvdata=cpl_image_get_data_float(v_img);
02816
02817 if ( y >= lly[i] && y < lly[i]+ily &&
02818 x >= llx[i] && x < llx[i]+ilx ) {
02819 posx = x - llx[i] ;
02820 posy = y - lly[i] ;
02821 if (!isnan(pvdata[posx+posy*ilx]) &&
02822 pvdata[posx+posy*ilx] != 0. &&
02823 (cpl_vector_get(msk,i) != 0)) {
02824 msk_sum+=pmdata[x+y*mlx];
02825 val_msk_sum+=pvdata[posx+posy*ilx]*
02826 pmdata[x+y*mlx];
02827 }
02828 }
02829 }
02830 podata[x+y*olx]=val_msk_sum/msk_sum;
02831 cpl_vector_delete(msk);
02832 }
02833 }
02834 }
02835 m++;
02836 }
02837
02838
02839
02840
02841
02842
02843
02844
02845
02846
02847
02848
02849
02850
02851
02852
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864
02865
02866
02867
02868
02869
02870
02871
02872
02873
02874
02875
02876
02877
02878
02879
02880
02881
02882
02883
02884
02885
02886
02887
02888
02889
02890
02891
02892 cpl_free(kernel) ;
02893 for( i = 0 ; i < n_cubes ; i++ )
02894 {
02895 cpl_imagelist_delete (tmpcubes[i]) ;
02896 }
02897
02898 cpl_free(tmpcubes);
02899 cpl_free(llx);
02900 cpl_free(lly) ;
02901 cpl_free(sub_offsetx) ;
02902 cpl_free(sub_offsety) ;
02903
02904
02905 return 0 ;
02906 }
02907
02908
02909
02910
02911
02912
02913
02914
02915
02916
02917
02918
02919
02920 cpl_imagelist *
02921 sinfo_new_interpol_cube_simple( cpl_imagelist * cube,
02922 cpl_imagelist * badcube,
02923 int maxdist )
02924 {
02925 cpl_imagelist * intercube ;
02926 float* goodNeighbors=NULL ;
02927 int z, row, col ;
02928 int nx, ny, nz ;
02929 int llx, lly, llz ;
02930 int zi, coli, rowi ;
02931 int n ;
02932
02933
02934
02935
02936 int clx=0;
02937 int cly=0;
02938 int blx=0;
02939 int bly=0;
02940
02941 int cnp=0;
02942
02943
02944 float* pbdata=NULL;
02945 float* pidata=NULL;
02946 float* pbzidata=NULL;
02947 float* pczidata=NULL;
02948
02949 cpl_image* c_img=NULL;
02950 cpl_image* b_img=NULL;
02951 cpl_image* i_img=NULL;
02952
02953 cpl_image* bzi_img=NULL;
02954 cpl_image* czi_img=NULL;
02955
02956
02957
02958 if ( cube == NULL || badcube == NULL )
02959 {
02960 sinfo_msg_error("no cube given!") ;
02961 return NULL ;
02962 }
02963 if ( maxdist < 1 )
02964 {
02965 sinfo_msg_error("wrong maxrad given!") ;
02966 return NULL ;
02967 }
02968 intercube = cpl_imagelist_duplicate(cube) ;
02969
02970 goodNeighbors=cpl_calloc((2*maxdist+1)*(2*maxdist+1)*(2*maxdist+1) -1,
02971 sizeof(float)) ;
02972
02973 cnp=cpl_imagelist_get_size(cube);
02974 for ( z = 0 ; z < cnp ; z++ )
02975 {
02976 b_img=cpl_imagelist_get(badcube,z);
02977 i_img=cpl_imagelist_get(intercube,z);
02978 pbdata=cpl_image_get_data_float(b_img);
02979 pidata=cpl_image_get_data_float(i_img);
02980 blx=cpl_image_get_size_x(b_img);
02981 bly=cpl_image_get_size_y(b_img);
02982
02983 c_img=cpl_imagelist_get(cube,z);
02984 clx=cpl_image_get_size_x(c_img);
02985 cly=cpl_image_get_size_y(c_img);
02986
02987 for ( row = 0 ; row < cly ; row++ )
02988 {
02989 for ( col = 0 ; col < clx ; col++ )
02990 {
02991 if ( pbdata[col+row*clx] == 0 )
02992 {
02993
02994 llx = col - maxdist ;
02995 nx = 2*maxdist +1 ;
02996 if (llx < 0)
02997 {
02998 nx += llx ;
02999 llx = 0 ;
03000 }
03001 if ( llx + nx > clx )
03002 {
03003 nx -= (llx + nx - clx) ;
03004 }
03005
03006 lly = row - maxdist ;
03007 ny = 2*maxdist +1 ;
03008 if (lly < 0)
03009 {
03010 ny += lly ;
03011 lly = 0 ;
03012 }
03013 if ( lly + ny > cly )
03014 {
03015 ny -= (lly + ny - cly) ;
03016 }
03017
03018 llz = z - maxdist ;
03019 nz = 2*maxdist +1 ;
03020 if (llz < 0)
03021 {
03022 nz += llz ;
03023 llz = 0 ;
03024 }
03025 if ( llz + nz > cnp )
03026 {
03027 nz -= (llz + nz - cnp) ;
03028 }
03029 n = 0 ;
03030 for ( zi = llz ; zi < llz+nz ; zi++ )
03031 {
03032 bzi_img=cpl_imagelist_get(badcube,zi);
03033 czi_img=cpl_imagelist_get(cube,zi);
03034 pbzidata=cpl_image_get_data_float(bzi_img);
03035 pczidata=cpl_image_get_data_float(czi_img);
03036
03037 for ( rowi = lly ; rowi < lly+ny ; rowi++ )
03038 {
03039 for ( coli = llx ; coli < llx+nx ; coli++ )
03040 {
03041 if ( pbzidata[coli+rowi*blx] == 1 )
03042 {
03043 goodNeighbors[n] = pczidata[coli+rowi*clx] ;
03044 n++ ;
03045 }
03046 }
03047 }
03048 }
03049 if ( n > 0 )
03050 {
03051 pidata[col+row*clx]=sinfo_new_median(goodNeighbors,n);
03052 pbdata[col+row*clx]=1 ;
03053 }
03054 else
03055 {
03056 continue ;
03057 }
03058 }
03059 }
03060 }
03061 }
03062 cpl_free(goodNeighbors) ;
03063 return intercube ;
03064 }
03065
03066
03067
03068
03069
03070
03071
03072
03073
03074
03075
03076
03077
03078
03079
03080
03081
03082
03083
03084
03085
03086
03087
03088
03089
03090
03091
03092
03093
03094
03095
03096
03097
03098
03099
03100
03101
03102
03103
03104 cpl_imagelist *
03105 sinfo_new_combine_cubes ( cpl_imagelist ** cubes,
03106 cpl_imagelist * mergedCube,
03107 int n_cubes,
03108 float * cumoffsetx,
03109 float * cumoffsety,
03110 float factor,
03111 char * kernel_type )
03112 {
03113 int i=0 ;
03114 int x=0;
03115 int y=0;
03116 int z=0;
03117 int llx0=0;
03118 int lly0=0;
03119 int posx=0;
03120 int posy=0;
03121 cpl_imagelist * mask=NULL ;
03122 double * kernel=NULL ;
03123 cpl_image * shiftedImage=NULL ;
03124 int n=0;
03125 int ns=0;
03126 double sum=0;
03127 double sum2=0;
03128 double mean=0;
03129 double sigma=0;
03130
03131 cpl_imagelist ** tmpcubes=NULL ;
03132
03133 int* llx=NULL ;
03134 int* lly=NULL ;
03135
03136 float* sub_offsetx=NULL ;
03137 float* sub_offsety=NULL ;
03138 float* cubedata=NULL ;
03139
03140 int mlx=0;
03141 int mly=0;
03142 int clx=0;
03143 int cly=0;
03144 int mnp=0;
03145 int cnp=0;
03146
03147
03148 float* ptdata=NULL;
03149 float* podata=NULL;
03150 float* pmdata=NULL;
03151
03152 cpl_image* tmp_img=NULL;
03153 cpl_image* o_img=NULL;
03154 cpl_image* m_img=NULL;
03155 cpl_image* c_img=NULL;
03156 cpl_image* t_img=NULL;
03157
03158
03159
03160
03161
03162 if ( cubes == NULL )
03163 {
03164 sinfo_msg_error ("no cube list given!") ;
03165 return NULL ;
03166 }
03167
03168 if ( n_cubes <= 0 )
03169 {
03170 sinfo_msg_error ("wrong number of data cubes in list!") ;
03171 return NULL ;
03172 }
03173
03174 if ( cumoffsetx == NULL || cumoffsety == NULL )
03175 {
03176 sinfo_msg_error ("no cumoffsetx/y given!") ;
03177 return NULL ;
03178 }
03179
03180 if ( factor <= 0. )
03181 {
03182 sinfo_msg_error ("wrong factor given!") ;
03183 return NULL ;
03184 }
03185
03186 m_img=cpl_imagelist_get(mergedCube,0);
03187 mlx=cpl_image_get_size_x(m_img);
03188 mly=cpl_image_get_size_y(m_img);
03189 cnp=cpl_imagelist_get_size(cubes[0]);
03190 c_img=cpl_imagelist_get(cubes[0],0);
03191 clx=cpl_image_get_size_x(c_img);
03192 cly=cpl_image_get_size_y(c_img);
03193
03194
03195 tmpcubes=(cpl_imagelist**)cpl_calloc(n_cubes,sizeof(cpl_imagelist*)) ;
03196
03197
03198
03199
03200
03201
03202
03203
03204
03205
03206
03207
03208
03209
03210
03211
03212
03213
03214
03215
03216
03217 tmpcubes[0]=cpl_imagelist_duplicate(cubes[0]);
03218
03219
03220
03221
03222
03223
03224
03225 llx0 = mlx/2 - clx/2 ;
03226 lly0 = mly/2 - cly/2 ;
03227
03228
03229
03230
03231
03232
03233
03234
03235 llx=cpl_calloc(n_cubes,sizeof(int)) ;
03236 lly=cpl_calloc(n_cubes,sizeof(int)) ;
03237 sub_offsetx=cpl_calloc(n_cubes,sizeof(float)) ;
03238 sub_offsety=cpl_calloc(n_cubes,sizeof(float)) ;
03239
03240 for ( i = 0 ; i < n_cubes ; i++ )
03241 {
03242 llx[i] = llx0 - sinfo_new_nint(cumoffsetx[i]) ;
03243 sub_offsetx[i] = (float)sinfo_new_nint(cumoffsetx[i]) - cumoffsetx[i] ;
03244 lly[i] = lly0 - sinfo_new_nint(cumoffsety[i]) ;
03245 sub_offsety[i] = (float)sinfo_new_nint(cumoffsety[i]) - cumoffsety[i] ;
03246 }
03247
03248
03249
03250
03251
03252
03253 if ( NULL == (kernel = sinfo_generate_interpolation_kernel(kernel_type)) )
03254 {
03255 sinfo_msg_warning ("could not generate desired interpolation kernel"
03256 " or no kernel_typ was given, the default kernel"
03257 " is used now!") ;
03258 }
03259
03260 for ( i = 0 ; i < n_cubes ; i++ )
03261 {
03262
03263
03264 for ( z = 0 ; z < cnp ; z++ )
03265 {
03266 tmp_img=cpl_imagelist_get(cubes[i],z);
03267 if ( NULL == (shiftedImage = sinfo_new_shift_image(tmp_img,
03268 sub_offsetx[i],
03269 sub_offsety[i],
03270 kernel ) ) )
03271 {
03272 sinfo_msg_error ("could not shift image plane no %d "
03273 "in cube no %d!", z, i) ;
03274 cpl_imagelist_delete(mergedCube) ;
03275 cpl_imagelist_delete(mask) ;
03276 cpl_free(kernel) ;
03277 return NULL ;
03278 }
03279 cpl_imagelist_set(tmpcubes[i],shiftedImage,z);
03280 }
03281 }
03282
03283 cubedata=cpl_calloc(n_cubes,sizeof(float)) ;
03284
03285 for ( y = 0 ; y < mly ; y++ )
03286 {
03287 for ( x = 0 ; x < mlx ; x++ )
03288 {
03289 for ( z = 0 ; z < mnp ; z++ )
03290 {
03291 sum = 0. ;
03292 sum2 = 0. ;
03293 n = 0 ;
03294 for ( i = 0 ; i < n_cubes ; i++ )
03295 {
03296 c_img=cpl_imagelist_get(cubes[i],z);
03297
03298 clx=cpl_image_get_size_x(c_img);
03299 cly=cpl_image_get_size_y(c_img);
03300
03301 t_img=cpl_imagelist_get(tmpcubes[i],z);
03302 ptdata=cpl_image_get_data_float(t_img);
03303
03304 m_img=cpl_imagelist_get(mergedCube,z);
03305 pmdata=cpl_image_get_data_float(m_img);
03306 o_img=cpl_imagelist_get(mask,z);
03307 podata=cpl_image_get_data_float(o_img);
03308
03309
03310
03311
03312 if ( y >= lly[i] && y < lly[i]+cly &&
03313 x >= llx[i] && x < llx[i]+clx )
03314 {
03315 posx = x - llx[i] ;
03316 posy = y - lly[i] ;
03317 if (!isnan(ptdata[posx+posy*clx]))
03318 {
03319 sum += ptdata[posx+posy*clx] ;
03320 sum2 += (ptdata[posx+posy*clx] *
03321 ptdata[posx+posy*clx]) ;
03322 cubedata[n] = ptdata[posx+posy*clx] ;
03323 n++ ;
03324 }
03325 }
03326 }
03327
03328 if ( n == 0 )
03329 {
03330 mean = 0. ;
03331 sigma = 0. ;
03332 pmdata[x+y*mlx] = 0. ;
03333 podata[x+y*mlx] = 0 ;
03334 }
03335 else if ( n == 1 )
03336 {
03337 mean = sum ;
03338 sigma = 0. ;
03339 pmdata[x+y*mlx] = mean ;
03340 podata[x+y*mlx] = 1 ;
03341 }
03342 else
03343 {
03344 mean = sum/(double)n ;
03345 sigma = sqrt( (sum2 - sum*mean) / (double)(n - 1) ) ;
03346 ns = 0 ;
03347 for ( i = 0 ; i < n ; i++ )
03348 {
03349 if ( cubedata[i] > mean+factor*sigma ||
03350 cubedata[i] < mean-factor*sigma )
03351 {
03352 continue ;
03353 }
03354 else
03355 {
03356 pmdata[x+y*mlx] += cubedata[i] ;
03357 ns++ ;
03358 }
03359 }
03360 if ( ns == 0 )
03361 {
03362 pmdata[x+y*mlx] = 0. ;
03363 }
03364 else
03365 {
03366 pmdata[x+y*mlx] /= (float)ns ;
03367 }
03368 podata[x+y*mlx] = (float)ns ;
03369 }
03370 }
03371 }
03372 }
03373
03374 for( i = 0 ; i < n_cubes ; i++ )
03375 {
03376 cpl_imagelist_delete (tmpcubes[i]) ;
03377 }
03378 cpl_free(tmpcubes);
03379 cpl_free(llx);
03380 cpl_free(lly);
03381 cpl_free(sub_offsetx);
03382 cpl_free(sub_offsety);
03383 cpl_free(cubedata);
03384
03385
03386 sinfo_new_convert_0_to_ZERO_for_cubes(mergedCube) ;
03387 cpl_free(kernel) ;
03388 return mask ;
03389 }
03390
03391 cpl_imagelist *
03392 sinfo_new_bin_cube(cpl_imagelist *cu,
03393 int xscale,
03394 int yscale,
03395 int xmin,
03396 int xmax,
03397 int ymin,
03398 int ymax)
03399 {
03400 int i,j,k;
03401 cpl_imagelist * cube;
03402 int ilx=0;
03403 int ily=0;
03404 int olx=0;
03405 int oly=0;
03406 int inp=0;
03407
03408 float* pidata=NULL;
03409 float* podata=NULL;
03410 cpl_image* i_img=NULL;
03411 cpl_image* o_img=NULL;
03412
03413
03414
03415
03416
03417
03418
03419
03420
03421 inp=cpl_imagelist_get_size(cu);
03422 i_img=cpl_imagelist_get(cu,0);
03423 ilx=cpl_image_get_size_x(i_img);
03424 ily=cpl_image_get_size_y(i_img);
03425 olx=xmax-xmin+1;
03426 oly=ymax-ymin+1;
03427
03428
03429 cube=cpl_imagelist_new();
03430 for ( i = 0 ; i < inp ; i++ ) {
03431 o_img = cpl_image_new(olx,oly,CPL_TYPE_FLOAT);
03432 cpl_imagelist_set(cube,o_img,i);
03433 }
03434
03435
03436 for (i=0;i<inp;i++){
03437 i_img=cpl_imagelist_get(cu,i);
03438 pidata=cpl_image_get_data_float(i_img);
03439 o_img=cpl_imagelist_get(cube,i);
03440 podata=cpl_image_get_data_float(o_img);
03441 for (j=0 ; j < olx ; j++) {
03442 for (k=0 ; k< oly ; k++) {
03443 podata[j+k*olx]=pidata[((int) (j+xmin)/xscale)+
03444 ((int) (k+ymin)/yscale)*ilx]/
03445 (xscale*yscale);
03446 }
03447 }
03448 }
03449
03450 return cube;
03451 }
03452
03453
03454 cpl_imagelist *
03455 sinfo_new_scale_cube(cpl_imagelist *cu,
03456 float xscale,
03457 float yscale,
03458 char * kernel_type)
03459 {
03460 cpl_imagelist * cube ;
03461 int i, j, k, l ;
03462 int lx_out, ly_out ;
03463 double cur ;
03464 double * invert_transform ;
03465 double neighbors[16] ;
03466 double rsc[8],
03467 sumrs ;
03468 double param[6];
03469 double x, y ;
03470 int px, py ;
03471 int pos ;
03472 int tabx, taby ;
03473 double * kernel ;
03474 int leaps[16] ;
03475 int ilx=0;
03476 int ily=0;
03477 int tlx=0;
03478 int tly=0;
03479 int inp;
03480 float* podata=0;
03481 cpl_image* in_img=NULL;
03482 cpl_image* ou_img=NULL;
03483
03484
03485 if (cu == NULL)
03486 {
03487 sinfo_msg_error ("null cube") ;
03488 return NULL ;
03489 }
03490
03491 param[0]=xscale;
03492 param[1]=0;
03493 param[2]=0;
03494 param[3]=0;
03495 param[4]=yscale;
03496 param[5]=0;
03497
03498
03499 invert_transform = sinfo_invert_linear_transform(param) ;
03500 if (invert_transform == NULL) {
03501 sinfo_msg_error("cannot compute sinfo_invert transform: "
03502 "aborting warping") ;
03503 return NULL ;
03504 }
03505
03506
03507 kernel = sinfo_generate_interpolation_kernel(kernel_type) ;
03508 if (kernel == NULL) {
03509 sinfo_msg_error("cannot generate kernel: aborting resampling") ;
03510 return NULL ;
03511 }
03512
03513
03514
03515 ilx=cpl_image_get_size_x(cpl_imagelist_get(cu,0));
03516 ily=cpl_image_get_size_y(cpl_imagelist_get(cu,0));
03517 inp=cpl_imagelist_get_size(cu);
03518
03519 lx_out = (int) ilx*xscale ;
03520 ly_out = (int) ily*yscale ;
03521
03522 cube=cpl_imagelist_new();
03523 for ( l = 0 ; l < inp ; i++ ) {
03524 in_img = cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
03525 cpl_imagelist_set(cube,in_img,l);
03526 }
03527
03528
03529
03530
03531
03532
03533
03534
03535
03536 for (l=0;l<inp;l++){
03537 in_img=cpl_imagelist_get(cu,l);
03538 ou_img=cpl_imagelist_get(cube,l);
03539 tlx=cpl_image_get_size_x(in_img);
03540 tly=cpl_image_get_size_y(in_img);
03541 podata=cpl_image_get_data_float(ou_img);
03542
03543 leaps[0] = -1 - tlx ;
03544 leaps[1] = - tlx ;
03545 leaps[2] = 1 - tlx ;
03546 leaps[3] = 2 - tlx ;
03547
03548 leaps[4] = -1 ;
03549 leaps[5] = 0 ;
03550 leaps[6] = 1 ;
03551 leaps[7] = 2 ;
03552
03553 leaps[8] = -1 + tlx ;
03554 leaps[9] = tlx ;
03555 leaps[10]= 1 + tlx ;
03556 leaps[11]= 2 + tlx ;
03557
03558 leaps[12]= -1 + 2*tlx ;
03559 leaps[13]= 2*tlx ;
03560 leaps[14]= 1 + 2*tlx ;
03561 leaps[15]= 2 + 2*tlx ;
03562
03563
03564 for (j=0 ; j < ly_out ; j++) {
03565 for (i=0 ; i< lx_out ; i++) {
03566
03567
03568 x = invert_transform[0] * (double)i +
03569 invert_transform[1] * (double)j +
03570 invert_transform[2] ;
03571
03572 y = invert_transform[3] * (double)i +
03573 invert_transform[4] * (double)j +
03574 invert_transform[5] ;
03575
03576
03577 px = (int)x ;
03578 py = (int)y ;
03579
03580 if ((px < 1) ||
03581 (px > (tlx-2)) ||
03582 (py < 1) ||
03583 (py > (tly-2)))
03584 podata[i+j*lx_out] = (pixelvalue)0.0 ;
03585 else {
03586
03587 pos = px + py * tlx ;
03588 for (k=0 ; k<16 ; k++){
03589 if(!isnan(podata[(int)(pos+leaps[k])])) neighbors[k] =
03590 (double)(podata[(int)(pos+leaps[k])]) ;
03591 else neighbors[k]=0;
03592 }
03593
03594
03595 tabx = (x - (double)px) * (double)(TABSPERPIX) ;
03596 taby = (y - (double)py) * (double)(TABSPERPIX) ;
03597
03598
03599
03600
03601 rsc[0] = kernel[TABSPERPIX + tabx] ;
03602 rsc[1] = kernel[tabx] ;
03603 rsc[2] = kernel[TABSPERPIX - tabx] ;
03604 rsc[3] = kernel[2 * TABSPERPIX - tabx] ;
03605 rsc[4] = kernel[TABSPERPIX + taby] ;
03606 rsc[5] = kernel[taby] ;
03607 rsc[6] = kernel[TABSPERPIX - taby] ;
03608 rsc[7] = kernel[2 * TABSPERPIX - taby] ;
03609
03610 sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) *
03611 (rsc[4]+rsc[5]+rsc[6]+rsc[7]) ;
03612
03613
03614 cur = rsc[4] * ( rsc[0]*neighbors[0] +
03615 rsc[1]*neighbors[1] +
03616 rsc[2]*neighbors[2] +
03617 rsc[3]*neighbors[3] ) +
03618 rsc[5] * ( rsc[0]*neighbors[4] +
03619 rsc[1]*neighbors[5] +
03620 rsc[2]*neighbors[6] +
03621 rsc[3]*neighbors[7] ) +
03622 rsc[6] * ( rsc[0]*neighbors[8] +
03623 rsc[1]*neighbors[9] +
03624 rsc[2]*neighbors[10] +
03625 rsc[3]*neighbors[11] ) +
03626 rsc[7] * ( rsc[0]*neighbors[12] +
03627 rsc[1]*neighbors[13] +
03628 rsc[2]*neighbors[14] +
03629 rsc[3]*neighbors[15] ) ;
03630
03631
03632 podata[i+j*lx_out] = (pixelvalue)(cur/sumrs) ;
03633
03634 }
03635 }
03636 }
03637 }
03638 cpl_free(kernel) ;
03639 cpl_free(invert_transform) ;
03640 return cube ;
03641 }
03642
03643
03653 cpl_imagelist *
03654 sinfo_cube_zshift(const cpl_imagelist * cube_inp,
03655 const double shift,
03656 double* sub_shift)
03657 {
03658
03659 cpl_imagelist * cube_out=NULL ;
03660 cpl_image* img_inp=NULL;
03661 cpl_image* img_out=NULL;
03662 int col, row,z ;
03663 int int_shift ;
03664 int ilx=0;
03665 int ily=0;
03666 int ilz=0;
03667
03668 int olx=0;
03669 int oly=0;
03670 int olz=0;
03671 int i=0;
03672 float* pidata=NULL;
03673 float* podata=NULL;
03674
03675 cknull(cube_inp,"no input cube given!") ;
03676 check_nomsg(img_inp=cpl_imagelist_get(cube_inp,0));
03677 check_nomsg(ilx=cpl_image_get_size_x(img_inp));
03678 check_nomsg(ily=cpl_image_get_size_y(img_inp));
03679 check_nomsg(ilz=cpl_imagelist_get_size(cube_inp));
03680
03681 olx=ilx;
03682 oly=ily;
03683 olz=ilz;
03684
03685 int_shift = sinfo_new_nint(shift) ;
03686 *sub_shift = shift - (double) int_shift ;
03687 if ( int_shift == 0 )
03688 {
03689 cube_out =cpl_imagelist_duplicate(cube_inp) ;
03690 return cube_out ;
03691 }
03692 else
03693 {
03694
03695 cknull(cube_out = cpl_imagelist_new(),"could not allocate memory!") ;
03696 for ( i = 0 ; i < olz ; i++ ) {
03697 check_nomsg(img_out=cpl_image_new(olx,oly,CPL_TYPE_FLOAT));
03698 check_nomsg(cpl_imagelist_set(cube_out,img_out,i));
03699 }
03700 }
03701
03702 for(z=0; z< ilz; z++) {
03703 if ( (z-int_shift >= 0 ) && (z - int_shift < olz) ) {
03704 check_nomsg(img_inp=cpl_imagelist_get(cube_inp,z));
03705 check_nomsg(img_out=cpl_imagelist_get(cube_out,z-int_shift));
03706 check_nomsg(pidata=cpl_image_get_data_float(img_inp));
03707 check_nomsg(podata=cpl_image_get_data_float(img_out));
03708 for ( col = 0 ; col < ilx ; col++ ) {
03709 for ( row = 0 ; row < ily ; row++ ) {
03710 podata[col+row*olx] = pidata[col+row*olx] ;
03711 }
03712 }
03713 }
03714 }
03715 return cube_out ;
03716
03717 cleanup:
03718 sinfo_free_imagelist(&cube_out);
03719 return NULL ;
03720 }
03721
03731 cpl_imagelist *
03732 sinfo_cube_zshift_poly(const cpl_imagelist * cube_inp,
03733 const double sub_shift,
03734 const int order)
03735 {
03736 cpl_imagelist * cube_out ;
03737
03738 float* spec=NULL ;
03739 float* corrected_spec=NULL ;
03740 float* xnum=NULL ;
03741
03742 float sum=0;
03743 float new_sum=0 ;
03744 float eval=0 ;
03745 float * imageptr=NULL ;
03746 int row=0;
03747 int col=0 ;
03748 int firstpos=0 ;
03749 int n_points=0 ;
03750 int i=0 ;
03751 int flag=0;
03752 int ilx=0;
03753 int ily=0;
03754 int ilz=0;
03755
03756 int olx=0;
03757 int oly=0;
03758 int olz=0;
03759 int z=0;
03760
03761 float* pidata=NULL;
03762 float* podata=NULL;
03763 cpl_image* img_inp=NULL;
03764 cpl_image* img_out=NULL;
03765
03766 if ( cube_inp == NULL ) {
03767 sinfo_msg_error("no imagelist given!") ;
03768 return NULL ;
03769 }
03770
03771 img_inp=cpl_imagelist_get(cube_inp,0);
03772
03773 ilx=cpl_image_get_size_x(img_inp);
03774 ily=cpl_image_get_size_y(img_inp);
03775 ilz=cpl_imagelist_get_size(cube_inp);
03776
03777 if ( order <= 0 ) {
03778 sinfo_msg_error("wrong order of interpolation polynom given!") ;
03779 return NULL ;
03780 }
03781
03782
03783 olx=ilx;
03784 oly=ily;
03785 olz=ilz;
03786
03787
03788 if ( NULL == (cube_out = cpl_imagelist_new()) ) {
03789 sinfo_msg_error ("could not allocate memory!") ;
03790 return NULL ;
03791 } else {
03792 for ( i = 0 ; i < ilz ; i++ ) {
03793 img_out=cpl_image_new(olx,oly,CPL_TYPE_FLOAT);
03794 cpl_imagelist_set(cube_out,img_out,i);
03795 }
03796 }
03797
03798
03799 n_points = order + 1 ;
03800 if ( n_points % 2 == 0 ) {
03801 firstpos = (int)(n_points/2) - 1 ;
03802 } else {
03803 firstpos = (int)(n_points/2) ;
03804 }
03805
03806 spec=cpl_calloc(ilz,sizeof(float)) ;
03807 corrected_spec=cpl_calloc(ilz,sizeof(float)) ;
03808 xnum=cpl_calloc(order+1,sizeof(float)) ;
03809
03810
03811
03812 for ( i = 0 ; i < n_points ; i++ ) {
03813 xnum[i] = i ;
03814 }
03815
03816 for ( col = 0 ; col < ilx ; col++ ) {
03817 for ( row = 0 ; row < ily ; row++ ) {
03818 for( z=0; z< ilz; z++) {
03819 corrected_spec[z] = 0. ;
03820 }
03821 sum = 0. ;
03822 for ( z = 0 ; z < ilz ; z++ ) {
03823 img_inp=cpl_imagelist_get(cube_inp,z);
03824 pidata=cpl_image_get_data_float(img_inp);
03825 spec[z] = pidata[col + row*ilx] ;
03826 if (isnan(spec[z]) ) {
03827 spec[z] = 0. ;
03828
03829 for ( i = z - firstpos ; i < z-firstpos+n_points ; i++ ) {
03830 if ( i < 0 ) continue ;
03831 if ( i >= ilz) continue ;
03832 corrected_spec[i] = ZERO ;
03833 }
03834 }
03835 if ( z != 0 && z != ilz - 1 ) {
03836 sum += spec[z] ;
03837 }
03838
03839 }
03840
03841 new_sum = 0. ;
03842 for ( z = 0 ; z < ilz ; z++ ) {
03843
03844
03845
03846
03847
03848
03849
03850 if (isnan(corrected_spec[z])) continue ;
03851 if ( z - firstpos < 0 ) {
03852 imageptr = &spec[0] ;
03853 eval = sub_shift + z ;
03854 } else if ( z - firstpos + n_points >= ilz ) {
03855 imageptr = &spec[ilz - n_points] ;
03856 eval = sub_shift + z + n_points - ilz ;
03857 } else {
03858 imageptr = &spec[z-firstpos] ;
03859 eval = sub_shift + firstpos ;
03860 }
03861
03862 flag=0;
03863 corrected_spec[z]=sinfo_new_nev_ille(xnum,imageptr,order,eval,&flag);
03864 if ( z != 0 && z != ilz - 1 ) {
03865 new_sum += corrected_spec[z] ;
03866 }
03867 }
03868
03869
03870 for (z = 0 ; z < ilz ; z++ )
03871 {
03872 img_out=cpl_imagelist_get(cube_out,z);
03873 podata=cpl_image_get_data_float(img_out);
03874 if ( new_sum == 0. ) {
03875 new_sum = 1. ;
03876 }
03877 if ( z == 0 ) {
03878 podata[col+row*olx] = ZERO ;
03879 } else if ( z == ilz - 1 ) {
03880 podata[col+row*olx] = ZERO ;
03881 } else if ( isnan(corrected_spec[z]) ) {
03882 podata[col+row*olx] = ZERO ;
03883 } else {
03884 corrected_spec[z] *= sum / new_sum ;
03885 podata[col+row*olx] = corrected_spec[z] ;
03886 }
03887 }
03888
03889 }
03890 }
03891
03892 cpl_free(spec) ;
03893 cpl_free(corrected_spec) ;
03894 cpl_free(xnum) ;
03895 return cube_out ;
03896
03897
03898 }
03899
03908 cpl_imagelist *
03909 sinfo_cube_zshift_spline3(const cpl_imagelist * cube_inp,
03910 const double sub_shift)
03911 {
03912
03913 cpl_imagelist * cube_out=NULL ;
03914 float* spec=NULL ;
03915 float* corrected_spec=NULL ;
03916 float* xnum=NULL ;
03917 float* eval=NULL ;
03918 float sum=0;
03919 float new_sum=0 ;
03920 int row=0;
03921 int col=0;
03922 int i=0;
03923 int z=0;
03924
03925 int ilx=0;
03926 int ily=0;
03927 int ilz=0;
03928 int olx=0;
03929 int oly=0;
03930 int olz=0;
03931
03932 float* pidata=NULL;
03933 float* podata=NULL;
03934 cpl_image* img_inp=NULL;
03935 cpl_image* img_out=NULL;
03936
03937 if ( cube_inp == NULL ) {
03938 sinfo_msg_error("no imagelist given!") ;
03939 return NULL ;
03940 }
03941
03942 img_inp=cpl_imagelist_get(cube_inp,0);
03943 ilx=cpl_image_get_size_x(img_inp);
03944 ily=cpl_image_get_size_y(img_inp);
03945 ilz=cpl_imagelist_get_size(cube_inp);
03946
03947
03948 olx=ilx;
03949 oly=ily;
03950 olz=ilz;
03951
03952 if ( NULL == (cube_out = cpl_imagelist_new()) ) {
03953 sinfo_msg_error ("could not allocate memory!") ;
03954 return NULL ;
03955 } else {
03956 for ( i = 0 ; i < ilz ; i++ ) {
03957 img_out=cpl_image_new(olx,oly,CPL_TYPE_FLOAT);
03958 cpl_imagelist_set(cube_out,img_out,i);
03959 }
03960 }
03961
03962 xnum=cpl_calloc(ilz,sizeof(float)) ;
03963
03964 for ( i = 0 ; i < ilz ; i++ ) {
03965 xnum[i] = i ;
03966 }
03967
03968 spec=cpl_calloc(ilz,sizeof(float)) ;
03969 corrected_spec=cpl_calloc(ilz,sizeof(float)) ;
03970 eval=cpl_calloc(ilz,sizeof(float)) ;
03971
03972 for ( col = 0 ; col < ilx ; col++ ) {
03973 for ( row = 0 ; row < ily ; row++ ) {
03974 sum = 0. ;
03975 for ( z = 0 ; z < ilz ; z++ ) {
03976 img_inp=cpl_imagelist_get(cube_inp,z);
03977 pidata=cpl_image_get_data_float(img_inp);
03978 spec[z] = pidata[col + row*ilx] ;
03979 if (isnan(spec[z]) ) {
03980 for ( i = z-1 ; i <= z+1 ; i++ ) {
03981 if ( i < 0 ) continue ;
03982 if ( i >= ilz) continue ;
03983 corrected_spec[i] = ZERO ;
03984 }
03985 spec[z] = 0. ;
03986 }
03987 sum += spec[z] ;
03988 eval[z] = (float)sub_shift+(float)z ;
03989 }
03990
03991 if ( -1 == sinfo_function1d_natural_spline( xnum, spec, ilz, eval,
03992 corrected_spec, ilz ) )
03993 {
03994 sinfo_msg_error("error in spline interpolation!") ;
03995 return NULL ;
03996 }
03997
03998 new_sum = 0. ;
03999 for ( z = 0 ; z < ilz ; z++ ) {
04000 if ( isnan(corrected_spec[z]) ) {
04001 continue ;
04002 }
04003 new_sum += corrected_spec[z] ;
04004 }
04005
04006 for ( z = 0 ; z < ilz ; z++ ) {
04007 img_out=cpl_imagelist_get(cube_out,z);
04008 podata=cpl_image_get_data_float(img_out);
04009 if ( new_sum == 0. ) new_sum =1. ;
04010 {
04011 if ( isnan(corrected_spec[z]) ) {
04012 podata[col+row*olx] = ZERO ;
04013 } else {
04014 corrected_spec[z] *= sum / new_sum ;
04015 podata[col+row*olx] = corrected_spec[z] ;
04016 }
04017 }
04018 }
04019 }
04020 }
04021 cpl_free(xnum);
04022 cpl_free(spec) ;
04023 cpl_free(corrected_spec) ;
04024 cpl_free(eval) ;
04025
04026 return cube_out ;
04027 }