00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076 #define POSIX_SOURCE 1
00077 #include "vltPort.h"
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087 #include "coltilt.h"
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113 float slopeOfSpectrum( OneImage * ns_image,
00114 int box_length,
00115 float fwhm,
00116 float minDiff )
00117 {
00118 int i, k, row, col ;
00119 int counter, iters ;
00120 int xdim, ndat, its, numpar ;
00121 float maxval ;
00122 float tol, lab ;
00123 float col_value[ns_image->ly] ;
00124 float column_value[ns_image->ly] ;
00125 pixelvalue col_position[ns_image->ly] ;
00126 int column_position[ns_image->ly] ;
00127 int col_median ;
00128 float x_position[ns_image->ly] ;
00129 float * xdat, * wdat ;
00130 int * mpar ;
00131 Vector * line ;
00132 FitParams ** dec_par ;
00133 FitParams * par ;
00134 float x[ns_image->lx], y[ns_image->lx] ;
00135 float sig[ns_image->lx] ;
00136 int position ;
00137 int ndata, mwt ;
00138 float a, b, siga, sigb, chi2, q ;
00139 int bad_ind ;
00140
00141
00142 if ( ns_image == NullImage )
00143 {
00144 cpl_msg_error("slopeOfSpectrum:"," sorry, no image given\n") ;
00145 return FLAG ;
00146 }
00147 if ( box_length <= 1 || box_length >= (int) sqrt(ns_image -> lx) )
00148 {
00149 cpl_msg_error("slopeOfSpectrum:"," wrong box length given\n") ;
00150 return FLAG ;
00151 }
00152 if ( fwhm < 1. || fwhm > 10. )
00153 {
00154 cpl_msg_error("slopeOfSpectrum","wrong full width at half maximum given\n") ;
00155 return FLAG ;
00156 }
00157 if ( minDiff < 1. )
00158 {
00159 cpl_msg_error("slopeOfSpectrum","wrong amplitude threshold given\n") ;
00160 return FLAG ;
00161 }
00162
00163
00164 for ( row = 0 ; row < ns_image -> ly ; row++ )
00165 {
00166 col_value[row] = -FLT_MAX ;
00167 col_position[row] = -1. ;
00168
00169 for ( col = 0 ; col < ns_image -> lx ; col++ )
00170 {
00171 if ( ns_image -> data[col+row*ns_image->lx] > col_value[row] )
00172 {
00173 col_value[row] = ns_image -> data[col+row*ns_image->lx] ;
00174 col_position[row] = (pixelvalue)col ;
00175 }
00176 }
00177 }
00178
00179
00180
00181 col_median = (int)median(col_position, ns_image->ly) ;
00182 cpl_msg_info ("slopeOfSpectrum","median column position of brightest spectrum %d\n", col_median) ;
00183
00184
00185 for ( row = 0 ; row < ns_image -> ly ; row++ )
00186 {
00187 x_position[row] = 0. ;
00188 column_value[row] = -FLT_MAX ;
00189 column_position[row] = -1 ;
00190 for ( col = col_median - box_length ; col <= col_median + box_length ; col++ )
00191 {
00192 if ( ns_image -> data[col+row*ns_image->lx] > column_value[row] )
00193 {
00194 column_value[row] = ns_image -> data[col+row*ns_image->lx] ;
00195 column_position[row] = col ;
00196 }
00197 }
00198
00199
00200 if ( NULL == (line = newVector (2*box_length + 1)) )
00201 {
00202 cpl_msg_error ("slopeOfSpectrum","cannot allocate new Vector in row %d\n", row) ;
00203 return FLAG ;
00204 }
00205
00206
00207 xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00208 wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00209 mpar = (int *) cpl_calloc( MAXPAR, sizeof (int) ) ;
00210 dec_par = newFitParams(1) ;
00211 par = dec_par[0];
00212
00213 counter = 0 ;
00214 bad_ind = 0 ;
00215
00216 for ( col = column_position[row] - box_length ; col <= column_position[row] + box_length ; col++ )
00217 {
00218 if ( col < 0 || col >= ns_image -> lx )
00219 {
00220 cpl_msg_error ("slopeOfSpectrum","wrong spectrum position or box_length given in row: %d\n", row) ;
00221 cpl_free (xdat) ;
00222 cpl_free (wdat) ;
00223 cpl_free (mpar) ;
00224 destroyFitParams(dec_par) ;
00225 destroyVector( line ) ;
00226 return FLAG ;
00227 }
00228 else if ( qfits_isnan(ns_image->data[col+row*ns_image->lx]) )
00229 {
00230 bad_ind = 1 ;
00231 }
00232 else
00233 {
00234 line -> data[counter] = ns_image -> data[col + row*ns_image->lx] ;
00235 counter++ ;
00236 }
00237 }
00238
00239
00240 if ( bad_ind == 1 )
00241 {
00242 cpl_msg_warning ("slopeOfSpectrum","sorry, bad pixel inside fitting box in row: %d\n", row) ;
00243 cpl_free (xdat) ;
00244 cpl_free (wdat) ;
00245 cpl_free (mpar) ;
00246 destroyFitParams(dec_par) ;
00247 destroyVector( line ) ;
00248 continue ;
00249 }
00250
00251
00252
00253
00254
00255 maxval = -FLT_MAX ;
00256 position = -INT32_MAX ;
00257 for ( i = 0 ; i < counter ; i++ )
00258 {
00259 xdat[i] = i ;
00260 wdat[i] = 1.0 ;
00261 if ( line -> data[i] >= maxval )
00262 {
00263 maxval = line -> data[i] ;
00264 position = i ;
00265 }
00266 }
00267
00268
00269 xdim = XDIM ;
00270 ndat = line -> n_elements ;
00271 numpar = MAXPAR ;
00272 tol = TOL ;
00273 lab = LAB ;
00274 its = ITS ;
00275 (*par).fit_par[1] = fwhm ;
00276 (*par).fit_par[2] = (float) position ;
00277 (*par).fit_par[3] = (float) (line -> data[0] + line -> data[line->n_elements - 1]) / 2.0 ;
00278
00279 (*par).fit_par[0] = maxval - ((*par).fit_par[3]) ;
00280
00281 if ( (*par).fit_par[0] < minDiff )
00282 {
00283 cpl_msg_warning ("slopeOfSpectrum","sorry, negative peak or signal of line too low to fit in row: %d\n", row) ;
00284 cpl_free (xdat) ;
00285 cpl_free (wdat) ;
00286 cpl_free (mpar) ;
00287 destroyFitParams(dec_par) ;
00288 destroyVector( line ) ;
00289 continue ;
00290 }
00291
00292 for ( k = 0 ; k < MAXPAR ; k++ )
00293 {
00294 (*par).derv_par[k] = 0.0 ;
00295 mpar[k] = 1 ;
00296 }
00297
00298
00299 if ( 0 > ( iters = lsqfit_c( xdat, &xdim, line -> data, wdat, &ndat, (*par).fit_par,
00300 (*par).derv_par, mpar, &numpar, &tol, &its, &lab )) )
00301 {
00302 cpl_msg_warning ("slopeOfSpectrum","lsqfit_c: least squares fit failed in row: %d, error no.: %d\n", row, iters) ;
00303 cpl_free (xdat) ;
00304 cpl_free (wdat) ;
00305 cpl_free (mpar) ;
00306 destroyFitParams(dec_par) ;
00307 destroyVector( line ) ;
00308 continue ;
00309 }
00310
00311
00312 if ( (*par).fit_par[0] <= 0. || (*par).fit_par[1] <= 0. ||
00313 (*par).fit_par[2] <= 0. )
00314 {
00315 cpl_msg_warning ("slopeOfSpectrum","negative parameters as fit result, not used! in row: %d\n", row) ;
00316 cpl_free (xdat) ;
00317 cpl_free (wdat) ;
00318 cpl_free (mpar) ;
00319 destroyFitParams(dec_par) ;
00320 destroyVector( line ) ;
00321 continue ;
00322 }
00323
00324
00325 x_position[row] = (float) (column_position[row] - box_length) + (*par).fit_par[2] ;
00326
00327
00328 sig[row] = (*par).derv_par[2] ;
00329
00330
00331 destroyFitParams(dec_par) ;
00332 destroyVector ( line ) ;
00333 cpl_free ( xdat ) ;
00334 cpl_free ( wdat ) ;
00335 cpl_free ( mpar ) ;
00336 }
00337
00338
00339
00340
00341
00342
00343 ndata = 0 ;
00344 for ( row = 0 ; row < ns_image -> lx ; row++ )
00345 {
00346 if ( x_position[row] == 0. || sig[row] == 0. )
00347 {
00348 continue ;
00349 }
00350 else
00351 {
00352 y[ndata] = x_position[row] ;
00353 x[ndata] = (float)row ;
00354 sig[ndata] = sig[row] ;
00355 ndata++ ;
00356 }
00357 }
00358 if ( ndata < 10 )
00359 {
00360 cpl_msg_error("slopeOfSpectrum","not enough positions to do the linear fit\n") ;
00361 return FLAG ;
00362 }
00363
00364
00365 mwt = 0 ;
00366 myfit(x, y, ndata, sig, mwt, &a, &b, &siga, &sigb, &chi2, &q) ;
00367 return b ;
00368 }
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381 OneImage * shiftRows( OneImage * image,
00382 float slope,
00383 int n_order )
00384 {
00385 OneImage * returnImage ;
00386 float xshift ;
00387 int intshift = 0 ;
00388 float sum, new_sum ;
00389 float xnum[n_order + 1] ;
00390 float row_data[image->lx] ;
00391 float corrected_row_data[image->lx] ;
00392 float eval ;
00393 float * imageptr ;
00394 int col, row ;
00395 int firstpos, n_points ;
00396 int i ;
00397 int flag;
00398
00399 if ( image == NullImage )
00400 {
00401 cpl_msg_error("shiftRows","sorry, no image given\n") ;
00402 return NullImage ;
00403 }
00404
00405 if ( slope == 0. )
00406 {
00407 cpl_msg_error("shiftRows","there is no need to shift the image rows!\n") ;
00408 return NullImage ;
00409 }
00410
00411 if ( n_order <= 0 )
00412 {
00413 cpl_msg_error("shiftRows","wrong order of interpolation polynom given!\n") ;
00414 return NullImage ;
00415 }
00416
00417 returnImage = copy_image( image ) ;
00418
00419 n_points = n_order + 1 ;
00420 if ( n_points % 2 == 0 )
00421 {
00422 firstpos = (int)(n_points/2) - 1 ;
00423 }
00424 else
00425 {
00426 firstpos = (int)(n_points/2) ;
00427 }
00428
00429
00430 for ( i = 0 ; i < n_points ; i++ )
00431 {
00432 xnum[i] = i ;
00433 }
00434
00435
00436 for ( row = 0 ; row < image -> ly ; row++ )
00437 {
00438
00439 xshift = slope * (float)( (image -> ly / 2) - row ) ;
00440
00441 intshift = nint(xshift) ;
00442 xshift = xshift - (float)intshift ;
00443
00444 for ( col = 0 ; col < image -> lx ; col++ )
00445 {
00446 corrected_row_data[col] = 0. ;
00447 }
00448 sum = 0. ;
00449
00450 for ( col = 0 ; col < image -> lx ; col++ )
00451 {
00452
00453 if ( intshift < 0 )
00454 {
00455 if ( col-intshift < image->lx )
00456 {
00457 row_data[col] = image->data[col-intshift+row*image->lx] ;
00458 }
00459 else
00460 {
00461 row_data[col] = 0. ;
00462 }
00463 }
00464 else if ( intshift > 0 )
00465 {
00466 if ( col - intshift >= 0 )
00467 {
00468 row_data[col] = image->data[col-intshift+row*image->lx] ;
00469 }
00470 else
00471 {
00472 row_data[col] = 0. ;
00473 }
00474 }
00475 else
00476 {
00477 row_data[col] = image -> data[col+row*image->lx] ;
00478 }
00479
00480
00481 if ( col != 0 && col != image->lx - 1 && !qfits_isnan(row_data[col]) )
00482 {
00483 sum += row_data[col] ;
00484 }
00485 if (qfits_isnan(row_data[col]))
00486 {
00487 row_data[col] = 0. ;
00488 for ( i = col - firstpos ; i < col - firstpos + n_points ; i++ )
00489 {
00490 if ( i < 0 ) continue ;
00491 if ( i >= image->lx) continue ;
00492 corrected_row_data[i] = ZERO ;
00493 }
00494 }
00495 }
00496
00497
00498
00499
00500
00501
00502 new_sum = 0. ;
00503 for ( col = 0 ; col < image -> lx ; col++ )
00504 {
00505
00506
00507
00508
00509
00510
00511 if ( qfits_isnan(corrected_row_data[col]) )
00512 {
00513 continue ;
00514 }
00515 if ( col - firstpos < 0 )
00516 {
00517 imageptr = &row_data[0] ;
00518 eval = (float)col - xshift ;
00519 }
00520 else if ( col - firstpos + n_points >= image -> lx )
00521 {
00522 imageptr = &row_data[image->lx - n_points] ;
00523 eval = (float)(col + n_points - image->lx) - xshift ;
00524 }
00525 else
00526 {
00527 imageptr = &row_data[col-firstpos] ;
00528 eval = (float)firstpos - xshift ;
00529 }
00530
00531 flag=0;
00532 corrected_row_data[col]=nev_ille(xnum,imageptr,n_order,eval,&flag);
00533
00534
00535
00536 if ( col != 0 && col != image->lx - 1 && !qfits_isnan(corrected_row_data[col]) )
00537 {
00538 new_sum += corrected_row_data[col] ;
00539 }
00540 }
00541
00542 if ( new_sum == 0. )
00543 {
00544 new_sum = 1. ;
00545 }
00546 for ( col = 0 ; col < image -> lx ; col++ )
00547 {
00548 if ( qfits_isnan(corrected_row_data[col]))
00549 {
00550 returnImage->data[col+row*image->lx] = ZERO ;
00551 }
00552 else
00553 {
00554
00555
00556
00557
00558
00559 returnImage->data[col+row*image->lx] = corrected_row_data[col] ;
00560 }
00561 }
00562 }
00563 return returnImage ;
00564 }
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576 void parameterToAscii ( float * parameter,
00577 int n,
00578 char * filename )
00579 {
00580 FILE * fp ;
00581 int i ;
00582
00583 if ( parameter == NULL || filename == NULL || n <= 0 )
00584 {
00585 cpl_msg_error ("parameterToAscii","input is missing or wrong!\n") ;
00586 return ;
00587 }
00588
00589 if ( NULL == (fp = fopen ( filename, "w" ) ) )
00590 {
00591 cpl_msg_error("parameterToAscii","cannot open %s\n", filename) ;
00592 return ;
00593 }
00594
00595 for ( i = 0 ; i < n ; i++ )
00596 {
00597 fprintf (fp, "%le\n", parameter[i] ) ;
00598 }
00599 fclose (fp ) ;
00600 }
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610 float * asciiToParameter ( char * filename,
00611 int * n )
00612 {
00613 FILE * fp ;
00614 float * parameter ;
00615 int i ;
00616
00617 if ( filename == NULL || n == NULL )
00618 {
00619 cpl_msg_error ("asciiToParameter","input is missing or wrong!\n") ;
00620 return NULL ;
00621 }
00622
00623 if ( NULL == (fp = fopen ( filename, "r" ) ) )
00624 {
00625 cpl_msg_error("asciiToParameter","cannot open %s\n", filename) ;
00626 return NULL ;
00627 }
00628
00629
00630
00631 if ( NULL == ( parameter = (float*) cpl_calloc (ESTIMATE, sizeof(float)) ) )
00632 {
00633 cpl_msg_error ("asciiToParameter","cannot allocate memory\n") ;
00634 return NULL ;
00635 }
00636
00637 i = 0 ;
00638 while ( fscanf(fp, "%g\n", ¶meter[i]) != EOF )
00639 {
00640 i++ ;
00641 }
00642 *n = i ;
00643
00644 fclose (fp ) ;
00645
00646 return parameter ;
00647 }
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672 double * curvatureOfSpectrum( OneImage * ns_image,
00673 int order,
00674 int box_length,
00675 int left_pos,
00676 int right_pos,
00677 float fwhm,
00678 float minDiff )
00679 {
00680 int i, k, row, col ;
00681 int counter, iters ;
00682 int xdim, ndat, its, numpar ;
00683 float maxval ;
00684 float tol, lab ;
00685 float col_value[ns_image->ly] ;
00686 float column_value[ns_image->ly] ;
00687 pixelvalue col_position[ns_image->ly] ;
00688 int column_position[ns_image->ly] ;
00689 int col_median ;
00690 float x_position[ns_image->ly] ;
00691 float * xdat, * wdat ;
00692 int * mpar ;
00693 Vector * line ;
00694 FitParams ** dec_par ;
00695 FitParams * par ;
00696 int position ;
00697 int ndata ;
00698 int bad_ind ;
00699 dpoint * list ;
00700 double * coeffs ;
00701 double offset ;
00702
00703
00704 if ( ns_image == NullImage )
00705 {
00706 cpl_msg_error("curvatureOfSpectrum","sorry, no image given\n") ;
00707 return NULL ;
00708 }
00709 if ( box_length <= 1 || box_length >= right_pos - left_pos )
00710 {
00711 cpl_msg_error("curvatureOfSpectrum","wrong box length given\n") ;
00712 return NULL ;
00713 }
00714 if ( fwhm < 1. || fwhm > 10. )
00715 {
00716 cpl_msg_error("curvatureOfSpectrum","wrong full width at half maximum given\n") ;
00717 return NULL ;
00718 }
00719 if ( left_pos < 0 || right_pos <= left_pos || right_pos > ns_image ->lx )
00720 {
00721 cpl_msg_error("curvatureOfSpectrum","wrong left and right positions\n") ;
00722 return NULL ;
00723 }
00724 if ( minDiff < 1. )
00725 {
00726 cpl_msg_error("curvatureOfSpectrum","wrong amplitude threshold given! \n") ;
00727 return NULL ;
00728 }
00729
00730
00731 for ( row = 0 ; row < ns_image -> ly ; row++ )
00732 {
00733 col_value[row] = -FLT_MAX ;
00734 col_position[row] = -1. ;
00735
00736 for ( col = left_pos ; col < right_pos ; col++ )
00737 {
00738 if ( ns_image -> data[col+row*ns_image->lx] > col_value[row] )
00739 {
00740 col_value[row] = ns_image -> data[col+row*ns_image->lx] ;
00741 col_position[row] = (pixelvalue)col ;
00742 }
00743 }
00744 }
00745
00746
00747
00748 col_median = (int)median(col_position, right_pos - left_pos) ;
00749
00750
00751 for ( row = 0 ; row < ns_image -> ly ; row++ )
00752 {
00753 x_position[row] = 0. ;
00754 column_value[row] = -FLT_MAX ;
00755 column_position[row] = -1 ;
00756 for ( col = col_median - box_length ; col <= col_median + box_length ; col++ )
00757 {
00758 if ( ns_image -> data[col+row*ns_image->lx] > column_value[row] )
00759 {
00760 column_value[row] = ns_image -> data[col+row*ns_image->lx] ;
00761 column_position[row] = col ;
00762 }
00763 }
00764
00765
00766 if ( NULL == (line = newVector (2*box_length + 1)) )
00767 {
00768 cpl_msg_error ("curvatureOfSpectrum","cannot allocate new Vector in row: %d\n", row) ;
00769 return NULL ;
00770 }
00771
00772
00773 xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00774 wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00775 mpar = (int *) cpl_calloc( MAXPAR, sizeof (int) ) ;
00776 dec_par = newFitParams(1) ;
00777 par = dec_par[0];
00778
00779 counter = 0 ;
00780 bad_ind = 0 ;
00781
00782 for ( col = column_position[row] - box_length ; col <= column_position[row] + box_length ; col++ )
00783 {
00784 if ( col < 0 || col >= ns_image -> lx )
00785 {
00786 cpl_msg_error ("curvatureOfSpectrum","wrong spectrum position or box_length given in row: %d\n", row) ;
00787 cpl_free (xdat) ;
00788 cpl_free (wdat) ;
00789 cpl_free (mpar) ;
00790 destroyFitParams(dec_par) ;
00791 destroyVector( line ) ;
00792 return NULL ;
00793 }
00794 else if ( qfits_isnan(ns_image->data[col+row*ns_image->lx]) )
00795 {
00796 bad_ind = 1 ;
00797 }
00798 else
00799 {
00800 line -> data[counter] = ns_image -> data[col + row*ns_image->lx] ;
00801 counter++ ;
00802 }
00803 }
00804
00805
00806 if ( bad_ind == 1 )
00807 {
00808 cpl_msg_warning ("curvatureOfSpectrum","sorry, bad pixel inside fitting box in row: %d\n", row) ;
00809 cpl_free (xdat) ;
00810 cpl_free (wdat) ;
00811 cpl_free (mpar) ;
00812 destroyFitParams(dec_par) ;
00813 destroyVector( line ) ;
00814 continue ;
00815 }
00816
00817
00818
00819
00820
00821 maxval = -FLT_MAX ;
00822 position = -INT32_MAX ;
00823 for ( i = 0 ; i < counter ; i++ )
00824 {
00825 xdat[i] = i ;
00826 wdat[i] = 1.0 ;
00827 if ( line -> data[i] >= maxval )
00828 {
00829 maxval = line -> data[i] ;
00830 position = i ;
00831 }
00832 }
00833
00834
00835 xdim = XDIM ;
00836 ndat = line -> n_elements ;
00837 numpar = MAXPAR ;
00838 tol = TOL ;
00839 lab = LAB ;
00840 its = ITS ;
00841 (*par).fit_par[1] = fwhm ;
00842 (*par).fit_par[2] = (float) position ;
00843 (*par).fit_par[3] = (float) (line -> data[0] + line -> data[line->n_elements - 1]) / 2.0 ;
00844
00845 (*par).fit_par[0] = maxval - ((*par).fit_par[3]) ;
00846
00847 if ( (*par).fit_par[0] < minDiff )
00848 {
00849 cpl_msg_warning ("curvatureOfSpectrum","sorry, negative peak or signal of line too low to fit in row: %d\n", row) ;
00850 cpl_free (xdat) ;
00851 cpl_free (wdat) ;
00852 cpl_free (mpar) ;
00853 destroyFitParams(dec_par) ;
00854 destroyVector( line ) ;
00855 continue ;
00856 }
00857
00858 for ( k = 0 ; k < MAXPAR ; k++ )
00859 {
00860 (*par).derv_par[k] = 0.0 ;
00861 mpar[k] = 1 ;
00862 }
00863
00864
00865 if ( 0 > ( iters = lsqfit_c( xdat, &xdim, line -> data, wdat, &ndat, (*par).fit_par,
00866 (*par).derv_par, mpar, &numpar, &tol, &its, &lab )) )
00867 {
00868 cpl_msg_warning ("curvatureOfSpectrum","lsqfit_c: least squares fit failed in row: %d, error no.: %d\n", row, iters) ;
00869 cpl_free (xdat) ;
00870 cpl_free (wdat) ;
00871 cpl_free (mpar) ;
00872 destroyFitParams(dec_par) ;
00873 destroyVector( line ) ;
00874 continue ;
00875 }
00876
00877
00878 if ( (*par).fit_par[0] <= 0. || (*par).fit_par[1] <= 1. ||
00879 (*par).fit_par[2] <= 0. )
00880 {
00881 cpl_msg_warning ("curvatureOfSpectrum","negative parameters as fit result, not used! in row: %d\n", row) ;
00882 cpl_free (xdat) ;
00883 cpl_free (wdat) ;
00884 cpl_free (mpar) ;
00885 destroyFitParams(dec_par) ;
00886 destroyVector( line ) ;
00887 continue ;
00888 }
00889
00890
00891 x_position[row] = (float) (column_position[row] - box_length) + (*par).fit_par[2] ;
00892 printf("%d %f %f\n",row, (*par).fit_par[1], x_position[row] ) ;
00893
00894
00895 destroyFitParams(dec_par) ;
00896 destroyVector ( line ) ;
00897 cpl_free ( xdat ) ;
00898 cpl_free ( wdat ) ;
00899 cpl_free ( mpar ) ;
00900 }
00901
00902 if ( NULL == ( list = (dpoint*) cpl_calloc ( ns_image->ly, sizeof (dpoint)) ) )
00903 {
00904 cpl_msg_error("curvatureOfSpectrum","could not allocate memory!\n") ;
00905 return NULL ;
00906 }
00907
00908
00909
00910
00911
00912
00913 offset = (double) ns_image->ly/2. ;
00914 ndata = 0 ;
00915 for ( row = 0 ; row < ns_image -> ly ; row++ )
00916 {
00917 if ( x_position[row] == 0. )
00918 {
00919 continue ;
00920 }
00921 else
00922 {
00923 list[ndata].y = (double)x_position[row] ;
00924 list[ndata].x = (double)row - offset ;
00925 ndata++ ;
00926 }
00927 }
00928
00929
00930 if ( NULL == (coeffs = fit_1d_poly(order, list, ndata, NULL)) )
00931 {
00932 cpl_msg_error("curvatureOfSpectrum","eclipse function fit_1d_poly() did not work!\n") ;
00933 return NULL ;
00934 }
00935 cpl_free ( list ) ;
00936
00937 return coeffs ;
00938 }
00939
00940
00941
00942
00943
00944
00945
00946 OneImage * image_warp( OneImage * image,
00947 char * kernel_type,
00948 char * poly_table )
00949 {
00950 OneImage * warped;
00951 char ch,temp[150];
00952 char poly_string[150];
00953 int count=0;
00954 FILE * poly_in;
00955
00956 poly2d * poly_u;
00957 poly2d * poly_v;
00958
00959 if((poly_in=fopen(poly_table,"r"))==NULL) {
00960 cpl_msg_error("Cannot open file %s", poly_table);
00961 return NULL;
00962 }
00963 fscanf(poly_in,"%s %s %s %s",temp,temp,temp,temp);
00964 ch=fgetc(poly_in);
00965 while(ch!='p'){
00966 ch=fgetc(poly_in);
00967 poly_string[count]=ch;
00968 count++;
00969 }
00970 poly_string[count-2]='\0';
00971
00972 cpl_msg_info("%s",poly_string);
00973
00974 fclose(poly_in);
00975
00976 poly_u = poly2d_build_from_string(poly_string);
00977
00978
00979 if (poly_u == NULL) {
00980 cpl_msg_error("image_warp","cannot read 2D poly from arc table") ;
00981 return NULL ;
00982 }
00983 poly_v=poly2d_build_from_string("0 1 1.0") ;
00984
00985 warped = warp_image_generic(image, kernel_type, poly_u, poly_v);
00986 if (poly_u!=NULL) poly2d_free(poly_u);
00987 if (poly_v!=NULL) poly2d_free(poly_v);
00988
00989 return warped;
00990 }
00991
00992 #include <cpl_table.h>
00993 OneImage * image_warp_fits( OneImage * image,
00994 char * kernel_type,
00995 char * poly_table )
00996 {
00997 const char* _id="image_warp_fits";
00998 OneImage * warped;
00999 char poly_string[150];
01000 char tmp_string[50];
01001 int nraw=0;
01002
01003 poly2d * poly_u;
01004 poly2d * poly_v;
01005 int* degx;
01006 int* degy;
01007 double* coef;
01008 int i=0;
01009 cpl_table* poly_tbl=NULL;
01010
01011 if(is_fits_file(poly_table) != 1) {
01012 cpl_msg_error(_id,"Input file %s is not FITS",poly_table);
01013 return NULL;
01014 }
01015 poly_tbl = cpl_table_load(poly_table,1,0);
01016 nraw = cpl_table_get_nrow(poly_tbl);
01017 degx = cpl_table_get_data_int(poly_tbl,"degx");
01018 degy = cpl_table_get_data_int(poly_tbl,"degy");
01019 coef = cpl_table_get_data_double(poly_tbl,"coeff");
01020
01021 sprintf(poly_string,"%2d%2d%s%g ",degx[0],degy[0]," ",coef[0]);
01022 for(i=1;i<nraw;i++) {
01023 sprintf(tmp_string,"%2d%2d%s%g ",degx[i],degy[i]," ", coef[i]);
01024 strcat(poly_string,tmp_string);
01025 }
01026
01027
01028
01029
01030 poly_u = poly2d_build_from_string(poly_string);
01031
01032
01033 if (poly_u == NULL) {
01034 cpl_msg_error(_id,"cannot read 2D poly from arc table") ;
01035 return NULL ;
01036 }
01037 poly_v=poly2d_build_from_string("0 1 1.0") ;
01038
01039
01040 warped = warp_image_generic(image, kernel_type, poly_u, poly_v);
01041 if (poly_u!=NULL) poly2d_free(poly_u);
01042 if (poly_v!=NULL) poly2d_free(poly_v);
01043 cpl_table_delete(poly_tbl);
01044 return warped;
01045 }
01046
01047
01048
01049
01050
01051