00001 /*******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 *
00005 *
00006 * who when what
00007 * -------- -------- ----------------------------------------------
00008 * schreib 30/08/00 created
00009 */
00010
00011 /************************************************************************
00012 * NAME
00013 * cube_construct.c -
00014 * some procedures to construct a data cube
00015 *
00016 * SYNOPSIS
00017 *
00018 * 1) OneImage * convolveNSImageByGauss( OneImage * lineImage,
00019 * int hw )
00020 *
00021 * 2) float * north_south_test( OneImage * ns_image,
00022 * int n_slitlets,
00023 * int halfWidth,
00024 * float fwhm,
00025 * float minDiff,
00026 * float estimated_dist,
00027 * float devtol )
00028 *
00029 * 3) OneCube * makeCube ( OneImage * calibImage,
00030 * float * distances,
00031 * float * correct_diff_dist )
00032 *
00033 * 4) OneCube * makeCubeSpi ( OneImage * calibImage,
00034 * float ** slit_edges,
00035 * float * shift )
00036 *
00037 * 5) OneCube * makeCubeDist ( OneImage * calibImage,
00038 * float firstCol,
00039 * float * distances,
00040 * float * shift )
00041 *
00042 * 6) OneCube * make3DCubeDist ( OneImage * calibImage,
00043 * float firstCol,
00044 * float * distances,
00045 * float * shift )
00046 *
00047 * 7) OneCube * make3DCube ( OneImage * calibImage,
00048 * int * kpixshift,
00049 * int kpixfirst )
00050 *
00051 * 8) OneCube * determineMaskCube ( OneCube * sourceMaskCube,
00052 * float lowLimit,
00053 * float highLimit )
00054 *
00055 * 9) OneCube * interpolCube ( OneCube * sourceCube,
00056 * OneCube * maskCube,
00057 * int n_neighbors,
00058 * int max_radius )
00059 *
00060 * 10) OneCube * fineTuneCube( OneCube * cube,
00061 * float * correct_diff_dist )
00062 *
00063 * 11) OneCube * fineTuneCubeByFFT( OneCube * cube,
00064 * float * correct_diff_dist )
00065 *
00066 * 12) OneCube * fineTuneCubeBySpline ( OneCube * cube,
00067 * float * correct_diff_dist )
00068 *
00069 * DESCRIPTION
00070 *
00071 * 1) convolves a north-south-test image with a gaussian
00072 * with user given integer half width by using the eclipse
00073 * routine function1d_filter_lowpass().
00074 * 2) determines the distances of the slitlets
00075 * 3) makes a data cube out of a resampled source image
00076 * this SPIFFI specific routine takes into account the
00077 * Spiffi slitlet order on the detector.
00078 * Also shifts the resulting image rows by one pixel if
00079 * necessary according to the distances array gained from
00080 * the north-south test routine.
00081 * Can do the same with the bad pixel map image to generate a
00082 * bad pixel mask cube.
00083 * 4) makes a data cube out of a resampled source image
00084 * this SPIFFI specific routine takes into account the
00085 * Spiffi slitlet order on the detector.
00086 * This routine takes fitted slitlet positions into account.
00087 * Can do the same with the bad pixel map image to generate a
00088 * bad pixel mask cube.
00089 * 5) makes a data cube out of a resampled source image
00090 * this SPIFFI specific routine takes into account the
00091 * Spiffi slitlet order on the detector.
00092 * Also shifts the resulting image rows by one pixel if
00093 * necessary according to the distances array gained from
00094 * the north-south test routine.
00095 * Can do the same with the bad pixel map image to generate a
00096 * bad pixel mask cube.
00097 * 6) makes a data cube out of a resampled source image
00098 * this 3D specific routine takes into account the
00099 * 3D slitlet order on the detector.
00100 * Also shifts the resulting image rows by one pixel if
00101 * necessary according to the distances array gained from
00102 * the north-south test routine.
00103 * Can do the same with the bad pixel map image to generate a
00104 * bad pixel mask cube.
00105 * 7) makes a data cube out of a resampled source image
00106 * this MPE 3D specific routine takes into account the
00107 * 3D slitlet order on the detector.
00108 * Also shifts the resulting image row by an integer pixel shift if
00109 * necessary according to the input kpixshift array
00110 * Can do the same with the bad pixel map image to generate a
00111 * bad pixel mask cube.
00112 * 8) converts resampled bad pixels to real bad pixels in data cubes.
00113 * 9) Bad pixel interpolation 3D like (saturated pixels exist):
00114 * interpolates the bad pixels of the source cube by
00115 * using the nearest neighbors.
00116 * first it is checked if the bad pixel is interpolatable:
00117 * it is only interpolatable if the number of good pixels
00118 * in its spectrum of length 2*n_neighbors+1 exceeds 3 and
00119 * if there is at least one good pixel on either side of the
00120 * central pixel.
00121 * Afterwards good neighboring pixels are searched within the
00122 * image plane of the bad pixel by using an increasing pixel radius.
00123 * Good pixels mean, the corresponding spectral pixels of the
00124 * bad pixel and its spatial neighboring pixel must have
00125 * at least 2 valid pixel pairs to be able to be used for
00126 * the interpolation. The search is stopped if 9 valid neighboring
00127 * pixels are found.
00128 * Now normalize the found spectral values, collect the valid pixels
00129 * (there must be at least 18) and take the median of the valid
00130 * pixels with which the bad pixel is replaced.
00131 * 10) fine tunes each row in the right position according
00132 * to the distances of the slitlets to each other
00133 * (output of the north-south test).
00134 * This means that the rows must be realigned by a
00135 * fraction of a pixel to accomodate non-integer slit
00136 * length. The fractional realignment is done by using
00137 * the polynomial interpolation algorithm of N.R.
00138 * Each row is rescaled so that the total flux is
00139 * conserved.
00140 * 11) fine tunes each row in the right position according
00141 * to the distances of the slitlets to each other
00142 * (output of the north-south test).
00143 * This means that the rows must be realigned by a
00144 * fraction of a pixel to accomodate non-integer slit
00145 * length. The fractional realignment is done by using
00146 * the FFT algorithm four1() of N.R.
00147 * 12) fine tunes each row in the right position according
00148 * to the distances of the slitlets to each other
00149 * (output of the north-south test).
00150 * This means that the rows must be realigned by a
00151 * fraction of a pixel to accomodate non-integer slit
00152 * length. The fractional realignment is done by using
00153 * the spline interpolation algorithm splint in connection
00154 * with the algorithm spline of N.R.
00155 * This algorithms assume that each row is a tabulated
00156 * function. The first derivatives of the interpolating
00157 * function at the first and last point must be given.
00158 * These are set higher than 1xe^30, so the routine
00159 * sets the corresponding boundary condition for a natural
00160 * spline, with zero second derivative on that boundary.
00161 * Each row is rescaled so that the total flux is
00162 * conserved.
00163 *
00164 * FILES
00165 *
00166 * ENVIRONMENT
00167 *
00168 * RETURN VALUES
00169 *
00170 * CAUTIONS
00171 *
00172 * EXAMPLES
00173 *
00174 * SEE ALSO
00175 *
00176 * BUGS
00177 *
00178 *------------------------------------------------------------------------
00179 */
00180
00181 #define POSIX_SOURCE 1
00182 #include "vltPort.h"
00183
00184 /*
00185 * System Headers
00186 */
00187
00188 /*
00189 * Local Headers
00190 */
00191
00192 #include "cube_construct.h"
00193
00194 /*----------------------------------------------------------------------------
00195 * Function codes
00196 *--------------------------------------------------------------------------*/
00197
00198 /*----------------------------------------------------------------------------
00199 Function : convolveNSImageByGauss()
00200 In : lineImage: North-south-test image
00201 hw: kernel half width of the gaussian
00202 response function
00203 Out : north-south-test image convolved with a gaussian
00204 Job : convolves a north-south-test image with a gaussian
00205 with user given integer half width by using the eclipse
00206 routine function1d_filter_lowpass().
00207 ---------------------------------------------------------------------------*/
00208
00209 OneImage * convolveNSImageByGauss( OneImage * lineImage,
00210 int hw )
00211 {
00212 OneImage * returnImage ;
00213 float row_buffer[lineImage->ly] ;
00214 float * filter ;
00215 int col, row ;
00216
00217 if ( lineImage == NullImage )
00218 {
00219 cpl_msg_error("convolveImageByGauss:","no input image given!\n") ;
00220 return NullImage ;
00221 }
00222
00223 if ( hw < 1 )
00224 {
00225 cpl_msg_error("convolveImageByGauss:"," wrong half width given!\n") ;
00226 return NullImage ;
00227 }
00228
00229 /* allocate memory for returned image */
00230 if ( NULL == ( returnImage = new_image( lineImage -> lx, lineImage -> ly ) ))
00231 {
00232 cpl_msg_error("convolveImageByGauss:","cannot allocate a new image\n");
00233 return NullImage ;
00234 }
00235
00236 /* go through the image rows and save them in a buffer */
00237 for ( row = 0 ; row < lineImage -> ly ; row++ )
00238 {
00239 for ( col = 0 ; col < lineImage -> lx ; col++ )
00240 {
00241 if ( qfits_isnan(lineImage->data[col+row*lineImage->lx]) )
00242 {
00243 row_buffer[col] = 0. ;
00244 }
00245 else
00246 {
00247 row_buffer[col] = lineImage -> data[col + row*lineImage->lx] ;
00248 }
00249 }
00250
00251 /*---------------------------------------------------------------------
00252 * now low pass filter the rows by the gaussian and fill the return
00253 * image.
00254 */
00255 filter = function1d_filter_lowpass( row_buffer,
00256 lineImage -> lx,
00257 LOW_PASS_GAUSSIAN,
00258 hw ) ;
00259 for ( col = 0 ; col < lineImage -> ly ; col++ )
00260 {
00261 returnImage -> data[col + row*lineImage->lx] = filter[col] ;
00262 }
00263 /* deallocate memory */
00264 function1d_del (filter) ;
00265 }
00266
00267 return returnImage ;
00268 }
00269
00270 /*----------------------------------------------------------------------------
00271 Function : north_south_test()
00272 In : ns_image: north-south image
00273 n_slitlets: number of slitlets
00274 halfWidth: half width of the box in which the lines
00275 are fit by a gaussian
00276 fwhm: first guess of the full width at half maximum
00277 minDiff: minimum amplitude below which the fit
00278 will not be carried through
00279 estimated_dist: estimated average distance of spectra
00280 devtol: maximal pixel deviation of slitlet distances
00281 Out : array of the distances of the slitlets
00282 from each other
00283 Job : determines the distances of the slitlets
00284 ---------------------------------------------------------------------------*/
00285
00286 float * north_south_test( OneImage * ns_image,
00287 int n_slitlets,
00288 int halfWidth,
00289 float fwhm,
00290 float minDiff,
00291 float estimated_dist,
00292 float devtol,
00293 int bottom,
00294 int top )
00295 {
00296 int i, j, k, m, row, col, n, ni, na ;
00297 int position, counter, iters ;
00298 int xdim, ndat, its, numpar ;
00299 pixelvalue row_buf[ns_image -> lx] ;
00300 float sum, mean, maxval ;
00301 float tol, lab ;
00302 float * distances ;
00303 float distances_buf[ns_image -> ly][n_slitlets-1] ;
00304 float x_position[n_slitlets] ;
00305 float * xdat, * wdat ;
00306 int * mpar ;
00307 int found[3*n_slitlets], found_clean[3*n_slitlets] ;
00308 int found_cleanit[3*n_slitlets] ;
00309 Vector * line ;
00310 FitParams ** par ;
00311 int foundit, begin, end ;
00312 int zeroindicator ;
00313
00314 if ( ns_image == NullImage )
00315 {
00316 cpl_msg_error("north_south_test:","sorry, no image given\n") ;
00317 return NULL ;
00318 }
00319 if ( n_slitlets < 1 )
00320 {
00321 cpl_msg_error("north_south_test:","wrong number of slitlets given\n") ;
00322 return NULL ;
00323 }
00324 if ( halfWidth < 0 || halfWidth >= estimated_dist )
00325 {
00326 cpl_msg_error("north_south_test:","wrong half width given\n") ;
00327 return NULL ;
00328 }
00329 if ( fwhm <= 0. )
00330 {
00331 cpl_msg_error("north_south_test:","wrong fwhm given\n") ;
00332 return NULL ;
00333 }
00334 if ( minDiff < 1. )
00335 {
00336 cpl_msg_error("north_south_test:","wrong minDiff given\n") ;
00337 return NULL ;
00338 }
00339
00340 /* allocate memory for output array */
00341 if (NULL == (distances = (float *) cpl_calloc ( n_slitlets - 1 , sizeof (float) )))
00342 {
00343 cpl_msg_error("north_south_test:","could not allocate memory\n") ;
00344 return NULL ;
00345 }
00346
00347 /* go through the image rows */
00348 for ( row = bottom ; row < top ; row++ )
00349 {
00350 zeroindicator = 0 ;
00351
00352 /* initialize the distance buffer */
00353 for ( i = 0 ; i < n_slitlets-1 ; i++ )
00354 {
00355 distances_buf[row][i] = ZERO ;
00356 }
00357
00358 /* fill the row buffer array with image data */
00359 for ( col = 0 ; col < ns_image -> lx ; col++ )
00360 {
00361 row_buf[col] = ns_image -> data[col + row*ns_image->lx] ;
00362 }
00363
00364 /* determine the mean of the row data */
00365 sum = 0. ;
00366 n = 0 ;
00367 for ( i = 0 ; i < ns_image -> lx ; i++ )
00368 {
00369 if ( qfits_isnan(row_buf[i]) )
00370 {
00371 continue ;
00372 }
00373 sum += row_buf[i] ;
00374 n++ ;
00375 }
00376 mean = sum / (float)n ;
00377
00378
00379 /* store the positions of image values greater than the mean */
00380 n = 0 ;
00381 for ( i = 0 ; i < ns_image -> lx ; i++ )
00382 {
00383 if (qfits_isnan(row_buf[i]))
00384 {
00385 continue ;
00386 }
00387 if ( row_buf[i] > sqrt(mean*mean*9) )
00388 {
00389 found[n] = i ;
00390 n++ ;
00391 }
00392 }
00393
00394 if ( n < n_slitlets )
00395 {
00396 cpl_msg_warning("north_south_test1:","t1 wrong number of intensity columns found in row: %d, found number: %d, mean: %g\n", row, n, mean) ;
00397 continue ;
00398 }
00399 else
00400 {
00401 /* find the maximum value position around the found columns */
00402 na = 0 ;
00403 for ( i = 1 ; i < n ; i ++ )
00404 {
00405 if ( found[i] - found[i-1] < halfWidth )
00406 {
00407 begin = found[i] - halfWidth ;
00408 if ( begin < 0 )
00409 {
00410 begin = 0 ;
00411 }
00412 end = found[i] + halfWidth ;
00413 if ( end >= ns_image->lx )
00414 {
00415 end = ns_image->lx - 1 ;
00416 }
00417 /* find the maximum value inside the box around the found positions*/
00418 maxval = -FLT_MAX ;
00419 foundit = 0 ;
00420 for ( j = begin ; j <= end ; j++ )
00421 {
00422 /* do not consider boxes that contain bad pixels */
00423 if (qfits_isnan(row_buf[j]))
00424 {
00425 continue ;
00426 }
00427 if (row_buf[j] >= maxval )
00428 {
00429 maxval = row_buf[j] ;
00430 foundit = j ;
00431 }
00432 }
00433 if (maxval == -FLT_MAX)
00434 {
00435 continue ;
00436 }
00437 for ( k = 0 ; k < na ; k++ )
00438 {
00439 if ( found_cleanit[k] >= begin && found_cleanit[k] < foundit )
00440 {
00441 na-- ;
00442 }
00443 }
00444 for ( k = 0 ; k < n ; k++ )
00445 {
00446 if ( found[k] == foundit)
00447 {
00448 if (na>0){
00449 if ( found_cleanit[na-1] != found[k] )
00450 {
00451 found_cleanit[na] = found[k] ;
00452 na++ ;
00453 }
00454 }
00455 else{
00456 found_cleanit[na] = found[k] ;
00457 na++ ;
00458 }
00459 }
00460 }
00461 }
00462 else
00463 {
00464 if ( i == 1 )
00465 {
00466 found_cleanit[na] = found[0] ;
00467 na++ ;
00468 found_cleanit[na] = found[1] ;
00469 na++ ;
00470 }
00471 else
00472 {
00473 if (na>0){
00474 if ( found_cleanit[na-1] != found[i-1])
00475 {
00476 found_cleanit[na] = found[i-1] ;
00477 na++ ;
00478 }
00479 if ( found_cleanit[na-1] != found[i])
00480 {
00481 found_cleanit[na] = found[i] ;
00482 na++ ;
00483 }
00484 }
00485 else
00486 {
00487 found_cleanit[na] = found[i] ;
00488 na++ ;
00489 }
00490 }
00491 }
00492 }
00493 /* determine only one pixel position for each slitlet intensity */
00494 j = 1 ;
00495 for ( i = 1 ; i < na ; i++ )
00496 {
00497 if ( (float)(found_cleanit[i] - found_cleanit[i-1]) < (estimated_dist - devtol) ||
00498 (float)(found_cleanit[i] - found_cleanit[i-1]) > (estimated_dist + devtol) )
00499 {
00500 continue ;
00501 }
00502 else
00503 {
00504 found_clean[j-1] = found_cleanit[i-1] ;
00505 found_clean[j] = found_cleanit[i] ;
00506 j++ ;
00507 }
00508 }
00509 }
00510 if ( j > n_slitlets )
00511 {
00512 /* check the distance again */
00513 ni = 1 ;
00514 for ( i = 1 ; i < j ; i++ )
00515 {
00516 if ( (float)(found_clean[i] - found_clean[i-1]) < (estimated_dist - devtol ) ||
00517 (float)(found_clean[i] - found_clean[i-1]) > (estimated_dist + devtol ) )
00518 {
00519 continue ;
00520 }
00521 else
00522 {
00523
00524 found_clean[ni-1] = found_clean[i-1] ;
00525 found_clean[ni] = found_clean[i] ;
00526 ni++ ;
00527 }
00528 }
00529 if ( ni != n_slitlets )
00530 {
00531 cpl_msg_warning("north_south_test2:","t2 wrong number of intensity columns found in row: %d, found number: %d\n", row, ni) ;
00532 continue ;
00533 }
00534 else
00535 {
00536 j = ni ;
00537 }
00538 }
00539 else if ( j < n_slitlets )
00540 {
00541 cpl_msg_debug ("north_south_test3:","t3 wrong number of intensity columns found in row: %d , found number: %d, mean: %g\n", row, j, mean) ;
00542 continue ;
00543 }
00544 counter = 0 ;
00545 /* go through the found intensity pixels in one row */
00546 for ( i = 0 ; i < j ; i++ )
00547 {
00548 /* allocate memory for the array where the line is fitted in */
00549 if ( NULL == (line = newVector (2*halfWidth + 1)) )
00550 {
00551 cpl_msg_error ("north_south_test:","cannot allocate new Vector \n") ;
00552 cpl_free(distances) ;
00553 return NULL ;
00554 }
00555
00556 /* allocate memory */
00557 xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00558 wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00559 mpar = (int *) cpl_calloc( MAXPAR, sizeof (int) ) ;
00560 par = newFitParams(1) ;
00561
00562 m = 0 ;
00563 for ( k = found_clean[i]-halfWidth ; k <= found_clean[i]+halfWidth ; k++ )
00564 {
00565 if ( k < 0 )
00566 {
00567 k = 0. ;
00568 }
00569 else if ( k >= ns_image -> lx )
00570 {
00571 k = ns_image->lx - 1 ;
00572 }
00573 else if ( qfits_isnan(row_buf[k]) )
00574 {
00575 zeroindicator = 1 ;
00576 break ;
00577 }
00578 else
00579 {
00580 line -> data[m] = row_buf[k] ;
00581 m++ ;
00582 }
00583 }
00584 if ( zeroindicator == 1 )
00585 {
00586 destroyVector(line) ;
00587 cpl_free(xdat) ;
00588 cpl_free(wdat) ;
00589 cpl_free(mpar) ;
00590 destroyFitParams(par) ;
00591 break ;
00592 }
00593
00594 /*--------------------------------------------------------------------
00595 * go through the spectral vector
00596 * determine the maximum pixel value in the spectral vector
00597 */
00598 maxval = -FLT_MAX ;
00599 position = -INT32_MAX ;
00600 for ( k = 0 ; k < m ; k++ )
00601 {
00602 xdat[k] = k ;
00603 wdat[k] = 1.0 ;
00604 if ( line -> data[k] >= maxval )
00605 {
00606 maxval = line -> data[k] ;
00607 position = k ;
00608 }
00609 }
00610
00611 /* set initial values for the fitting routine */
00612 xdim = XDIM ;
00613 ndat = line -> n_elements ;
00614 numpar = MAXPAR ;
00615 tol = TOL ;
00616 lab = LAB ;
00617 its = ITS ;
00618 (*par) -> fit_par[1] = fwhm ;
00619 (*par) -> fit_par[2] = (float) position ;
00620 (*par) -> fit_par[3] = (float) (line -> data[0] + line -> data[line->n_elements - 1]) / 2.0 ;
00621 (*par) -> fit_par[0] = maxval - ((*par) -> fit_par[3]) ;
00622
00623
00624 /* exclude negative peaks and low signal cases */
00625 if ( (*par) -> fit_par[0] < minDiff )
00626 {
00627 cpl_msg_warning ("north_south_test:","sorry, signal of line too low to fit in row: %d in slitlet %d\n", row, i) ;
00628 destroyVector(line) ;
00629 cpl_free(xdat) ;
00630 cpl_free(wdat) ;
00631 cpl_free(mpar) ;
00632 destroyFitParams(par) ;
00633 continue ;
00634 }
00635
00636 for ( k = 0 ; k < MAXPAR ; k++ )
00637 {
00638 (*par) -> derv_par[k] = 0.0 ;
00639 mpar[k] = 1 ;
00640 }
00641 /* finally, do the least square fit using a gaussian */
00642 if ( 0 > ( iters = lsqfit_c( xdat, &xdim, line -> data, wdat, &ndat, (*par) -> fit_par,
00643 (*par) -> derv_par, mpar, &numpar, &tol, &its, &lab )) )
00644 {
00645 /*
00646 cpl_msg_debug ("north_south_test:","lsqfit_c: least squares fit failed, error no.: %d in row: %d in slitlet %d\n", iters, row, i) ;
00647 */
00648 destroyVector(line) ;
00649 cpl_free(xdat) ;
00650 cpl_free(wdat) ;
00651 cpl_free(mpar) ;
00652 destroyFitParams(par) ;
00653 continue ;
00654 }
00655
00656 /* check for negative fit results */
00657 if ( (*par) -> fit_par[0] <= 0. || (*par) -> fit_par[1] <= 0. ||
00658 (*par) -> fit_par[2] < 0. )
00659 {
00660 cpl_msg_warning ("north_south_test:","negative parameters as fit result, not used! in row %d in slitlet %d\n", row, i) ;
00661 destroyVector(line) ;
00662 cpl_free(xdat) ;
00663 cpl_free(wdat) ;
00664 cpl_free(mpar) ;
00665 destroyFitParams(par) ;
00666 continue ;
00667 }
00668
00669 /* correct the fitted position for the given row of the line in image coordinates */
00670 (*par) -> fit_par[2] = (float) (found_clean[i] - halfWidth) + (*par) -> fit_par[2] ;
00671 x_position[counter] = (*par) -> fit_par[2] ;
00672 counter ++ ;
00673
00674 /* free memory */
00675 destroyFitParams(par) ;
00676 destroyVector ( line ) ;
00677 cpl_free ( xdat ) ;
00678 cpl_free ( wdat ) ;
00679 cpl_free ( mpar ) ;
00680 }
00681 if (zeroindicator == 1)
00682 {
00683 cpl_msg_debug ("north_south_test:","bad pixel in fitting box in row: %d\n", row) ;
00684 continue ;
00685 }
00686
00687 if ( counter != n_slitlets )
00688 {
00689 continue ;
00690 cpl_msg_warning("north_south_test:","wrong number of slitlets found in row: %d\n", row) ;
00691 }
00692 /* store the distances between the sources in a buffer */
00693 for ( i = 1 ; i < n_slitlets ; i++ )
00694 {
00695 distances_buf[row][i-1] = x_position[i] - x_position[i-1] ;
00696 }
00697 }
00698
00699 /* ----------------------------------------------------------------
00700 * go through the rows again and take the mean of the distances,
00701 * throw away the runaways
00702 */
00703 for ( i = 0 ; i < n_slitlets-1 ; i++ )
00704 {
00705 n = 0 ;
00706 sum = 0. ;
00707 for ( row = bottom ; row < top ; row++ )
00708 {
00709 if ( fabs( distances_buf[row][i] - estimated_dist ) > devtol ||
00710 qfits_isnan(distances_buf[row][i]) )
00711 {
00712 /*
00713 printf("check1=%g ref1=%g check3=%d\n",
00714 fabs( distances_buf[row][i] - estimated_dist ),
00715 devtol,
00716 qfits_isnan(distances_buf[row][i]));
00717 */
00718 continue ;
00719 }
00720 sum += distances_buf[row][i] ;
00721 n++ ;
00722 }
00723 if ( n < 2 )
00724 {
00725 cpl_msg_error( "north_south_test:","distances array could not be determined completely!, \
00726 deviations of distances from number of slitlets too big\n" ) ;
00727 cpl_free(distances) ;
00728 return NULL ;
00729 }
00730 else
00731 {
00732 distances[i] = sum / (float)n ;
00733 }
00734 }
00735 return distances ;
00736 }
00737
00738
00739 /*----------------------------------------------------------------------------
00740 Function : makeCube()
00741 In : calibImage: resampled source image
00742 distances: distances of the slitlets from each other
00743 output of function ns_test
00744 correct_diff_dist: dummy array with 32 elements
00745 Out : resulting source data cube
00746 correct_diff_dist: differences of the slitlets from
00747 distance 32 given in the correct
00748 Spiffi row sequence. The first slitlet
00749 is the reference, therefore element
00750 23 is set 0.
00751 Job : makes a data cube out of a resampled source image
00752 this SPIFFI specific routine takes into account the
00753 Spiffi slitlet order on the detector.
00754 Also shifts the resulting image rows by one pixel if
00755 necessary according to the distances array gained from
00756 the north-south test routine.
00757 Can do the same with the bad pixel map image to generate a
00758 bad pixel mask cube.
00759 ---------------------------------------------------------------------------*/
00760
00761 OneCube * makeCube ( OneImage * calibImage,
00762 float * distances,
00763 float * correct_diff_dist )
00764 {
00765 OneCube * returnCube ;
00766 int imsize, kslit, kpix ;
00767 int slit_index ;
00768 int z, col, recol ;
00769
00770 if ( NullImage == calibImage )
00771 {
00772 cpl_msg_error("makeCube:","no resampled image given!\n") ;
00773 return NullCube ;
00774 }
00775
00776 if ( NULL == distances )
00777 {
00778 cpl_msg_error("makeCube:","no distances array from ns_test given!/n") ;
00779 return NullCube ;
00780 }
00781
00782 if ( NULL == correct_diff_dist )
00783 {
00784 cpl_msg_error("makeCube:","correct_diff_dist array is not allocated!/n") ;
00785 return NullCube ;
00786 }
00787
00788 if ( N_SLITLETS != 32 )
00789 {
00790 cpl_msg_error ("makeCube:","wrong number of slitlets given \n" ) ;
00791 return NullCube ;
00792 }
00793 imsize = calibImage->lx / N_SLITLETS ;
00794
00795 /* allocate memory */
00796 if ( NULL == (returnCube = newCube(imsize, N_SLITLETS, calibImage->ly)) )
00797 {
00798 cpl_msg_error ("makeCube:","cannot allocate new cube \n" ) ;
00799 return NullCube ;
00800 }
00801
00802 /* now build the data cube out of the resampled image */
00803 for ( z = 0 ; z < calibImage -> ly ; z++ ) /* go through the z-axis */
00804 {
00805 kpix = 0 ;
00806 kslit = 0 ;
00807 slit_index = -1 ;
00808 recol = -1 ;
00809 for ( col = 0 ; col < calibImage->lx ; col++ ) /* go through the image columns */
00810 {
00811 if ( col % imsize == 0 )
00812 {
00813 recol = 0 ;
00814 kslit = col/imsize ;
00815 /* sort the slitlets in the right spiffi specific way */
00816 switch (kslit)
00817 {
00818 case 0:
00819 slit_index = 8 ;
00820 break ;
00821 case 1:
00822 slit_index = 7 ;
00823 break ;
00824 case 2:
00825 slit_index = 9 ;
00826 break ;
00827 case 3:
00828 slit_index = 6 ;
00829 break ;
00830 case 4:
00831 slit_index = 10 ;
00832 break ;
00833 case 5:
00834 slit_index = 5 ;
00835 break ;
00836 case 6:
00837 slit_index = 11 ;
00838 break ;
00839 case 7:
00840 slit_index = 4 ;
00841 break ;
00842 case 8:
00843 slit_index = 12 ;
00844 break ;
00845 case 9:
00846 slit_index = 3 ;
00847 break ;
00848 case 10:
00849 slit_index = 13 ;
00850 break ;
00851 case 11:
00852 slit_index = 2 ;
00853 break ;
00854 case 12:
00855 slit_index = 14 ;
00856 break ;
00857 case 13:
00858 slit_index = 1 ;
00859 break ;
00860 case 14:
00861 slit_index = 15 ;
00862 break ;
00863 case 15:
00864 slit_index = 0 ;
00865 break ;
00866 case 16:
00867 slit_index = 31 ;
00868 break ;
00869 case 17:
00870 slit_index = 16 ;
00871 break ;
00872 case 18:
00873 slit_index = 30 ;
00874 break ;
00875 case 19:
00876 slit_index = 17 ;
00877 break ;
00878 case 20:
00879 slit_index = 29 ;
00880 break ;
00881 case 21:
00882 slit_index = 18 ;
00883 break ;
00884 case 22:
00885 slit_index = 28 ;
00886 break ;
00887 case 23:
00888 slit_index = 19 ;
00889 break ;
00890 case 24:
00891 slit_index = 27 ;
00892 break ;
00893 case 25:
00894 slit_index = 20 ;
00895 break ;
00896 case 26:
00897 slit_index = 26 ;
00898 break ;
00899 case 27:
00900 slit_index = 21 ;
00901 break ;
00902 case 28:
00903 slit_index = 25 ;
00904 break ;
00905 case 29:
00906 slit_index = 22 ;
00907 break ;
00908 case 30:
00909 slit_index = 24 ;
00910 break ;
00911 case 31:
00912 slit_index = 23 ;
00913 break ;
00914 default:
00915 cpl_msg_error("makeCube:","wrong slitlet index: couldn't be a spiffi image, \
00916 there must be 32 slitlets!\n") ;
00917 destroy_cube(returnCube) ;
00918 return NullCube ;
00919 break ;
00920 }
00921
00922 if ( kslit != 0 )
00923 {
00924 /*------------------------------------------------------------------
00925 * shift the first pixel by an integer if the absolute amount of distances[]
00926 * is bigger than 0.5
00927 */
00928 kpix = nint(distances[kslit-1]) ;
00929
00930 /*-----------------------------------------------------------------
00931 * now sort the distances array according to the row order and
00932 * add a 0 value for the first (reference) slitlet that means row 8
00933 */
00934 correct_diff_dist[slit_index] = distances[kslit-1] - (float)kpix ;
00935 }
00936 /* refer all distances to the first slitlet */
00937 else
00938 {
00939 correct_diff_dist[slit_index] = 0. ;
00940 }
00941 }
00942
00943 /* fill each cube plane with one image row */
00944 returnCube->plane[z]->data[recol+slit_index*imsize] =
00945 calibImage->data[col+kpix+z*calibImage->lx] ;
00946 recol++ ;
00947
00948 if ( recol > imsize )
00949 {
00950 cpl_msg_error("makeCube:","wrong column of reconstructed image, shouldn't happen!\n") ;
00951 destroy_cube(returnCube) ;
00952 return NullCube ;
00953 }
00954 }
00955 }
00956 return returnCube ;
00957 }
00958
00959 /*----------------------------------------------------------------------------
00960 Function : makeCubeSpi()
00961 In : calibImage: resampled source image
00962 slit_edges: absolute beginning and ending positions of
00963 slitlet, output of fitSlits().
00964 shift: sub_pixel shifts referred to the reference slit
00965 edge
00966 Out : resulting source data cube
00967 shift: sub_pixel shifts referred to the reference slit
00968 edge
00969 Job : makes a data cube out of a resampled source image
00970 this SPIFFI specific routine takes into account the
00971 Spiffi slitlet order on the detector.
00972 This routine takes fitted slitlet positions into account.
00973 Can do the same with the bad pixel map image to generate a
00974 bad pixel mask cube.
00975 ---------------------------------------------------------------------------*/
00976
00977 OneCube * makeCubeSpi ( OneImage * calibImage,
00978 float ** slit_edges,
00979 float * shift )
00980 {
00981 OneCube * returnCube ;
00982 float diff, start ;
00983 float * center ;
00984 int * row_index ;
00985 int slit ;
00986 int col, z ;
00987 int imsize ;
00988 int * beginCol ;
00989 int col_counter ;
00990
00991 if ( NullImage == calibImage )
00992 {
00993 cpl_msg_error("makeCubeSpi:","no resampled image given!\n") ;
00994 return NullCube ;
00995 }
00996
00997 if ( NULL == slit_edges )
00998 {
00999 cpl_msg_error("makeCubeSpi:","no slit_edges array given from fitSlits()!/n") ;
01000 return NullCube ;
01001 }
01002
01003 if ( N_SLITLETS != 32 )
01004 {
01005 cpl_msg_error ("makeCube:","wrong number of slitlets given \n" ) ;
01006 return NullCube ;
01007 }
01008 imsize = calibImage->lx / N_SLITLETS ;
01009
01010 /* allocate memory */
01011 if ( NULL == (row_index = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01012 {
01013 cpl_msg_error ("makeCubeSpi:","cannot allocate memory \n" ) ;
01014 return NullCube ;
01015 }
01016 if ( NULL == (beginCol = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01017 {
01018 cpl_msg_error ("makeCubeSpi:","cannot allocate memory \n" ) ;
01019 cpl_free(row_index) ;
01020 return NullCube ;
01021 }
01022 if ( NULL == (center = (float*) cpl_calloc(N_SLITLETS, sizeof(float)) ) )
01023 {
01024 cpl_msg_error ("makeCubeSpi:","cannot allocate memory \n" ) ;
01025 cpl_free (row_index) ;
01026 cpl_free (beginCol) ;
01027 return NullCube ;
01028 }
01029 if ( NULL == (returnCube = newCube(imsize, N_SLITLETS, calibImage->ly)) )
01030 {
01031 cpl_msg_error ("makeCubeSpi:","cannot allocate new cube \n" ) ;
01032 cpl_free (row_index) ;
01033 cpl_free (beginCol) ;
01034 cpl_free (center) ;
01035 return NullCube ;
01036 }
01037 /* determine the absolute center of the slitlets and the distances inside the image*/
01038 for ( slit = 0 ; slit < N_SLITLETS ; slit++ ) /* go through the slitlets of each
01039 row of the resampled image */
01040 {
01041 center[slit] = (slit_edges[slit][1] + slit_edges[slit][0]) / 2. ;
01042 /* -------------------------------------------------------------------------
01043 * sort the slitlets in the right spiffi specific way
01044 * the row_index describes the row index of the current slitlet
01045 * in the resulting cube images.
01046 */
01047 switch (slit)
01048 {
01049 case 0:
01050 row_index[0] = 8 ;
01051 break ;
01052 case 1:
01053 row_index[1] = 7 ;
01054 break ;
01055 case 2:
01056 row_index[2] = 9 ;
01057 break ;
01058 case 3:
01059 row_index[3] = 6 ;
01060 break ;
01061 case 4:
01062 row_index[4] = 10 ;
01063 break ;
01064 case 5:
01065 row_index[5] = 5 ;
01066 break ;
01067 case 6:
01068 row_index[6] = 11 ;
01069 break ;
01070 case 7:
01071 row_index[7] = 4 ;
01072 break ;
01073 case 8:
01074 row_index[8] = 12 ;
01075 break ;
01076 case 9:
01077 row_index[9] = 3 ;
01078 break ;
01079 case 10:
01080 row_index[10] = 13 ;
01081 break ;
01082 case 11:
01083 row_index[11] = 2 ;
01084 break ;
01085 case 12:
01086 row_index[12] = 14 ;
01087 break ;
01088 case 13:
01089 row_index[13] = 1 ;
01090 break ;
01091 case 14:
01092 row_index[14] = 15 ;
01093 break ;
01094 case 15:
01095 row_index[15] = 0 ;
01096 break ;
01097 case 16:
01098 row_index[16] = 31 ;
01099 break ;
01100 case 17:
01101 row_index[17] = 16 ;
01102 break ;
01103 case 18:
01104 row_index[18] = 30 ;
01105 break ;
01106 case 19:
01107 row_index[19] = 17 ;
01108 break ;
01109 case 20:
01110 row_index[20] = 29 ;
01111 break ;
01112 case 21:
01113 row_index[21] = 18 ;
01114 break ;
01115 case 22:
01116 row_index[22] = 28 ;
01117 break ;
01118 case 23:
01119 row_index[23] = 19 ;
01120 break ;
01121 case 24:
01122 row_index[24] = 27 ;
01123 break ;
01124 case 25:
01125 row_index[25] = 20 ;
01126 break ;
01127 case 26:
01128 row_index[26] = 26 ;
01129 break ;
01130 case 27:
01131 row_index[27] = 21 ;
01132 break ;
01133 case 28:
01134 row_index[28] = 25 ;
01135 break ;
01136 case 29:
01137 row_index[29] = 22 ;
01138 break ;
01139 case 30:
01140 row_index[30] = 24 ;
01141 break ;
01142 case 31:
01143 row_index[31] = 23 ;
01144 break ;
01145 default:
01146 cpl_msg_error("makeCubeSpi:","wrong slitlet index: couldn't be a spiffi image, there must be 32 slitlets!\n") ;
01147 destroy_cube(returnCube) ;
01148 cpl_free (row_index) ;
01149 cpl_free (beginCol) ;
01150 cpl_free (center) ;
01151 return NullCube ;
01152 }
01153 /* determine the integer column on which the slitlet starts, center the
01154 slitlet on the image row */
01155 start = center[slit] - (float) (imsize - 1)/2. ;
01156 beginCol[slit] = nint (start) ;
01157 /* determine the error of using integer pixels */
01158 diff = start - (float)beginCol[slit] ;
01159
01160 /*-----------------------------------------------------------------------------------
01161 * determine the output shift values by which the rows are finally shifted,
01162 * consider the integer pixel errors
01163 * resort shift array to get the row index
01164 */
01165 shift[row_index[slit]] = diff ;
01166 }
01167
01168 /* now build the data cube out of the resampled image */
01169 for ( z = 0 ; z < calibImage -> ly ; z++ ) /* go through the z-axis */
01170 {
01171 for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01172 {
01173 col_counter = beginCol[slit] ;
01174 /* each slitlet is centered on the final image row */
01175 for ( col = 0 ; col < imsize ; col++ )
01176 {
01177 if ( col_counter > calibImage->lx-1 )
01178 {
01179 col_counter-- ;
01180 }
01181 if ( col_counter + z*calibImage->lx < 0 )
01182 {
01183 returnCube->plane[z]->data[col+row_index[slit]*imsize] =
01184 calibImage->data[0] ;
01185 }
01186 else
01187 {
01188 returnCube->plane[z]->data[col+row_index[slit]*imsize] =
01189 calibImage->data[col_counter+z*calibImage->lx] ;
01190 }
01191
01192 col_counter++ ;
01193 }
01194 }
01195 }
01196 cpl_free (row_index) ;
01197 cpl_free (beginCol) ;
01198 cpl_free (center) ;
01199
01200 return returnCube ;
01201 }
01202
01203 /*----------------------------------------------------------------------------
01204 Function : makeCubeDist()
01205 In : calibImage: resampled source image
01206 firstCol: floating point value of the first column of
01207 the first slitlet in the resampled image,
01208 determined "by hand"
01209 distances: distances of the slitlets from each other
01210 output of function ns_test
01211 shift: dummy array with 32 elements
01212 Out : resulting source data cube
01213 shift: differences of the slitlets from
01214 distance 32 given in the correct
01215 Spiffi row sequence. The first slitlet
01216 is the reference, therefore element
01217 23 is set 0.
01218 Job : makes a data cube out of a resampled source image
01219 this SPIFFI specific routine takes into account the
01220 Spiffi slitlet order on the detector.
01221 Also shifts the resulting image rows by one pixel if
01222 necessary according to the distances array gained from
01223 the north-south test routine.
01224 Can do the same with the bad pixel map image to generate a
01225 bad pixel mask cube.
01226 ---------------------------------------------------------------------------*/
01227
01228 OneCube * makeCubeDist ( OneImage * calibImage,
01229 float firstCol,
01230 float * distances,
01231 float * shift )
01232 {
01233 OneCube * returnCube ;
01234 float di ;
01235 float diff, start ;
01236 int * row_index ;
01237 int slit ;
01238 int col, z ;
01239 int imsize ;
01240 int * beginCol ;
01241 int col_counter ;
01242
01243 if ( NullImage == calibImage )
01244 {
01245 cpl_msg_error("makeCubeDist:"," no resampled image given!\n") ;
01246 return NullCube ;
01247 }
01248 if ( NULL == distances )
01249 {
01250 cpl_msg_error("makeCubeDist:","no distances array given from north_south_test()!\n") ;
01251 return NullCube ;
01252 }
01253
01254 if ( N_SLITLETS != 32 )
01255 {
01256 cpl_msg_error ("makeCube:","wrong number of slitlets given \n" ) ;
01257 return NullCube ;
01258 }
01259 imsize = calibImage->lx / N_SLITLETS ;
01260
01261 /* allocate memory */
01262 if ( NULL == (row_index = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01263 {
01264 cpl_msg_error ("makeCubeDist:","cannot allocate memory \n" ) ;
01265 return NullCube ;
01266 }
01267 if ( NULL == (beginCol = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01268 {
01269 cpl_msg_error ("makeCubeDist:","cannot allocate memory \n" ) ;
01270 cpl_free(row_index) ;
01271 return NullCube ;
01272 }
01273 if ( NULL == (returnCube = newCube(imsize, N_SLITLETS, calibImage->ly)) )
01274 {
01275 cpl_msg_error ("makeCubeDist:","cannot allocate new cube \n" ) ;
01276 cpl_free(row_index) ;
01277 cpl_free(beginCol) ;
01278 return NullCube ;
01279 }
01280
01281 di = 0. ;
01282 /* determine the absolute beginning of the slitlets and the distances inside the image*/
01283 for ( slit = 0 ; slit < N_SLITLETS ; slit++ ) /* go through the slitlets of
01284 each row of the resampled image */
01285 {
01286
01287 /* -------------------------------------------------------------------------
01288 * sort the slitlets in the right spiffi specific way
01289 * the row_index describes the row index of the current slitlet
01290 * in the resulting cube images.
01291 */
01292 switch (slit)
01293 {
01294 case 0:
01295 row_index[0] = 8 ;
01296 break ;
01297 case 1:
01298 row_index[1] = 7 ;
01299 break ;
01300 case 2:
01301 row_index[2] = 9 ;
01302 break ;
01303 case 3:
01304 row_index[3] = 6 ;
01305 break ;
01306 case 4:
01307 row_index[4] = 10 ;
01308 break ;
01309 case 5:
01310 row_index[5] = 5 ;
01311 break ;
01312 case 6:
01313 row_index[6] = 11 ;
01314 break ;
01315 case 7:
01316 row_index[7] = 4 ;
01317 break ;
01318 case 8:
01319 row_index[8] = 12 ;
01320 break ;
01321 case 9:
01322 row_index[9] = 3 ;
01323 break ;
01324 case 10:
01325 row_index[10] = 13 ;
01326 break ;
01327 case 11:
01328 row_index[11] = 2 ;
01329 break ;
01330 case 12:
01331 row_index[12] = 14 ;
01332 break ;
01333 case 13:
01334 row_index[13] = 1 ;
01335 break ;
01336 case 14:
01337 row_index[14] = 15 ;
01338 break ;
01339 case 15:
01340 row_index[15] = 0 ;
01341 break ;
01342 case 16:
01343 row_index[16] = 31 ;
01344 break ;
01345 case 17:
01346 row_index[17] = 16 ;
01347 break ;
01348 case 18:
01349 row_index[18] = 30 ;
01350 break ;
01351 case 19:
01352 row_index[19] = 17 ;
01353 break ;
01354 case 20:
01355 row_index[20] = 29 ;
01356 break ;
01357 case 21:
01358 row_index[21] = 18 ;
01359 break ;
01360 case 22:
01361 row_index[22] = 28 ;
01362 break ;
01363 case 23:
01364 row_index[23] = 19 ;
01365 break ;
01366 case 24:
01367 row_index[24] = 27 ;
01368 break ;
01369 case 25:
01370 row_index[25] = 20 ;
01371 break ;
01372 case 26:
01373 row_index[26] = 26 ;
01374 break ;
01375 case 27:
01376 row_index[27] = 21 ;
01377 break ;
01378 case 28:
01379 row_index[28] = 25 ;
01380 break ;
01381 case 29:
01382 row_index[29] = 22 ;
01383 break ;
01384 case 30:
01385 row_index[30] = 24 ;
01386 break ;
01387 case 31:
01388 row_index[31] = 23 ;
01389 break ;
01390 default:
01391 cpl_msg_error("makeCubeDist:","wrong slitlet index: couldn't be a spiffi image, there must be 32 slitlets!\n") ;
01392 destroy_cube(returnCube) ;
01393 cpl_free(row_index) ;
01394 cpl_free(beginCol) ;
01395 return NullCube ;
01396 }
01397
01398 /* determine the integer column on which the slitlet starts */
01399 if ( slit == 0 )
01400 {
01401 start = firstCol ;
01402 }
01403 else
01404 {
01405 di += distances[slit-1] ;
01406 start = firstCol + di ;
01407 }
01408 beginCol[slit] = nint(start) ;
01409
01410 /* determine the error of using integer pixels, its always smaller than 1 */
01411 diff = start - (float)beginCol[slit] ;
01412
01413 /*-----------------------------------------------------------------------------------
01414 * determine the output shift values by which the rows are finally shifted,
01415 * consider the integer pixel errors and resort shift array to get the row index
01416 */
01417 shift[row_index[slit]] = diff ;
01418 }
01419
01420 /* now build the data cube out of the resampled image */
01421 for ( z = 0 ; z < calibImage -> ly ; z++ ) /* go through the z-axis */
01422 {
01423 for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01424 {
01425 col_counter = beginCol[slit] ;
01426 /* each slitlet is centered on the final image row */
01427 for ( col = 0 ; col < imsize ; col++ )
01428 {
01429 if ( col_counter > calibImage->lx-1 )
01430 {
01431 col_counter-- ;
01432 }
01433 if ( col_counter + z*calibImage->lx < 0 )
01434 {
01435 returnCube->plane[z]->data[col+row_index[slit]*imsize] =
01436 calibImage->data[0] ;
01437 }
01438 else
01439 {
01440
01441 returnCube->plane[z]->data[col+row_index[slit]*imsize] =
01442 calibImage->data[col_counter+z*calibImage->lx] ;
01443 }
01444
01445 col_counter++ ;
01446 }
01447 }
01448 }
01449 cpl_free (row_index) ;
01450 cpl_free (beginCol) ;
01451
01452 return returnCube ;
01453 }
01454
01455 /*----------------------------------------------------------------------------
01456 Function : make3DCubeDist()
01457 In : calibImage: resampled source image
01458 firstCol: floating point value of the first column of
01459 the first slitlet in the resampled image,
01460 determined "by hand"
01461 distances: distances of the slitlets from each other
01462 output of function ns_test
01463 shift: dummy array with 32 elements
01464 Out : resulting source data cube
01465 shift: differences of the slitlets from
01466 distance 32 given in the correct
01467 Spiffi row sequence. The first slitlet
01468 is the reference, therefore element
01469 23 is set 0.
01470 Job : makes a data cube out of a resampled source image
01471 this 3D specific routine takes into account the
01472 3D slitlet order on the detector.
01473 Also shifts the resulting image rows by one pixel if
01474 necessary according to the distances array gained from
01475 the north-south test routine.
01476 Can do the same with the bad pixel map image to generate a
01477 bad pixel mask cube.
01478 ---------------------------------------------------------------------------*/
01479
01480 OneCube * make3DCubeDist ( OneImage * calibImage,
01481 float firstCol,
01482 float * distances,
01483 float * shift )
01484 {
01485 OneCube * returnCube ;
01486 float di ;
01487 float diff, start ;
01488 int * row_index ;
01489 int slit ;
01490 int col, z ;
01491 int imsize ;
01492 int * beginCol ;
01493 int col_counter ;
01494
01495 if ( NullImage == calibImage )
01496 {
01497 cpl_msg_error("makeCubeDist:"," no resampled image given!\n") ;
01498 return NullCube ;
01499 }
01500 if ( NULL == distances )
01501 {
01502 cpl_msg_error("makeCubeDist:","no distances array given from north_south_test()!\n") ;
01503 return NullCube ;
01504 }
01505
01506 if ( N_SLITLETS != 16 )
01507 {
01508 cpl_msg_error ("makeCube:","wrong number of slitlets given \n" ) ;
01509 return NullCube ;
01510 }
01511 imsize = calibImage->lx / N_SLITLETS ;
01512
01513 /* allocate memory */
01514 if ( NULL == (row_index = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01515 {
01516 cpl_msg_error ("makeCubeDist:","cannot allocate memory \n" ) ;
01517 return NullCube ;
01518 }
01519 if ( NULL == (beginCol = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01520 {
01521 cpl_msg_error ("makeCubeDist:","cannot allocate memory \n" ) ;
01522 cpl_free(row_index) ;
01523 return NullCube ;
01524 }
01525 if ( NULL == (returnCube = newCube(imsize, N_SLITLETS, calibImage->ly)) )
01526 {
01527 cpl_msg_error ("makeCubeDist:","cannot allocate new cube \n" ) ;
01528 cpl_free(row_index) ;
01529 cpl_free(beginCol) ;
01530 return NullCube ;
01531 }
01532
01533 di = 0. ;
01534 /* determine the absolute beginning of the slitlets and the distances inside the image*/
01535 for ( slit = 0 ; slit < N_SLITLETS ; slit++ ) /* go through the slitlets of
01536 each row of the resampled image */
01537 {
01538
01539 /* -------------------------------------------------------------------------
01540 * sort the slitlets in the right 3D specific way
01541 * the row_index describes the row index of the current slitlet
01542 * in the resulting cube images.
01543 */
01544 row_index[slit] = slit ;
01545
01546 /* determine the integer column on which the slitlet starts */
01547 if ( slit == 0 )
01548 {
01549 start = firstCol ;
01550 }
01551 else
01552 {
01553 di += distances[slit-1] ;
01554 start = firstCol + di ;
01555 }
01556 beginCol[slit] = nint(start) ;
01557
01558 /* determine the error of using integer pixels, its always smaller than 1 */
01559 diff = start - (float)beginCol[slit] ;
01560
01561 /*-----------------------------------------------------------------------------------
01562 * determine the output shift values by which the rows are finally shifted,
01563 * consider the integer pixel errors and resort shift array to get the row index
01564 */
01565 shift[row_index[slit]] = diff ;
01566 }
01567
01568 /* now build the data cube out of the resampled image */
01569 for ( z = 0 ; z < calibImage -> ly ; z++ ) /* go through the z-axis */
01570 {
01571 for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01572 {
01573 col_counter = beginCol[slit] ;
01574 /* each slitlet is centered on the final image row */
01575 for ( col = 0 ; col < imsize ; col++ )
01576 {
01577 if ( col_counter > calibImage->lx-1 )
01578 {
01579 col_counter-- ;
01580 }
01581 returnCube->plane[z]->data[col+row_index[slit]*imsize] =
01582 calibImage->data[col_counter+z*calibImage->lx] ;
01583 col_counter++ ;
01584 }
01585 }
01586 }
01587 cpl_free (row_index) ;
01588 cpl_free (beginCol) ;
01589
01590 return returnCube ;
01591 }
01592
01593
01594 /*----------------------------------------------------------------------------
01595 Function : make3DCube()
01596 In : calibImage: resampled source image
01597 kpixshift: integer pixel shifts of the resulting image
01598 plane rows.
01599 kpixfirst: first valid pixel
01600 Out : resulting source data cube
01601 Job : makes a data cube out of a resampled source image
01602 this MPE 3D specific routine takes into account the
01603 3D slitlet order on the detector.
01604 Also shifts the resulting image row by an integer pixel shift if
01605 necessary according to the input kpixshift array
01606 Can do the same with the bad pixel map image to generate a
01607 bad pixel mask cube.
01608 ---------------------------------------------------------------------------*/
01609
01610 OneCube * make3DCube ( OneImage * calibImage,
01611 int * kpixshift,
01612 int kpixfirst )
01613 {
01614 OneCube * returnCube ;
01615 int imsize, kslit, kpix ;
01616 int z, col, recol ;
01617
01618 if ( NullImage == calibImage )
01619 {
01620 cpl_msg_error("make3DCube:","no resampled image given!\n") ;
01621 return NullCube ;
01622 }
01623
01624 if ( NULL == kpixshift )
01625 {
01626 cpl_msg_error("make3DCube:","no shift array given!/n") ;
01627 return NullCube ;
01628 }
01629
01630 if ( kpixfirst < 0 )
01631 {
01632 cpl_msg_error("make3DCube:","wrong first valid pixel given!/n") ;
01633 return NullCube ;
01634 }
01635
01636 if ( N_SLITLETS != 16 )
01637 {
01638 cpl_msg_error ("makeCube:","wrong number of slitlets given \n" ) ;
01639 return NullCube ;
01640 }
01641 imsize = calibImage->lx / N_SLITLETS ;
01642
01643 if ( NULL == (returnCube = newCube(imsize, N_SLITLETS, calibImage->ly)) )
01644 {
01645 cpl_msg_error ("make3DCube:","cannot allocate new cube \n" ) ;
01646 return NullCube ;
01647 }
01648
01649 /* now build the data cube out of the resampled image */
01650 for ( z = 0 ; z < calibImage -> ly ; z++ ) /* go through the z-axis */
01651 {
01652 kpix = 0 ;
01653 kslit = 0 ;
01654 recol = -1 ;
01655 for ( col = 0 ; col < calibImage->lx ; col++ ) /* go through the image columns */
01656 {
01657 if ( col % imsize == 0 )
01658 {
01659 recol = 0 ;
01660 kslit = col/imsize ;
01661 kpix = kpixfirst + kpixshift[kslit] ;
01662 }
01663
01664 /* fill each cube plane with one image row */
01665 returnCube->plane[z]->data[recol+kslit*imsize] =
01666 calibImage->data[col+kpix+z*calibImage->lx] ;
01667 recol++ ;
01668 if ( recol > imsize )
01669 {
01670 cpl_msg_error("make3DCube:","wrong column of reconstructed image, shouldn't happen!\n") ;
01671 destroy_cube(returnCube) ;
01672 return NullCube ;
01673 }
01674 }
01675 }
01676 return returnCube ;
01677 }
01678
01679 /*----------------------------------------------------------------------------
01680 Function : determineMaskCube()
01681 In : sourceMaskCube: bad pixel mask cube generated by using
01682 the bad pixel mask frame (0: bad, 1: good)
01683 and going through
01684 the same reduction steps as with the
01685 observation frame.
01686 lowLimit: low limit of pixel value (about -0.7)
01687 highLimit: high limit of bad pixel value (about 0.7)
01688 Out : resulting bad pixel data cube (bad pixels: 0, good pixels: 1).
01689 Job : converts resampled bad pixels to real bad pixels in data cubes.
01690 ---------------------------------------------------------------------------*/
01691
01692 OneCube * determineMaskCube ( OneCube * sourceMaskCube,
01693 float lowLimit,
01694 float highLimit )
01695 {
01696 OneCube * retCube ;
01697 int z, n ;
01698
01699 if ( sourceMaskCube == NullCube )
01700 {
01701 cpl_msg_error("determineMaskCube:","no cube given!\n") ;
01702 return NullCube ;
01703 }
01704 if ( lowLimit > 0. )
01705 {
01706 cpl_msg_error("determineMaskCube:","lowLimit wrong!\n") ;
01707 return NullCube ;
01708 }
01709 if ( highLimit >= 1. || highLimit < 0. )
01710 {
01711 cpl_msg_error("determineMaskCube:","highLimit wrong!\n") ;
01712 return NullCube ;
01713 }
01714
01715 retCube = copy_cube (sourceMaskCube) ;
01716
01717 for ( z = 0 ; z < retCube -> np ; z++ )
01718 {
01719 for ( n = 0 ; n < (int) retCube -> plane[z] -> nbpix; n++ )
01720 {
01721 if ( retCube -> plane[z] -> data[n] == 0. )
01722 {
01723 continue ;
01724 }
01725 if ( retCube -> plane[z] -> data[n] == 1. )
01726 {
01727 continue ;
01728 }
01729 if ( retCube -> plane[z] -> data[n] >= lowLimit &&
01730 retCube -> plane[z] -> data[n] <= highLimit )
01731 {
01732 retCube -> plane[z] -> data[n] = 0. ;
01733 }
01734 else
01735 {
01736 retCube -> plane[z] -> data[n] = 1. ;
01737 }
01738 }
01739 }
01740 return retCube ;
01741 }
01742
01743 /*----------------------------------------------------------------------------
01744 Function : interpolCube()
01745 In : sourceCube: reconstructed source cube from makeCube
01746 without fine tuning of rows
01747 maskCube: bad pixel mask cube, bad pixel are marked
01748 with 0., good and interpolated pixels with 1.
01749 this maskCube is changed within the routine
01750 if a bad pixel was interpolated.
01751 n_neighbors: number of neighbors in one spectral direction
01752 with which the bad pixel will be interpolated (7)
01753 max_radius: maximal pixel radius within an image plane
01754 inside which valid pixels are searched to
01755 be used for interpolation. If there aren't
01756 found 9 good neighboring pixels within this
01757 radius the loop is left. (5)
01758 Out : resulting interpolated data cube.
01759 maskCube changed at the positions of the interpolated pixels
01760 Job : Bad pixel interpolation 3D like (saturated pixels exist):
01761 interpolates the bad pixels of the source cube by
01762 using the nearest neighbors.
01763 first it is checked if the bad pixel is interpolatable:
01764 it is only interpolatable if the number of good pixels
01765 in its spectrum of length 2*n_neighbors+1 exceeds 3 and
01766 if there is at least one good pixel on either side of the
01767 central pixel.
01768 Afterwards good neighboring pixels are searched within the
01769 image plane of the bad pixel by using an increasing pixel radius.
01770 Good pixels mean, the corresponding spectral pixels of the
01771 bad pixel and its spatial neighboring pixel must have
01772 at least 2 valid pixel pairs to be able to be used for
01773 the interpolation. The search is stopped if 9 valid neighboring
01774 pixels are found.
01775 Now normalize the found spectral values, collect the valid pixels
01776 (there must be at least 18) and take the median of the valid
01777 pixels with which the bad pixel is replaced.
01778 ---------------------------------------------------------------------------*/
01779
01780 OneCube * interpolCube ( OneCube * sourceCube,
01781 OneCube * maskCube,
01782 int n_neighbors, /* 7 */
01783 int max_radius ) /* 5 */
01784 {
01785 OneCube * returnCube ;
01786 float spec[100][2*n_neighbors+1] ;
01787 float spec1[300] ;
01788 int n_im, n_bad, n_bad1, n_bad2 ;
01789 int n_planes, specn, nspec1 ;
01790 int i, m, n, z, ni, kk, p ;
01791 int dis, dismin, dismax ;
01792 int agreed ;
01793 int xcordi, ycordi, xcordm, ycordm ;
01794
01795 if ( NullCube == sourceCube )
01796 {
01797 cpl_msg_error("interpolCube:"," no source cube given!\n") ;
01798 return NullCube ;
01799 }
01800
01801 if ( NullCube == maskCube )
01802 {
01803 cpl_msg_error("interpolCube:","no bad pixel mask cube given!\n") ;
01804 return NullCube ;
01805 }
01806
01807 if ( n_neighbors <= 0 )
01808 {
01809 cpl_msg_error("interpolCube:","wrong number of neighbors in the spectral direction given!\n") ;
01810 return NullCube ;
01811 }
01812
01813 if ( max_radius <= 0 )
01814 {
01815 cpl_msg_error("interpolCube:","wrong maximal radius for interpolation inside an image plane given!\n") ;
01816 return NullCube ;
01817 }
01818
01819 returnCube = copy_cube(sourceCube) ;
01820
01821 n_im = sourceCube->lx * sourceCube->ly ;
01822 n_planes = sourceCube->np ;
01823
01824 /* loop over the image planes and look for bad pixels and correct them */
01825 for ( z = 0 ; z < n_planes ; z++ ) /* go through image planes */
01826 {
01827
01828 /*-------------------------------------------------------------------
01829 * determine n, the length of one wing in one spectrum with which the
01830 * bad pixel will be interpolated. The length of one wing is n_neighbors
01831 * but less at the edges of the cube.
01832 */
01833 if ( z < n_neighbors )
01834 {
01835 n = z ;
01836 }
01837 else if ( n_planes - z <= n_neighbors)
01838 {
01839 n = n_planes - z -1 ;
01840 }
01841 else
01842 {
01843 n = n_neighbors ;
01844 }
01845
01846 for ( i = 0 ; i < n_im ; i ++ ) /* go through one image */
01847 {
01848 /* continue if the pixel is a good one */
01849 if ( maskCube -> plane[z] -> data[i] != 0. )
01850 {
01851 continue ;
01852 }
01853
01854 /*-----------------------------------------------------------------
01855 * exclude pixels with too many bad neighbors in the spectrum.
01856 * exit if: too few good pixels in the neighboring spectrum or
01857 * good pixels are only on one side of the spectrum.
01858 */
01859 n_bad = 0 ;
01860 n_bad1 = 0 ;
01861 n_bad2 = 0 ;
01862 for ( ni = z-n ; ni <= z+n ; ni++ ) /* go through the neighbor spectral pixels */
01863 {
01864 if ( maskCube -> plane[ni] -> data[i] == 0. )
01865 {
01866 n_bad++ ;
01867 /* count bad pixels on either spectral side of the bad pixel to be interpolated */
01868 if ( ni < z )
01869 {
01870 n_bad1++ ;
01871 }
01872 if ( ni > z )
01873 {
01874 n_bad2++ ;
01875 }
01876 }
01877 }
01878
01879 /*----------------------------------------------------------------------
01880 * now the criteria are checked which the neighborhood in the spectral
01881 * dimension has to match if the pixel is interpolatable.
01882 * The total number of the good pixel in the spectrum must be more than
01883 * 3 and there must be at least one good pixel on either side of the
01884 * central pixel.
01885 */
01886 if ( (2*n+1 - n_bad) < 3 || (n - n_bad1) < 1 || (n - n_bad2) < 1 )
01887 {
01888 continue ;
01889 }
01890
01891 /* read the master spectrum into the first row of the array spec */
01892 kk = 0 ;
01893 for ( ni = z-n ; ni <= z+n ; ni++ )
01894 {
01895 spec[1][kk] = maskCube -> plane[ni] -> data[i] != 0. ?
01896 sourceCube -> plane[ni] -> data[i] : ZERO ;
01897 kk++ ; /* length of spectrum */
01898 }
01899
01900 /* look for appropriate neighbors in the x-y neighborhood */
01901 agreed = 1 ; /* loop guard */
01902 specn = 2 ; /* number of spectra in spec. First is master spectrum */
01903 dismin = 0 ; /* x+y minimal distance to bad pixel */
01904 dismax = 1 ; /* x+y maximal distance to bad pixel */
01905 do
01906 {
01907 for ( m = 0 ; m < n_im ; m++ )
01908 {
01909 if ( maskCube -> plane[z] -> data[m] == 0. )
01910 {
01911 continue ;
01912 }
01913
01914 /* --------------------------------------------------------
01915 * determine the x and y coordinates of the bad pixel (i)
01916 * and the pixels used to interpolate (m)
01917 */
01918 xcordi = i % sourceCube->lx ;
01919 xcordm = m % sourceCube->lx ;
01920 ycordi = i / sourceCube->lx ;
01921 ycordm = m / sourceCube->lx ;
01922 /*-----------------------------------------------------
01923 * check the distance: take only close pixels
01924 * extension 'i' is coordinate of the bad pixel to be interpolated
01925 */
01926 dis = abs(xcordi-xcordm) + abs(ycordi-ycordm) ;
01927 if ( dis <= dismin || dis > dismax )
01928 {
01929 continue ;
01930 }
01931 /*------------------------------------------------------------
01932 * check on number of bad pixels in the spectrum of a neighbor
01933 * pixel; reject it if it contains less than 2 usable pixel pairs.
01934 * a bit more explanation:
01935 * let this be a 15 pixel spectrum with the pixel to be interpolated
01936 * denoted by '0' and other bad pixels marked with 'b'.
01937 * Good pixels are marked with '1'. Below a neighbor spectrum
01938 * is drawn containing bad pixels as well. The third line
01939 * shows the position of the usable pixel pairs, spectral
01940 * positions, where both spectra have valid pixels.
01941 *
01942 * 1 1 1 b b 1 1 0 b 1 b b 1 b b
01943 * b 1 1 1 b b 1 1 1 1 1 1 b b 1
01944 * ^ ^ ^ ^ 4 good pixel pairs
01945 */
01946
01947 n_bad = 0 ;
01948 for ( ni = z-n ; ni <= z+n ; ni++ )
01949 {
01950 if ( maskCube->plane[ni]->data[i] == 0. ||
01951 maskCube->plane[ni]->data[m] == 0. )
01952 {
01953 n_bad++ ;
01954 }
01955 }
01956 if ( n_bad > 2*n-1 ) /* we need at least 2 usable pixel pairs */
01957 {
01958 continue ;
01959 }
01960
01961 /* transfer the spectrum to the next position of array spec */
01962 kk = 0 ;
01963 for ( ni = z-n ; ni <= z+n ; ni++ )
01964 {
01965 spec[specn][kk] = maskCube->plane[ni]->data[m] != 0. ?
01966 sourceCube->plane[ni]->data[m] : ZERO ;
01967 kk++ ;
01968 }
01969 specn++ ;
01970 if ( specn > 10 ) /* if we have 9 neighbors then break */
01971 {
01972 agreed = 0 ;
01973 break ;
01974 }
01975 }
01976 /* if no break, increase search radius and continue */
01977 dismin++ ;
01978 dismax++ ;
01979 /* if search radius is too big, exit with fewer good neighbors */
01980 if ( dismax > max_radius )
01981 {
01982 agreed = 0 ;
01983 }
01984 } while(agreed) ;
01985
01986 specn-- ;
01987 dismax -= 2 ;
01988
01989 /* ---------------------------------------------------------------------------
01990 * Take the master spectrum with the bad pixel in the middle and divide it by
01991 * each of the neighbor spectra and normalize the division with the value in
01992 * the center position.
01993 */
01994 for ( kk = 0 ; kk < 2*n+1 ; kk++ )
01995 {
01996 if ( kk == n ) /* do not divide the master bad pixel */
01997 {
01998 continue ;
01999 }
02000
02001 /* do not divide bad pixels in the master spectrum */
02002 if ( qfits_isnan(spec[1][kk]) )
02003 {
02004 for ( p = 2 ; p <= specn ; p++ )
02005 {
02006 spec[p][kk] = ZERO ;
02007 }
02008 }
02009 else /* all is well, now divide */
02010 {
02011 for ( p = 2 ; p <= specn ; p++ )
02012 {
02013 if ( !qfits_isnan(spec[p][kk]) && spec[p][kk] != 0. &&
02014 !qfits_isnan(spec[p][n]) )
02015 {
02016 spec[p][kk] = spec[1][kk] / spec[p][kk] * spec[p][n] ;
02017 }
02018 else
02019 {
02020 spec[p][kk] = ZERO ;
02021 }
02022 }
02023 }
02024 }
02025
02026 /*----------------------------------------------------------------------------
02027 * determine the median of all values. With 9 good neighbors and at least
02028 * 2 good values per neighbor we have between 18 and 9*14 values for the
02029 * statistics. If there are not enough good neighbors available, only
02030 * continue if we have collected at least 18 values.
02031 */
02032 nspec1 = 0 ;
02033 /* collect the good values in the array spec1 */
02034 for ( p = 2 ; p <= specn ; p++ )
02035 {
02036 for ( kk = 0 ; kk < 2*n+1 ; kk++ )
02037 {
02038 if ( !qfits_isnan(spec[p][kk]) && kk != n )
02039 {
02040 spec1[nspec1] = spec[p][kk] ;
02041 nspec1++ ;
02042 }
02043 }
02044 }
02045
02046 /* now test if we have at least 18 values */
02047 if ( nspec1 < 18 )
02048 {
02049 continue ;
02050 }
02051
02052 /* interpolate the bad pixel by the median of spec1 */
02053 returnCube->plane[z]->data[i] = median(spec1, nspec1) ;
02054 maskCube ->plane[z]->data[i] = 1 ;
02055 }
02056 }
02057
02058 return returnCube ;
02059 }
02060
02061
02062 /*----------------------------------------------------------------------------
02063 Function : fineTuneCube()
02064 In : cube: cube, output of makeCube
02065 correct_diff_dist: differences of the slitlets from
02066 distance 32 given in the correct
02067 Spiffi row sequence. The first slitlet
02068 is the reference
02069 Output of makeCube*!
02070 Out : resulting data cube having the exact row positions
02071 Job : fine tunes each row in the right position according
02072 to the distances of the slitlets to each other
02073 (output of the north-south test).
02074 This means that the rows must be realigned by a
02075 fraction of a pixel to accomodate non-integer slit
02076 length. The fractional realignment is done by using
02077 tanh interpolation.
02078 Each row is rescaled so that the total flux is
02079 conserved.
02080 ---------------------------------------------------------------------------*/
02081
02082 OneCube * fineTuneCube( OneCube * cube,
02083 float * correct_diff_dist,
02084 int n_order )
02085 {
02086 OneCube * returnCube ;
02087 float row_data[cube->lx] ;
02088 float corrected_row_data[cube->lx] ;
02089 float xnum[n_order+1] ;
02090 float sum, new_sum ;
02091 float eval/*, dy*/ ;
02092 float * imageptr ;
02093 int row, col ;
02094 int i, z ;
02095 int imsize, n_points ;
02096 int firstpos ;
02097 int flag;
02098
02099 if ( NullCube == cube )
02100 {
02101 cpl_msg_error("fineTuneCube:","no input cube given!\n") ;
02102 return NullCube ;
02103 }
02104
02105 if ( NULL == correct_diff_dist )
02106 {
02107 cpl_msg_error("fineTuneCube:","no distances array from ns_test given!n") ;
02108 return NullCube ;
02109 }
02110
02111 if ( n_order <= 0 )
02112 {
02113 cpl_msg_error("fineTuneCube:","wrong order of interpolation polynom given!") ;
02114 returnCube = copy_cube(cube);
02115 return returnCube ;
02116 }
02117
02118 returnCube = copy_cube(cube);
02119
02120 imsize = cube -> ly ;
02121 if ( imsize != N_SLITLETS )
02122 {
02123 cpl_msg_error ("makeCube:","wrong image size\n" ) ;
02124 return NullCube ;
02125 }
02126
02127 n_points = n_order + 1 ;
02128 if ( n_points % 2 == 0 )
02129 {
02130 firstpos = (int)(n_points/2) - 1 ;
02131 }
02132 else
02133 {
02134 firstpos = (int)(n_points/2) ;
02135 }
02136
02137 for ( i = 0 ; i < n_points ; i++ )
02138 {
02139 xnum[i] = i ;
02140 }
02141
02142 for ( z = 0 ; z < cube->np ; z++ )
02143 {
02144 for ( row = 0 ; row < imsize ; row++ )
02145 {
02146 for ( col = 0 ; col < cube->lx ; col++ )
02147 {
02148 corrected_row_data[col] = 0. ;
02149 }
02150 sum = 0. ;
02151 for ( col = 0 ; col < cube->lx ; col++ )
02152 {
02153 row_data[col] = cube -> plane[z] -> data[col+row*cube->lx] ;
02154 if ( qfits_isnan(row_data[col]) )
02155 {
02156 row_data[col] = 0. ;
02157 for ( i = col - firstpos ; i < col -firstpos+n_points ; i++ )
02158 {
02159 if ( i < 0 ) continue ;
02160 if ( i >= cube->lx) continue ;
02161 corrected_row_data[i] = ZERO ;
02162 }
02163 }
02164 if ( col != 0 && col != cube->lx - 1 )
02165 {
02166 sum += row_data[col] ;
02167 }
02168 }
02169
02170
02171 new_sum = 0. ;
02172 for ( col = 0 ; col < cube->lx ; col++ )
02173 {
02174
02175 if ( qfits_isnan(corrected_row_data[col]) )
02176 {
02177 continue ;
02178 }
02179 if ( col - firstpos < 0 )
02180 {
02181 imageptr = &row_data[0] ;
02182 eval = correct_diff_dist[row] + col ;
02183 }
02184 else if ( col - firstpos + n_points >= cube->lx )
02185 {
02186 imageptr = &row_data[cube->lx - n_points] ;
02187 eval = correct_diff_dist[row] + col + n_points - cube->lx ;
02188 }
02189 else
02190 {
02191 imageptr = &row_data[col-firstpos] ;
02192 eval = correct_diff_dist[row] + firstpos ;
02193 }
02194
02195
02196 flag = 0;
02197 corrected_row_data[col]=nev_ille(xnum, imageptr, n_order, eval, &flag);
02198
02199
02200 if ( col != 0 && col != cube->lx - 1 )
02201 {
02202 new_sum += corrected_row_data[col] ;
02203 }
02204 }
02205 for ( col = 0 ; col < cube->lx ; col++ )
02206 {
02207
02208 if ( col == 0 )
02209 {
02210 returnCube->plane[z]->data[col+row*cube->lx] = ZERO ;
02211 }
02212 else if ( col == cube->lx - 1 )
02213 {
02214 returnCube->plane[z]->data[col+row*cube->lx] = ZERO ;
02215 }
02216 else
02217 {
02218 if ( qfits_isnan(corrected_row_data[col]) )
02219 {
02220 returnCube->plane[z]->data[col+row*cube->lx] = ZERO ;
02221 }
02222 else
02223 {
02224 if ( new_sum == 0. ) new_sum = 1. ;
02225
02226 returnCube->plane[z]->data[col+row*cube->lx] = corrected_row_data[col] ;
02227 }
02228 }
02229 }
02230 }
02231 }
02232 return returnCube ;
02233 }
02234
02235
02236
02237 /*----------------------------------------------------------------------------
02238 Function : fineTuneCubeByFFT()
02239 In : cube: cube, output of makeCube
02240 correct_diff_dist: differences of the slitlets from
02241 distance 32 given in the correct
02242 Spiffi row sequence. The first slitlet
02243 is the reference, therefore element
02244 23 is set 0.
02245 Output of makeCube!
02246 Out : resulting data cube having the exact row positions
02247 Job : fine tunes each row in the right position according
02248 to the distances of the slitlets to each other
02249 (output of the north-south test).
02250 This means that the rows must be realigned by a
02251 fraction of a pixel to accomodate non-integer slit
02252 length. The fractional realignment is done by using
02253 the FFT algorithm four1() of N.R.
02254 ---------------------------------------------------------------------------*/
02255
02256 OneCube * fineTuneCubeByFFT( OneCube * cube,
02257 float * correct_diff_dist )
02258 {
02259 OneCube * returnCube ;
02260 float row_data[cube->lx] ;
02261 dcomplex data[cube->lx] ;
02262 dcomplex corrected_data[cube->lx] ;
02263 unsigned nn[1];
02264 /*float corrected_row_data[cube->lx] ;*/
02265 float phi, pphi ;
02266 float coph, siph ;
02267 int row, col ;
02268 int i, z ;
02269 int imsize ;
02270 int blank_indicator ;
02271 nn[1] = cube->lx ;
02272
02273 if ( NullCube == cube )
02274 {
02275 cpl_msg_error("fineTuneCubeByFFT:"," no input cube given!\n") ;
02276 return NullCube ;
02277 }
02278
02279 if ( NULL == correct_diff_dist )
02280 {
02281 cpl_msg_error("fineTuneCubeByFFT:","no distances array from ns_test given!/n") ;
02282 return NullCube ;
02283 }
02284
02285 returnCube = copy_cube( cube ) ;
02286
02287 imsize = cube -> ly ;
02288 if ( imsize != N_SLITLETS )
02289 {
02290 cpl_msg_error ("makeCube:","wrong image size\n" ) ;
02291 return NullCube ;
02292 }
02293
02294 /* loop over the image planes */
02295 for ( z = 0 ; z < cube->np ; z++ )
02296 {
02297 /* consider one row at a time */
02298 for ( row = 0 ; row < imsize ; row++ )
02299 {
02300 blank_indicator = 1 ;
02301 for ( col = 0 ; col < cube->lx ; col++ )
02302 {
02303 /* transfer the row data to a double sized array */
02304 row_data[col] = cube -> plane[z] -> data[col+row*cube->lx] ;
02305 data[col].x = row_data[col] ;
02306 data[col].y = 0. ;
02307 /* if row contains a blank pixel proceed */
02308 if ( qfits_isnan(row_data[col]) )
02309 {
02310 blank_indicator = 0 ;
02311 }
02312 }
02313
02314 /* if row contains a blank don't apply FFT but proceed */
02315 if ( blank_indicator == 0 )
02316 {
02317 for ( col = 0 ; col < cube->lx ; col++ )
02318 {
02319 returnCube->plane[z]->data[col+row*cube->lx] = ZERO ;
02320 }
02321 continue ;
02322 }
02323
02324 /* FFT algorithm of eclipse */
02325 fftn( data, nn, 1, 1 ) ;
02326
02327 /* calculate the corrected phase shift for each frequency */
02328 phi = 2*PI_NUMB/(float)cube->lx * correct_diff_dist[row] ;
02329 for ( i = 0 ; i < cube->lx ; i++ )
02330 {
02331 /* positive frequencies */
02332 if ( i <= cube->lx/2 )
02333 {
02334 /* phase shift */
02335 pphi = phi * (float)(i) ;
02336 /* Euler factor */
02337 coph = cos ( pphi ) ;
02338 siph = sin ( pphi ) ;
02339 }
02340 else /* negative frequencies */
02341 {
02342 /* phase shift */
02343 pphi = phi * (float)(i - cube->lx/2) ;
02344 /* Euler factor */
02345 coph = cos ( pphi ) ;
02346 siph = sin ( pphi ) ;
02347 }
02348
02349 /* ----------------------------------------------------------------
02350 * now calculate the shift in the pixel space by multiplying
02351 * the fourier transform by the Euler factor of the phase shift
02352 * and inverse fourier transforming.
02353 * used Fourier pair: h(x-x0) <==> H(k)*exp(2*pi*i*k*x0)
02354 */
02355 /* calculate real part */
02356 corrected_data[i].x = data[i].x * coph - data[i].y * siph ;
02357 /* calculate imaginary part */
02358 corrected_data[i].y = data[i].x * siph + data[i].y * coph ;
02359 }
02360
02361 /* transform back: inverse FFT */
02362 fftn( corrected_data, nn, 1, -1 ) ;
02363
02364 /* normalize */
02365 for ( i = 0 ; i < cube->lx ; i++ )
02366 {
02367 corrected_data[i].x /= cube->lx ;
02368 corrected_data[i].y /= cube->lx ;
02369 }
02370
02371 /* now transfer row to output, leave the left-most and right-most pixel column */
02372 for ( col = 0 ; col < cube->lx ; col++ )
02373 {
02374 if ( col == 0 )
02375 {
02376 returnCube->plane[z]->data[col+row*cube->lx] = ZERO ;
02377 }
02378 else if ( col == cube->lx - 1 )
02379 {
02380 returnCube->plane[z]->data[col+row*cube->lx] = ZERO ;
02381 }
02382 else
02383 {
02384 returnCube->plane[z]->data[col+row*cube->lx] = corrected_data[col].x ;
02385 }
02386 }
02387 }
02388 }
02389 return returnCube ;
02390 }
02391
02392 /*----------------------------------------------------------------------------
02393 Function : fineTuneCubeBySpline()
02394 In : cube: cube, output of makeCube
02395 correct_diff_dist: differences of the slitlets from
02396 distance 32 given in the correct
02397 Spiffi row sequence. The first slitlet
02398 is the reference, therefore element
02399 23 is set 0.
02400 Output of makeCube!
02401 Out : resulting data cube having the exact row positions
02402 Job : fine tunes each row in the right position according
02403 to the distances of the slitlets to each other
02404 (output of the north-south test).
02405 This means that the rows must be realigned by a
02406 fraction of a pixel to accomodate non-integer slit
02407 length. The fractional realignment is done by using
02408 the spline interpolation algorithm splint in connection
02409 with the algorithm spline of N.R.
02410 This algorithms assume that each row is a tabulated
02411 function. The first derivatives of the interpolating
02412 function at the first and last point must be given.
02413 These are set higher than 1xe^30, so the routine
02414 sets the corresponding boundary condition for a natural
02415 spline, with zero second derivative on that boundary.
02416 Each row is rescaled so that the total flux is
02417 conserved.
02418 ---------------------------------------------------------------------------*/
02419
02420 OneCube * fineTuneCubeBySpline ( OneCube * cube,
02421 float * correct_diff_dist )
02422 {
02423 OneCube * returnCube ;
02424 float row_data[cube->lx] ;
02425 float corrected_row_data[cube->lx] ;
02426 float xnum[cube->lx] ;
02427 float sum, new_sum ;
02428 float eval[cube->lx] ;
02429 int row, col ;
02430 int i, z ;
02431 int imsize ;
02432
02433 if ( NullCube == cube )
02434 {
02435 cpl_msg_error("fineTuneCubeBySpline:","no input cube given!\n") ;
02436 return NullCube ;
02437 }
02438
02439 if ( NULL == correct_diff_dist )
02440 {
02441 cpl_msg_error("fineTuneCubeBySpline:","no distances array from ns_test given!/n") ;
02442 return NullCube ;
02443 }
02444
02445 imsize = cube -> ly ;
02446 if ( imsize != N_SLITLETS )
02447 {
02448 cpl_msg_error ("makeCube:","wrong image size\n" ) ;
02449 return NullCube ;
02450 }
02451
02452 returnCube = copy_cube( cube ) ;
02453
02454 /* fill the xa[] array for a polynomial interpolation */
02455 for ( i = 0 ; i < cube->lx ; i++ )
02456 {
02457 xnum[i] = i ;
02458 }
02459
02460 /* loop over the image planes */
02461 for ( z = 0 ; z < cube->np ; z++ )
02462 {
02463 /* consider 1 row at a time */
02464 for ( row = 0 ; row < imsize ; row++ )
02465 {
02466 for ( col = 0 ; col < cube->lx ; col++ )
02467 {
02468 corrected_row_data[col] = 0. ;
02469 }
02470 sum = 0. ; /* initialize flux for later rescaling */
02471 /* go through the columns and compute the flux for each row (leave the edge points) */
02472 for ( col = 0 ; col < cube->lx ; col++ )
02473 {
02474 eval[col] = correct_diff_dist[row] + (float)col ;
02475 row_data[col] = cube -> plane[z] -> data[col+row*cube->lx] ;
02476 if (col != 0 && col != cube->lx - 1 && !qfits_isnan(row_data[col]) )
02477 {
02478 sum += row_data[col] ;
02479 }
02480 if (qfits_isnan(row_data[col]) )
02481 {
02482 for ( i = col -1 ; i <= col+1 ; i++ )
02483 {
02484 if ( i < 0 ) continue ;
02485 if ( i >= cube->lx ) continue ;
02486 corrected_row_data[i] = ZERO ;
02487 }
02488 row_data[col] = 0. ;
02489 }
02490 }
02491
02492
02493 /* --------------------------------------------------------------------
02494 * now we do the cubic spline interpolation to achieve the fractional
02495 * (see eclipse).
02496 */
02497 if ( -1 == function1d_natural_spline( xnum, row_data, cube->lx, eval, corrected_row_data, cube->lx ) )
02498 {
02499 cpl_msg_error("fineTuneCubeBySpline:","error in spline interpolation\n") ;
02500 destroy_cube(returnCube) ;
02501 return NullCube ;
02502 }
02503
02504 new_sum = 0. ;
02505 for ( col = 0 ; col < cube->lx ; col++ )
02506 {
02507 if (qfits_isnan(corrected_row_data[col])) continue ;
02508 /* don't take the edge points to calculate the scaling factor */
02509 if ( col != 0 && col != cube->lx - 1 )
02510 {
02511 new_sum += corrected_row_data[col] ;
02512 }
02513 }
02514 for ( col = 0 ; col < cube->lx ; col++ )
02515 {
02516 /* --------------------------------------------------------------------
02517 * rescale the row data and fill the returned cube,
02518 * leave the left-most and right-most
02519 * pixel column
02520 */
02521 if ( col == 0 )
02522 {
02523 returnCube->plane[z]->data[col+row*cube->lx] = ZERO ;
02524 }
02525 else if ( col == cube->lx - 1 )
02526 {
02527 returnCube->plane[z]->data[col+row*cube->lx] = ZERO ;
02528 }
02529 else
02530 {
02531 if ( qfits_isnan(corrected_row_data[col]) )
02532 {
02533 returnCube->plane[z]->data[col+row*cube->lx] = ZERO ;
02534 }
02535 else
02536 {
02537 if (new_sum == 0.) new_sum = 1. ;
02538 /* rescaling is commented out because it delivers wrong results
02539 in case of appearance of blanks or bad pixels */
02540 /* corrected_row_data[col] *= sum / new_sum ; */
02541 returnCube->plane[z]->data[col+row*cube->lx] = corrected_row_data[col] ;
02542 }
02543 }
02544 }
02545 }
02546 }
02547 return returnCube ;
02548 }
02549
02550 /*----------------------------------------------------------------------------
02551 Function : calibrate_ns_test()
02552 In : ns_image: north-south image
02553 n_slitlets: number of slitlets
02554 halfWidth: half width of the box in which the lines
02555 are fit by a gaussian
02556 fwhm: first guess of the full width at half maximum
02557 minDiff: minimum amplitude below which the fit
02558 will not be carried through
02559 estimated_dist: estimated average distance of spectra
02560 devtol: maximal pixel deviation of the distances from
02561 slitlet center
02562 Out : array of the distances of the slitlets
02563 from each other
02564 Job : determines the distances of the spectra from the center
02565 of the slitlets to be able to find the correct angle of
02566 the north-south entrance slit.
02567 ---------------------------------------------------------------------------*/
02568
02569 float * calibrate_ns_test( OneImage * ns_image,
02570 int n_slitlets,
02571 int halfWidth,
02572 float fwhm,
02573 float minDiff,
02574 float estimated_dist,
02575 float devtol,
02576 int bottom,
02577 int top )
02578 {
02579 int i, j, k, m, row, col, n, ni, na ;
02580 int position, counter, iters ;
02581 int xdim, ndat, its, numpar ;
02582 pixelvalue row_buf[ns_image -> lx] ;
02583 float sum, mean, maxval ;
02584 float tol, lab ;
02585 float * distances ;
02586 float * ret_distances ;
02587 float distances_buf[ns_image -> ly][n_slitlets] ;
02588 float x_position[n_slitlets] ;
02589 float * xdat, * wdat ;
02590 int * mpar ;
02591 int found[3*n_slitlets], found_clean[3*n_slitlets] ;
02592 int found_cleanit[3*n_slitlets] ;
02593 Vector * line ;
02594 FitParams ** par ;
02595 int foundit, begin, end ;
02596 int zeroindicator ;
02597 int row_index ;
02598
02599 if ( ns_image == NullImage )
02600 {
02601 cpl_msg_error("calibrate_ns_test:","sorry, no image given\n") ;
02602 return NULL ;
02603 }
02604 if ( n_slitlets < 1 )
02605 {
02606 cpl_msg_error("calibrate_ns_test:","wrong number of slitlets given\n") ;
02607 return NULL ;
02608 }
02609 if ( halfWidth < 0 || halfWidth >= estimated_dist )
02610 {
02611 cpl_msg_error("calibrate_ns_test:","wrong half width given\n") ;
02612 return NULL ;
02613 }
02614 if ( fwhm <= 0. )
02615 {
02616 cpl_msg_error("calibrate_ns_test:","wrong fwhm given\n") ;
02617 return NULL ;
02618 }
02619 if ( minDiff < 1. )
02620 {
02621 cpl_msg_error("calibrate_ns_test:","wrong minDiff given\n") ;
02622 return NULL ;
02623 }
02624
02625 /* allocate memory for output array */
02626 if (NULL == (distances = (float *) cpl_calloc ( n_slitlets , sizeof (float) )))
02627 {
02628 cpl_msg_error("calibrate_ns_test:","could not allocate memory\n") ;
02629 return NULL ;
02630 }
02631 /* allocate memory for output array */
02632 if (NULL == (ret_distances = (float *) cpl_calloc ( n_slitlets , sizeof (float) )))
02633 {
02634 cpl_msg_error("calibrate_ns_test:","could not allocate memory\n") ;
02635 return NULL ;
02636 }
02637
02638 /* go through the image rows */
02639 for ( row = 0 ; row < ns_image -> ly ; row++ )
02640 {
02641 zeroindicator = 0 ;
02642
02643 /* initialize the distance buffer */
02644 for ( i = 0 ; i < n_slitlets ; i++ )
02645 {
02646 distances_buf[row][i] = ZERO ;
02647 }
02648
02649 /* fill the row buffer array with image data */
02650 for ( col = 0 ; col < ns_image -> lx ; col++ )
02651 {
02652 row_buf[col] = ns_image -> data[col + row*ns_image->lx] ;
02653 }
02654
02655 /* determine the mean of the row data */
02656 sum = 0. ;
02657 n = 0 ;
02658 for ( i = 0 ; i < ns_image -> lx ; i++ )
02659 {
02660 if ( qfits_isnan(row_buf[i]) )
02661 {
02662 continue ;
02663 }
02664 sum += row_buf[i] ;
02665 n++ ;
02666 }
02667 mean = sum / (float)n ;
02668
02669 /* store the positions of image values greater than the mean */
02670 n = 0 ;
02671 for ( i = 0 ; i < ns_image -> lx ; i++ )
02672 {
02673 if (qfits_isnan(row_buf[i]))
02674 {
02675 continue ;
02676 }
02677 if ( row_buf[i] > mean + ESTIMATE )
02678 {
02679 found[n] = i ;
02680 n++ ;
02681 }
02682 }
02683
02684 if ( n < n_slitlets )
02685 {
02686 cpl_msg_warning("calibrate_ns_test4:","t4 wrong number of intensity columns found in row: %d, found number: %d\n", row, n) ;
02687 continue ;
02688 }
02689 else
02690 {
02691 /* find the maximum value position around the found columns */
02692 na = 0 ;
02693 for ( i = 1 ; i < n ; i ++ )
02694 {
02695 if ( found[i] - found[i-1] < halfWidth )
02696 {
02697 begin = found[i] - halfWidth ;
02698 if ( begin < 0 )
02699 {
02700 begin = 0 ;
02701 }
02702 end = found[i] + halfWidth ;
02703 if ( end >= ns_image->lx )
02704 {
02705 end = ns_image->lx - 1 ;
02706 }
02707 /* find the maximum value inside the box around the found positions*/
02708 maxval = -FLT_MAX ;
02709 foundit = 0 ;
02710 for ( j = begin ; j <= end ; j++ )
02711 {
02712 /* do not consider boxes that contain bad pixels */
02713 if (qfits_isnan(row_buf[j]))
02714 {
02715 continue ;
02716 }
02717 if (row_buf[j] >= maxval )
02718 {
02719 maxval = row_buf[j] ;
02720 foundit = j ;
02721 }
02722 }
02723 if (maxval == -FLT_MAX)
02724 {
02725 continue ;
02726 }
02727 for ( k = 0 ; k < na ; k++ )
02728 {
02729 if ( found_cleanit[k] >= begin && found_cleanit[k] < foundit )
02730 {
02731 na-- ;
02732 }
02733 }
02734 for ( k = 0 ; k < n ; k++ )
02735 {
02736 if ( found[k] == foundit)
02737 {
02738 if ( found_cleanit[na-1] != found[k] )
02739 {
02740 found_cleanit[na] = found[k] ;
02741 na++ ;
02742 }
02743 }
02744 }
02745 }
02746 else
02747 {
02748 if ( i == 1 )
02749 {
02750 found_cleanit[na] = found[0] ;
02751 na++ ;
02752 found_cleanit[na] = found[1] ;
02753 na++ ;
02754 }
02755 else
02756 {
02757 if ( found_cleanit[na-1] != found[i-1])
02758 {
02759 found_cleanit[na] = found[i-1] ;
02760 na++ ;
02761 }
02762 if ( found_cleanit[na-1] != found[i])
02763 {
02764 found_cleanit[na] = found[i] ;
02765 na++ ;
02766 }
02767 }
02768 }
02769 }
02770
02771 /* determine only one pixel position for each slitlet intensity */
02772 j = 1 ;
02773 for ( i = 1 ; i < na ; i++ )
02774 {
02775 if ( (float)(found_cleanit[i] - found_cleanit[i-1]) < (estimated_dist - devtol) ||
02776 (float)(found_cleanit[i] - found_cleanit[i-1]) > (estimated_dist + devtol) )
02777 {
02778 continue ;
02779 }
02780 else
02781 {
02782 found_clean[j-1] = found_cleanit[i-1] ;
02783 found_clean[j] = found_cleanit[i] ;
02784 j++ ;
02785 }
02786 }
02787 }
02788 if ( j > n_slitlets )
02789 {
02790 /* check the distance again */
02791 ni = 1 ;
02792 for ( i = 1 ; i < j ; i++ )
02793 {
02794 if ( (float)(found_clean[i] - found_clean[i-1]) < (estimated_dist - devtol ) ||
02795 (float)(found_clean[i] - found_clean[i-1]) > (estimated_dist + devtol ) )
02796 {
02797 continue ;
02798 }
02799 else
02800 {
02801 found_clean[ni-1] = found_clean[i-1] ;
02802 found_clean[ni] = found_clean[i] ;
02803 ni++ ;
02804 }
02805 }
02806 if ( ni != n_slitlets )
02807 {
02808 cpl_msg_warning("calibrate_ns_test5:","t5 wrong number of intensity columns found in row: %d, found number: %d\n", row, ni) ;
02809 continue ;
02810 }
02811 else
02812 {
02813 j = ni ;
02814 }
02815 }
02816 else if ( j < n_slitlets )
02817 {
02818 cpl_msg_warning("calibrate_ns_test6:","t6 wrong number of intensity columns found in row: %d , found number: %d\n", row, j) ;
02819 continue ;
02820 }
02821 counter = 0 ;
02822 /* go through the found intensity pixels in one row */
02823 for ( i = 0 ; i < j ; i++ )
02824 {
02825 /* allocate memory for the array where the line is fitted in */
02826 if ( NULL == (line = newVector (2*halfWidth + 1)) )
02827 {
02828 cpl_msg_error ("calibrate_ns_test:","cannot allocate new Vector \n") ;
02829 cpl_free(distances) ;
02830 return NULL ;
02831 }
02832
02833 /* allocate memory */
02834 xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
02835 wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
02836 mpar = (int *) cpl_calloc( MAXPAR, sizeof (int) ) ;
02837 par = newFitParams(1) ;
02838
02839 m = 0 ;
02840 for ( k = found_clean[i]-halfWidth ; k <= found_clean[i]+halfWidth ; k++ )
02841 {
02842 if ( k < 0 )
02843 {
02844 k = 0. ;
02845 }
02846 else if ( k >= ns_image -> lx )
02847 {
02848 k = ns_image->lx - 1 ;
02849 }
02850 else if ( qfits_isnan(row_buf[k]) )
02851 {
02852 zeroindicator = 1 ;
02853 break ;
02854 }
02855 else
02856 {
02857 line -> data[m] = row_buf[k] ;
02858 m++ ;
02859 }
02860 }
02861 if ( zeroindicator == 1 )
02862 {
02863 destroyVector(line) ;
02864 cpl_free(xdat) ;
02865 cpl_free(wdat) ;
02866 cpl_free(mpar) ;
02867 destroyFitParams(par) ;
02868 break ;
02869 }
02870
02871 /*--------------------------------------------------------------------
02872 * go through the spectral vector
02873 * determine the maximum pixel value in the spectral vector
02874 */
02875 maxval = -FLT_MAX ;
02876 position = -INT32_MAX ;
02877 for ( k = 0 ; k < m ; k++ )
02878 {
02879 xdat[k] = k ;
02880 wdat[k] = 1.0 ;
02881 if ( line -> data[k] >= maxval )
02882 {
02883 maxval = line -> data[k] ;
02884 position = k ;
02885 }
02886 }
02887
02888 /* set initial values for the fitting routine */
02889 xdim = XDIM ;
02890 ndat = line -> n_elements ;
02891 numpar = MAXPAR ;
02892 tol = TOL ;
02893 lab = LAB ;
02894 its = ITS ;
02895 (*par) -> fit_par[1] = fwhm ;
02896 (*par) -> fit_par[2] = (float) position ;
02897 (*par) -> fit_par[3] = (float) (line -> data[0] + line -> data[line->n_elements - 1]) / 2.0 ;
02898 (*par) -> fit_par[0] = maxval - ((*par) -> fit_par[3]) ;
02899
02900
02901 /* exclude negative peaks and low signal cases */
02902 if ( (*par) -> fit_par[0] < minDiff )
02903 {
02904 cpl_msg_warning ("calibrate_ns_test:","sorry, signal of line too low to fit in row: %d in slitlet %d\n", row, i) ;
02905 destroyVector(line) ;
02906 cpl_free(xdat) ;
02907 cpl_free(wdat) ;
02908 cpl_free(mpar) ;
02909 destroyFitParams(par) ;
02910 continue ;
02911 }
02912
02913 for ( k = 0 ; k < MAXPAR ; k++ )
02914 {
02915 (*par) -> derv_par[k] = 0.0 ;
02916 mpar[k] = 1 ;
02917 }
02918 /* finally, do the least square fit using a gaussian */
02919 if ( 0 > ( iters = lsqfit_c( xdat, &xdim, line -> data, wdat, &ndat, (*par) -> fit_par,
02920 (*par) -> derv_par, mpar, &numpar, &tol, &its, &lab )) )
02921 {
02922 /*
02923 cpl_msg_debug ("calibrate_ns_test:","lsqfit_c: least squares fit failed, error no.: %d in row: %d in slitlet %d\n", iters, row, i) ;
02924 */
02925 destroyVector(line) ;
02926 cpl_free(xdat) ;
02927 cpl_free(wdat) ;
02928 cpl_free(mpar) ;
02929 destroyFitParams(par) ;
02930 continue ;
02931 }
02932
02933 /* check for negative fit results */
02934 if ( (*par) -> fit_par[0] <= 0. || (*par) -> fit_par[1] <= 0. ||
02935 (*par) -> fit_par[2] < 0. )
02936 {
02937 cpl_msg_warning ("calibrate_ns_test:","negative parameters as fit result, not used! in row %d in slitlet %d\n", row, i) ;
02938 destroyVector(line) ;
02939 cpl_free(xdat) ;
02940 cpl_free(wdat) ;
02941 cpl_free(mpar) ;
02942 destroyFitParams(par) ;
02943 continue ;
02944 }
02945
02946 /* correct the fitted position for the given row of the line in image coordinates */
02947 (*par) -> fit_par[2] = (float) (found_clean[i] - halfWidth) + (*par) -> fit_par[2] ;
02948 x_position[counter] = (*par) -> fit_par[2] ;
02949 counter ++ ;
02950
02951 /* free memory */
02952 destroyFitParams(par) ;
02953 destroyVector ( line ) ;
02954 cpl_free ( xdat ) ;
02955 cpl_free ( wdat ) ;
02956 cpl_free ( mpar ) ;
02957 }
02958 if (zeroindicator == 1)
02959 {
02960 cpl_msg_debug ("calibrate_ns_test:","bad pixel in fitting box in row: %d\n", row) ;
02961 continue ;
02962 }
02963
02964 if ( counter != n_slitlets )
02965 {
02966 cpl_msg_warning("calibrate_ns_test:","wrong number of slitlets found in row: %d\n", row) ;
02967 continue ;
02968 }
02969 /* store the distances between the sources and the slitlet centers */
02970 for ( i = 0 ; i < n_slitlets ; i++ )
02971 {
02972 distances_buf[row][i] = x_position[i] - (15.5 + 32.*(float)i) ;
02973 }
02974 }
02975
02976 /* ----------------------------------------------------------------
02977 * go through the rows again and take the mean of the distances,
02978 * throw away the runaways
02979 */
02980 for ( i = 0 ; i < n_slitlets ; i++ )
02981 {
02982 n = 0 ;
02983 sum = 0. ;
02984 for ( row = bottom ; row < top ; row++ )
02985 {
02986 if ( fabs( distances_buf[row][i] ) > devtol ||
02987 qfits_isnan(distances_buf[row][i]) )
02988 {
02989 /*
02990 printf("A check1=%g ref1=%g check3=%d\n",
02991 fabs( distances_buf[row][i] ),
02992 devtol,
02993 qfits_isnan(distances_buf[row][i]));
02994 */
02995 continue ;
02996 }
02997 sum += distances_buf[row][i] ;
02998 n++ ;
02999 }
03000 if ( n < 2 )
03001 {
03002 cpl_msg_error( "calibrate_ns_test:","distances array could not be determined completely!, deviations of distances from devtol too big\n" ) ;
03003 cpl_free(distances) ;
03004 return NULL ;
03005 }
03006 else
03007 {
03008 distances[i] = sum / (float)n ;
03009 }
03010 }
03011
03012 /* now sort the result according to the row sequence in the reconstructed image*/
03013 for ( i = 0 ; i < n_slitlets ; i++ )
03014 {
03015 switch (i)
03016 {
03017 case 0:
03018 row_index = 8 ;
03019 break ;
03020 case 1:
03021 row_index = 7 ;
03022 break ;
03023 case 2:
03024 row_index = 9 ;
03025 break ;
03026 case 3:
03027 row_index = 6 ;
03028 break ;
03029 case 4:
03030 row_index = 10 ;
03031 break ;
03032 case 5:
03033 row_index = 5 ;
03034 break ;
03035 case 6:
03036 row_index = 11 ;
03037 break ;
03038 case 7:
03039 row_index = 4 ;
03040 break ;
03041 case 8:
03042 row_index = 12 ;
03043 break ;
03044 case 9:
03045 row_index = 3 ;
03046 break ;
03047 case 10:
03048 row_index = 13 ;
03049 break ;
03050 case 11:
03051 row_index = 2 ;
03052 break ;
03053 case 12:
03054 row_index = 14 ;
03055 break ;
03056 case 13:
03057 row_index = 1 ;
03058 break ;
03059 case 14:
03060 row_index = 15 ;
03061 break ;
03062 case 15:
03063 row_index = 0 ;
03064 break ;
03065 case 16:
03066 row_index = 31 ;
03067 break ;
03068 case 17:
03069 row_index = 16 ;
03070 break ;
03071 case 18:
03072 row_index = 30 ;
03073 break ;
03074 case 19:
03075 row_index = 17 ;
03076 break ;
03077 case 20:
03078 row_index = 29 ;
03079 break ;
03080 case 21:
03081 row_index = 18 ;
03082 break ;
03083 case 22:
03084 row_index = 28 ;
03085 break ;
03086 case 23:
03087 row_index = 19 ;
03088 break ;
03089 case 24:
03090 row_index = 27 ;
03091 break ;
03092 case 25:
03093 row_index = 20 ;
03094 break ;
03095 case 26:
03096 row_index = 26 ;
03097 break ;
03098 case 27:
03099 row_index = 21 ;
03100 break ;
03101 case 28:
03102 row_index = 25 ;
03103 break ;
03104 case 29:
03105 row_index = 22 ;
03106 break ;
03107 case 30:
03108 row_index = 24 ;
03109 break ;
03110 case 31:
03111 row_index = 23 ;
03112 break ;
03113 default:
03114 cpl_msg_error("calibrate_ns_test:","wrong number of a slitlet\n") ;
03115 cpl_free (distances) ;
03116 return NULL ;
03117 }
03118 ret_distances[row_index] = distances[i] ;
03119 }
03120 cpl_free(distances) ;
03121 return ret_distances ;
03122 }
03123 /*----------------------------------------------------------------------------
03124 Function : makeCube()
03125 In : calibImage: resampled source image
03126 distances: distances of the slitlets from each other
03127 output of function ns_test
03128 correct_diff_dist: dummy array with 32 elements
03129 Out : resulting source data cube
03130 correct_diff_dist: differences of the slitlets from
03131 distance 32 given in the correct
03132 Spiffi row sequence. The first slitlet
03133 is the reference, therefore element
03134 23 is set 0.
03135 Job : makes a data cube out of a resampled source image
03136 this SPIFFI specific routine takes into account the
03137 Spiffi slitlet order on the detector.
03138 Also shifts the resulting image rows by one pixel if
03139 necessary according to the distances array gained from
03140 the north-south test routine.
03141 Can do the same with the bad pixel map image to generate a
03142 bad pixel mask cube.
03143 ---------------------------------------------------------------------------*/
03144
03145 OneImage * maketrueResamp ( OneImage * calibImage, float * distances , OneImage * wavemap)
03146 {
03147 OneImage * returnImage ;
03148 float edges[33] ;
03149 int imsize, kslit,i,j ;
03150 int slit_index ;
03151 int z, col, recol ;
03152
03153 /*edges[0]=0;
03154 for(i=0;i<31;i++)
03155 edges[i+1]=distances[i]+edges[i];
03156 edges[32]=2048;
03157
03158 for(i=0;i<33;i++)
03159 cpl_msg_warning("%g",edges[i]);*/
03160
03161 edges[0]=0;
03162 j=1;
03163 for(i=0;i<wavemap->lx-1;i++)
03164 {
03165 if((wavemap->data[i]-wavemap->data[i+1])>0.0025 || (wavemap->data[i]-wavemap->data[i+1])<-0.0025)
03166 {
03167 cpl_msg_error("maketrueResamp","wavemap edge %d", i+1);
03168 edges[j]=i+1;
03169 j++;
03170 }
03171 }
03172 edges[32]=2048;
03173
03174 imsize = calibImage->lx / N_SLITLETS ;
03175
03176 /* allocate memory */
03177 returnImage = new_image(calibImage->lx,calibImage->ly);
03178 for ( z = 0 ; z < calibImage -> ly ; z++ ) /* go through the z-axis */
03179 {
03180 for ( col = 0 ; col < calibImage->lx ; col++ ) /* go through the image columns */
03181 returnImage->data[col+z*calibImage->lx]=ZERO;
03182 }
03183
03184
03185 /* now build the data cube out of the resampled image */
03186 for ( z = 0 ; z < calibImage -> ly ; z++ ) /* go through the z-axis */
03187 {
03188 kslit = 0 ;
03189 slit_index = -1 ;
03190 recol = -1 ;
03191 for ( col = 0 ; col < calibImage->lx ; col++ ) /* go through the image columns */
03192 {
03193 /*if ( col % imsize == 0 )
03194 {*/
03195 recol = 0 ;
03196 /*kslit = col/imsize ;*/
03197 for(i=0;i<32;i++)
03198 {
03199 if(col>=nint(edges[i]) && col<nint(edges[i+1]))
03200 kslit=i;
03201 }
03202 /* sort the slitlets in the right spiffi specific way */
03203 switch (kslit)
03204 {
03205 case 0:
03206 slit_index = 8 ;
03207 break ;
03208 case 1:
03209 slit_index = 7 ;
03210 break ;
03211 case 2:
03212 slit_index = 9 ;
03213 break ;
03214 case 3:
03215 slit_index = 6 ;
03216 break ;
03217 case 4:
03218 slit_index = 10 ;
03219 break ;
03220 case 5:
03221 slit_index = 5 ;
03222 break ;
03223 case 6:
03224 slit_index = 11 ;
03225 break ;
03226 case 7:
03227 slit_index = 4 ;
03228 break ;
03229 case 8:
03230 slit_index = 12 ;
03231 break ;
03232 case 9:
03233 slit_index = 3 ;
03234 break ;
03235 case 10:
03236 slit_index = 13 ;
03237 break ;
03238 case 11:
03239 slit_index = 2 ;
03240 break ;
03241 case 12:
03242 slit_index = 14 ;
03243 break ;
03244 case 13:
03245 slit_index = 1 ;
03246 break ;
03247 case 14:
03248 slit_index = 15 ;
03249 break ;
03250 case 15:
03251 slit_index = 0 ;
03252 break ;
03253 case 16:
03254 slit_index = 31 ;
03255 break ;
03256 case 17:
03257 slit_index = 16 ;
03258 break ;
03259 case 18:
03260 slit_index = 30 ;
03261 break ;
03262 case 19:
03263 slit_index = 17 ;
03264 break ;
03265 case 20:
03266 slit_index = 29 ;
03267 break ;
03268 case 21:
03269 slit_index = 18 ;
03270 break ;
03271 case 22:
03272 slit_index = 28 ;
03273 break ;
03274 case 23:
03275 slit_index = 19 ;
03276 break ;
03277 case 24:
03278 slit_index = 27 ;
03279 break ;
03280 case 25:
03281 slit_index = 20 ;
03282 break ;
03283 case 26:
03284 slit_index = 26 ;
03285 break ;
03286 case 27:
03287 slit_index = 21 ;
03288 break ;
03289 case 28:
03290 slit_index = 25 ;
03291 break ;
03292 case 29:
03293 slit_index = 22 ;
03294 break ;
03295 case 30:
03296 slit_index = 24 ;
03297 break ;
03298 case 31:
03299 slit_index = 23 ;
03300 break ;
03301 default:
03302 cpl_msg_error("makeCube:","wrong slitlet index: couldn't be a spiffi image, \
03303 there must be 32 slitlets!\n") ;
03304 break ;
03305 }
03306
03307 /*}*/
03308
03309 /* fill each cube plane with one image row */
03310 if((col-nint(edges[kslit]))>0 && (col-nint(edges[kslit]))<imsize-1 )
03311 returnImage->data[(col-nint(edges[kslit]))+slit_index*imsize+z*calibImage->lx] =
03312 calibImage->data[col+z*calibImage->lx] ;
03313 else
03314 returnImage->data[(col-nint(edges[kslit]))+slit_index*imsize+z*calibImage->lx] = ZERO;
03315 /*recol++ ;*/
03316
03317 }
03318 }
03319 return returnImage ;
03320 }
03321
03322 /*The old slitlet order*/
03323 /*switch (kslit)
03324 {
03325 case 0:
03326 slit_index = 23 ;
03327 break ;
03328 case 1:
03329 slit_index = 24 ;
03330 break ;
03331 case 2:
03332 slit_index = 22 ;
03333 break ;
03334 case 3:
03335 slit_index = 25 ;
03336 break ;
03337 case 4:
03338 slit_index = 21 ;
03339 break ;
03340 case 5:
03341 slit_index = 26 ;
03342 break ;
03343 case 6:
03344 slit_index = 20 ;
03345 break ;
03346 case 7:
03347 slit_index = 27 ;
03348 break ;
03349 case 8:
03350 slit_index = 19 ;
03351 break ;
03352 case 9:
03353 slit_index = 28 ;
03354 break ;
03355 case 10:
03356 slit_index = 18 ;
03357 break ;
03358 case 11:
03359 slit_index = 29 ;
03360 break ;
03361 case 12:
03362 slit_index = 17 ;
03363 break ;
03364 case 13:
03365 slit_index = 30 ;
03366 break ;
03367 case 14:
03368 slit_index = 16 ;
03369 break ;
03370 case 15:
03371 slit_index = 31 ;
03372 break ;
03373 case 16:
03374 slit_index = 0 ;
03375 break ;
03376 case 17:
03377 slit_index = 15 ;
03378 break ;
03379 case 18:
03380 slit_index = 1 ;
03381 break ;
03382 case 19:
03383 slit_index = 14 ;
03384 break ;
03385 case 20:
03386 slit_index = 2 ;
03387 break ;
03388 case 21:
03389 slit_index = 13 ;
03390 break ;
03391 case 22:
03392 slit_index = 3 ;
03393 break ;
03394 case 23:
03395 slit_index = 12 ;
03396 break ;
03397 case 24:
03398 slit_index = 4 ;
03399 break ;
03400 case 25:
03401 slit_index = 11 ;
03402 break ;
03403 case 26:
03404 slit_index = 5 ;
03405 break ;
03406 case 27:
03407 slit_index = 10 ;
03408 break ;
03409 case 28:
03410 slit_index = 6 ;
03411 break ;
03412 case 29:
03413 slit_index = 9 ;
03414 break ;
03415 case 30:
03416 slit_index = 7 ;
03417 break ;
03418 case 31:
03419 slit_index = 8 ;
03420 break ;
03421 default:
03422 cpl_msg_error("makeCube:","wrong slitlet index: couldn't be a spiffi image, \
03423 there must be 32 slitlets!\n") ;
03424 destroy_cube(returnCube) ;
03425 return NullCube ;
03426 break ;
03427 }*/
03428
03429
03430 /*--------------------------------------------------------------------------*/
1.2.13.1 written by Dimitri van Heesch,
© 1997-2001