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 #ifdef HAVE_CONFIG_H
00201 # include <config.h>
00202 #endif
00203 #define POSIX_SOURCE 1
00204 #include "sinfo_vltPort.h"
00205
00206
00207
00208
00209
00210
00211
00212
00213 #include "sinfo_function_1d.h"
00214 #include "sinfo_cube_construct.h"
00215 #include "sinfo_spectrum_ops.h"
00216 #include "sinfo_wave_calibration.h"
00217 #include "sinfo_utilities.h"
00218 #include "sinfo_local_types.h"
00219 #include "sinfo_fft_base.h"
00228
00229
00230
00241 cpl_image *
00242 sinfo_new_convolve_ns_image_by_gauss( cpl_image * lineImage,
00243 int hw )
00244 {
00245 cpl_image * returnImage ;
00246 float* row_buffer=NULL ;
00247 float * filter ;
00248 int col, row ;
00249 int ilx=0;
00250 int ily=0;
00251
00252 float* pidata=NULL;
00253 float* podata=NULL;
00254
00255 if ( lineImage == NULL )
00256 {
00257 sinfo_msg_error("no input image given!\n") ;
00258 return NULL ;
00259 }
00260 ilx=cpl_image_get_size_x(lineImage);
00261 ily=cpl_image_get_size_y(lineImage);
00262 pidata=cpl_image_get_data_float(lineImage);
00263 if ( hw < 1 )
00264 {
00265 sinfo_msg_error(" wrong half width given!\n") ;
00266 return NULL ;
00267 }
00268
00269
00270 if ( NULL == ( returnImage = cpl_image_new(ilx,ily,CPL_TYPE_FLOAT ) ))
00271 {
00272 sinfo_msg_error("cannot allocate a new image\n");
00273 return NULL ;
00274 }
00275 podata=cpl_image_get_data_float(returnImage);
00276
00277
00278 row_buffer=cpl_calloc(ily,sizeof(float)) ;
00279
00280 for ( row = 0 ; row < ily ; row++ )
00281 {
00282 for ( col = 0 ; col < ilx ; col++ )
00283 {
00284 if ( isnan(pidata[col+row*ilx]) )
00285 {
00286 row_buffer[col] = 0. ;
00287 }
00288 else
00289 {
00290 row_buffer[col] = pidata[col + row*ilx] ;
00291 }
00292 }
00293
00294
00295
00296
00297
00298 filter = sinfo_function1d_filter_lowpass( row_buffer,
00299 ilx,
00300 LOW_PASS_GAUSSIAN,
00301 hw ) ;
00302 for ( col = 0 ; col < ily ; col++ )
00303 {
00304 podata[col + row*ilx] = filter[col] ;
00305 }
00306
00307 sinfo_function1d_del (filter) ;
00308 }
00309 cpl_free(row_buffer);
00310 return returnImage ;
00311 }
00312
00329 float *
00330 sinfo_north_south_test( cpl_image * ns_image,
00331 int n_slitlets,
00332 int halfWidth,
00333 float fwhm,
00334 float minDiff,
00335 float estimated_dist,
00336 float devtol,
00337 int bottom,
00338 int top )
00339 {
00340 int i, j, k, m, row, col, n, ni, na ;
00341 int position, counter, iters ;
00342 int xdim, ndat, its, numpar ;
00343 pixelvalue row_buf[cpl_image_get_size_x(ns_image)] ;
00344 float sum, mean, maxval ;
00345 float tol, lab ;
00346 float * distances ;
00347 float distances_buf[cpl_image_get_size_y(ns_image)][n_slitlets-1] ;
00348 float x_position[n_slitlets] ;
00349 float * xdat, * wdat ;
00350 int * mpar ;
00351 int found[3*n_slitlets], found_clean[3*n_slitlets] ;
00352 int found_cleanit[3*n_slitlets] ;
00353 Vector * line ;
00354 FitParams ** par ;
00355 int foundit, begin, end ;
00356 int zeroindicator ;
00357 int ilx=0;
00358 int ily=0;
00359
00360 float* pidata=NULL;
00361
00362 if ( ns_image == NULL )
00363 {
00364 sinfo_msg_error("sorry, no image given\n") ;
00365 return NULL ;
00366 }
00367 ilx=cpl_image_get_size_x(ns_image);
00368 ily=cpl_image_get_size_y(ns_image);
00369 pidata=cpl_image_get_data_float(ns_image);
00370
00371
00372 if ( n_slitlets < 1 )
00373 {
00374 sinfo_msg_error("wrong number of slitlets given\n") ;
00375 return NULL ;
00376 }
00377 if ( halfWidth < 0 || halfWidth >= estimated_dist )
00378 {
00379 sinfo_msg_error("wrong half width given\n") ;
00380 return NULL ;
00381 }
00382 if ( fwhm <= 0. )
00383 {
00384 sinfo_msg_error("wrong fwhm given\n") ;
00385 return NULL ;
00386 }
00387 if ( minDiff < 1. )
00388 {
00389 sinfo_msg_error("wrong minDiff given\n") ;
00390 return NULL ;
00391 }
00392
00393
00394 if (NULL == (distances = (float *) cpl_calloc ( n_slitlets - 1 ,
00395 sizeof (float) )))
00396 {
00397 sinfo_msg_error("could not allocate memory\n") ;
00398 return NULL ;
00399 }
00400
00401
00402 for ( row = bottom ; row < top ; row++ )
00403 {
00404 zeroindicator = 0 ;
00405
00406
00407 for ( i = 0 ; i < n_slitlets-1 ; i++ )
00408 {
00409 distances_buf[row][i] = ZERO ;
00410 }
00411
00412
00413 for ( col = 0 ; col < ilx ; col++ )
00414 {
00415 row_buf[col] = pidata[col + row*ilx] ;
00416 }
00417
00418
00419 sum = 0. ;
00420 n = 0 ;
00421 for ( i = 0 ; i < ilx ; i++ )
00422 {
00423 if ( isnan(row_buf[i]) )
00424 {
00425 continue ;
00426 }
00427 sum += row_buf[i] ;
00428 n++ ;
00429 }
00430 mean = sum / (float)n ;
00431
00432
00433
00434 n = 0 ;
00435 for ( i = 0 ; i < ilx ; i++ )
00436 {
00437 if (isnan(row_buf[i]))
00438 {
00439 continue ;
00440 }
00441 if ( row_buf[i] > sqrt(mean*mean*9) )
00442 {
00443 found[n] = i ;
00444 n++ ;
00445 }
00446 }
00447
00448 if ( n < n_slitlets )
00449 {
00450 sinfo_msg_warning("t1 wrong number of intensity columns found "
00451 "in row: %d, found number: %d, mean: %g",
00452 row, n, mean) ;
00453 continue ;
00454 }
00455 else
00456 {
00457
00458 na = 0 ;
00459 for ( i = 1 ; i < n ; i ++ )
00460 {
00461 if ( found[i] - found[i-1] < halfWidth )
00462 {
00463 begin = found[i] - halfWidth ;
00464 if ( begin < 0 )
00465 {
00466 begin = 0 ;
00467 }
00468 end = found[i] + halfWidth ;
00469 if ( end >= ilx )
00470 {
00471 end = ilx - 1 ;
00472 }
00473
00474
00475 maxval = -FLT_MAX ;
00476 foundit = 0 ;
00477 for ( j = begin ; j <= end ; j++ )
00478 {
00479
00480 if (isnan(row_buf[j]))
00481 {
00482 continue ;
00483 }
00484 if (row_buf[j] >= maxval )
00485 {
00486 maxval = row_buf[j] ;
00487 foundit = j ;
00488 }
00489 }
00490 if (maxval == -FLT_MAX)
00491 {
00492 continue ;
00493 }
00494 for ( k = 0 ; k < na ; k++ )
00495 {
00496 if ( found_cleanit[k] >= begin &&
00497 found_cleanit[k] < foundit )
00498 {
00499 na-- ;
00500 }
00501 }
00502 for ( k = 0 ; k < n ; k++ )
00503 {
00504 if ( found[k] == foundit)
00505 {
00506 if (na>0){
00507 if ( found_cleanit[na-1] != found[k] )
00508 {
00509 found_cleanit[na] = found[k] ;
00510 na++ ;
00511 }
00512 }
00513 else{
00514 found_cleanit[na] = found[k] ;
00515 na++ ;
00516 }
00517 }
00518 }
00519 }
00520 else
00521 {
00522 if ( i == 1 )
00523 {
00524 found_cleanit[na] = found[0] ;
00525 na++ ;
00526 found_cleanit[na] = found[1] ;
00527 na++ ;
00528 }
00529 else
00530 {
00531 if (na>0){
00532 if ( found_cleanit[na-1] != found[i-1])
00533 {
00534 found_cleanit[na] = found[i-1] ;
00535 na++ ;
00536 }
00537 if ( found_cleanit[na-1] != found[i])
00538 {
00539 found_cleanit[na] = found[i] ;
00540 na++ ;
00541 }
00542 }
00543 else
00544 {
00545 found_cleanit[na] = found[i] ;
00546 na++ ;
00547 }
00548 }
00549 }
00550 }
00551
00552 j = 1 ;
00553 for ( i = 1 ; i < na ; i++ )
00554 {
00555 if ( (float)(found_cleanit[i] - found_cleanit[i-1]) <
00556 (estimated_dist - devtol) ||
00557 (float)(found_cleanit[i] - found_cleanit[i-1]) >
00558 (estimated_dist + devtol) )
00559 {
00560 continue ;
00561 }
00562 else
00563 {
00564 found_clean[j-1] = found_cleanit[i-1] ;
00565 found_clean[j] = found_cleanit[i] ;
00566 j++ ;
00567 }
00568 }
00569 }
00570 if ( j > n_slitlets )
00571 {
00572
00573 ni = 1 ;
00574 for ( i = 1 ; i < j ; i++ )
00575 {
00576 if ( (float)(found_clean[i] - found_clean[i-1]) <
00577 (estimated_dist - devtol ) ||
00578 (float)(found_clean[i] - found_clean[i-1]) >
00579 (estimated_dist + devtol ) )
00580 {
00581 continue ;
00582 }
00583 else
00584 {
00585
00586 found_clean[ni-1] = found_clean[i-1] ;
00587 found_clean[ni] = found_clean[i] ;
00588 ni++ ;
00589 }
00590 }
00591 if ( ni != n_slitlets )
00592 {
00593 sinfo_msg_warning("t2 wrong number of intensity columns"
00594 " found in row: %d, found number: %d",
00595 row, ni) ;
00596 continue ;
00597 }
00598 else
00599 {
00600 j = ni ;
00601 }
00602 }
00603 else if ( j < n_slitlets )
00604 {
00605 cpl_msg_debug ("north_south_test3:",
00606 "t3 wrong number of intensity columns "
00607 "found in row: %d , found number: %d, mean: %g\n",
00608 row, j, mean) ;
00609 continue ;
00610 }
00611 counter = 0 ;
00612
00613 for ( i = 0 ; i < j ; i++ )
00614 {
00615
00616 if ( NULL == (line = sinfo_new_vector (2*halfWidth + 1)) )
00617 {
00618 sinfo_msg_error ("cannot allocate new Vector \n") ;
00619 cpl_free(distances) ;
00620 return NULL ;
00621 }
00622
00623
00624 xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00625 wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00626 mpar = (int *) cpl_calloc( MAXPAR, sizeof (int) ) ;
00627 par = sinfo_new_fit_params(1) ;
00628
00629 m = 0 ;
00630 for ( k = found_clean[i]-halfWidth ;
00631 k <= found_clean[i]+halfWidth ; k++ )
00632 {
00633 if ( k < 0 )
00634 {
00635 k = 0. ;
00636 }
00637 else if ( k >= ilx )
00638 {
00639 k = ilx - 1 ;
00640 }
00641 else if ( isnan(row_buf[k]) )
00642 {
00643 zeroindicator = 1 ;
00644 break ;
00645 }
00646 else
00647 {
00648 line -> data[m] = row_buf[k] ;
00649 m++ ;
00650 }
00651 }
00652 if ( zeroindicator == 1 )
00653 {
00654 sinfo_new_destroy_vector(line) ;
00655 cpl_free(xdat) ;
00656 cpl_free(wdat) ;
00657 cpl_free(mpar) ;
00658 sinfo_new_destroy_fit_params(&par) ;
00659 break ;
00660 }
00661
00662
00663
00664
00665
00666 maxval = -FLT_MAX ;
00667 position = -INT32_MAX ;
00668 for ( k = 0 ; k < m ; k++ )
00669 {
00670 xdat[k] = k ;
00671 wdat[k] = 1.0 ;
00672 if ( line -> data[k] >= maxval )
00673 {
00674 maxval = line -> data[k] ;
00675 position = k ;
00676 }
00677 }
00678
00679
00680 xdim = XDIM ;
00681 ndat = line -> n_elements ;
00682 numpar = MAXPAR ;
00683 tol = TOL ;
00684 lab = LAB ;
00685 its = ITS ;
00686 (*par) -> fit_par[1] = fwhm ;
00687 (*par) -> fit_par[2] = (float) position ;
00688 (*par) -> fit_par[3] = (float) (line -> data[0] +
00689 line -> data[line->n_elements - 1]) / 2.0 ;
00690 (*par) -> fit_par[0] = maxval - ((*par) -> fit_par[3]) ;
00691
00692
00693
00694 if ( (*par) -> fit_par[0] < minDiff )
00695 {
00696 sinfo_msg_warning ("sorry, signal of line too low to fit "
00697 "in row: %d in slitlet %d\n", row, i) ;
00698 sinfo_new_destroy_vector(line) ;
00699 cpl_free(xdat) ;
00700 cpl_free(wdat) ;
00701 cpl_free(mpar) ;
00702 sinfo_new_destroy_fit_params(&par) ;
00703 continue ;
00704 }
00705
00706 for ( k = 0 ; k < MAXPAR ; k++ )
00707 {
00708 (*par) -> derv_par[k] = 0.0 ;
00709 mpar[k] = 1 ;
00710 }
00711
00712 if ( 0 > ( iters = sinfo_new_lsqfit_c( xdat, &xdim,
00713 line -> data, wdat, &ndat,
00714 (*par) -> fit_par,
00715 (*par) -> derv_par, mpar,
00716 &numpar, &tol, &its, &lab)) )
00717 {
00718
00719
00720
00721
00722
00723
00724 sinfo_new_destroy_vector(line) ;
00725 cpl_free(xdat) ;
00726 cpl_free(wdat) ;
00727 cpl_free(mpar) ;
00728 sinfo_new_destroy_fit_params(&par) ;
00729 continue ;
00730 }
00731
00732
00733 if ( (*par) -> fit_par[0] <= 0. ||
00734 (*par) -> fit_par[1] <= 0. ||
00735 (*par) -> fit_par[2] < 0. )
00736 {
00737 sinfo_msg_warning ("negative parameters as fit result, "
00738 "not used! in row %d in slitlet %d",
00739 row, i) ;
00740 sinfo_new_destroy_vector(line) ;
00741 cpl_free(xdat) ;
00742 cpl_free(wdat) ;
00743 cpl_free(mpar) ;
00744 sinfo_new_destroy_fit_params(&par) ;
00745 continue ;
00746 }
00747
00748
00749
00750 (*par) -> fit_par[2] = (float) (found_clean[i] - halfWidth) +
00751 (*par) -> fit_par[2] ;
00752 x_position[counter] = (*par) -> fit_par[2] ;
00753 counter ++ ;
00754
00755
00756 sinfo_new_destroy_fit_params(&par) ;
00757 sinfo_new_destroy_vector ( line ) ;
00758 cpl_free ( xdat ) ;
00759 cpl_free ( wdat ) ;
00760 cpl_free ( mpar ) ;
00761 }
00762 if (zeroindicator == 1)
00763 {
00764 sinfo_msg_debug ("bad pixel in fitting box in row: %d\n", row) ;
00765 continue ;
00766 }
00767
00768 if ( counter != n_slitlets )
00769 {
00770 continue ;
00771 sinfo_msg_warning("wrong number of slitlets found in row: %d",row);
00772 }
00773
00774 for ( i = 1 ; i < n_slitlets ; i++ )
00775 {
00776 distances_buf[row][i-1] = x_position[i] - x_position[i-1] ;
00777 }
00778 }
00779
00780
00781
00782
00783
00784 for ( i = 0 ; i < n_slitlets-1 ; i++ )
00785 {
00786 n = 0 ;
00787 sum = 0. ;
00788 for ( row = bottom ; row < top ; row++ )
00789 {
00790 if ( fabs( distances_buf[row][i] - estimated_dist ) > devtol ||
00791 isnan(distances_buf[row][i]) )
00792 {
00793
00794
00795
00796
00797
00798
00799 continue ;
00800 }
00801 sum += distances_buf[row][i] ;
00802 n++ ;
00803 }
00804 if ( n < 2 )
00805 {
00806 sinfo_msg_error("distances array could not be determined "
00807 "completely!, deviations of distances from number "
00808 "of slitlets too big\n" ) ;
00809 cpl_free(distances) ;
00810 return NULL ;
00811 }
00812 else
00813 {
00814 distances[i] = sum / (float)n ;
00815 }
00816 }
00817 return distances ;
00818 }
00819
00843 cpl_imagelist *
00844 sinfo_new_make_cube ( cpl_image * calibImage,
00845 float * distances,
00846 float * correct_diff_dist )
00847 {
00848 cpl_imagelist * returnCube ;
00849 int imsize, kslit, kpix ;
00850 int slit_index ;
00851 int z, col, recol ;
00852 int ilx=0;
00853 int ily=0;
00854
00855 float* podata=NULL;
00856 float* pidata=NULL;
00857 cpl_image* o_img;
00858
00859 if ( NULL == calibImage )
00860 {
00861 sinfo_msg_error("no resampled image given!\n") ;
00862 return NULL ;
00863 }
00864 ilx=cpl_image_get_size_x(calibImage);
00865 ily=cpl_image_get_size_y(calibImage);
00866 pidata=cpl_image_get_data_float(calibImage);
00867
00868 if ( NULL == distances )
00869 {
00870 sinfo_msg_error("no distances array from ns_test given!/n") ;
00871 return NULL ;
00872 }
00873
00874 if ( NULL == correct_diff_dist )
00875 {
00876 sinfo_msg_error("correct_diff_dist array is not allocated!/n") ;
00877 return NULL ;
00878 }
00879
00880 if ( N_SLITLETS != 32 )
00881 {
00882 sinfo_msg_error ("wrong number of slitlets given \n" ) ;
00883 return NULL ;
00884 }
00885 imsize = ilx / N_SLITLETS ;
00886
00887
00888 if ( NULL == (returnCube = cpl_imagelist_new()) )
00889 {
00890 sinfo_msg_error ("cannot allocate new cube \n" ) ;
00891 return NULL ;
00892 }
00893
00894
00895 for ( z = 0 ; z < ily ; z++ )
00896 {
00897
00898 o_img=cpl_image_new(imsize,N_SLITLETS,CPL_TYPE_FLOAT);
00899 podata=cpl_image_get_data_float(o_img);
00900 kpix = 0 ;
00901 kslit = 0 ;
00902 slit_index = -1 ;
00903 recol = -1 ;
00904 for ( col = 0 ; col < ilx ; col++ )
00905 {
00906 if ( col % imsize == 0 )
00907 {
00908 recol = 0 ;
00909 kslit = col/imsize ;
00910
00911 switch (kslit)
00912 {
00913 case 0:
00914 slit_index = 8 ;
00915 break ;
00916 case 1:
00917 slit_index = 7 ;
00918 break ;
00919 case 2:
00920 slit_index = 9 ;
00921 break ;
00922 case 3:
00923 slit_index = 6 ;
00924 break ;
00925 case 4:
00926 slit_index = 10 ;
00927 break ;
00928 case 5:
00929 slit_index = 5 ;
00930 break ;
00931 case 6:
00932 slit_index = 11 ;
00933 break ;
00934 case 7:
00935 slit_index = 4 ;
00936 break ;
00937 case 8:
00938 slit_index = 12 ;
00939 break ;
00940 case 9:
00941 slit_index = 3 ;
00942 break ;
00943 case 10:
00944 slit_index = 13 ;
00945 break ;
00946 case 11:
00947 slit_index = 2 ;
00948 break ;
00949 case 12:
00950 slit_index = 14 ;
00951 break ;
00952 case 13:
00953 slit_index = 1 ;
00954 break ;
00955 case 14:
00956 slit_index = 15 ;
00957 break ;
00958 case 15:
00959 slit_index = 0 ;
00960 break ;
00961 case 16:
00962 slit_index = 31 ;
00963 break ;
00964 case 17:
00965 slit_index = 16 ;
00966 break ;
00967 case 18:
00968 slit_index = 30 ;
00969 break ;
00970 case 19:
00971 slit_index = 17 ;
00972 break ;
00973 case 20:
00974 slit_index = 29 ;
00975 break ;
00976 case 21:
00977 slit_index = 18 ;
00978 break ;
00979 case 22:
00980 slit_index = 28 ;
00981 break ;
00982 case 23:
00983 slit_index = 19 ;
00984 break ;
00985 case 24:
00986 slit_index = 27 ;
00987 break ;
00988 case 25:
00989 slit_index = 20 ;
00990 break ;
00991 case 26:
00992 slit_index = 26 ;
00993 break ;
00994 case 27:
00995 slit_index = 21 ;
00996 break ;
00997 case 28:
00998 slit_index = 25 ;
00999 break ;
01000 case 29:
01001 slit_index = 22 ;
01002 break ;
01003 case 30:
01004 slit_index = 24 ;
01005 break ;
01006 case 31:
01007 slit_index = 23 ;
01008 break ;
01009 default:
01010 sinfo_msg_error("wrong slitlet index: couldn't be a "
01011 "spiffi image, there must be 32 slitlets!") ;
01012 cpl_imagelist_delete(returnCube) ;
01013 return NULL ;
01014 break ;
01015 }
01016
01017 if ( kslit != 0 )
01018 {
01019
01020
01021
01022
01023
01024 kpix = sinfo_new_nint(distances[kslit-1]) ;
01025
01026
01027
01028
01029
01030
01031 correct_diff_dist[slit_index] = distances[kslit-1] -
01032 (float)kpix ;
01033 }
01034
01035 else
01036 {
01037 correct_diff_dist[slit_index] = 0. ;
01038 }
01039 }
01040
01041
01042 podata[recol+slit_index*imsize] = pidata[col+kpix+z*ilx];
01043 recol++ ;
01044
01045 if ( recol > imsize )
01046 {
01047 sinfo_msg_error("wrong column of reconstructed "
01048 "image, shouldn't happen!\n") ;
01049 cpl_imagelist_delete(returnCube) ;
01050 return NULL ;
01051 }
01052 }
01053 }
01054 return returnCube ;
01055 }
01072 cpl_imagelist *
01073 sinfo_new_make_cube_spi ( cpl_image * calibImage,
01074 float ** slit_edges,
01075 float * shift )
01076 {
01077 cpl_imagelist * returnCube ;
01078 float diff, start ;
01079 float * center ;
01080 int * row_index ;
01081 int slit ;
01082 int col, z ;
01083 int imsize ;
01084 int * beginCol ;
01085 int col_counter ;
01086 int ilx=0;
01087 int ily=0;
01088
01089 float* podata=NULL;
01090 float* pidata=NULL;
01091 cpl_image* o_img;
01092
01093
01094 if ( NULL == calibImage )
01095 {
01096 sinfo_msg_error("no resampled image given!\n") ;
01097 return NULL ;
01098 }
01099 ilx=cpl_image_get_size_x(calibImage);
01100 ily=cpl_image_get_size_y(calibImage);
01101 pidata=cpl_image_get_data_float(calibImage);
01102
01103 if ( NULL == slit_edges )
01104 {
01105 sinfo_msg_error("no slit_edges array given from sinfo_fitSlits()!/n") ;
01106 return NULL ;
01107 }
01108
01109 if ( N_SLITLETS != 32 )
01110 {
01111 sinfo_msg_error ("wrong number of slitlets given \n" ) ;
01112 return NULL ;
01113 }
01114 imsize = ilx / N_SLITLETS ;
01115
01116
01117 if ( NULL == (row_index = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01118 {
01119 sinfo_msg_error ("cannot allocate memory \n" ) ;
01120 return NULL ;
01121 }
01122 if ( NULL == (beginCol = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01123 {
01124 sinfo_msg_error ("cannot allocate memory \n" ) ;
01125 cpl_free(row_index) ;
01126 return NULL ;
01127 }
01128 if ( NULL == (center = (float*) cpl_calloc(N_SLITLETS, sizeof(float)) ) )
01129 {
01130 sinfo_msg_error ("cannot allocate memory \n" ) ;
01131 cpl_free (row_index) ;
01132 cpl_free (beginCol) ;
01133 return NULL ;
01134 }
01135 if ( NULL == (returnCube = cpl_imagelist_new()) )
01136 {
01137 sinfo_msg_error ("cannot allocate new cube \n" ) ;
01138 cpl_free (row_index) ;
01139 cpl_free (beginCol) ;
01140 cpl_free (center) ;
01141 return NULL ;
01142 }
01143
01144
01145 for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01146
01147 {
01148 center[slit] = (slit_edges[slit][1] + slit_edges[slit][0]) / 2. ;
01149
01150
01151
01152
01153
01154 switch (slit)
01155 {
01156 case 0:
01157 row_index[0] = 8 ;
01158 break ;
01159 case 1:
01160 row_index[1] = 7 ;
01161 break ;
01162 case 2:
01163 row_index[2] = 9 ;
01164 break ;
01165 case 3:
01166 row_index[3] = 6 ;
01167 break ;
01168 case 4:
01169 row_index[4] = 10 ;
01170 break ;
01171 case 5:
01172 row_index[5] = 5 ;
01173 break ;
01174 case 6:
01175 row_index[6] = 11 ;
01176 break ;
01177 case 7:
01178 row_index[7] = 4 ;
01179 break ;
01180 case 8:
01181 row_index[8] = 12 ;
01182 break ;
01183 case 9:
01184 row_index[9] = 3 ;
01185 break ;
01186 case 10:
01187 row_index[10] = 13 ;
01188 break ;
01189 case 11:
01190 row_index[11] = 2 ;
01191 break ;
01192 case 12:
01193 row_index[12] = 14 ;
01194 break ;
01195 case 13:
01196 row_index[13] = 1 ;
01197 break ;
01198 case 14:
01199 row_index[14] = 15 ;
01200 break ;
01201 case 15:
01202 row_index[15] = 0 ;
01203 break ;
01204 case 16:
01205 row_index[16] = 31 ;
01206 break ;
01207 case 17:
01208 row_index[17] = 16 ;
01209 break ;
01210 case 18:
01211 row_index[18] = 30 ;
01212 break ;
01213 case 19:
01214 row_index[19] = 17 ;
01215 break ;
01216 case 20:
01217 row_index[20] = 29 ;
01218 break ;
01219 case 21:
01220 row_index[21] = 18 ;
01221 break ;
01222 case 22:
01223 row_index[22] = 28 ;
01224 break ;
01225 case 23:
01226 row_index[23] = 19 ;
01227 break ;
01228 case 24:
01229 row_index[24] = 27 ;
01230 break ;
01231 case 25:
01232 row_index[25] = 20 ;
01233 break ;
01234 case 26:
01235 row_index[26] = 26 ;
01236 break ;
01237 case 27:
01238 row_index[27] = 21 ;
01239 break ;
01240 case 28:
01241 row_index[28] = 25 ;
01242 break ;
01243 case 29:
01244 row_index[29] = 22 ;
01245 break ;
01246 case 30:
01247 row_index[30] = 24 ;
01248 break ;
01249 case 31:
01250 row_index[31] = 23 ;
01251 break ;
01252 default:
01253 sinfo_msg_error("wrong slitlet index: couldn't be a spiffi "
01254 "image, there must be 32 slitlets!\n") ;
01255 cpl_imagelist_delete(returnCube) ;
01256 cpl_free (row_index) ;
01257 cpl_free (beginCol) ;
01258 cpl_free (center) ;
01259 return NULL ;
01260 }
01261
01262
01263 start = center[slit] - (float) (imsize - 1)/2. ;
01264 beginCol[slit] = sinfo_new_nint (start) ;
01265
01266 diff = start - (float)beginCol[slit] ;
01267
01268
01269
01270
01271
01272
01273 shift[row_index[slit]] = diff ;
01274 }
01275
01276
01277 for ( z = 0 ; z < ily ; z++ )
01278 {
01279 o_img=cpl_image_new(imsize,N_SLITLETS,CPL_TYPE_FLOAT);
01280 podata=cpl_image_get_data_float(o_img);
01281 for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01282 {
01283 col_counter = beginCol[slit] ;
01284
01285 for ( col = 0 ; col < imsize ; col++ )
01286 {
01287 if ( col_counter > ilx-1 )
01288 {
01289 col_counter-- ;
01290 }
01291 if ( col_counter + z*ilx < 0 )
01292 {
01293 podata[col+row_index[slit]*imsize] = pidata[0] ;
01294 }
01295 else
01296 {
01297 podata[col+row_index[slit]*imsize]=pidata[col_counter+z*ilx];
01298 }
01299
01300 col_counter++ ;
01301 }
01302 }
01303 cpl_imagelist_set(returnCube,o_img,z);
01304 }
01305 cpl_free (row_index) ;
01306 cpl_free (beginCol) ;
01307 cpl_free (center) ;
01308
01309 return returnCube ;
01310 }
01337 cpl_imagelist *
01338 sinfo_new_make_cube_dist ( cpl_image * calibImage,
01339 float firstCol,
01340 float * distances,
01341 float * shift )
01342 {
01343 cpl_imagelist * returnCube ;
01344 float di ;
01345 float diff, start ;
01346 int * row_index ;
01347 int slit ;
01348 int col, z ;
01349 int imsize ;
01350 int * beginCol ;
01351 int col_counter ;
01352 int ilx=0;
01353 int ily=0;
01354
01355 float* podata=NULL;
01356 float* pidata=NULL;
01357 cpl_image* o_img;
01358
01359 if ( NULL == calibImage )
01360 {
01361 sinfo_msg_error(" no resampled image given!\n") ;
01362 return NULL ;
01363 }
01364 ilx=cpl_image_get_size_x(calibImage);
01365 ily=cpl_image_get_size_y(calibImage);
01366 pidata=cpl_image_get_data_float(calibImage);
01367
01368 if ( NULL == distances )
01369 {
01370 sinfo_msg_error("no distances array given from north_south_test()!") ;
01371 return NULL ;
01372 }
01373
01374 if ( N_SLITLETS != 32 )
01375 {
01376 sinfo_msg_error ("wrong number of slitlets given \n" ) ;
01377 return NULL ;
01378 }
01379 imsize = ilx / N_SLITLETS ;
01380
01381
01382 if ( NULL == (row_index = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01383 {
01384 sinfo_msg_error ("cannot allocate memory \n" ) ;
01385 return NULL ;
01386 }
01387 if ( NULL == (beginCol = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01388 {
01389 sinfo_msg_error ("cannot allocate memory \n" ) ;
01390 cpl_free(row_index) ;
01391 return NULL ;
01392 }
01393 if ( NULL == (returnCube = cpl_imagelist_new()) )
01394 {
01395 sinfo_msg_error ("cannot allocate new cube \n" ) ;
01396 cpl_free(row_index) ;
01397 cpl_free(beginCol) ;
01398 return NULL ;
01399 }
01400
01401 di = 0. ;
01402
01403
01404 for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01405
01406 {
01407
01408
01409
01410
01411
01412
01413 switch (slit)
01414 {
01415 case 0:
01416 row_index[0] = 8 ;
01417 break ;
01418 case 1:
01419 row_index[1] = 7 ;
01420 break ;
01421 case 2:
01422 row_index[2] = 9 ;
01423 break ;
01424 case 3:
01425 row_index[3] = 6 ;
01426 break ;
01427 case 4:
01428 row_index[4] = 10 ;
01429 break ;
01430 case 5:
01431 row_index[5] = 5 ;
01432 break ;
01433 case 6:
01434 row_index[6] = 11 ;
01435 break ;
01436 case 7:
01437 row_index[7] = 4 ;
01438 break ;
01439 case 8:
01440 row_index[8] = 12 ;
01441 break ;
01442 case 9:
01443 row_index[9] = 3 ;
01444 break ;
01445 case 10:
01446 row_index[10] = 13 ;
01447 break ;
01448 case 11:
01449 row_index[11] = 2 ;
01450 break ;
01451 case 12:
01452 row_index[12] = 14 ;
01453 break ;
01454 case 13:
01455 row_index[13] = 1 ;
01456 break ;
01457 case 14:
01458 row_index[14] = 15 ;
01459 break ;
01460 case 15:
01461 row_index[15] = 0 ;
01462 break ;
01463 case 16:
01464 row_index[16] = 31 ;
01465 break ;
01466 case 17:
01467 row_index[17] = 16 ;
01468 break ;
01469 case 18:
01470 row_index[18] = 30 ;
01471 break ;
01472 case 19:
01473 row_index[19] = 17 ;
01474 break ;
01475 case 20:
01476 row_index[20] = 29 ;
01477 break ;
01478 case 21:
01479 row_index[21] = 18 ;
01480 break ;
01481 case 22:
01482 row_index[22] = 28 ;
01483 break ;
01484 case 23:
01485 row_index[23] = 19 ;
01486 break ;
01487 case 24:
01488 row_index[24] = 27 ;
01489 break ;
01490 case 25:
01491 row_index[25] = 20 ;
01492 break ;
01493 case 26:
01494 row_index[26] = 26 ;
01495 break ;
01496 case 27:
01497 row_index[27] = 21 ;
01498 break ;
01499 case 28:
01500 row_index[28] = 25 ;
01501 break ;
01502 case 29:
01503 row_index[29] = 22 ;
01504 break ;
01505 case 30:
01506 row_index[30] = 24 ;
01507 break ;
01508 case 31:
01509 row_index[31] = 23 ;
01510 break ;
01511 default:
01512 sinfo_msg_error("wrong slitlet index: couldn't be a "
01513 "spiffi image, there must be 32 slitlets!") ;
01514 cpl_imagelist_delete(returnCube) ;
01515 cpl_free(row_index) ;
01516 cpl_free(beginCol) ;
01517 return NULL ;
01518 }
01519
01520
01521 if ( slit == 0 )
01522 {
01523 start = firstCol ;
01524 }
01525 else
01526 {
01527 di += distances[slit-1] ;
01528 start = firstCol + di ;
01529 }
01530 beginCol[slit] = sinfo_new_nint(start) ;
01531
01532
01533
01534 diff = start - (float)beginCol[slit] ;
01535
01536
01537
01538
01539
01540
01541 shift[row_index[slit]] = diff ;
01542 }
01543
01544
01545 for ( z = 0 ; z < ily ; z++ )
01546 {
01547 o_img=cpl_image_new(imsize,N_SLITLETS,CPL_TYPE_FLOAT);
01548 podata=cpl_image_get_data_float(o_img);
01549 for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01550 {
01551 col_counter = beginCol[slit] ;
01552
01553 for ( col = 0 ; col < imsize ; col++ )
01554 {
01555 if ( col_counter > ilx-1 )
01556 {
01557 col_counter-- ;
01558 }
01559 if ( col_counter + z*ilx < 0 )
01560 {
01561 podata[col+row_index[slit]*imsize] = podata[0] ;
01562 }
01563 else
01564 {
01565 podata[col+row_index[slit]*imsize]=pidata[col_counter+z*ilx];
01566 }
01567
01568 col_counter++ ;
01569 }
01570 }
01571 cpl_imagelist_set(returnCube,o_img,z);
01572 }
01573 cpl_free (row_index) ;
01574 cpl_free (beginCol) ;
01575
01576 return returnCube ;
01577 }
01604 cpl_imagelist *
01605 sinfo_new_make_3D_cube_dist ( cpl_image * calibImage,
01606 float firstCol,
01607 float * distances,
01608 float * shift )
01609 {
01610 cpl_imagelist * returnCube ;
01611 float di ;
01612 float diff, start ;
01613 int * row_index ;
01614 int slit ;
01615 int col, z ;
01616 int imsize ;
01617 int * beginCol ;
01618 int col_counter ;
01619 int ilx=0;
01620 int ily=0;
01621
01622 float* podata=NULL;
01623 float* pidata=NULL;
01624 cpl_image* o_img;
01625
01626 if ( NULL == calibImage )
01627 {
01628 sinfo_msg_error(" no resampled image given!\n") ;
01629 return NULL ;
01630 }
01631 ilx=cpl_image_get_size_x(calibImage);
01632 ily=cpl_image_get_size_y(calibImage);
01633 pidata=cpl_image_get_data_float(calibImage);
01634
01635 if ( NULL == distances )
01636 {
01637 sinfo_msg_error("no distances array given from north_south_test()!") ;
01638 return NULL ;
01639 }
01640
01641 if ( N_SLITLETS != 16 )
01642 {
01643 sinfo_msg_error ("wrong number of slitlets given \n" ) ;
01644 return NULL ;
01645 }
01646 imsize = ilx / N_SLITLETS ;
01647
01648
01649 if ( NULL == (row_index = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01650 {
01651 sinfo_msg_error ("cannot allocate memory \n" ) ;
01652 return NULL ;
01653 }
01654 if ( NULL == (beginCol = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01655 {
01656 sinfo_msg_error ("cannot allocate memory \n" ) ;
01657 cpl_free(row_index) ;
01658 return NULL ;
01659 }
01660 if ( NULL == (returnCube = cpl_imagelist_new()) )
01661 {
01662 sinfo_msg_error ("cannot allocate new cube \n" ) ;
01663 cpl_free(row_index) ;
01664 cpl_free(beginCol) ;
01665 return NULL ;
01666 }
01667
01668 di = 0. ;
01669
01670
01671 for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01672
01673 {
01674
01675
01676
01677
01678
01679
01680 row_index[slit] = slit ;
01681
01682
01683 if ( slit == 0 )
01684 {
01685 start = firstCol ;
01686 }
01687 else
01688 {
01689 di += distances[slit-1] ;
01690 start = firstCol + di ;
01691 }
01692 beginCol[slit] = sinfo_new_nint(start) ;
01693
01694
01695
01696 diff = start - (float)beginCol[slit] ;
01697
01698
01699
01700
01701
01702
01703 shift[row_index[slit]] = diff ;
01704 }
01705
01706
01707 for ( z = 0 ; z < ily ; z++ )
01708 {
01709 o_img=cpl_image_new(imsize,N_SLITLETS,CPL_TYPE_FLOAT);
01710 podata=cpl_image_get_data_float(o_img);
01711 for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01712 {
01713 col_counter = beginCol[slit] ;
01714
01715 for ( col = 0 ; col < imsize ; col++ )
01716 {
01717 if ( col_counter > ilx-1 )
01718 {
01719 col_counter-- ;
01720 }
01721 podata[col+row_index[slit]*imsize]=pidata[col_counter+z*ilx];
01722 col_counter++ ;
01723 }
01724 }
01725 cpl_imagelist_set(returnCube,o_img,z);
01726 }
01727 cpl_free (row_index) ;
01728 cpl_free (beginCol) ;
01729
01730 return returnCube ;
01731 }
01732
01750 cpl_imagelist *
01751 sinfo_new_make_3D_cube ( cpl_image * calibImage,
01752 int * kpixshift,
01753 int kpixfirst )
01754 {
01755 cpl_imagelist * returnCube ;
01756 int imsize, kslit, kpix ;
01757 int z, col, recol ;
01758 int ilx=0;
01759 int ily=0;
01760
01761 float* podata=NULL;
01762 float* pidata=NULL;
01763 cpl_image* o_img;
01764
01765 if ( NULL == calibImage )
01766 {
01767 sinfo_msg_error("no resampled image given!\n") ;
01768 return NULL ;
01769 }
01770 ilx=cpl_image_get_size_x(calibImage);
01771 ily=cpl_image_get_size_y(calibImage);
01772 pidata=cpl_image_get_data_float(calibImage);
01773
01774 if ( NULL == kpixshift )
01775 {
01776 sinfo_msg_error("no shift array given!/n") ;
01777 return NULL ;
01778 }
01779
01780 if ( kpixfirst < 0 )
01781 {
01782 sinfo_msg_error("wrong first valid pixel given!/n") ;
01783 return NULL ;
01784 }
01785
01786 if ( N_SLITLETS != 16 )
01787 {
01788 sinfo_msg_error ("wrong number of slitlets given \n" ) ;
01789 return NULL ;
01790 }
01791 imsize = ilx / N_SLITLETS ;
01792
01793 if ( NULL == (returnCube = cpl_imagelist_new()) )
01794 {
01795 sinfo_msg_error ("cannot allocate new cube \n" ) ;
01796 return NULL ;
01797 }
01798
01799
01800 for ( z = 0 ; z < ily ; z++ )
01801 {
01802 o_img=cpl_image_new(imsize,N_SLITLETS,CPL_TYPE_FLOAT);
01803 podata=cpl_image_get_data_float(o_img);
01804 kpix = 0 ;
01805 kslit = 0 ;
01806 recol = -1 ;
01807 for ( col = 0 ; col < ilx ; col++ )
01808 {
01809 if ( col % imsize == 0 )
01810 {
01811 recol = 0 ;
01812 kslit = col/imsize ;
01813 kpix = kpixfirst + kpixshift[kslit] ;
01814 }
01815
01816
01817 podata[recol+kslit*imsize] = pidata[col+kpix+z*ilx] ;
01818 recol++ ;
01819 if ( recol > imsize )
01820 {
01821 sinfo_msg_error("wrong column of reconstructed image, i"
01822 "shouldn't happen!\n") ;
01823 cpl_imagelist_delete(returnCube) ;
01824 return NULL ;
01825 }
01826 }
01827 cpl_imagelist_set(returnCube,o_img,z);
01828 }
01829 return returnCube ;
01830 }
01831
01844 cpl_imagelist *
01845 sinfo_new_determine_mask_cube ( cpl_imagelist * sourceMaskCube,
01846 float lowLimit,
01847 float highLimit )
01848 {
01849 cpl_imagelist * retCube ;
01850 int z, n ;
01851 int ilx=0;
01852 int ily=0;
01853 int inp=0;
01854 int olx=0;
01855 int oly=0;
01856 int onp=0;
01857 float* podata=NULL;
01858 cpl_image* o_img;
01859
01860 if ( sourceMaskCube == NULL )
01861 {
01862 sinfo_msg_error("no cube given!\n") ;
01863 return NULL ;
01864 }
01865 ilx=cpl_image_get_size_x(cpl_imagelist_get(sourceMaskCube,0));
01866 ily=cpl_image_get_size_y(cpl_imagelist_get(sourceMaskCube,0));
01867 inp=cpl_imagelist_get_size(sourceMaskCube);
01868
01869
01870 if ( lowLimit > 0. )
01871 {
01872 sinfo_msg_error("lowLimit wrong!\n") ;
01873 return NULL ;
01874 }
01875 if ( highLimit >= 1. || highLimit < 0. )
01876 {
01877 sinfo_msg_error("highLimit wrong!\n") ;
01878 return NULL ;
01879 }
01880
01881 retCube = cpl_imagelist_duplicate (sourceMaskCube) ;
01882 onp=inp;
01883 olx=ilx;
01884 oly=ily;
01885
01886 for ( z = 0 ; z < onp ; z++ )
01887 {
01888 o_img=cpl_imagelist_get(retCube,0);
01889 podata=cpl_image_get_data_float(o_img);
01890 for ( n = 0 ; n < (int) olx*oly; n++ )
01891 {
01892 if ( podata[n] == 0. )
01893 {
01894 continue ;
01895 }
01896 if ( podata[n] == 1. )
01897 {
01898 continue ;
01899 }
01900 if ( podata[n] >= lowLimit &&
01901 podata[n] <= highLimit )
01902 {
01903 podata[n] = 0. ;
01904 }
01905 else
01906 {
01907 podata[n] = 1. ;
01908 }
01909 }
01910 }
01911 return retCube ;
01912 }
01951 cpl_imagelist *
01952 sinfo_new_interpol_cube ( cpl_imagelist * sourceCube,
01953 cpl_imagelist * maskCube,
01954 int n_neighbors,
01955 int max_radius )
01956 {
01957 cpl_imagelist * returnCube ;
01958 float** spec=NULL ;
01959 float* spec1=NULL ;
01960 int n_im, n_bad, n_bad1, n_bad2 ;
01961 int n_planes, specn, nspec1 ;
01962 int i, m, n, z, ni, kk, p ;
01963 int dis, dismin, dismax ;
01964 int agreed ;
01965 int xcordi, ycordi, xcordm, ycordm ;
01966
01967
01968
01969 int ilx=0;
01970 int ily=0;
01971 int inp=0;
01972
01973 float* pidata=NULL;
01974 float* pmdata=NULL;
01975 float* podata=NULL;
01976 cpl_image* i_img=NULL;
01977 cpl_image* m_img=NULL;
01978 cpl_image* o_img=NULL;
01979
01980 if ( NULL == sourceCube )
01981 {
01982 sinfo_msg_error(" no source cube given!\n") ;
01983 return NULL ;
01984 }
01985
01986
01987 ilx=cpl_image_get_size_x(cpl_imagelist_get(sourceCube,0));
01988 ily=cpl_image_get_size_y(cpl_imagelist_get(sourceCube,0));
01989 inp=cpl_imagelist_get_size(sourceCube);
01990
01991 if ( NULL == maskCube )
01992 {
01993 sinfo_msg_error("no bad pixel mask cube given!\n") ;
01994 return NULL ;
01995 }
01996
01997 if ( n_neighbors <= 0 )
01998 {
01999 sinfo_msg_error("wrong number of neighbors in the spectral "
02000 "direction given!") ;
02001 return NULL ;
02002 }
02003
02004 if ( max_radius <= 0 )
02005 {
02006 sinfo_msg_error("wrong maximal radius for interpolation inside "
02007 "an image plane given!") ;
02008 return NULL ;
02009 }
02010
02011 returnCube = cpl_imagelist_duplicate(sourceCube) ;
02012
02013 n_im = ilx * ily ;
02014 n_planes = inp ;
02015
02016 spec1=cpl_calloc(300,sizeof(float)) ;
02017 spec=sinfo_new_2Dfloatarray(100,2*n_neighbors+1) ;
02018
02019
02020 for ( z = 0 ; z < n_planes ; z++ )
02021 {
02022 m_img=cpl_imagelist_get(maskCube,z);
02023 pmdata=cpl_image_get_data_float(m_img);
02024 o_img=cpl_imagelist_get(returnCube,z);
02025 podata=cpl_image_get_data_float(o_img);
02026
02027
02028
02029
02030
02031
02032 if ( z < n_neighbors )
02033 {
02034 n = z ;
02035 }
02036 else if ( n_planes - z <= n_neighbors)
02037 {
02038 n = n_planes - z -1 ;
02039 }
02040 else
02041 {
02042 n = n_neighbors ;
02043 }
02044
02045 for ( i = 0 ; i < n_im ; i ++ )
02046 {
02047
02048 if ( pmdata[i] != 0. )
02049 {
02050 continue ;
02051 }
02052
02053
02054
02055
02056
02057
02058 n_bad = 0 ;
02059 n_bad1 = 0 ;
02060 n_bad2 = 0 ;
02061
02062 for ( ni = z-n ; ni <= z+n ; ni++ )
02063 {
02064 if ( pmdata[i] == 0. )
02065 {
02066 n_bad++ ;
02067
02068
02069 if ( ni < z )
02070 {
02071 n_bad1++ ;
02072 }
02073 if ( ni > z )
02074 {
02075 n_bad2++ ;
02076 }
02077 }
02078 }
02079
02080
02081
02082
02083
02084
02085
02086
02087 if ( (2*n+1 - n_bad) < 3 || (n - n_bad1) < 1 || (n - n_bad2) < 1 )
02088 {
02089 continue ;
02090 }
02091
02092
02093 kk = 0 ;
02094 for ( ni = z-n ; ni <= z+n ; ni++ )
02095 {
02096 i_img=cpl_imagelist_get(sourceCube,ni);
02097 pidata=cpl_image_get_data_float(i_img);
02098 spec[1][kk] = pmdata[i] != 0. ? pidata[i] : ZERO ;
02099 kk++ ;
02100 }
02101
02102
02103 agreed = 1 ;
02104 specn = 2 ;
02105
02106 dismin = 0 ;
02107 dismax = 1 ;
02108 do
02109 {
02110 for ( m = 0 ; m < n_im ; m++ )
02111 {
02112 if ( pmdata[m] == 0. )
02113 {
02114 continue ;
02115 }
02116
02117
02118
02119
02120
02121 xcordi = i % ilx ;
02122 xcordm = m % ilx ;
02123 ycordi = i / ilx ;
02124 ycordm = m / ilx ;
02125
02126
02127
02128
02129
02130 dis = abs(xcordi-xcordm) + abs(ycordi-ycordm) ;
02131 if ( dis <= dismin || dis > dismax )
02132 {
02133 continue ;
02134 }
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152 n_bad = 0 ;
02153 for ( ni = z-n ; ni <= z+n ; ni++ )
02154 {
02155 if ( pmdata[i] == 0. || pmdata[m] == 0. )
02156 {
02157 n_bad++ ;
02158 }
02159 }
02160 if ( n_bad > 2*n-1 )
02161
02162 {
02163 continue ;
02164 }
02165
02166
02167
02168 kk = 0 ;
02169 for ( ni = z-n ; ni <= z+n ; ni++ )
02170 {
02171 i_img=cpl_imagelist_get(sourceCube,ni);
02172 pidata=cpl_image_get_data_float(i_img);
02173 spec[specn][kk] = pmdata[m] != 0. ? pidata[m] : ZERO ;
02174 kk++ ;
02175 }
02176 specn++ ;
02177 if ( specn > 10 )
02178 {
02179 agreed = 0 ;
02180 break ;
02181 }
02182 }
02183
02184 dismin++ ;
02185 dismax++ ;
02186
02187 if ( dismax > max_radius )
02188 {
02189 agreed = 0 ;
02190 }
02191 } while(agreed) ;
02192
02193 specn-- ;
02194 dismax -= 2 ;
02195
02196
02197
02198
02199
02200
02201 for ( kk = 0 ; kk < 2*n+1 ; kk++ )
02202 {
02203 if ( kk == n )
02204 {
02205 continue ;
02206 }
02207
02208
02209 if ( isnan(spec[1][kk]) )
02210 {
02211 for ( p = 2 ; p <= specn ; p++ )
02212 {
02213 spec[p][kk] = ZERO ;
02214 }
02215 }
02216 else
02217 {
02218 for ( p = 2 ; p <= specn ; p++ )
02219 {
02220 if ( !isnan(spec[p][kk]) && spec[p][kk] != 0. &&
02221 !isnan(spec[p][n]) )
02222 {
02223 spec[p][kk] = spec[1][kk] /
02224 spec[p][kk] * spec[p][n] ;
02225 }
02226 else
02227 {
02228 spec[p][kk] = ZERO ;
02229 }
02230 }
02231 }
02232 }
02233
02234
02235
02236
02237
02238
02239
02240
02241 nspec1 = 0 ;
02242
02243 for ( p = 2 ; p <= specn ; p++ )
02244 {
02245 for ( kk = 0 ; kk < 2*n+1 ; kk++ )
02246 {
02247 if ( !isnan(spec[p][kk]) && kk != n )
02248 {
02249 spec1[nspec1] = spec[p][kk] ;
02250 nspec1++ ;
02251 }
02252 }
02253 }
02254
02255
02256 if ( nspec1 < 18 )
02257 {
02258 continue ;
02259 }
02260
02261
02262 podata[i] = sinfo_new_median(spec1, nspec1) ;
02263 pmdata[i] = 1 ;
02264 }
02265 }
02266 sinfo_new_destroy_2Dfloatarray(&spec,2*n_neighbors+1) ;
02267 cpl_free(spec1);
02268 return returnCube ;
02269 }
02290 cpl_imagelist *
02291 sinfo_new_fine_tune_cube( cpl_imagelist * cube,
02292 float * correct_diff_dist,
02293 int n_order )
02294 {
02295 cpl_imagelist * returnCube ;
02296 float* row_data=NULL ;
02297 float* corrected_row_data=NULL ;
02298 float* xnum=NULL ;
02299 float sum, new_sum ;
02300 float eval ;
02301 float * imageptr ;
02302 int row, col ;
02303 int i, z ;
02304 int imsize, n_points ;
02305 int firstpos ;
02306 int flag;
02307 int ilx=0;
02308 int ily=0;
02309 int inp=0;
02310
02311 float* pidata=NULL;
02312 float* podata=NULL;
02313 cpl_image* i_img=NULL;
02314 cpl_image* o_img=NULL;
02315
02316
02317 if ( NULL == cube )
02318 {
02319 sinfo_msg_error("no input cube given!\n") ;
02320 return NULL ;
02321 }
02322 ilx=cpl_image_get_size_x(cpl_imagelist_get(cube,0));
02323 ily=cpl_image_get_size_y(cpl_imagelist_get(cube,0));
02324 inp=cpl_imagelist_get_size(cube);
02325
02326 if ( NULL == correct_diff_dist )
02327 {
02328 sinfo_msg_error("no distances array from ns_test given!n") ;
02329 return NULL ;
02330 }
02331
02332 if ( n_order <= 0 )
02333 {
02334 sinfo_msg_error("wrong order of interpolation polynom given!") ;
02335 returnCube = cpl_imagelist_duplicate(cube);
02336 return returnCube ;
02337 }
02338
02339 returnCube = cpl_imagelist_duplicate(cube);
02340
02341 imsize = ily ;
02342 if ( imsize != N_SLITLETS )
02343 {
02344 sinfo_msg_error ("wrong image size\n" ) ;
02345 return NULL ;
02346 }
02347
02348 n_points = n_order + 1 ;
02349 if ( n_points % 2 == 0 )
02350 {
02351 firstpos = (int)(n_points/2) - 1 ;
02352 }
02353 else
02354 {
02355 firstpos = (int)(n_points/2) ;
02356 }
02357 xnum=cpl_calloc(n_order+1,sizeof(float)) ;
02358
02359 for ( i = 0 ; i < n_points ; i++ )
02360 {
02361 xnum[i] = i ;
02362 }
02363
02364 row_data=cpl_calloc(ilx,sizeof(float)) ;
02365 corrected_row_data=cpl_calloc(ilx,sizeof(float)) ;
02366
02367 for ( z = 0 ; z < inp ; z++ )
02368 {
02369 i_img=cpl_imagelist_get(cube,z);
02370 pidata=cpl_image_get_data_float(i_img);
02371 o_img=cpl_imagelist_get(returnCube,z);
02372 podata=cpl_image_get_data_float(o_img);
02373
02374
02375 for ( row = 0 ; row < imsize ; row++ )
02376 {
02377 for ( col = 0 ; col < ilx ; col++ )
02378 {
02379 corrected_row_data[col] = 0. ;
02380 }
02381 sum = 0. ;
02382 for ( col = 0 ; col < ilx ; col++ )
02383 {
02384 row_data[col] = pidata[col+row*ilx] ;
02385 if ( isnan(row_data[col]) )
02386 {
02387 row_data[col] = 0. ;
02388 for ( i = col - firstpos ;
02389 i < col -firstpos+n_points ; i++ )
02390 {
02391 if ( i < 0 ) continue ;
02392 if ( i >= ilx) continue ;
02393 corrected_row_data[i] = ZERO ;
02394 }
02395 }
02396 if ( col != 0 && col != ilx - 1 )
02397 {
02398 sum += row_data[col] ;
02399 }
02400 }
02401
02402
02403 new_sum = 0. ;
02404 for ( col = 0 ; col < ilx ; col++ )
02405 {
02406
02407 if ( isnan(corrected_row_data[col]) )
02408 {
02409 continue ;
02410 }
02411 if ( col - firstpos < 0 )
02412 {
02413 imageptr = &row_data[0] ;
02414 eval = correct_diff_dist[row] + col ;
02415 }
02416 else if ( col - firstpos + n_points >= ilx )
02417 {
02418 imageptr = &row_data[ilx - n_points] ;
02419 eval = correct_diff_dist[row] + col + n_points - ilx ;
02420 }
02421 else
02422 {
02423 imageptr = &row_data[col-firstpos] ;
02424 eval = correct_diff_dist[row] + firstpos ;
02425 }
02426
02427
02428 flag = 0;
02429 corrected_row_data[col]=sinfo_new_nev_ille(xnum, imageptr,
02430 n_order, eval, &flag);
02431
02432
02433 if ( col != 0 && col != ilx - 1 )
02434 {
02435 new_sum += corrected_row_data[col] ;
02436 }
02437 }
02438 for ( col = 0 ; col < ilx ; col++ )
02439 {
02440
02441 if ( col == 0 )
02442 {
02443 podata[col+row*ilx] = ZERO ;
02444 }
02445 else if ( col == ilx - 1 )
02446 {
02447 podata[col+row*ilx] = ZERO ;
02448 }
02449 else
02450 {
02451 if ( isnan(corrected_row_data[col]) )
02452 {
02453 podata[col+row*ilx] = ZERO ;
02454 }
02455 else
02456 {
02457 if ( new_sum == 0. ) new_sum = 1. ;
02458
02459 podata[col+row*ilx] = corrected_row_data[col] ;
02460 }
02461 }
02462 }
02463 }
02464 }
02465
02466 cpl_free(xnum) ;
02467 cpl_free(row_data) ;
02468 cpl_free(corrected_row_data) ;
02469
02470 return returnCube ;
02471 }
02472
02492 cpl_imagelist *
02493 sinfo_new_fine_tune_cube_by_FFT( cpl_imagelist * cube,
02494 float * correct_diff_dist )
02495 {
02496 cpl_imagelist * returnCube ;
02497
02498 float* row_data=NULL ;
02499 dcomplex* data=NULL ;
02500 dcomplex* corrected_data=NULL ;
02501
02502 unsigned nn[1];
02503
02504 float phi, pphi ;
02505 float coph, siph ;
02506 int row, col ;
02507 int i, z ;
02508 int imsize ;
02509 int blank_indicator ;
02510
02511
02512 int ilx=0;
02513 int ily=0;
02514 int inp=0;
02515
02516 float* pidata=NULL;
02517 float* podata=NULL;
02518 cpl_image* i_img=NULL;
02519 cpl_image* o_img=NULL;
02520
02521
02522
02523 if ( NULL == cube )
02524 {
02525 sinfo_msg_error(" no input cube given!\n") ;
02526 return NULL ;
02527 }
02528 ilx=cpl_image_get_size_x(cpl_imagelist_get(cube,0));
02529 ily=cpl_image_get_size_y(cpl_imagelist_get(cube,0));
02530 inp=cpl_imagelist_get_size(cube);
02531
02532 nn[1] = ilx ;
02533 if ( NULL == correct_diff_dist )
02534 {
02535 sinfo_msg_error("no distances array from ns_test given!") ;
02536 return NULL ;
02537 }
02538
02539 returnCube = cpl_imagelist_duplicate( cube ) ;
02540
02541 imsize = ily ;
02542 if ( imsize != N_SLITLETS )
02543 {
02544 sinfo_msg_error ("wrong image size\n" ) ;
02545 return NULL ;
02546 }
02547
02548 data=cpl_calloc(ilx,sizeof(dcomplex)) ;
02549 corrected_data=cpl_calloc(ilx,sizeof(dcomplex)) ;
02550
02551 row_data=cpl_calloc(ilx,sizeof(float)) ;
02552
02553 for ( z = 0 ; z < inp ; z++ )
02554 {
02555 i_img=cpl_imagelist_get(cube,z);
02556 pidata=cpl_image_get_data_float(i_img);
02557 o_img=cpl_imagelist_get(returnCube,z);
02558 podata=cpl_image_get_data_float(o_img);
02559
02560 for ( row = 0 ; row < imsize ; row++ )
02561 {
02562 blank_indicator = 1 ;
02563 for ( col = 0 ; col < ilx ; col++ )
02564 {
02565
02566 row_data[col] = pidata[col+row*ilx] ;
02567 data[col].x = row_data[col] ;
02568 data[col].y = 0. ;
02569
02570 if ( isnan(row_data[col]) )
02571 {
02572 blank_indicator = 0 ;
02573 }
02574 }
02575
02576
02577 if ( blank_indicator == 0 )
02578 {
02579 for ( col = 0 ; col < ilx ; col++ )
02580 {
02581 podata[col+row*ilx] = ZERO ;
02582 }
02583 continue ;
02584 }
02585
02586
02587 sinfo_fftn( data, nn, 1, 1 ) ;
02588
02589
02590 phi = 2*PI_NUMB/(float)ilx * correct_diff_dist[row] ;
02591 for ( i = 0 ; i < ilx ; i++ )
02592 {
02593
02594 if ( i <= ilx/2 )
02595 {
02596
02597 pphi = phi * (float)(i) ;
02598
02599 coph = cos ( pphi ) ;
02600 siph = sin ( pphi ) ;
02601 }
02602 else
02603 {
02604
02605 pphi = phi * (float)(i - ilx/2) ;
02606
02607 coph = cos ( pphi ) ;
02608 siph = sin ( pphi ) ;
02609 }
02610
02611
02612
02613
02614
02615
02616
02617
02618 corrected_data[i].x = data[i].x * coph - data[i].y * siph ;
02619
02620 corrected_data[i].y = data[i].x * siph + data[i].y * coph ;
02621 }
02622
02623
02624 sinfo_fftn( corrected_data, nn, 1, -1 ) ;
02625
02626
02627 for ( i = 0 ; i < ilx ; i++ )
02628 {
02629 corrected_data[i].x /= ilx ;
02630 corrected_data[i].y /= ilx ;
02631 }
02632
02633
02634
02635 for ( col = 0 ; col < ilx ; col++ )
02636 {
02637 if ( col == 0 )
02638 {
02639 podata[col+row*ilx] = ZERO ;
02640 }
02641 else if ( col == ilx - 1 )
02642 {
02643 podata[col+row*ilx] = ZERO ;
02644 }
02645 else
02646 {
02647 podata[col+row*ilx] = corrected_data[col].x ;
02648 }
02649 }
02650 }
02651 }
02652
02653 cpl_free(data) ;
02654 cpl_free(corrected_data) ;
02655
02656
02657 cpl_free(row_data);
02658 return returnCube ;
02659 }
02687 cpl_imagelist * sinfo_new_fine_tune_cube_by_spline ( cpl_imagelist * cube,
02688 float * correct_diff_dist )
02689 {
02690 cpl_imagelist * returnCube ;
02691
02692 float* row_data=NULL ;
02693 float* corrected_row_data=NULL ;
02694 float* xnum=NULL ;
02695 float* eval=NULL ;
02696
02697 float sum, new_sum ;
02698 int row, col ;
02699 int i, z ;
02700 int imsize ;
02701 int ilx=0;
02702 int ily=0;
02703 int inp=0;
02704
02705 float* pidata=NULL;
02706 float* podata=NULL;
02707 cpl_image* i_img=NULL;
02708 cpl_image* o_img=NULL;
02709
02710
02711 if ( NULL == cube )
02712 {
02713 sinfo_msg_error("no input cube given!\n") ;
02714 return NULL ;
02715 }
02716 ilx=cpl_image_get_size_x(cpl_imagelist_get(cube,0));
02717 ily=cpl_image_get_size_y(cpl_imagelist_get(cube,0));
02718 inp=cpl_imagelist_get_size(cube);
02719
02720 if ( NULL == correct_diff_dist )
02721 {
02722 sinfo_msg_error("no distances array from ns_test given!/n") ;
02723 return NULL ;
02724 }
02725
02726 imsize = ily ;
02727 if ( imsize != N_SLITLETS )
02728 {
02729 sinfo_msg_error ("wrong image size\n" ) ;
02730 return NULL ;
02731 }
02732
02733 returnCube = cpl_imagelist_duplicate( cube ) ;
02734
02735 row_data=cpl_calloc(ilx,sizeof(float)) ;
02736 corrected_row_data=cpl_calloc(ilx,sizeof(float)) ;
02737 xnum=cpl_calloc(ilx,sizeof(float)) ;
02738 eval=cpl_calloc(ilx,sizeof(float)) ;
02739
02740
02741 for ( i = 0 ; i < ilx ; i++ )
02742 {
02743 xnum[i] = i ;
02744 }
02745
02746
02747 for ( z = 0 ; z < inp ; z++ )
02748 {
02749 i_img=cpl_imagelist_get(cube,z);
02750 pidata=cpl_image_get_data_float(i_img);
02751 o_img=cpl_imagelist_get(returnCube,z);
02752 podata=cpl_image_get_data_float(o_img);
02753
02754 for ( row = 0 ; row < imsize ; row++ )
02755 {
02756 for ( col = 0 ; col < ilx ; col++ )
02757 {
02758 corrected_row_data[col] = 0. ;
02759 }
02760 sum = 0. ;
02761
02762
02763 for ( col = 0 ; col < ilx ; col++ )
02764 {
02765 eval[col] = correct_diff_dist[row] + (float)col ;
02766 row_data[col] = pidata[col+row*ilx] ;
02767 if (col != 0 && col != ilx - 1 && !isnan(row_data[col]) )
02768 {
02769 sum += row_data[col] ;
02770 }
02771 if (isnan(row_data[col]) )
02772 {
02773 for ( i = col -1 ; i <= col+1 ; i++ )
02774 {
02775 if ( i < 0 ) continue ;
02776 if ( i >= ilx ) continue ;
02777 corrected_row_data[i] = ZERO ;
02778 }
02779 row_data[col] = 0. ;
02780 }
02781 }
02782
02783
02784
02785
02786
02787
02788 if ( -1 == sinfo_function1d_natural_spline(xnum,row_data, ilx,
02789 eval,corrected_row_data,
02790 ilx ) )
02791 {
02792 sinfo_msg_error("error in spline interpolation\n") ;
02793 cpl_imagelist_delete(returnCube) ;
02794 return NULL ;
02795 }
02796
02797 new_sum = 0. ;
02798 for ( col = 0 ; col < ilx ; col++ )
02799 {
02800 if (isnan(corrected_row_data[col])) continue ;
02801
02802
02803 if ( col != 0 && col != ilx - 1 )
02804 {
02805 new_sum += corrected_row_data[col] ;
02806 }
02807 }
02808 for ( col = 0 ; col < ilx ; col++ )
02809 {
02810
02811
02812
02813
02814
02815 if ( col == 0 )
02816 {
02817 podata[col+row*ilx] = ZERO ;
02818 }
02819 else if ( col == ilx - 1 )
02820 {
02821 podata[col+row*ilx] = ZERO ;
02822 }
02823 else
02824 {
02825 if ( isnan(corrected_row_data[col]) )
02826 {
02827 podata[col+row*ilx] = ZERO ;
02828 }
02829 else
02830 {
02831 if (new_sum == 0.) new_sum = 1. ;
02832
02833
02834
02835
02836 podata[col+row*ilx] = corrected_row_data[col] ;
02837 }
02838 }
02839 }
02840 }
02841 }
02842
02843 cpl_free(row_data) ;
02844 cpl_free(corrected_row_data) ;
02845 cpl_free(xnum) ;
02846 cpl_free(eval) ;
02847
02848 return returnCube ;
02849 }
02850
02868 float *
02869 sinfo_new_calibrate_ns_test( cpl_image * ns_image,
02870 int n_slitlets,
02871 int halfWidth,
02872 float fwhm,
02873 float minDiff,
02874 float estimated_dist,
02875 float devtol,
02876 int bottom,
02877 int top )
02878 {
02879 int i, j, k, m, row, col, n, ni, na ;
02880 int position, counter, iters ;
02881 int xdim, ndat, its, numpar ;
02882 float sum, mean, maxval ;
02883 float tol, lab ;
02884 float * distances ;
02885 float * ret_distances ;
02886
02887 float * xdat, * wdat ;
02888 int * mpar ;
02889
02890 pixelvalue* row_buf=NULL ;
02891 float** distances_buf=NULL ;
02892 float* x_position=NULL ;
02893 int* found=NULL;
02894 int* found_clean=NULL ;
02895 int* found_cleanit=NULL ;
02896
02897 Vector * line ;
02898 FitParams ** par ;
02899 int foundit, begin, end ;
02900 int zeroindicator ;
02901 int row_index ;
02902
02903 int ilx=0;
02904 int ily=0;
02905 float* pidata=NULL;
02906
02907 if ( ns_image == NULL )
02908 {
02909 sinfo_msg_error("sorry, no image given\n") ;
02910 return NULL ;
02911 }
02912 if ( n_slitlets < 1 )
02913 {
02914 sinfo_msg_error("wrong number of slitlets given\n") ;
02915 return NULL ;
02916 }
02917 if ( halfWidth < 0 || halfWidth >= estimated_dist )
02918 {
02919 sinfo_msg_error("wrong half width given\n") ;
02920 return NULL ;
02921 }
02922 if ( fwhm <= 0. )
02923 {
02924 sinfo_msg_error("wrong fwhm given\n") ;
02925 return NULL ;
02926 }
02927 if ( minDiff < 1. )
02928 {
02929 sinfo_msg_error("wrong minDiff given\n") ;
02930 return NULL ;
02931 }
02932
02933
02934 if (NULL==(distances=(float *)cpl_calloc( n_slitlets , sizeof (float) )))
02935 {
02936 sinfo_msg_error("could not allocate memory\n") ;
02937 return NULL ;
02938 }
02939
02940 if (NULL == (ret_distances = (float *) cpl_calloc ( n_slitlets ,
02941 sizeof (float) )))
02942 {
02943 sinfo_msg_error("could not allocate memory\n") ;
02944 return NULL ;
02945 }
02946
02947 ilx=cpl_image_get_size_x(ns_image);
02948 ily=cpl_image_get_size_y(ns_image);
02949 pidata=cpl_image_get_data_float(ns_image);
02950
02951 row_buf=(pixelvalue*)cpl_calloc(ilx,sizeof(pixelvalue)) ;
02952 x_position=cpl_calloc(n_slitlets,sizeof(float)) ;
02953 found=cpl_calloc(3*n_slitlets,sizeof(int));
02954 found_clean=cpl_calloc(3*n_slitlets,sizeof(int)) ;
02955 found_cleanit=cpl_calloc(3*n_slitlets,sizeof(int)) ;
02956 distances_buf=sinfo_new_2Dfloatarray(ily,n_slitlets) ;
02957
02958
02959 for ( row = 0 ; row < ily ; row++ )
02960 {
02961 zeroindicator = 0 ;
02962
02963
02964 for ( i = 0 ; i < n_slitlets ; i++ )
02965 {
02966 distances_buf[row][i] = ZERO ;
02967 }
02968
02969
02970 for ( col = 0 ; col < ilx ; col++ )
02971 {
02972 row_buf[col] = pidata[col + row*ilx] ;
02973 }
02974
02975
02976 sum = 0. ;
02977 n = 0 ;
02978 for ( i = 0 ; i < ilx ; i++ )
02979 {
02980 if ( isnan(row_buf[i]) )
02981 {
02982 continue ;
02983 }
02984 sum += row_buf[i] ;
02985 n++ ;
02986 }
02987 mean = sum / (float)n ;
02988
02989
02990 n = 0 ;
02991 for ( i = 0 ; i < ilx ; i++ )
02992 {
02993 if (isnan(row_buf[i]))
02994 {
02995 continue ;
02996 }
02997 if ( row_buf[i] > mean + ESTIMATE )
02998 {
02999 found[n] = i ;
03000 n++ ;
03001 }
03002 }
03003
03004 if ( n < n_slitlets )
03005 {
03006 sinfo_msg_warning("t4 wrong number of intensity columns "
03007 "found in row: %d, found number: %d", row, n) ;
03008 continue ;
03009 }
03010 else
03011 {
03012
03013 na = 0 ;
03014 for ( i = 1 ; i < n ; i ++ )
03015 {
03016 if ( found[i] - found[i-1] < halfWidth )
03017 {
03018 begin = found[i] - halfWidth ;
03019 if ( begin < 0 )
03020 {
03021 begin = 0 ;
03022 }
03023 end = found[i] + halfWidth ;
03024 if ( end >= ilx )
03025 {
03026 end = ilx - 1 ;
03027 }
03028
03029
03030 maxval = -FLT_MAX ;
03031 foundit = 0 ;
03032 for ( j = begin ; j <= end ; j++ )
03033 {
03034
03035 if (isnan(row_buf[j]))
03036 {
03037 continue ;
03038 }
03039 if (row_buf[j] >= maxval )
03040 {
03041 maxval = row_buf[j] ;
03042 foundit = j ;
03043 }
03044 }
03045 if (maxval == -FLT_MAX)
03046 {
03047 continue ;
03048 }
03049 for ( k = 0 ; k < na ; k++ )
03050 {
03051 if ( found_cleanit[k] >= begin &&
03052 found_cleanit[k] < foundit )
03053 {
03054 na-- ;
03055 }
03056 }
03057 for ( k = 0 ; k < n ; k++ )
03058 {
03059 if ( found[k] == foundit)
03060 {
03061 if ( found_cleanit[na-1] != found[k] )
03062 {
03063 found_cleanit[na] = found[k] ;
03064 na++ ;
03065 }
03066 }
03067 }
03068 }
03069 else
03070 {
03071 if ( i == 1 )
03072 {
03073 found_cleanit[na] = found[0] ;
03074 na++ ;
03075 found_cleanit[na] = found[1] ;
03076 na++ ;
03077 }
03078 else
03079 {
03080 if ( found_cleanit[na-1] != found[i-1])
03081 {
03082 found_cleanit[na] = found[i-1] ;
03083 na++ ;
03084 }
03085 if ( found_cleanit[na-1] != found[i])
03086 {
03087 found_cleanit[na] = found[i] ;
03088 na++ ;
03089 }
03090 }
03091 }
03092 }
03093
03094
03095 j = 1 ;
03096 for ( i = 1 ; i < na ; i++ )
03097 {
03098 if ( (float)(found_cleanit[i] - found_cleanit[i-1]) <
03099 (estimated_dist - devtol) ||
03100 (float)(found_cleanit[i] - found_cleanit[i-1]) >
03101 (estimated_dist + devtol) )
03102 {
03103 continue ;
03104 }
03105 else
03106 {
03107 found_clean[j-1] = found_cleanit[i-1] ;
03108 found_clean[j] = found_cleanit[i] ;
03109 j++ ;
03110 }
03111 }
03112 }
03113 if ( j > n_slitlets )
03114 {
03115
03116 ni = 1 ;
03117 for ( i = 1 ; i < j ; i++ )
03118 {
03119 if ( (float)(found_clean[i] - found_clean[i-1]) <
03120 (estimated_dist - devtol ) ||
03121 (float)(found_clean[i] - found_clean[i-1]) >
03122 (estimated_dist + devtol ) )
03123 {
03124 continue ;
03125 }
03126 else
03127 {
03128 found_clean[ni-1] = found_clean[i-1] ;
03129 found_clean[ni] = found_clean[i] ;
03130 ni++ ;
03131 }
03132 }
03133 if ( ni != n_slitlets )
03134 {
03135 sinfo_msg_warning("t5 wrong number of intensity columns "
03136 "found in row: %d, found number: %d",
03137 row,ni) ;
03138 continue ;
03139 }
03140 else
03141 {
03142 j = ni ;
03143 }
03144 }
03145 else if ( j < n_slitlets )
03146 {
03147 sinfo_msg_warning("t6 wrong number of intensity columns found "
03148 "in row: %d , found number: %d\n", row, j) ;
03149 continue ;
03150 }
03151 counter = 0 ;
03152
03153 for ( i = 0 ; i < j ; i++ )
03154 {
03155
03156 if ( NULL == (line = sinfo_new_vector (2*halfWidth + 1)) )
03157 {
03158 sinfo_msg_error ("cannot allocate new Vector \n") ;
03159 cpl_free(distances) ;
03160 return NULL ;
03161 }
03162
03163
03164 xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
03165 wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
03166 mpar = (int *) cpl_calloc( MAXPAR, sizeof (int) ) ;
03167 par = sinfo_new_fit_params(1) ;
03168
03169 m = 0 ;
03170 for ( k = found_clean[i]-halfWidth ;
03171 k <= found_clean[i]+halfWidth ; k++ )
03172 {
03173 if ( k < 0 )
03174 {
03175 k = 0. ;
03176 }
03177 else if ( k >= ilx )
03178 {
03179 k = ilx - 1 ;
03180 }
03181 else if ( isnan(row_buf[k]) )
03182 {
03183 zeroindicator = 1 ;
03184 break ;
03185 }
03186 else
03187 {
03188 line -> data[m] = row_buf[k] ;
03189 m++ ;
03190 }
03191 }
03192 if ( zeroindicator == 1 )
03193 {
03194 sinfo_new_destroy_vector(line) ;
03195 cpl_free(xdat) ;
03196 cpl_free(wdat) ;
03197 cpl_free(mpar) ;
03198 sinfo_new_destroy_fit_params(&par) ;
03199 break ;
03200 }
03201
03202
03203
03204
03205
03206 maxval = -FLT_MAX ;
03207 position = -INT32_MAX ;
03208 for ( k = 0 ; k < m ; k++ )
03209 {
03210 xdat[k] = k ;
03211 wdat[k] = 1.0 ;
03212 if ( line -> data[k] >= maxval )
03213 {
03214 maxval = line -> data[k] ;
03215 position = k ;
03216 }
03217 }
03218
03219
03220 xdim = XDIM ;
03221 ndat = line -> n_elements ;
03222 numpar = MAXPAR ;
03223 tol = TOL ;
03224 lab = LAB ;
03225 its = ITS ;
03226 (*par) -> fit_par[1] = fwhm ;
03227 (*par) -> fit_par[2] = (float) position ;
03228 (*par) -> fit_par[3] = (float) (line -> data[0] +
03229 line -> data[line->n_elements - 1]) / 2.0 ;
03230 (*par) -> fit_par[0] = maxval - ((*par) -> fit_par[3]) ;
03231
03232
03233
03234 if ( (*par) -> fit_par[0] < minDiff )
03235 {
03236 sinfo_msg_warning ("sorry, signal of line too low to fit "
03237 "in row: %d in slitlet %d\n", row, i) ;
03238 sinfo_new_destroy_vector(line) ;
03239 cpl_free(xdat) ;
03240 cpl_free(wdat) ;
03241 cpl_free(mpar) ;
03242 sinfo_new_destroy_fit_params(&par) ;
03243 continue ;
03244 }
03245
03246 for ( k = 0 ; k < MAXPAR ; k++ )
03247 {
03248 (*par) -> derv_par[k] = 0.0 ;
03249 mpar[k] = 1 ;
03250 }
03251
03252 if ( 0 > ( iters = sinfo_new_lsqfit_c(xdat, &xdim,
03253 line -> data, wdat,
03254 &ndat, (*par) -> fit_par,
03255 (*par) -> derv_par, mpar,
03256 &numpar, &tol,
03257 &its, &lab )) )
03258 {
03259
03260
03261
03262
03263
03264
03265 sinfo_new_destroy_vector(line) ;
03266 cpl_free(xdat) ;
03267 cpl_free(wdat) ;
03268 cpl_free(mpar) ;
03269 sinfo_new_destroy_fit_params(&par) ;
03270 continue ;
03271 }
03272
03273
03274 if ( (*par) -> fit_par[0] <= 0. || (*par) -> fit_par[1] <= 0. ||
03275 (*par) -> fit_par[2] < 0. )
03276 {
03277 sinfo_msg_warning ("negative parameters as fit result, not "
03278 "used! in row %d in slitlet %d", row, i) ;
03279 sinfo_new_destroy_vector(line) ;
03280 cpl_free(xdat) ;
03281 cpl_free(wdat) ;
03282 cpl_free(mpar) ;
03283 sinfo_new_destroy_fit_params(&par) ;
03284 continue ;
03285 }
03286
03287
03288
03289 (*par) -> fit_par[2] = (float) (found_clean[i] - halfWidth) +
03290 (*par) -> fit_par[2] ;
03291 x_position[counter] = (*par) -> fit_par[2] ;
03292 counter ++ ;
03293
03294
03295 sinfo_new_destroy_fit_params(&par) ;
03296 sinfo_new_destroy_vector ( line ) ;
03297 cpl_free ( xdat ) ;
03298 cpl_free ( wdat ) ;
03299 cpl_free ( mpar ) ;
03300 }
03301 if (zeroindicator == 1)
03302 {
03303 sinfo_msg_debug ("bad pixel in fitting box in row: %d\n", row) ;
03304 continue ;
03305 }
03306
03307 if ( counter != n_slitlets )
03308 {
03309 sinfo_msg_warning("wrong number of slitlets found "
03310 "in row: %d", row) ;
03311 continue ;
03312 }
03313
03314 for ( i = 0 ; i < n_slitlets ; i++ )
03315 {
03316 distances_buf[row][i] = x_position[i] - (15.5 + 32.*(float)i) ;
03317 }
03318 }
03319
03320
03321
03322
03323
03324 for ( i = 0 ; i < n_slitlets ; i++ )
03325 {
03326 n = 0 ;
03327 sum = 0. ;
03328 for ( row = bottom ; row < top ; row++ )
03329 {
03330 if ( fabs( distances_buf[row][i] ) > devtol ||
03331 isnan(distances_buf[row][i]) )
03332 {
03333
03334
03335
03336
03337
03338
03339 continue ;
03340 }
03341 sum += distances_buf[row][i] ;
03342 n++ ;
03343 }
03344 if ( n < 2 )
03345 {
03346 sinfo_msg_error("distances array could not be determined"
03347 " completely!, deviations of distances from"
03348 " devtol too big" ) ;
03349 cpl_free(distances) ;
03350 return NULL ;
03351 }
03352 else
03353 {
03354 distances[i] = sum / (float)n ;
03355 }
03356 }
03357
03358
03359
03360 for ( i = 0 ; i < n_slitlets ; i++ )
03361 {
03362 switch (i)
03363 {
03364 case 0:
03365 row_index = 8 ;
03366 break ;
03367 case 1:
03368 row_index = 7 ;
03369 break ;
03370 case 2:
03371 row_index = 9 ;
03372 break ;
03373 case 3:
03374 row_index = 6 ;
03375 break ;
03376 case 4:
03377 row_index = 10 ;
03378 break ;
03379 case 5:
03380 row_index = 5 ;
03381 break ;
03382 case 6:
03383 row_index = 11 ;
03384 break ;
03385 case 7:
03386 row_index = 4 ;
03387 break ;
03388 case 8:
03389 row_index = 12 ;
03390 break ;
03391 case 9:
03392 row_index = 3 ;
03393 break ;
03394 case 10:
03395 row_index = 13 ;
03396 break ;
03397 case 11:
03398 row_index = 2 ;
03399 break ;
03400 case 12:
03401 row_index = 14 ;
03402 break ;
03403 case 13:
03404 row_index = 1 ;
03405 break ;
03406 case 14:
03407 row_index = 15 ;
03408 break ;
03409 case 15:
03410 row_index = 0 ;
03411 break ;
03412 case 16:
03413 row_index = 31 ;
03414 break ;
03415 case 17:
03416 row_index = 16 ;
03417 break ;
03418 case 18:
03419 row_index = 30 ;
03420 break ;
03421 case 19:
03422 row_index = 17 ;
03423 break ;
03424 case 20:
03425 row_index = 29 ;
03426 break ;
03427 case 21:
03428 row_index = 18 ;
03429 break ;
03430 case 22:
03431 row_index = 28 ;
03432 break ;
03433 case 23:
03434 row_index = 19 ;
03435 break ;
03436 case 24:
03437 row_index = 27 ;
03438 break ;
03439 case 25:
03440 row_index = 20 ;
03441 break ;
03442 case 26:
03443 row_index = 26 ;
03444 break ;
03445 case 27:
03446 row_index = 21 ;
03447 break ;
03448 case 28:
03449 row_index = 25 ;
03450 break ;
03451 case 29:
03452 row_index = 22 ;
03453 break ;
03454 case 30:
03455 row_index = 24 ;
03456 break ;
03457 case 31:
03458 row_index = 23 ;
03459 break ;
03460 default:
03461 sinfo_msg_error("wrong number of a slitlet\n") ;
03462 cpl_free (distances) ;
03463 return NULL ;
03464 }
03465 ret_distances[row_index] = distances[i] ;
03466 }
03467 cpl_free(distances) ;
03468
03469 cpl_free(row_buf) ;
03470 cpl_free(x_position) ;
03471 cpl_free(found);
03472 cpl_free(found_clean) ;
03473 cpl_free(found_cleanit) ;
03474 sinfo_new_destroy_2Dfloatarray(&distances_buf,n_slitlets) ;
03475
03476
03477 return ret_distances ;
03478 }
03503 cpl_image *
03504 sinfo_new_make_true_resamp(cpl_image * calibImage,
03505 cpl_image * wavemap)
03506 {
03507 cpl_image * returnImage ;
03508 float edges[33] ;
03509 int imsize, kslit,i,j ;
03510 int slit_index ;
03511 int z, col, recol ;
03512 int wlx=0;
03513 int wly=0;
03514 int clx=0;
03515 int cly=0;
03516
03517 float* pcdata=NULL;
03518 float* pwdata=NULL;
03519 float* prdata=NULL;
03520
03521
03522 wlx=cpl_image_get_size_x(wavemap);
03523 wly=cpl_image_get_size_y(wavemap);
03524 pwdata=cpl_image_get_data_float(wavemap);
03525
03526 edges[0]=0;
03527 j=1;
03528 for(i=0;i<wlx-1;i++)
03529 {
03530 if((pwdata[i]-pwdata[i+1])>0.0025 || (pwdata[i]-pwdata[i+1])<-0.0025)
03531 {
03532 sinfo_msg_error("wavemap sinfo_edge %d", i+1);
03533 edges[j]=i+1;
03534 j++;
03535 }
03536 }
03537 edges[32]=2048;
03538
03539 clx=cpl_image_get_size_x(calibImage);
03540 cly=cpl_image_get_size_y(calibImage);
03541 pcdata=cpl_image_get_data_float(calibImage);
03542
03543 imsize = clx / N_SLITLETS ;
03544
03545
03546 returnImage = cpl_image_new(clx,cly,CPL_TYPE_FLOAT);
03547 prdata=cpl_image_get_data_float(returnImage);
03548 for ( z = 0 ; z < cly ; z++ )
03549 {
03550 for ( col = 0 ; col < clx ; col++ )
03551 prdata[col+z*clx]=ZERO;
03552 }
03553
03554
03555
03556 for ( z = 0 ; z < cly ; z++ )
03557 {
03558 kslit = 0 ;
03559 slit_index = -1 ;
03560 recol = -1 ;
03561 for ( col = 0 ; col < clx ; col++ )
03562 {
03563
03564
03565 recol = 0 ;
03566
03567 for(i=0;i<32;i++)
03568 {
03569 if(col>=sinfo_new_nint(edges[i]) &&
03570 col<sinfo_new_nint(edges[i+1]))
03571 kslit=i;
03572 }
03573
03574 switch (kslit)
03575 {
03576 case 0:
03577 slit_index = 8 ;
03578 break ;
03579 case 1:
03580 slit_index = 7 ;
03581 break ;
03582 case 2:
03583 slit_index = 9 ;
03584 break ;
03585 case 3:
03586 slit_index = 6 ;
03587 break ;
03588 case 4:
03589 slit_index = 10 ;
03590 break ;
03591 case 5:
03592 slit_index = 5 ;
03593 break ;
03594 case 6:
03595 slit_index = 11 ;
03596 break ;
03597 case 7:
03598 slit_index = 4 ;
03599 break ;
03600 case 8:
03601 slit_index = 12 ;
03602 break ;
03603 case 9:
03604 slit_index = 3 ;
03605 break ;
03606 case 10:
03607 slit_index = 13 ;
03608 break ;
03609 case 11:
03610 slit_index = 2 ;
03611 break ;
03612 case 12:
03613 slit_index = 14 ;
03614 break ;
03615 case 13:
03616 slit_index = 1 ;
03617 break ;
03618 case 14:
03619 slit_index = 15 ;
03620 break ;
03621 case 15:
03622 slit_index = 0 ;
03623 break ;
03624 case 16:
03625 slit_index = 31 ;
03626 break ;
03627 case 17:
03628 slit_index = 16 ;
03629 break ;
03630 case 18:
03631 slit_index = 30 ;
03632 break ;
03633 case 19:
03634 slit_index = 17 ;
03635 break ;
03636 case 20:
03637 slit_index = 29 ;
03638 break ;
03639 case 21:
03640 slit_index = 18 ;
03641 break ;
03642 case 22:
03643 slit_index = 28 ;
03644 break ;
03645 case 23:
03646 slit_index = 19 ;
03647 break ;
03648 case 24:
03649 slit_index = 27 ;
03650 break ;
03651 case 25:
03652 slit_index = 20 ;
03653 break ;
03654 case 26:
03655 slit_index = 26 ;
03656 break ;
03657 case 27:
03658 slit_index = 21 ;
03659 break ;
03660 case 28:
03661 slit_index = 25 ;
03662 break ;
03663 case 29:
03664 slit_index = 22 ;
03665 break ;
03666 case 30:
03667 slit_index = 24 ;
03668 break ;
03669 case 31:
03670 slit_index = 23 ;
03671 break ;
03672 default:
03673 sinfo_msg_error("wrong slitlet index: couldn't be a "
03674 "spiffi image, there must be 32 "
03675 "slitlets!") ;
03676 break ;
03677 }
03678
03679
03680
03681
03682 if((col-sinfo_new_nint(edges[kslit]))>0 &&
03683 (col-sinfo_new_nint(edges[kslit]))<imsize-1 )
03684 prdata[(col-sinfo_new_nint(edges[kslit]))+
03685 slit_index*imsize+z*clx] =
03686 pcdata[col+z*clx] ;
03687 else
03688 prdata[(col-sinfo_new_nint(edges[kslit]))+
03689 slit_index*imsize+z*clx] = ZERO;
03690
03691
03692 }
03693 }
03694 return returnImage ;
03695 }
03696
03697
03698
03699
03700
03701
03702
03703
03704
03705
03706
03707
03708
03709
03710
03711
03712
03713
03714
03715
03716
03717
03718
03719
03720
03721
03722
03723
03724
03725
03726
03727
03728
03729
03730
03731
03732
03733
03734
03735
03736
03737
03738
03739
03740
03741
03742
03743
03744
03745
03746
03747
03748
03749
03750
03751
03752
03753
03754
03755
03756
03757
03758
03759
03760
03761
03762
03763
03764
03765
03766
03767
03768
03769
03770
03771
03772
03773
03774
03775
03776
03777
03778
03779
03780
03781
03782
03783
03784
03785
03786
03787
03788
03789
03790
03791
03792
03793
03794
03795
03796
03797
03798
03799
03800
03801
03802
03803
03804
03805
03806