00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "wave_cal.h"
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 int wave_cal (cpl_parameterlist* config, cpl_frameset* sof) {
00070
00071 const char* _id = "wave_cal";
00072
00073 wave_config * cfg ;
00074 OneImage * map ;
00075
00076 OneImage * im ;
00077 FitParams** par;
00078 float* wave;
00079 float* intens;
00080 float** bcoefs;
00081
00082 int check = 0;
00083 int lx = 0;
00084 int ly = 0;
00085 int n_lines=0;
00086 int i = 0;
00087 int j = 0;
00088 int n = 0;
00089 int* n_found_lines;
00090 int* sum_pointer;
00091 int** row_clean;
00092 float** wavelength_clean;
00093
00094 int sum=0;
00095 int er=0;
00096 float a=0.;
00097
00098 float val1=0.;
00099 float val2=0.;
00100 float val3=0.;
00101 float val4=0.;
00102 float val5=0.;
00103
00104 float** val_x;
00105 float** val_y;
00106
00107 float* vec_x;
00108 float* vec_y;
00109 int cnt = 0 ;
00110
00111 float** slit_pos;
00112 char* tmp_str;
00113 int readsum=0;
00114 int fit=0;
00115 float** val;
00116 FILE * fp ;
00117 FILE * bcoefs_file;
00118 FILE * bcoefs_list;
00119 FILE * par_list;
00120 FILE * slit_list;
00121 FILE * slitpos_file;
00122 cpl_frameset* raw=NULL;
00123
00124
00125
00126
00127
00128
00129
00130
00131 cfg = parse_cpl_input_wave(config,sof,&raw) ;
00132 if (cfg == NULL)
00133 {
00134 cpl_msg_error ("wave_cal:","Error could not parse cpl input!") ;
00135 return -1 ;
00136 }
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 if(is_fits_file (cfg->inFrame) != 1) {
00164 cpl_msg_error(_id,"Input file %s is not FITS",cfg->inFrame);
00165 return -1;
00166 }
00167
00168 if (cfg->slitposIndicator == 1 && cfg->estimateIndicator == 1) {
00169 if (is_fits_file(cfg->slitposGuessName) != 1) {
00170 cpl_msg_error("wave_cal:","slitlet position guess list not given!");
00171 exit(-1);
00172 }
00173 }
00174
00175 if (cfg->calibIndicator == 0 && cfg->wavemapInd == 1) {
00176 if (is_fits_file(cfg->coeffsName) != 1) {
00177 cpl_msg_error("wave_cal:","coefficients list not given!");
00178 exit(-1);
00179
00180 }
00181 }
00182 if (cfg->slitposIndicator == 1) {
00183 if (cfg->calibIndicator != 1 && cfg->estimateIndicator != 1) {
00184 if (is_fits_file(cfg->paramsList) != 1) {
00185 cpl_msg_error("wave_cal:","parameter list not given!");
00186 exit(-1);
00187 }
00188 }
00189 }
00190
00191
00192 im = load_image(cfg->inFrame);
00193 if (im == NULL) {
00194 cpl_msg_error ("wave_cal","could not load image");
00195 exit(-1);
00196 }
00197
00198 lx = im->lx;
00199 ly = im->ly;
00200
00201 if (cfg->calibIndicator == 1 || cfg->wavemapInd == 1) {
00202 if (is_fits_file(cfg->lineList) != 1) {
00203 cpl_msg_error("wave_cal:","linelist not given!");
00204 exit(-1);
00205 }
00206 }
00207
00208
00209
00210 if ( NULL == (fp = fopen (cfg->lineList, "r" ) ) )
00211 {
00212 cpl_msg_error("wave_cal:","cannot open %s", cfg->lineList) ;
00213 return -1 ;
00214 }
00215
00216 cnt = 0 ;
00217 while ( fscanf( fp, "%f %f", &vec_x[cnt], &vec_y[cnt] ) != EOF )
00218 {
00219 cnt ++ ;
00220 }
00221 n= cnt ;
00222
00223 vec_x = new_float_array (n);
00224 vec_y = new_float_array (n);
00225
00226 cnt=0;
00227 while ( fscanf( fp, "%f %f", &vec_x[cnt], &vec_y[cnt] ) != EOF )
00228 {}
00229
00230
00231
00232 wave = new_float_array (n);
00233 intens = new_float_array (n);
00234 if (wave == NULL || intens == NULL) {
00235 cpl_msg_error ( "wave_cal","could not allocate memory for the line list values" );
00236 exit(-1);
00237 }
00238
00239
00240 n_lines = readList(cfg->lineList, wave, intens);
00241 if (-1 == n_lines) {
00242 cpl_msg_error("wave_cal","readList failed!");
00243 exit(-1);
00244 }
00245
00246
00247 fclose(fp);
00248 bcoefs = new_2Dfloat_array(cfg->nrDispCoefficients, cfg->nrCoefCoefficients);
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 if (cfg->calibIndicator == 1) {
00260
00261 n_found_lines = new_int_array(lx);
00262 row_clean = new_2Dint_array(lx, n_lines);
00263 wavelength_clean = new_2Dfloat_array(lx, n_lines);
00264 sum_pointer = new_int_array(1);
00265
00266
00267 intarray_set_value(sum_pointer, 0, 0);
00268 check = findLines(im, wave, intens, n_lines, row_clean,
00269 wavelength_clean, cfg->guessBeginWavelength,
00270 cfg->guessDispersion1, cfg->guessDispersion2,
00271 cfg->mindiff, cfg->halfWidth,
00272 n_found_lines, cfg->sigma, sum_pointer );
00273 if (-1 == check) {
00274 cpl_msg_error ("wave_cal","findLines failed!");
00275 exit(-1);
00276 }
00277
00278
00279
00280
00281
00282
00283
00284 sum = intarray_get_value(sum_pointer,0);
00285
00286
00287 par = newFitParams( sum );
00288 if (par == NULL) {
00289 cpl_msg_error("wave_cal","newFitParams failed!");
00290 exit(-1);
00291 }
00292
00293
00294
00295
00296
00297
00298 er = wavelengthCalibration(im, par, bcoefs, wave, n_lines, row_clean,
00299 wavelength_clean, n_found_lines,
00300 cfg->guessDispersion1,
00301 cfg->halfWidth, cfg->minAmplitude,
00302 cfg->maxResidual, cfg->fwhm,
00303 cfg->nrDispCoefficients,
00304 cfg->nrCoefCoefficients,
00305 cfg->sigmaFactor,
00306 cfg->pixel_tolerance);
00307 if (-1 == er) {
00308 cpl_msg_error ("wave_cal","wavelengthCalibration failed!");
00309 exit(-1);
00310 }
00311
00312
00313
00314
00315
00316 tmp_str = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char));
00317 if (cfg->writeCoeffsInd == 1) {
00318
00319 if ( NULL == (fp = fopen ( cfg->coeffsName, "w" ) ) )
00320 {
00321 cpl_msg_error("wave_cal:","cannot open %s", cfg->coeffsName) ;
00322 return -1 ;
00323 }
00324 for (i = 0; i< cfg->nrDispCoefficients; i++) {
00325 for (j=0; j< cfg->nrCoefCoefficients; j++) {
00326 a = array2D_get_value(bcoefs, i, j);
00327 sprintf(tmp_str,"%f",a);
00328 strcat(tmp_str," \n");
00329 fprintf(bcoefs_file,"%s",tmp_str);
00330 }
00331
00332 if ( 0 == fclose ( bcoefs_file ) )
00333 {
00334 cpl_msg_error("wave_cal:","cannot close %s", cfg->coeffsName) ;
00335 return -1 ;
00336 }
00337 }
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350 }
00351
00352
00353
00354
00355 if (cfg->writeParInd == 1) {
00356 dumpFitParamsToAscii(par, cfg->paramsList);
00357 }
00358
00359 destroy_2Darray ( wavelength_clean, lx );
00360 destroy_2Dintarray (row_clean, lx);
00361 destroy_intarray ( n_found_lines );
00362 destroy_intarray ( sum_pointer );
00363
00364 }
00365
00366
00367
00368
00369
00370
00371
00372 if (cfg->wavemapInd == 1) {
00373
00374
00375
00376 if ( NULL == (bcoefs_list = fopen (cfg->coeffsName, "r" ) ) )
00377 {
00378 cpl_msg_error("wave_cal:","cannot open %s", cfg->coeffsName) ;
00379 return -1 ;
00380 }
00381
00382 for (i = 0; i< cfg->nrDispCoefficients; i++) {
00383 for (j = 0; j< cfg->nrCoefCoefficients; j++) {
00384 fscanf( bcoefs_list, "%f", &val[i][j]);
00385 array2D_set_value (bcoefs, val[i][j], i, j);
00386 }
00387 }
00388 fclose(bcoefs_list);
00389
00390 map = waveMap ( im, bcoefs, cfg->nrDispCoefficients,
00391 cfg->nrCoefCoefficients,
00392 wave, intens, n_lines, cfg->magFactor );
00393 if ( map == NULL ) {
00394 cpl_msg_error("wave_cal:","waveMap failed!");
00395 exit(-1);
00396 }
00397 save_image_to_fits( map, cfg->outName, -32 );
00398
00399
00400 destroy_image(map);
00401
00402 }
00403
00404
00405
00406
00407
00408
00409 if (cfg->slitposIndicator == 1) {
00410 if (cfg->calibIndicator != 1 && cfg->estimateIndicator != 1) {
00411
00412
00413
00414
00415
00416 if ( NULL == (par_list = fopen (cfg->paramsList, "r" ) ) )
00417 {
00418 cpl_msg_error("wave_cal:","cannot open %s", cfg->paramsList) ;
00419 return -1 ;
00420 }
00421
00422 cnt = 0 ;
00423 while ( fscanf( par_list, "%f %f %f %f %f",
00424 &val1,&val2,&val3,&val4,&val5 ) != EOF )
00425 {
00426 cnt ++ ;
00427 }
00428 n= cnt ;
00429 readsum = val5;
00430 fclose(par_list);
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440 par = newFitParams( readsum );
00441 if (par == NULL ) {
00442 cpl_msg_error("wave_cal","newFitParams failed!");
00443 exit(-1);
00444 }
00445 dumpAsciiToFitParams(par, cfg->paramsList);
00446
00447 }
00448
00449 slit_pos = new_2Dfloat_array(32,2);
00450
00451
00452
00453 if (cfg->fitBoltzIndicator == 1) {
00454 if (cfg->estimateIndicator == 1) {
00455
00456
00457
00458 if ( NULL == (slit_list = fopen (cfg->slitposGuessName, "r" ) ) )
00459 {
00460 cpl_msg_error("wave_cal:","cannot open %s", cfg->slitposGuessName) ;
00461 return -1 ;
00462 }
00463
00464
00465 for (i =0; i< 32; i++){
00466 fscanf(slit_list, "%f %f",
00467 &val_x[i][0],&val_y[i][1] );
00468
00469 array2D_set_value(slit_pos, val_x[i][0], i, 0);
00470 array2D_set_value(slit_pos, val_y[i][1], i, 1);
00471
00472 }
00473 fclose(slit_list);
00474
00475
00476 fit = fitSlitsBoltzWithEstimate( im, slit_pos, cfg->boxLength,
00477 cfg->yBox, cfg->diffTol,
00478 cfg->loPos, cfg->hiPos );
00479 if (fit < 0){
00480 cpl_msg_error ( "wave_cal:","fitSlitsBoltz failed" );
00481 exit(-1);
00482 }
00483
00484 } else {
00485
00486 fit = fitSlitsBoltz( im, par, slit_pos,
00487 cfg->boxLength, cfg->yBox,
00488 cfg->diffTol );
00489
00490 }
00491 } else if ( cfg->fitEdgeIndicator == 1) {
00492 if (cfg->estimateIndicator == 1) {
00493
00494
00495 if (NULL == (slit_list = fopen (cfg->slitposGuessName, "r" )))
00496 {
00497 cpl_msg_error("wave_cal:","cannot open %s", cfg->slitposGuessName) ;
00498 return -1 ;
00499 }
00500
00501
00502 for (i =0; i< 32; i++){
00503 fscanf(slit_list, "%f %f",
00504 &val_x[i][0],&val_y[i][1] );
00505
00506 array2D_set_value(slit_pos, val_x[i][0], i, 0);
00507 array2D_set_value(slit_pos, val_y[i][1], i, 1);
00508
00509 }
00510 fclose(slit_list);
00511
00512
00513 fit = fitSlitsEdgeWithEstimate( im, slit_pos, cfg->boxLength,
00514 cfg->yBox, cfg->diffTol,
00515 cfg->loPos, cfg->hiPos );
00516 if (fit < 0) {
00517 cpl_msg_error ( "wave_cal:","fitSlitsBoltz failed" );
00518 exit(-1);
00519 }
00520 } else {
00521 fit = fitSlitsEdge( im, par, slit_pos,
00522 cfg->boxLength, cfg->yBox,
00523 cfg->diffTol );
00524 if (fit < 0) {
00525 cpl_msg_error ( "wave_cal:","fitSlitsEdge failed" );
00526 exit(-1);
00527 }
00528
00529 }
00530 } else {
00531 cpl_msg_error ( "wave_cal:","no indication of desired fit function given, will exit" );
00532 exit(-1);
00533
00534 }
00535
00536
00537
00538 if (NULL == (slitpos_file = fopen (cfg->slitposName, "w" )))
00539 {
00540 cpl_msg_error("wave_cal:","cannot open %s", cfg->slitposName) ;
00541 return -1 ;
00542 }
00543
00544
00545 for (i =0; i< 32; i++){
00546 fprintf(slitpos_file, "%f %s %f\n",
00547 &slit_pos[i][0]," ", &slit_pos[i][1] );
00548 }
00549 fclose(slitpos_file);
00550
00551 destroy_2Darray ( slit_pos, 32 );
00552 }
00553
00554
00555
00556
00557
00558
00559 if ((cfg->slitposIndicator == 1 && cfg->estimateIndicator != 1) ||
00560 cfg->calibIndicator == 1){
00561 destroyFitParams(par);
00562 }
00563 if (cfg->calibIndicator == 1 || cfg->wavemapInd == 1) {
00564 destroy_array(intens);
00565 destroy_array(wave);
00566 destroy_2Darray ( bcoefs, cfg->nrDispCoefficients );
00567 }
00568 destroy_image ( im );
00569 wavecal_free(cfg);
00570 return 0;
00571
00572 }