00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "focus_det.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 int focus_det (cpl_parameterlist* config,cpl_frameset* sof)
00044 {
00045 const char* _id = "focus";
00046 focus_config * cfg ;
00047
00048 fits_header * head;
00049
00050 OneImage * im ;
00051 OneImage * img ;
00052 OneImage * clean_im ;
00053 OneImage * out_image ;
00054 OneCube * cube ;
00055 OneCube * out_cube ;
00056
00057 int npar =7;
00058 char* name;
00059 char* tbl_name=NULL;
00060
00061 int i = 0;
00062 int fra=0;
00063 int nslits=0;
00064
00065 float** slit_edges;
00066 int* mpar;
00067 double* derv_par;
00068 double* fit_par;
00069 float* correct_dist;
00070 float* vec_x;
00071 float* vec_y;
00072
00073 double par=0.;
00074 int iters=0;
00075 float* distances;
00076 float fc=0;
00077 char* out_Name_No;
00078 char* gaussplotNameNo;
00079 char* string_par;
00080 char* tmp_string_par;
00081
00082 FILE * fitpar_file ;
00083 FILE * slit_list ;
00084 FILE * dist_list ;
00085
00086 int cnt = 0 ;
00087 float edge_x=0;
00088 float edge_y=0;
00089 float tmp_float=0;
00090 cpl_table* tbl_focus=NULL;
00091 cpl_propertylist* tbl_focus_phead=NULL;
00092 cpl_propertylist* tbl_focus_head=NULL;
00093 char* col_name=NULL;
00094 char* tbl_slitpos_name=NULL;
00095 cpl_table* tbl_slitpos=NULL;
00096 char* tbl_distances_name=NULL;
00097 cpl_table* tbl_distances=NULL;
00098 char* tbl_firstcol_name=NULL;
00099 cpl_table* tbl_firstcol=NULL;
00100 int* status=NULL;
00101
00102 cpl_image* o_img=NULL;
00103 cpl_image* g_img=NULL;
00104 cpl_frameset* stk=NULL;
00105 int check=0;
00106
00107
00108
00109
00110
00111
00112
00113 cfg = parse_cpl_input_focus(config,sof,&stk) ;
00114
00115 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00116 cpl_msg_error(_id,"error parsing cpl input");
00117 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00118 return -1;
00119 }
00120 if (cfg == NULL)
00121 {
00122 cpl_msg_error (_id,"could not parse cpl input!\n") ;
00123 return -1 ;
00124 }
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142 if (is_fits_file(cfg->poslist) != 1) {
00143 cpl_msg_error(_id,"poslist not given!");
00144 return -1;
00145 }
00146 out_Name_No = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00147 gaussplotNameNo = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00148
00149
00150
00151 if ( NULL == (fitpar_file = fopen (cfg->fitlist, "w" ) ) )
00152 {
00153 cpl_msg_error(_id,"cannot open %s\n", cfg->fitlist) ;
00154 return -1 ;
00155 }
00156
00157 col_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00158 tbl_focus = cpl_table_new(cfg->nframes);
00159 for (i=0; i< npar; i++) {
00160 sprintf(col_name,"%s%d","par",i);
00161 cpl_table_new_column(tbl_focus,col_name, CPL_TYPE_DOUBLE);
00162 cpl_table_set_column_format(tbl_focus,col_name,"%14.10f");
00163
00164 }
00165
00166 printf("nframes=%d\n",cfg->nframes);
00167 for (fra=0; fra < cfg->nframes; fra++) {
00168
00169 name = cfg->inFrameList[fra];
00170 if(is_fits_file(name) != 1 ) {
00171 cpl_msg_error(_id,"Input file %s is not FITS",cfg->inFrameList[fra]);
00172 return -1;
00173 }
00174 im = load_image(name);
00175
00176 if (im == NULL) {
00177 cpl_msg_error(_id, " could not load inFrame\n" );
00178 return -1;
00179 }
00180
00181
00182
00183 head = fits_read_header(name);
00184
00185
00186
00187
00188
00189
00190 clean_im = cleanMeanOfColumns( im, cfg->lo_reject, cfg->hi_reject );
00191 if (clean_im == NULL) {
00192 cpl_msg_error(_id," could not take a clean mean of the spectra!\n");
00193 return -1;
00194 }
00195 destroy_image (im);
00196
00197
00198 correct_dist = new_float_array(cfg->nslits) ;
00199
00200 if (cfg->northsouthInd == 0) {
00201
00202
00203 tbl_slitpos_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00204 strcpy(tbl_slitpos_name,cfg->poslist);
00205 tbl_slitpos = cpl_table_load(tbl_slitpos_name,1,0);
00206
00207
00208 for (i =0 ; i< cfg->nslits; i++){
00209
00210 edge_x=cpl_table_get_float(tbl_slitpos,"start",i,status);
00211 edge_y=cpl_table_get_float(tbl_slitpos,"end",i,status);
00212 array2D_set_value(slit_edges,edge_x,i,0);
00213 array2D_set_value(slit_edges,edge_y,i,1);
00214 }
00215
00216 cpl_table_delete(tbl_slitpos);
00217 cpl_free(tbl_slitpos_name);
00218
00219 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00220 cpl_msg_error(_id,"error writing tbl %s",tbl_slitpos_name);
00221 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00222 return -1;
00223 }
00224 slit_edges = new_2Dfloat_array(cfg->nslits, 2) ;
00225
00226 if ( NULL == (slit_list = fopen (cfg->poslist, "r" ) ) )
00227 {
00228 cpl_msg_error(_id,"cannot open %s\n", cfg->poslist) ;
00229 return -1 ;
00230 }
00231
00232 cnt = 0 ;
00233 while ( fscanf( slit_list, "%f %f", &vec_x[cnt], &vec_y[cnt] ) != EOF )
00234 {
00235 cnt ++ ;
00236 }
00237
00238 vec_x=new_float_array(cnt);
00239 vec_y=new_float_array(cnt);
00240
00241 cnt=0;
00242 while ( fscanf( slit_list, "%f %f", &vec_x[cnt], &vec_y[cnt] ) != EOF )
00243 {
00244 cnt ++ ;
00245 }
00246
00247
00248 if (cnt != cfg->nslits) {
00249 cpl_msg_error(_id," wrong file format!\n");
00250 fclose(slit_list);
00251 return -1;
00252 }
00253 fclose(slit_list);
00254
00255
00256
00257 for (i=0; i< nslits; i++) {
00258 array2D_set_value(slit_edges, vec_x[i], i, 0);
00259 array2D_set_value(slit_edges, vec_y[i], i, 1);
00260
00261 }
00262 fclose(slit_list);
00263
00264
00265 cube = makeCubeSpi( clean_im, slit_edges, correct_dist );
00266 if (cube == NULL) {
00267 cpl_msg_error(_id," could not construct data cube!\n");
00268 return -1;
00269 }
00270 destroy_2Darray(slit_edges, cfg->nslits);
00271 } else {
00272 distances = new_float_array(cfg->nslits-1) ;
00273
00274
00275
00276 if ( NULL == (dist_list = fopen (cfg->poslist, "r" ) ) )
00277 {
00278 cpl_msg_error(_id,"cannot open %s\n", cfg->poslist) ;
00279 return -1 ;
00280 }
00281
00282
00283 tbl_distances_name=cfg->poslist;
00284
00285 tbl_distances = cpl_table_load(tbl_distances_name,1,0);
00286 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00287 cpl_msg_error(_id,"error loadting tbl %s",tbl_distances_name);
00288 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00289 return -1;
00290 }
00291 for (i =0 ; i< cfg->nslits-1; i++){
00292 tmp_float=cpl_table_get_float(tbl_distances,
00293 "slitlet_distance",i,status);
00294 array_set_value(distances,tmp_float,i);
00295
00296 }
00297 cpl_table_delete(tbl_distances);
00298 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00299 cpl_msg_error(_id,"error writing tbl %s",tbl_distances_name);
00300 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00301 return -1;
00302 }
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339 tbl_firstcol_name=cfg->firstCol;
00340 if ((tbl_firstcol=cpl_table_load(tbl_firstcol_name,1,0)) != NULL) {
00341 fc=cpl_table_get_float(tbl_firstcol,"first_col",0,status);
00342 cpl_table_delete(tbl_firstcol);
00343 } else {
00344 cpl_msg_error(_id,"table %s not found! Exit!",tbl_firstcol_name);
00345 return -1 ;
00346 }
00347
00348
00349 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00350 cpl_msg_error(_id,"error loading tbl %s",tbl_firstcol_name);
00351 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00352 return -1;
00353 }
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367 cube = makeCubeDist( clean_im, fc, distances, correct_dist );
00368 if (cube == NULL) {
00369 cpl_msg_error(_id," could not construct a data cube!\n");
00370 return -1;
00371 }
00372 destroy_array (distances);
00373
00374 }
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385 if (strcmp(cfg->method,"P") == 0) {
00386 out_cube = fineTuneCube( cube, correct_dist, cfg->order );
00387 if (out_cube == NULL) {
00388 cpl_msg_error(_id," could not fine tune the data cube\n");
00389 return -1;
00390 }
00391 } else if (strcmp(cfg->method, "F") == 0) {
00392 for (i=0; i< cfg->nslits; i++) {
00393 array_set_value(correct_dist, -array_get_value(correct_dist,i),
00394 i);
00395 }
00396 out_cube = fineTuneCubeByFFT( cube, correct_dist );
00397 if (out_cube == NULL) {
00398 cpl_msg_error(_id," could not fine tune the data cube\n");
00399 return -1;
00400 }
00401 } else if ( strcmp(cfg->method, "S") == 0) {
00402 out_cube = fineTuneCubeBySpline( cube, correct_dist );
00403 if (out_cube == NULL) {
00404 cpl_msg_error(_id," could not fine tune the data cube\n");
00405 return -1;
00406 }
00407 } else {
00408 cpl_msg_error(_id," wrong method indicator given!");
00409 return -1;
00410 }
00411
00412 out_image = medianCube(out_cube);
00413 if (out_image == NULL) {
00414 cpl_msg_error(_id," could not do medianCube()\n");
00415 return -1;
00416 }
00417
00418 if (strstr(cfg->outName, "." ) != NULL ) {
00419 sprintf(out_Name_No, "%s%d%s",
00420 get_rootname(cfg->outName), fra,strstr(cfg->outName,"."));
00421 } else {
00422 sprintf(out_Name_No, "%s%d", cfg->outName, fra);
00423 }
00424
00425 o_img=cpl_image_wrap_float(out_image->lx,out_image->ly,out_image->data);
00426
00427
00428
00429 if(-1 == sinfoni_pro_dump_ima(o_img,stk,sof,out_Name_No,
00430 PRO_FOCUS,_id,config)) {
00431 cpl_msg_error(_id,"cannot dump ima %s", out_Name_No);
00432 destroy_array(correct_dist);
00433 destroy_image(clean_im);
00434 destroy_cube(cube);
00435 destroy_cube (out_cube);
00436 destroy_image (out_image);
00437 fits_header_destroy(head);
00438 cpl_free(col_name);
00439 fclose(fitpar_file);
00440 cpl_free(out_Name_No);
00441 cpl_free(gaussplotNameNo);
00442 free_focus(cfg);
00443 return -1;
00444 }
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454 fit_par = new_double_array(npar);
00455 derv_par = new_double_array(npar);
00456 mpar = new_int_array(npar);
00457
00458
00459 intarray_set_value(mpar, cfg->mpar0, 0);
00460 intarray_set_value(mpar, cfg->mpar1, 1);
00461 intarray_set_value(mpar, cfg->mpar2, 2);
00462 intarray_set_value(mpar, cfg->mpar3, 3);
00463 intarray_set_value(mpar, cfg->mpar4, 4);
00464 intarray_set_value(mpar, cfg->mpar5, 5);
00465 intarray_set_value(mpar, cfg->mpar6, 6);
00466
00467 iters = fit2dGaussian( out_image, fit_par, derv_par,
00468 mpar, cfg->llx, cfg->lly,
00469 cfg->halfbox_x, cfg->halfbox_y,&check );
00470 if (iters < 0) {
00471 cpl_msg_error(_id," could not fit 2d-Gaussian\n");
00472 return -1;
00473 }
00474
00475
00476
00477 string_par = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00478 tmp_string_par = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00479
00480 sprintf(string_par,"%s","");
00481 for (i=0; i< npar; i++) {
00482 sprintf(col_name,"%s%d","par",i);
00483 par = doublearray_get_value(fit_par, i);
00484 sprintf(tmp_string_par,"%14.10f%s",par," ");
00485
00486 strcat(string_par,tmp_string_par);
00487 cpl_table_set_double(tbl_focus,col_name,fra,par);
00488 }
00489 fprintf(fitpar_file, "%s \n", string_par);
00490 cpl_free(string_par);
00491 cpl_free(tmp_string_par);
00492
00493 if (cfg->plotGaussInd == 1) {
00494 img = plotGaussian( out_image, fit_par );
00495 if (img == NULL) {
00496 cpl_msg_error(_id," could not do plotGaussian()\n");
00497 return -1;
00498 }
00499
00500 if (strstr(cfg->gaussplotName, "." ) != NULL ) {
00501 sprintf(gaussplotNameNo, "%s%d%s",
00502 get_rootname(cfg->gaussplotName), fra,
00503 strstr(cfg->gaussplotName,"."));
00504
00505 } else {
00506 sprintf(gaussplotNameNo, "%s%d", cfg->gaussplotName, fra);
00507 }
00508
00509 g_img=cpl_image_wrap_float(img->lx,img->ly,img->data);
00510
00511
00512 if(-1 == sinfoni_pro_dump_ima(g_img,stk,sof,gaussplotNameNo,
00513 PRO_FOCUS_GAUSS,_id,config)) {
00514 cpl_msg_error(_id,"cannot dump ima %s", gaussplotNameNo);
00515
00516 destroy_intarray(mpar);
00517 destroy_doublearray(fit_par);
00518 destroy_doublearray(derv_par);
00519 destroy_image(img);
00520
00521 destroy_array(correct_dist);
00522 destroy_image(clean_im);
00523 destroy_cube(cube);
00524 destroy_cube (out_cube);
00525 destroy_image (out_image);
00526 fits_header_destroy(head);
00527 cpl_free(col_name);
00528 fclose(fitpar_file);
00529 cpl_free(out_Name_No);
00530 cpl_free(gaussplotNameNo);
00531 free_focus(cfg);
00532 return -1;
00533 }
00534
00535
00536 destroy_image(img);
00537
00538 }
00539
00540
00541 destroy_intarray(mpar);
00542 destroy_doublearray(fit_par);
00543 destroy_doublearray(derv_par);
00544 destroy_array(correct_dist);
00545 destroy_image(clean_im);
00546 destroy_cube(cube);
00547 destroy_cube (out_cube);
00548 destroy_image (out_image);
00549 fits_header_destroy(head);
00550
00551 }
00552 tbl_name=cfg->fitlist;
00553 cpl_table_save(tbl_focus,tbl_focus_phead,tbl_focus_head,tbl_name,0);
00554 cpl_table_delete(tbl_focus);
00555 cpl_free(col_name);
00556 fclose(fitpar_file);
00557 cpl_free(out_Name_No);
00558 cpl_free(gaussplotNameNo);
00559 free_focus(cfg);
00560
00561
00562 return 0 ;
00563 }
00564