00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151 #define POSIX_SOURCE 1
00152 #include "vltPort.h"
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 #include "spectrum_ops.h"
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175 Vector * newVector( ulong32 n_elements )
00176 {
00177 Vector * new_vector ;
00178
00179 if ( n_elements <= 0 )
00180 {
00181 cpl_msg_error ("newVector:"," wrong number of elements\n") ;
00182 return NullVector ;
00183 }
00184
00185
00186 new_vector = (Vector *) cpl_malloc (sizeof (Vector)) ;
00187 new_vector -> n_elements = n_elements ;
00188 new_vector -> data = (pixelvalue *) cpl_calloc (n_elements, sizeof (pixelvalue)) ;
00189
00190 return new_vector ;
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 void destroyVector( Vector *vector )
00202 {
00203 if ( vector == NULL )
00204 {
00205 cpl_msg_error("destroyVector:"," NULL Vector given!\n") ;
00206 return ;
00207 }
00208
00209 cpl_free ( vector -> data ) ;
00210 cpl_free ( vector ) ;
00211 }
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221 OneImage * vectorToImage( Vector * spectrum )
00222 {
00223 OneImage * returnIm ;
00224 int i ;
00225
00226 if ( spectrum == NULL )
00227 {
00228 cpl_msg_error("vectorToImage:"," no spectrum given!\n") ;
00229 return NullImage ;
00230 }
00231
00232
00233 if ( NullImage == (returnIm = new_image(1, spectrum->n_elements)) )
00234 {
00235 cpl_msg_error("vectorToImage:"," no spectrum given!\n") ;
00236 destroyVector(spectrum) ;
00237 return NullImage ;
00238 }
00239
00240 for ( i = 0 ; i < spectrum->n_elements ; i++ )
00241 {
00242 returnIm -> data[i] = spectrum -> data[i] ;
00243 }
00244
00245 destroyVector (spectrum) ;
00246 return returnIm ;
00247 }
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257 Vector * imageToVector( OneImage * spectrum )
00258 {
00259 Vector * returnVector ;
00260 int i ;
00261
00262 if ( spectrum == NULL )
00263 {
00264 cpl_msg_error("imageToVector:"," no spectrum given!\n") ;
00265 return NULL ;
00266 }
00267
00268
00269 if ( NULL == (returnVector = newVector(spectrum->nbpix)) )
00270 {
00271 cpl_msg_error("imageToVector:"," cannot allocate memory!\n") ;
00272 destroy_image(spectrum) ;
00273 return NULL ;
00274 }
00275
00276 for ( i = 0 ; i < (int) spectrum->nbpix ; i++ )
00277 {
00278 returnVector -> data[i] = spectrum -> data[i] ;
00279 }
00280
00281 destroy_image (spectrum) ;
00282 return returnVector ;
00283 }
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295 OneImage * extractSpectrumFromResampledFlat( OneImage * resflat,
00296 float loreject,
00297 float hireject )
00298 {
00299 OneImage * retIm ;
00300 int col, row ;
00301 int n ;
00302 float array[resflat->ly] ;
00303 float cleanMean ;
00304 Vector * spectrum ;
00305
00306 if ( resflat == NullImage )
00307 {
00308 cpl_msg_error("extractSpectrumFromResampledFlat:"," no flatfield given!\n") ;
00309 return NullImage ;
00310 }
00311
00312
00313 if ( NullVector == (spectrum = newVector(resflat->ly) ) )
00314 {
00315 cpl_msg_error("extractSpectrumFromResampledFlat:"," could not allocate memory!\n") ;
00316 return NullImage ;
00317 }
00318
00319
00320 for ( row = 0 ; row < resflat -> ly ; row++ )
00321 {
00322 n = 0 ;
00323 for ( col = 0 ; col < resflat -> lx ; col++ )
00324 {
00325 if ( !qfits_isnan(resflat -> data[col + row*resflat->lx]) )
00326 {
00327 array[n] = resflat->data[col+row*resflat->lx] ;
00328 n++ ;
00329 }
00330 }
00331 if ( n == 0 )
00332 {
00333 cpl_msg_warning("extractSpectrumFromResampledFlat:"," only bad pixels in row: %d!\n", row) ;
00334 cleanMean = ZERO ;
00335 }
00336 else
00337 {
00338 if ( FLT_MAX == (cleanMean = clean_mean(array, n, loreject, hireject)) )
00339 {
00340 cpl_msg_error("extractSpectrumFromResampledFlat:"," could not do clean_mean!\n") ;
00341 destroyVector(spectrum) ;
00342 return NullImage ;
00343 }
00344 }
00345 spectrum->data[row] = cleanMean ;
00346 }
00347 if ( NullImage == ( retIm = vectorToImage( spectrum ) ) )
00348 {
00349 cpl_msg_error("extractSpectrumFromResampledFlat:"," could not do vectorToImage!\n") ;
00350 destroyVector(spectrum) ;
00351 return NullImage ;
00352 }
00353 return retIm ;
00354 }
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366 OneImage * multiplyImageWithSpectrum( OneImage * image, OneImage * spectrum )
00367 {
00368 int col, row ;
00369 OneImage * retImage ;
00370
00371 if ( image == NullImage )
00372 {
00373 cpl_msg_error("multiplyImageWithSpectrum:"," no image given!\n") ;
00374 return NullImage ;
00375 }
00376
00377 if ( spectrum == NullImage )
00378 {
00379 cpl_msg_error("multiplyImageWithSpectrum:"," no spectrum image given!\n") ;
00380 return NullImage ;
00381 }
00382 if ( spectrum -> ly != image -> ly )
00383 {
00384 cpl_msg_error("multiplyImageWithSpectrum:"," images are not compatible in pixel length!\n") ;
00385 return NullImage ;
00386 }
00387
00388 if ( NullImage == (retImage = copy_image(image)) )
00389 {
00390 cpl_msg_error("multiplyImageWithSpectrum:"," could not copy original image!\n") ;
00391 return NullImage ;
00392 }
00393
00394 for ( col = 0 ; col < image -> lx ; col++ )
00395 {
00396 for ( row = 0 ; row < image -> ly ; row++ )
00397 {
00398 if ( !qfits_isnan(image->data[col+row*image->lx]) &&
00399 !qfits_isnan(spectrum->data[col+row*image->lx]))
00400 {
00401 retImage->data[col+row*image->lx] =
00402 image->data[col+row*image->lx] * spectrum->data[row] ;
00403
00404 }
00405 }
00406 }
00407 return retImage ;
00408 }
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441 OneImage * optimalExtractionFromCube( OneCube * cube,
00442 int llx,
00443 int lly,
00444 int halfbox_x,
00445 int halfbox_y,
00446 float fwhm_factor,
00447 float backvariance,
00448 float sky,
00449 float gain,
00450 float exptime,
00451 const char* name,
00452 cpl_table** spectrum,
00453 int qc_info,
00454 int* check2)
00455 {
00456 int col, row, z ;
00457 OneImage * averagedIm ;
00458 OneImage * retIm ;
00459 double fit_par[7] ;
00460 double derv_par[7] ;
00461 int mpar[7] ;
00462 double gfit_par[7] ;
00463 double gderv_par[7] ;
00464 int gmpar[7] ;
00465 int fitInd ;
00466 int i ;
00467 double sum ;
00468 double weight[cube->lx][cube->ly] ;
00469 double psf[cube->lx][cube->ly] ;
00470
00471 double variance ;
00472 double xdat[2] ;
00473 float weighted_sum ;
00474 float counts_tot ;
00475 float counts_bkg ;
00476 float bkg_tot ;
00477
00478
00479 int first_col, last_col ;
00480 int first_row, last_row ;
00481 float norm ;
00482 float sum_psf=0;
00483 float sum_wgt=0;
00484 float cenpix = 0;
00485 float cenLambda = 0;
00486 float dispersion = 0;
00487 float lambda=0;
00488 float lambda_start=0;
00489
00490
00491 const char* _id="optimalExtractionFromCube";
00492
00493 if ( NullCube == cube )
00494 {
00495 cpl_msg_error("optimalExtractionFromCube:"," no cube given!\n") ;
00496 return NullImage ;
00497 }
00498 if ( llx < 0 || llx + 2*halfbox_x >= cube->lx || lly < 0 || lly+2*halfbox_y >= cube->ly )
00499 {
00500 cpl_msg_error("optimalExtractionFromCube:"," lower left edge points wrong position!\n") ;
00501 return NullImage ;
00502 }
00503 if ( halfbox_x <= 0 || halfbox_y <= 0 || 2*halfbox_x > cube->lx || 2*halfbox_y > cube->ly )
00504 {
00505 cpl_msg_error("optimalExtractionFromCube:"," wrong halfbox width given!\n") ;
00506 return NullImage ;
00507 }
00508 if ( fwhm_factor <= 0. )
00509 {
00510 cpl_msg_error("optimalExtractionFromCube:"," wrong fwhm_factor given!\n") ;
00511 return NullImage ;
00512 }
00513 if ( backvariance < 0. )
00514 {
00515 cpl_msg_error("optimalExtractionFromCube:"," wrong backvariance given!\n") ;
00516 return NullImage ;
00517 }
00518 if ( exptime <= 0. || exptime == FLAG )
00519 {
00520 cpl_msg_error("optimalExtractionFromCube:"," wrong exposure time given!\n") ;
00521 return NullImage ;
00522 }
00523
00524
00525 if ( NULL == (retIm = new_image(1, cube->np)) )
00526 {
00527 cpl_msg_error("optimalExtractionFromCube:"," memory allocation failed!\n") ;
00528 return NullImage ;
00529 }
00530
00531 if ( NullImage == (averagedIm = averageCubeToImage(cube)) )
00532 {
00533 cpl_msg_error("optimalExtractionFromCube:"," averageCubeToImage failed!\n") ;
00534 destroy_image(retIm) ;
00535 return NullImage ;
00536 }
00537
00538
00539 for ( i = 0 ; i < 7 ; i++ )
00540 {
00541 mpar[i] = 1 ;
00542 }
00543
00544 if ( -1 == (fitInd = fit2dGaussian(averagedIm, fit_par, derv_par, mpar, llx, lly, halfbox_x, halfbox_y, check2 )) )
00545 {
00546 cpl_msg_error("optimalExtractionFromCube:"," fit2dGaussian failed!\n") ;
00547 destroy_image(retIm) ;
00548 destroy_image(averagedIm) ;
00549 return NullImage ;
00550 }
00551
00552
00553 sum = 0 ;
00554 for ( row = 0 ; row < cube->ly ; row++ )
00555 {
00556 for ( col = 0 ; col < cube->lx ; col++ )
00557 {
00558 xdat[0] = (double) col ;
00559 xdat[1] = (double) row ;
00560 psf[col][row] = gaussianEllipse(xdat , fit_par) - fit_par[3] ;
00561 sum += psf[col][row] ;
00562 }
00563 }
00564
00565 norm = 0. ;
00566 variance = 0. ;
00567 sum_psf=0;
00568 for ( row = 0 ; row < cube->ly ; row++ )
00569 {
00570 for ( col = 0 ; col < cube->lx ; col++ )
00571 {
00572 psf[col][row] /= sum ;
00573 sum_psf += psf[col][row];
00574 if ( !qfits_isnan(averagedIm->data[col+row*cube->lx]) )
00575 {
00576
00577
00578
00579 variance = averagedIm -> data[col+row*cube->lx] / gain ;
00580
00581 }
00582 else
00583 {
00584 weight[col][row] = 0. ;
00585 }
00586 if (variance == 0.)
00587 {
00588 weight[col][row] = 0. ;
00589 }
00590 else
00591 {
00592
00593 weight[col][row] = psf[col][row]/variance ;
00594
00595 norm += weight[col][row] * weight[col][row] * variance ;
00596
00597 }
00598
00599 }
00600 }
00601
00602 sum_wgt=0;
00603 for ( row = 0 ; row < cube->ly ; row++ )
00604 {
00605 for ( col = 0 ; col < cube->lx ; col++ )
00606 {
00607 weight[col][row] /= norm;
00608 sum_wgt += weight[col][row]*psf[col][row];
00609 }
00610 }
00611 cpl_msg_info(_id,"sum_psf=%f sum_wgt=%f norm=%f",sum_psf,sum_wgt,norm);
00612 destroy_image(averagedIm) ;
00613 if ( norm == 0. )
00614 {
00615 cpl_msg_error("optimalExtractionFromCube:"," normalization sum is zero\n") ;
00616 destroy_image(retIm) ;
00617 return NullImage ;
00618 }
00619
00620
00621 first_col = (int) (fit_par[0] - fwhm_factor*fit_par[4]*cos(fit_par[6])) ;
00622 first_col = (first_col > 2 ) ? first_col : 2;
00623 last_col = (int) (fit_par[0] + fwhm_factor*fit_par[4]*cos(fit_par[6])) ;
00624 last_col = (last_col < 63 ) ? last_col : 63;
00625
00626 first_row = (int) (fit_par[1] - fwhm_factor*fit_par[5]*cos(fit_par[6])) ;
00627 first_row = (first_row > 2 ) ? first_row : 2;
00628 last_row = (int) (fit_par[1] + fwhm_factor*fit_par[5]*cos(fit_par[6])) ;
00629 last_row = (last_row < 63 ) ? last_row : 63;
00630
00631 cpl_msg_info(_id,"fit_par: %f %f %f %f %f %f %f", fit_par[0],fit_par[1],fit_par[2],fit_par[3],
00632 fit_par[4],fit_par[5],fit_par[6]);
00633
00634 if ( first_col < 0 || first_row < 0 || last_col >= cube->lx || last_row >= cube->ly )
00635 {
00636 cpl_msg_error("optimalExtractionFromCube:"," star badly centered in FOV!\n") ;
00637 destroy_image(retIm) ;
00638 return NullImage ;
00639 }
00640
00641
00642 cpl_table_new_column(*spectrum,"wavelength", CPL_TYPE_FLOAT);
00643
00644 cpl_table_new_column(*spectrum,"counts_tot" , CPL_TYPE_FLOAT);
00645 cpl_table_new_column(*spectrum,"counts_bkg" , CPL_TYPE_FLOAT);
00646 cpl_table_new_column(*spectrum,"bkg_tot" , CPL_TYPE_FLOAT);
00647
00648 if(qc_info==1) {
00649 cpl_table_new_column(*spectrum,"AMP" , CPL_TYPE_FLOAT);
00650 cpl_table_new_column(*spectrum,"XC" , CPL_TYPE_FLOAT);
00651 cpl_table_new_column(*spectrum,"YC" , CPL_TYPE_FLOAT);
00652 cpl_table_new_column(*spectrum,"BKG" , CPL_TYPE_FLOAT);
00653 cpl_table_new_column(*spectrum,"FWHMX" , CPL_TYPE_FLOAT);
00654 cpl_table_new_column(*spectrum,"FWHMY" , CPL_TYPE_FLOAT);
00655 cpl_table_new_column(*spectrum,"ANGLE" , CPL_TYPE_FLOAT);
00656 }
00657
00658 cenpix = spiffi_get_cenpix((char*)name);
00659 cenLambda = spiffi_get_cenLambda((char*)name);
00660 dispersion = spiffi_get_dispersion((char*)name);
00661 lambda_start=cenLambda-cenpix*dispersion;
00662
00663 cpl_msg_info(_id,"frow %d lrow %d fcol %d lcol %d",first_row, last_row, first_col, last_col);
00664
00665 for ( z = 0 ; z < cube->np ; z++ )
00666 {
00667 weighted_sum = 0. ;
00668 counts_tot=0.;
00669 counts_bkg=0.;
00670
00671 bkg_tot=0.;
00672
00673 if(qc_info==1) {
00674 fit2dGaussian(cube->plane[z],gfit_par,gderv_par,gmpar,llx,lly,halfbox_x,halfbox_y,check2);
00675 }
00676
00677 for ( row = first_row ; row <= last_row ; row++ )
00678 {
00679 for ( col = first_col ; col < last_col ; col++ )
00680 {
00681 if ( !qfits_isnan(cube->plane[z]->data[col+row*cube->lx]) )
00682 {
00683
00684 weighted_sum += weight[col][row] * (cube->plane[z]->data[col+row*cube->lx] - fit_par[3]);
00685
00686 counts_bkg += (cube->plane[z]->data[col+row*cube->lx] - fit_par[3]);
00687 counts_tot += (cube->plane[z]->data[col+row*cube->lx]);
00688 bkg_tot += fit_par[3];
00689
00690 }
00691 }
00692 }
00693
00694 if (weighted_sum == 0.)
00695 {
00696 weighted_sum = ZERO ;
00697 counts_tot = ZERO;
00698 counts_bkg = ZERO;
00699 bkg_tot = ZERO;
00700
00701 }
00702 else
00703 {
00704
00705
00706
00707
00708
00709 }
00710
00711 retIm->data[z] = weighted_sum ;
00712 lambda=lambda_start+z*dispersion;
00713 cpl_table_set_float(*spectrum,"wavelength" ,z,lambda);
00714
00715 cpl_table_set_float(*spectrum,"counts_tot" ,z,counts_tot);
00716 cpl_table_set_float(*spectrum,"counts_bkg" ,z,counts_bkg);
00717 cpl_table_set_float(*spectrum,"bkg_tot" ,z,bkg_tot);
00718
00719 if(qc_info==1) {
00720 cpl_table_set_float(*spectrum,"AMP" ,z,gfit_par[0]);
00721 cpl_table_set_float(*spectrum,"XC" ,z,gfit_par[1]);
00722 cpl_table_set_float(*spectrum,"YC" ,z,gfit_par[2]);
00723 cpl_table_set_float(*spectrum,"BKG" ,z,gfit_par[3]);
00724 cpl_table_set_float(*spectrum,"FWHMX" ,z,gfit_par[4]);
00725 cpl_table_set_float(*spectrum,"FWHMY" ,z,gfit_par[5]);
00726 cpl_table_set_float(*spectrum,"ANGLE" ,z,gfit_par[6]);
00727 }
00728
00729 }
00730
00731 return retIm ;
00732 }
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754 Vector * extractSkyFromCube( OneCube * cube,
00755 float loReject,
00756 float hiReject,
00757 int * position,
00758 int tolerance,
00759 int posindicator )
00760 {
00761 Vector * spectrum ;
00762 int x, y, z ;
00763 int n ;
00764 int n_sky ;
00765 int x_low , x_high ;
00766 int y_low , y_high ;
00767 int hi_x, lo_x ;
00768 float * to_average ;
00769 float cleanMean ;
00770
00771 if ( NullCube == cube )
00772 {
00773 cpl_msg_error("extractSkyFromCube:"," no cube given!\n") ;
00774 return NullVector ;
00775 }
00776 if ( loReject < 0. || hiReject < 0. || loReject + hiReject >= 90. )
00777 {
00778 cpl_msg_error("extractSkyFromCube:"," wrong or unrealistic loReject and hiReject values!\n") ;
00779 return NullVector ;
00780 }
00781 if ( position == NULL)
00782 {
00783 cpl_msg_error("extractSkyFromCube:"," no position array given!\n") ;
00784 return NullVector ;
00785 }
00786 if ( position[0] < 0 || position[1] < 0 || position[0] > cube->lx || position[1] > cube->ly )
00787 {
00788 cpl_msg_error("extractSkyFromCube:"," wrong position of sky spider!\n") ;
00789 return NullVector ;
00790 }
00791 if ( tolerance < 0 || tolerance >= cube->lx )
00792 {
00793 cpl_msg_error("extractSkyFromCube:"," wrong tolerance given!\n") ;
00794 return NullVector ;
00795 }
00796 if ( posindicator == 0 )
00797 {
00798 cpl_msg_error("extractSkyFromCube:"," no edge indicator given!\n") ;
00799 return NullVector ;
00800 }
00801
00802
00803 switch(posindicator)
00804 {
00805
00806 case 1:
00807 x_low = position[0] + tolerance ;
00808 x_high = cube->lx ;
00809 y_low = 0 ;
00810 y_high = position[1] - tolerance ;
00811 break ;
00812
00813 case 2:
00814 x_low = position[0] + tolerance ;
00815 x_high = cube->lx ;
00816 y_low = position[1] + tolerance ;
00817 y_high = cube->ly ;
00818 break ;
00819
00820 case 3:
00821 x_low = 0 ;
00822 x_high = position[0] - tolerance ;
00823 y_low = position [1] + tolerance ;
00824 y_high = cube->ly ;
00825 break ;
00826 default:
00827 cpl_msg_error("extractSkyFromCube:"," wrong position indicator index!\n") ;
00828 return NullVector ;
00829 break ;
00830 }
00831 if ( x_low >= cube -> lx || x_high < 1 || y_low >= cube->ly || y_high < 1 )
00832 {
00833 cpl_msg_error("extractSkyFromCube:"," tolerance too high!\n") ;
00834 return NullVector ;
00835 }
00836 if ( x_high - x_low != y_high - y_low )
00837 {
00838 cpl_msg_error("extractSkyFromCube:"," sky edge is not a diagonal line!\n") ;
00839 return NullVector ;
00840 }
00841
00842
00843
00844 n_sky = (x_high - x_low) * (x_high - x_low - 1) / 2 ;
00845 if ( n_sky <= 0 )
00846 {
00847 cpl_msg_error("extractSkyFromCube:"," no sky spectrum in found in cube!\n") ;
00848 return NullVector ;
00849 }
00850 if ( n_sky == 1 )
00851 {
00852 cpl_msg_warning("extractSkyFromCube:"," only one sky spectrum is taken, no averaging!\n") ;
00853 }
00854
00855
00856 if ( NullVector == (spectrum = newVector(cube->np)) )
00857 {
00858 cpl_msg_error("extractSkyFromCube:"," could not allocate memory!\n") ;
00859 return NullVector ;
00860 }
00861
00862
00863 for ( z = 0 ; z < cube->np ; z++ )
00864 {
00865
00866 if ( NULL == ( to_average = (float*) cpl_calloc(n_sky, sizeof (float)) ) )
00867 {
00868 cpl_msg_error("extractSkyFromCube:"," could not allocate memory!\n") ;
00869 destroyVector(spectrum) ;
00870 return NullVector ;
00871 }
00872 n = 0 ;
00873 switch(posindicator)
00874 {
00875
00876 case 1:
00877 lo_x = x_low ;
00878 for ( y = y_low ; y < y_high - 1 ; y++ )
00879 {
00880 lo_x++ ;
00881 for ( x = lo_x ; x < x_high ; x++ )
00882 {
00883 to_average[n] = cube->plane[z]->data[x+y*cube->lx] ;
00884 n++ ;
00885 }
00886 }
00887 break ;
00888
00889 case 2:
00890 hi_x = x_high ;
00891 for ( y = y_low ; y < y_high - 1 ; y++ )
00892 {
00893 hi_x-- ;
00894 for ( x = x_low ; x < hi_x ; x++ )
00895 {
00896 to_average[n] = cube->plane[z]->data[x+y*cube->lx] ;
00897 n++ ;
00898 }
00899 }
00900 break ;
00901
00902 case 3:
00903 lo_x = x_high ;
00904 for ( y = y_low+1 ; y < y_high ; y++ )
00905 {
00906 lo_x-- ;
00907 for ( x = lo_x ; x < x_high ; x++ )
00908 {
00909 to_average[n] = cube->plane[z]->data[x+y*cube->lx] ;
00910 n++ ;
00911 }
00912 }
00913 break ;
00914
00915 case 4:
00916 hi_x = x_low ;
00917 for ( y = y_low+1 ; y < y_high ; y++ )
00918 {
00919 hi_x++ ;
00920 for ( x = x_low ; x < hi_x ; x++ )
00921 {
00922 to_average[n] = cube->plane[z]->data[x+y*cube->lx] ;
00923 n++ ;
00924 }
00925 }
00926 break ;
00927 default:
00928 cpl_msg_error("extractSkyFromCube:"," wrong position indicator index!\n") ;
00929 return NullVector ;
00930 break ;
00931 }
00932 if ( n != n_sky )
00933 {
00934 cpl_msg_warning("extractSkyFromCube:","number of stored sky image pixels does not equal number of computed sky pixels!\n") ;
00935 }
00936
00937
00938 cleanMean = clean_mean (to_average, n, loReject, hiReject) ;
00939 if (cleanMean == FLT_MAX)
00940 {
00941 cpl_msg_error("extractSkyFromCube:"," could not take a clean mean!\n") ;
00942 destroyVector(spectrum) ;
00943 cpl_free(to_average) ;
00944 return NullVector ;
00945 }
00946 spectrum->data[z] = cleanMean ;
00947 cpl_free (to_average) ;
00948 }
00949
00950 return spectrum ;
00951 }
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964 Vector * sumRectangleOfCubeSpectra( OneCube * cube,
00965 int llx,
00966 int lly,
00967 int urx,
00968 int ury )
00969 {
00970 Vector * sum ;
00971 pixelvalue *rectangle ;
00972 int i, j, k, l, m ;
00973 int recsize ;
00974
00975 if ( cube == NULL || cube -> np < 1 )
00976 {
00977 cpl_msg_error ("sumRectangleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
00978 return NullVector ;
00979 }
00980
00981 if ((llx<0) || (llx>=cube->lx) ||
00982 (urx<0) || (urx>=cube->lx) ||
00983 (lly<0) || (lly>=cube->ly) ||
00984 (ury<0) || (ury>=cube->ly) ||
00985 (llx>=urx) || (lly>=ury))
00986 {
00987 cpl_msg_error("sumRectangleOfCubeSpectra:"," invalid rectangle coordinates:") ;
00988 cpl_msg_error("sumRectangleOfCubeSpectra:","lower left is [%d %d] upper right is [%d %d]", llx, lly, urx, ury) ;
00989 return NullVector ;
00990 }
00991
00992 recsize = (urx - llx + 1) * (ury - lly + 1) ;
00993
00994
00995 if (NULL == (sum = newVector (cube -> np)) )
00996 {
00997 cpl_msg_error ("sumRectangleOfCubeSpectra:"," cannot allocate a new vector \n") ;
00998 return NullVector ;
00999 }
01000
01001
01002
01003
01004
01005 for ( i = 0 ; i < cube->np ; i++ )
01006 {
01007 m = 0 ;
01008 rectangle = (pixelvalue *) cpl_calloc (recsize, sizeof (pixelvalue*));
01009
01010 for ( j = lly ; j <= ury ; j++ )
01011 {
01012 for ( k = llx ; k <= urx ; k++ )
01013 {
01014 rectangle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01015 m ++ ;
01016 }
01017 }
01018 for ( l = 0 ; l < recsize ; l++ )
01019 {
01020 if ( qfits_isnan(rectangle[l]) )
01021 {
01022 continue ;
01023 }
01024 sum -> data[i] += rectangle[l] ;
01025 }
01026 cpl_free ( rectangle ) ;
01027 }
01028 return sum ;
01029 }
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042 Vector * sumCircleOfCubeSpectra( OneCube * cube,
01043 int centerx,
01044 int centery,
01045 int radius )
01046 {
01047 Vector * sum ;
01048 pixelvalue * circle ;
01049 int i, j, k, l, m, n ;
01050 int circsize ;
01051
01052 if ( cube == NULL || cube -> np < 1 )
01053 {
01054 cpl_msg_error ("sumCircleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01055 return NullVector ;
01056 }
01057
01058 if ((centerx+radius>=cube->lx) ||
01059 (centery+radius>=cube->ly) ||
01060 (centerx-radius<0) ||
01061 (centery-radius<0))
01062 {
01063 cpl_msg_error("sumCircleOfCubeSpectra:"," invalid circular coordinates") ;
01064 return NullVector ;
01065 }
01066
01067 n = 0 ;
01068 for ( j = centery - radius ; j <= centery + radius ; j++ )
01069 {
01070 for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01071 {
01072 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01073 {
01074 n ++ ;
01075 }
01076 }
01077 }
01078 if (n == 0)
01079 {
01080 cpl_msg_error ("sumCircleOfCubeSpectra:"," no data points found!\n") ;
01081 return NullVector ;
01082 }
01083 circsize = n ;
01084
01085
01086 if (NULL == (sum = newVector (cube -> np)) )
01087 {
01088 cpl_msg_error ("sumCircleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01089 return NullVector ;
01090 }
01091
01092
01093
01094
01095
01096 for ( i = 0 ; i < cube->np ; i++ )
01097 {
01098 m = 0 ;
01099 circle = (pixelvalue *) cpl_calloc (circsize, sizeof (pixelvalue*));
01100
01101 for ( j = centery - radius ; j <= centery + radius ; j++ )
01102 {
01103 for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01104 {
01105 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01106 {
01107 circle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01108 m ++ ;
01109 }
01110 }
01111 }
01112
01113 for ( l = 0 ; l < circsize ; l++ )
01114 {
01115 if ( qfits_isnan(circle[l]) )
01116 {
01117 continue ;
01118 }
01119 sum -> data[i] += circle[l] ;
01120 }
01121 cpl_free (circle) ;
01122 }
01123 return sum ;
01124 }
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137 Vector * meanRectangleOfCubeSpectra( OneCube * cube,
01138 int llx,
01139 int lly,
01140 int urx,
01141 int ury )
01142 {
01143 Vector * mean ;
01144 pixelvalue *rectangle ;
01145 int i, j, k, l, m ;
01146 int recsize, nv ;
01147
01148 if ( cube == NULL || cube -> np < 1 )
01149 {
01150 cpl_msg_error ("meanRectangleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01151 return NullVector ;
01152 }
01153
01154 if ((llx<0) || (llx>=cube->lx) ||
01155 (urx<0) || (urx>=cube->lx) ||
01156 (lly<0) || (lly>=cube->ly) ||
01157 (ury<0) || (ury>=cube->ly) ||
01158 (llx>=urx) || (lly>=ury))
01159 {
01160 cpl_msg_error("meanRectangleOfCubeSpectra:"," invalid rectangle coordinates:") ;
01161 cpl_msg_error("meanRectangleOfCubeSpectra:","lower left is [%d %d] upper right is [%d %d]",
01162 llx, lly, urx, ury) ;
01163 return NullVector ;
01164 }
01165
01166 recsize = (urx - llx + 1) * (ury - lly + 1) ;
01167
01168
01169 if (NULL == (mean = newVector (cube -> np)) )
01170 {
01171 cpl_msg_error ("meanRectangleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01172 return NullVector ;
01173 }
01174
01175
01176
01177
01178
01179 for ( i = 0 ; i < cube->np ; i++ )
01180 {
01181 m = 0 ;
01182 rectangle = (pixelvalue *) cpl_calloc (recsize, sizeof (pixelvalue*));
01183
01184 for ( j = lly ; j <= ury ; j++ )
01185 {
01186 for ( k = llx ; k <= urx ; k++ )
01187 {
01188 rectangle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01189 m ++ ;
01190 }
01191 }
01192 nv = 0 ;
01193 for ( l = 0 ; l < recsize ; l++ )
01194 {
01195 if ( qfits_isnan(rectangle[l]) )
01196 {
01197 continue ;
01198 }
01199 mean -> data[i] += rectangle[l] ;
01200 nv ++;
01201 }
01202 if ( nv == 0 )
01203 {
01204 mean -> data[i] = ZERO ;
01205 }
01206 else
01207 {
01208 mean -> data[i] /= nv ;
01209 }
01210 cpl_free ( rectangle ) ;
01211 }
01212 return mean ;
01213 }
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226 Vector * meanCircleOfCubeSpectra( OneCube * cube,
01227 int centerx,
01228 int centery,
01229 int radius )
01230 {
01231 Vector * mean ;
01232 pixelvalue * circle ;
01233 int i, j, k, l, m, n ;
01234 int circsize, nv ;
01235
01236 if ( cube == NULL || cube -> np < 1 )
01237 {
01238 cpl_msg_error ("meanCircleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01239 return NullVector ;
01240 }
01241
01242 if ((centerx+radius>=cube->lx) ||
01243 (centery+radius>=cube->ly) ||
01244 (centerx-radius<0) ||
01245 (centery-radius<0))
01246 {
01247 cpl_msg_error("meanCircleOfCubeSpectra:"," invalid circular coordinates") ;
01248 return NullVector ;
01249 }
01250
01251 n = 0 ;
01252 for ( j = centery - radius ; j <= centery + radius ; j++ )
01253 {
01254 for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01255 {
01256 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01257 {
01258 n ++ ;
01259 }
01260 }
01261 }
01262 if (n == 0)
01263 {
01264 cpl_msg_error ("meanCircleOfCubeSpectra:"," no data points found!\n") ;
01265 return NullVector ;
01266 }
01267 circsize = n ;
01268
01269
01270 if (NULL == (mean = newVector (cube -> np)) )
01271 {
01272 cpl_msg_error ("meanCircleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01273 return NullVector ;
01274 }
01275
01276
01277
01278
01279
01280 for ( i = 0 ; i < cube->np ; i++ )
01281 {
01282 m = 0 ;
01283 circle = (pixelvalue *) cpl_calloc (circsize, sizeof (pixelvalue*));
01284
01285 for ( j = centery - radius ; j <= centery + radius ; j++ )
01286 {
01287 for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01288 {
01289 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01290 {
01291 circle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01292 m ++ ;
01293 }
01294 }
01295 }
01296
01297 nv = 0 ;
01298 for ( l = 0 ; l < circsize ; l++ )
01299 {
01300 if ( qfits_isnan(circle[l]) )
01301 {
01302 continue ;
01303 }
01304 mean -> data[i] += circle[l] ;
01305 nv ++;
01306 }
01307 if ( nv == 0 )
01308 {
01309 mean -> data[i] = ZERO ;
01310 }
01311 else
01312 {
01313 mean -> data[i] /= nv ;
01314 }
01315
01316 cpl_free (circle) ;
01317 }
01318 return mean ;
01319 }
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331 Vector * blackbodySpectrum( char * templateSpec, double temp )
01332 {
01333 Vector * retSpec ;
01334 fits_header * hdr ;
01335 int n ;
01336 double cenpix ;
01337 int npix ;
01338 double cenLambda ;
01339 double firstLambda ;
01340 double disp ;
01341 double lambda ;
01342 double intens ;
01343 double denom ;
01344 double norm ;
01345
01346 if ( NULL == templateSpec )
01347 {
01348 cpl_msg_error ("blackbodySpectrum:"," now input image given!\n") ;
01349 return NULL ;
01350 }
01351 if ( temp < 0. )
01352 {
01353 cpl_msg_error ("blackbodySpectrum:"," wrong temperature given!\n") ;
01354 return NULL ;
01355 }
01356
01357
01358 if ( NULL == (hdr = fits_read_header( templateSpec )) )
01359 {
01360 cpl_msg_error ("blackbodySpectrum:"," cannot read fits header and allocate data structure\n") ;
01361 return NULL ;
01362 }
01363
01364 if ( -1. == (cenpix = fits_header_getdouble(hdr, "CRPIX2", -1.)) )
01365 {
01366 cpl_msg_error ("blackbodySpectrum:"," cannot get CRPIX2\n") ;
01367 fits_header_destroy(hdr) ;
01368 return NULL ;
01369 }
01370 if ( -1. == (cenLambda = fits_header_getdouble(hdr, "CRVAL2", -1.)) )
01371 {
01372 cpl_msg_error ("blackbodySpectrum:"," cannot get CRVAL2\n") ;
01373 fits_header_destroy(hdr) ;
01374 return NULL ;
01375 }
01376 if ( -1. == (disp = fits_header_getdouble(hdr, "CDELT2", -1.)) )
01377 {
01378 cpl_msg_error ("blackbodySpectrum:"," cannot get CDELT2\n") ;
01379 fits_header_destroy(hdr) ;
01380 return NULL ;
01381 }
01382 if ( -1. == (npix = fits_header_getint(hdr, "NAXIS2", -1.)) )
01383 {
01384 cpl_msg_error ("blackbodySpectrum:"," cannot get NAXIS2\n") ;
01385 fits_header_destroy(hdr) ;
01386 return NULL ;
01387 }
01388 fits_header_destroy(hdr) ;
01389
01390 if (NULL == (retSpec = newVector (npix)))
01391 {
01392 cpl_msg_error ("blackbodySpectrum:"," could not allocate memory!\n") ;
01393 return NULL ;
01394 }
01395
01396
01397 cenpix-- ;
01398
01399 firstLambda = cenLambda - cenpix * disp ;
01400 for ( n = 0 ; n < npix ; n++ )
01401 {
01402 lambda = firstLambda + disp * (double)n ;
01403
01404
01405 lambda /= 1.0e6 ;
01406 denom = 1./(exp(PLANCK*SPEED_OF_LIGHT/(lambda*BOLTZMANN*temp)) - 1.) ;
01407 intens = 2.*PI_NUMB*PLANCK*SPEED_OF_LIGHT*SPEED_OF_LIGHT / pow(lambda, 5) * denom ;
01408 retSpec->data[n] = intens ;
01409 }
01410 norm = retSpec->data[npix/2] ;
01411 for ( n = 0 ; n < npix ; n++ )
01412 {
01413 retSpec->data[n] /= norm ;
01414 }
01415
01416 return retSpec ;
01417 }
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431 Vector * medianRectangleOfCubeSpectra( OneCube * cube,
01432 int llx,
01433 int lly,
01434 int urx,
01435 int ury )
01436 {
01437 Vector * med ;
01438 pixelvalue *rectangle ;
01439 int i, j, k, m ;
01440 int recsize ;
01441
01442 if ( cube == NULL || cube -> np < 1 )
01443 {
01444 cpl_msg_error ("medianRectangleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01445 return NullVector ;
01446 }
01447
01448 if ((llx<0) || (llx>=cube->lx) ||
01449 (urx<0) || (urx>=cube->lx) ||
01450 (lly<0) || (lly>=cube->ly) ||
01451 (ury<0) || (ury>=cube->ly) ||
01452 (llx>=urx) || (lly>=ury))
01453 {
01454 cpl_msg_error("medianRectangleOfCubeSpectra:"," invalid rectangle coordinates:") ;
01455 cpl_msg_error("medianRectangleOfCubeSpectra:","lower left is [%d %d] upper right is [%d %d]", llx, lly, urx, ury) ;
01456 return NullVector ;
01457 }
01458
01459 recsize = (urx - llx + 1) * (ury - lly + 1) ;
01460
01461
01462 if (NULL == (med = newVector (cube -> np)) )
01463 {
01464 cpl_msg_error ("medianRectangleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01465 return NullVector ;
01466 }
01467
01468
01469
01470
01471
01472 for ( i = 0 ; i < cube->np ; i++ )
01473 {
01474 m = 0 ;
01475 rectangle = (pixelvalue *) cpl_calloc (recsize, sizeof (pixelvalue*));
01476
01477 for ( j = lly ; j <= ury ; j++ )
01478 {
01479 for ( k = llx ; k <= urx ; k++ )
01480 {
01481 if ( qfits_isnan(cube->plane[i]->data[k+j*cube->lx]) )
01482 {
01483 continue ;
01484 }
01485 else
01486 {
01487 rectangle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01488 m ++ ;
01489 }
01490 }
01491 }
01492 if ( m == 0 )
01493 {
01494 med->data[i] = 0. ;
01495 }
01496 else
01497 {
01498 med->data[i] = median(rectangle, m) ;
01499 }
01500 cpl_free ( rectangle ) ;
01501 }
01502 return med ;
01503 }
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516 Vector * medianCircleOfCubeSpectra( OneCube * cube,
01517 int centerx,
01518 int centery,
01519 int radius )
01520 {
01521 Vector * med ;
01522 pixelvalue * circle ;
01523 int i, j, k, l, m, n ;
01524 int circsize, nv ;
01525
01526 if ( cube == NULL || cube -> np < 1 )
01527 {
01528 cpl_msg_error ("medianCircleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01529 return NullVector ;
01530 }
01531
01532 if ((centerx+radius>=cube->lx) ||
01533 (centery+radius>=cube->ly) ||
01534 (centerx-radius<0) ||
01535 (centery-radius<0))
01536 {
01537 cpl_msg_error("medianCircleOfCubeSpectra:"," invalid circular coordinates") ;
01538 return NullVector ;
01539 }
01540
01541 n = 0 ;
01542 for ( j = centery - radius ; j <= centery + radius ; j++ )
01543 {
01544 for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01545 {
01546 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01547 {
01548 n ++ ;
01549 }
01550 }
01551 }
01552 if (n == 0)
01553 {
01554 cpl_msg_error ("medianCircleOfCubeSpectra:"," no data points found!\n") ;
01555 return NullVector ;
01556 }
01557 circsize = n ;
01558
01559
01560 if (NULL == (med = newVector (cube -> np)) )
01561 {
01562 cpl_msg_error ("medianCircleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01563 return NullVector ;
01564 }
01565
01566
01567
01568
01569
01570 for ( i = 0 ; i < cube->np ; i++ )
01571 {
01572 m = 0 ;
01573 circle = (pixelvalue *) cpl_calloc (circsize, sizeof (pixelvalue*));
01574
01575 for ( j = centery - radius ; j <= centery + radius ; j++ )
01576 {
01577 for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01578 {
01579 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01580 {
01581 circle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01582 m ++ ;
01583 }
01584 }
01585 }
01586
01587 nv = 0 ;
01588 for ( l = 0 ; l < circsize ; l++ )
01589 {
01590 if ( qfits_isnan(circle[l]) )
01591 {
01592 continue ;
01593 }
01594 med -> data[i] += circle[l] ;
01595 nv ++;
01596 }
01597 if ( nv == 0 )
01598 {
01599 med->data[i] = 0. ;
01600 }
01601 else
01602 {
01603 med->data[i] = median(circle, nv) ;
01604 }
01605 cpl_free (circle) ;
01606 }
01607 return med ;
01608 }
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621 Vector * cleanmeanRectangleOfCubeSpectra( OneCube * cube,
01622 int llx,
01623 int lly,
01624 int urx,
01625 int ury,
01626 float lo_reject,
01627 float hi_reject )
01628 {
01629 Vector * clean ;
01630 pixelvalue *rectangle ;
01631 int i, j, k, m ;
01632 int recsize ;
01633
01634 if ( cube == NULL || cube -> np < 1 )
01635 {
01636 cpl_msg_error ("cleanmeanRectangleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01637 return NullVector ;
01638 }
01639
01640 if ((llx<0) || (llx>=cube->lx) ||
01641 (urx<0) || (urx>=cube->lx) ||
01642 (lly<0) || (lly>=cube->ly) ||
01643 (ury<0) || (ury>=cube->ly) ||
01644 (llx>=urx) || (lly>=ury))
01645 {
01646 cpl_msg_error("cleanmeanRectangleOfCubeSpectra:"," invalid rectangle coordinates:") ;
01647 cpl_msg_error("cleanmeanRectangleOfCubeSpectra:",
01648 "lower left is [%d %d] upper right is [%d %d]",
01649 llx, lly, urx, ury) ;
01650 return NullVector ;
01651 }
01652
01653 recsize = (urx - llx + 1) * (ury - lly + 1) ;
01654
01655
01656 if (NULL == (clean = newVector (cube -> np)) )
01657 {
01658 cpl_msg_error ("cleanmeanRectangleOfCubeSpectra:"," cannot allocate a new vector") ;
01659 return NullVector ;
01660 }
01661
01662
01663
01664
01665
01666 for ( i = 0 ; i < cube->np ; i++ )
01667 {
01668 m = 0 ;
01669 rectangle = (pixelvalue *) cpl_calloc (recsize, sizeof (pixelvalue*));
01670
01671 for ( j = lly ; j <= ury ; j++ )
01672 {
01673 for ( k = llx ; k <= urx ; k++ )
01674 {
01675 if ( qfits_isnan(cube->plane[i]->data[k+j*cube->lx]) )
01676 {
01677 continue ;
01678 }
01679 else
01680 {
01681 rectangle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01682 m ++ ;
01683 }
01684 }
01685 }
01686 if ( m == 0 )
01687 {
01688 clean->data[i] = 0. ;
01689 }
01690 else
01691 {
01692 clean->data[i] = clean_mean(rectangle, m, lo_reject, hi_reject) ;
01693 }
01694 cpl_free ( rectangle ) ;
01695 }
01696 return clean ;
01697 }
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710 Vector * cleanmeanCircleOfCubeSpectra( OneCube * cube,
01711 int centerx,
01712 int centery,
01713 int radius,
01714 float lo_reject,
01715 float hi_reject )
01716 {
01717 Vector * clean ;
01718 pixelvalue * circle ;
01719 int i, j, k, l, m, n ;
01720 int circsize, nv ;
01721
01722 if ( cube == NULL || cube -> np < 1 )
01723 {
01724 cpl_msg_error ("cleanmeanCircleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01725 return NullVector ;
01726 }
01727
01728 if ((centerx+radius>=cube->lx) ||
01729 (centery+radius>=cube->ly) ||
01730 (centerx-radius<0) ||
01731 (centery-radius<0))
01732 {
01733 cpl_msg_error("cleanmeanCircleOfCubeSpectra:"," invalid circular coordinates") ;
01734 return NullVector ;
01735 }
01736
01737 n = 0 ;
01738 for ( j = centery - radius ; j <= centery + radius ; j++ )
01739 {
01740 for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01741 {
01742 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01743 {
01744 n ++ ;
01745 }
01746 }
01747 }
01748 if (n == 0)
01749 {
01750 cpl_msg_error ("cleanmeanCircleOfCubeSpectra:"," no data points found!\n") ;
01751 return NullVector ;
01752 }
01753 circsize = n ;
01754
01755
01756 if (NULL == (clean = newVector (cube -> np)) )
01757 {
01758 cpl_msg_error ("cleanmeanCircleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01759 return NullVector ;
01760 }
01761
01762
01763
01764
01765
01766 for ( i = 0 ; i < cube->np ; i++ )
01767 {
01768 m = 0 ;
01769 circle = (pixelvalue *) cpl_calloc (circsize, sizeof (pixelvalue*));
01770
01771 for ( j = centery - radius ; j <= centery + radius ; j++ )
01772 {
01773 for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01774 {
01775 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01776 {
01777 circle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01778 m ++ ;
01779 }
01780 }
01781 }
01782
01783 nv = 0 ;
01784 for ( l = 0 ; l < circsize ; l++ )
01785 {
01786 if ( qfits_isnan(circle[l]) )
01787 {
01788 continue ;
01789 }
01790 clean -> data[i] += circle[l] ;
01791 nv ++;
01792 }
01793 if ( nv == 0 )
01794 {
01795 clean->data[i] = 0. ;
01796 }
01797 else
01798 {
01799 clean->data[i] = clean_mean(circle, nv, lo_reject, hi_reject) ;
01800 }
01801 cpl_free (circle) ;
01802 }
01803 return clean ;
01804 }
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817 float * shiftArray ( float * input, int n_elements, float shift, double * ker )
01818 {
01819 float * shifted ;
01820 int samples = KERNEL_SAMPLES ;
01821 int i ;
01822 float fx ;
01823 float rx ;
01824 int px ;
01825 int tabx ;
01826 float value ;
01827
01828 register float * pix ;
01829 int mid;
01830 float norm ;
01831
01832
01833 if (input==NULL)
01834 {
01835 cpl_msg_error("shiftArray:"," no input array given!\n") ;
01836 return NULL ;
01837 }
01838 if (n_elements<=0)
01839 {
01840 cpl_msg_error("shiftArray:"," wrong number of elements in input array given!\n") ;
01841 return NULL ;
01842 }
01843
01844 shifted = (float*) cpl_calloc(n_elements, sizeof(float)) ;
01845
01846
01847 if ((fabs(shift)<1e-2))
01848 {
01849 for (i = 0 ; i < n_elements ; i++ )
01850 {
01851 shifted[i] = input[i] ;
01852 }
01853 return shifted ;
01854 }
01855
01856 mid = (int)samples/(int)2 ;
01857
01858 for (i=1 ; i< n_elements-2 ; i++)
01859 {
01860 fx = (float)i+shift ;
01861 px = nint(fx) ;
01862 rx = fx - (float)px ;
01863 pix = input ;
01864
01865 if ((px>=1) && (px<(n_elements-2)))
01866 {
01867 tabx = (int)(fabs((float)mid * rx)) ;
01868
01869 if (qfits_isnan(pix[i]))
01870 {
01871 value = ZERO ;
01872 }
01873 else
01874 {
01875 if (qfits_isnan(pix[i-1]))
01876 {
01877 pix[i-1] = 0. ;
01878 }
01879 if (qfits_isnan(pix[i+1]))
01880 {
01881 pix[i+1] = 0. ;
01882 }
01883 if (qfits_isnan(pix[i+2]))
01884 {
01885 pix[i+2] = 0. ;
01886 }
01887
01888
01889
01890
01891
01892 value = pix[i-1] * ker[mid+tabx] +
01893 pix[i] * ker[tabx] +
01894 pix[i+1] * ker[mid-tabx] +
01895 pix[i+2] * ker[samples-tabx-1] ;
01896
01897
01898
01899
01900 norm = ker[mid+tabx] +
01901 ker[tabx] +
01902 ker[mid-tabx] +
01903 ker[samples-tabx-1] ;
01904 if (fabs(norm) > 1e-4)
01905 {
01906 value /= norm ;
01907 }
01908 }
01909 }
01910 else
01911 {
01912 value = 0.0 ;
01913 }
01914 if ( qfits_isnan(value) )
01915 {
01916 shifted[i] = ZERO ;
01917 }
01918 else
01919 {
01920 shifted[i] = value ;
01921 }
01922 }
01923 return shifted ;
01924 }
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938 OneImage * divImageBySpectrum( OneImage * image, OneImage * spectrum )
01939 {
01940 int col, row ;
01941 OneImage * retImage ;
01942 if ( image == NullImage )
01943 {
01944 cpl_msg_error("divImageBySpectrum:","no image given!") ;
01945 return NullImage ;
01946 }
01947 if ( spectrum == NullImage )
01948 {
01949 cpl_msg_error("divImageBySpectrum:","no spectrum image given!") ;
01950 return NullImage ;
01951 }
01952 if ( spectrum -> ly != image -> ly )
01953 {
01954 cpl_msg_error("divImageBySpectrum:","images are not compatible in pixel length!") ;
01955 return NullImage ;
01956 }
01957 if ( NullImage == (retImage = copy_image(image)) )
01958 {
01959 cpl_msg_error("divImageBySpectrum:","could not copy original image!") ;
01960 return NullImage ;
01961 }
01962 for ( col = 0 ; col < image -> lx ; col++ )
01963 {
01964 for ( row = 0 ; row < image -> ly ; row++ )
01965 {
01966 if ( !qfits_isnan(image->data[col+row*image->lx]) &&
01967 !qfits_isnan(spectrum->data[col+row*image->lx]))
01968 {
01969 retImage->data[col+row*image->lx] =
01970 image->data[col+row*image->lx] / spectrum->data[row] ;
01971 }
01972 }
01973 }
01974 return retImage ;
01975 }
01976
01977