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 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032 #include <fors_image.h>
00033
00034 #include <fors_dfs.h>
00035 #include <fors_utils.h>
00036 #include <fors_pfits.h>
00037 #include <fors_double.h>
00038
00039 #include <moses.h>
00040
00041 #include <cpl.h>
00042
00043 #include <math.h>
00044 #include <stdbool.h>
00045 #include <stdio.h>
00046
00062 struct _fors_image
00063 {
00064 cpl_image *data;
00065 cpl_image *variance;
00066
00067
00068
00069
00070
00071
00072
00073 };
00074
00078 const cpl_type FORS_IMAGE_TYPE = CPL_TYPE_FLOAT;
00079 #define FORS_IMAGE_TYPE_MAX FLT_MAX
00080
00081
00082 #undef cleanup
00083
00084
00085
00086
00087
00088
00089
00090
00091 static int
00092 max_filter(const float *ibuffer, float *obuffer, int length, int size)
00093 {
00094 float max;
00095 int start = size / 2;
00096 int end = length - size / 2;
00097 int i, j;
00098
00099
00100 for (i = start; i < end; i++) {
00101 max = ibuffer[i-start];
00102 for (j = i - start + 1; j <= i + start; j++)
00103 if (max < ibuffer[j])
00104 max = ibuffer[j];
00105 obuffer[i] = max;
00106 }
00107
00108 for (i = 0; i < start; i++)
00109 obuffer[i] = obuffer[start];
00110
00111 for (i = end; i < length; i++)
00112 obuffer[i] = obuffer[end-1];
00113
00114 return 0;
00115 }
00116
00117 #define cleanup
00118
00127 fors_image *
00128 fors_image_new(cpl_image *data, cpl_image *variance)
00129 {
00130 fors_image *image = NULL;
00131
00132 assure( data != NULL, return NULL, NULL );
00133 assure( variance != NULL, return NULL, NULL );
00134
00135 assure( cpl_image_get_type(data) == FORS_IMAGE_TYPE, return NULL,
00136 "Provided data image type is %s, must be %s",
00137 fors_type_get_string(cpl_image_get_type(data)),
00138 fors_type_get_string(FORS_IMAGE_TYPE) );
00139
00140 assure( cpl_image_get_type(variance) == FORS_IMAGE_TYPE, return NULL,
00141 "Provided weight image type is %s, must be %s",
00142 fors_type_get_string(cpl_image_get_type(variance)),
00143 fors_type_get_string(FORS_IMAGE_TYPE) );
00144
00145 assure( cpl_image_get_size_x(data) == cpl_image_get_size_x(variance) &&
00146 cpl_image_get_size_y(data) == cpl_image_get_size_y(variance),
00147 return NULL,
00148 "Incompatible data and weight image sizes: %dx%d and %dx%d",
00149 cpl_image_get_size_x(data), cpl_image_get_size_y(data),
00150 cpl_image_get_size_x(variance), cpl_image_get_size_y(variance));
00151
00152 assure( cpl_image_get_min(variance) >= 0, return NULL,
00153 "Variances must be non-negative, minimum is %f",
00154 cpl_image_get_min(variance));
00155
00156 image = cpl_malloc(sizeof(*image));
00157
00158 image->data = data;
00159 image->variance = variance;
00160
00161 return image;
00162 }
00163
00164 #undef cleanup
00165 #define cleanup
00166
00171 fors_image *
00172 fors_image_duplicate(const fors_image *image)
00173 {
00174 assure( image != NULL, return NULL, NULL );
00175
00176 return fors_image_new(cpl_image_duplicate(image->data),
00177 cpl_image_duplicate(image->variance));
00178 }
00179
00184 void
00185 fors_image_delete(fors_image **image)
00186 {
00187 if (image && *image) {
00188 cpl_image_delete((*image)->data);
00189 cpl_image_delete((*image)->variance);
00190 cpl_free(*image); *image = NULL;
00191 }
00192 return;
00193 }
00194
00199 void
00200 fors_image_delete_const(const fors_image **image)
00201 {
00202 fors_image_delete((fors_image **)image);
00203
00204 return;
00205 }
00206
00207
00208 #if 0
00209
00213 static void
00214 fors_image_dump(const fors_image *image, FILE *file)
00215 {
00216 if (image == NULL) {
00217 fprintf(file, "Null image\n");
00218 }
00219 else {
00220 cpl_stats *stats;
00221
00222 fprintf(file, "Data:\n");
00223 stats = cpl_stats_new_from_image(image->data, CPL_STATS_ALL);
00224 cpl_stats_dump(stats, CPL_STATS_ALL, file);
00225 cpl_stats_delete(stats);
00226
00227 fprintf(file, "Variance:\n");
00228 stats = cpl_stats_new_from_image(image->variance, CPL_STATS_ALL);
00229 cpl_stats_dump(stats, CPL_STATS_ALL, file);
00230 cpl_stats_delete(stats);
00231 }
00232
00233 return;
00234 }
00235 #endif
00236
00237 #undef cleanup
00238 #define cleanup \
00239 do { \
00240 double_list_delete(&sat_percent, double_delete); \
00241 } while (0)
00242
00256 fors_image_list *
00257 fors_image_load_list(const cpl_frameset *frames, const fors_image *bias,
00258 const fors_setting *setting,
00259 double *saturated)
00260 {
00261 fors_image_list *ilist = fors_image_list_new();
00262 double_list *sat_percent = double_list_new();
00263
00264 assure( frames != NULL, return ilist, NULL );
00265 assure( !cpl_frameset_is_empty(frames), return ilist, "Empty frameset");
00266
00267 {
00268 const cpl_frame *f;
00269
00270 for (f = cpl_frameset_get_first_const(frames);
00271 f != NULL;
00272 f = cpl_frameset_get_next_const(frames)) {
00273
00274 double saturated_one;
00275
00276 fors_image *i = fors_image_load(f, bias, setting,
00277 saturated != NULL ?
00278 &saturated_one : NULL);
00279 assure( !cpl_error_get_code(), return ilist, NULL );
00280
00281 fors_image_list_insert(ilist, i);
00282 if (saturated != NULL) {
00283 double_list_insert(sat_percent,
00284 double_duplicate(&saturated_one));
00285 }
00286 }
00287 }
00288
00289
00290
00291
00292
00293 if (saturated != NULL) {
00294 *saturated = double_list_mean(sat_percent, double_eval, NULL);
00295 }
00296
00297 cleanup;
00298 return ilist;
00299 }
00300
00312 const fors_image_list *
00313 fors_image_load_list_const(const cpl_frameset *frames, const fors_image *bias,
00314 const fors_setting *setting,
00315 double *saturated)
00316 {
00317 return (const fors_image_list *)
00318 fors_image_load_list(frames, bias, setting, saturated);
00319 }
00320
00321 #undef cleanup
00322 #define cleanup \
00323 do { \
00324 cpl_propertylist_delete(header); \
00325 cpl_table_delete(overscans); \
00326 cpl_image_delete(temp); \
00327 cpl_mask_delete(non_saturated); \
00328 fors_setting_delete(&frame_setting); \
00329 } while (0)
00330
00348 fors_image *
00349 fors_image_load(const cpl_frame *frame, const fors_image *bias,
00350 const fors_setting *setting,
00351 double *saturated)
00352 {
00353 fors_image *image = NULL;
00354 fors_setting *frame_setting = NULL;
00355 cpl_image *data = NULL;
00356 cpl_image *variance = NULL;
00357 cpl_image *temp = NULL;
00358 cpl_propertylist *header = NULL;
00359 cpl_table *overscans = NULL;
00360 cpl_mask *non_saturated = NULL;
00361 const char *filename;
00362 int extension= 0;
00363 const int plane = 0;
00364
00365 assure( frame != NULL, return image, NULL );
00366
00367 assure( setting != NULL, return image, NULL );
00368 filename = cpl_frame_get_filename(frame);
00369 assure( filename != NULL, return image,
00370 "NULL filename received");
00371
00372 cpl_msg_info(cpl_func, "Loading %s: %s",
00373
00374 (cpl_frame_get_tag(frame) != NULL) ?
00375 cpl_frame_get_tag(frame) : "NULL",
00376 filename);
00377
00378
00379 fors_setting_verify(setting, frame, &frame_setting);
00380 assure( !cpl_error_get_code(), return image,
00381 "Could not verify %s setting",
00382 filename);
00383
00384
00385
00386 header = cpl_propertylist_load(filename, extension);
00387 assure( !cpl_error_get_code(), return image,
00388 "Could not load %s extension %d header",
00389 filename, extension);
00390
00391
00392 data = cpl_image_load(filename,
00393 FORS_IMAGE_TYPE, plane, extension);
00394
00395 assure( !cpl_error_get_code(), return image,
00396 "Could not load image from %s extension %d",
00397 filename, extension);
00398
00399
00400 {
00401 int xlow, ylow, xhig, yhig;
00402 bool check_consistency = false;
00403
00404 overscans = mos_load_overscans_vimos(header, check_consistency);
00405
00406
00407
00408 assure( !cpl_error_get_code(), return image,
00409 "Could not read overscan information from %s extension %d",
00410 filename, extension);
00411
00412 xlow = cpl_table_get_int(overscans, "xlow", 0, NULL);
00413 ylow = cpl_table_get_int(overscans, "ylow", 0, NULL);
00414 xhig = cpl_table_get_int(overscans, "xhig", 0, NULL);
00415 yhig = cpl_table_get_int(overscans, "yhig", 0, NULL);
00416
00417 cpl_table_delete(overscans); overscans = NULL;
00418
00419 temp = cpl_image_duplicate(data);
00420 cpl_image_delete(data);
00421 data = cpl_image_extract(temp, xlow+1, ylow+1, xhig, yhig);
00422
00423 cpl_image_delete(temp); temp = NULL;
00424 }
00425
00426 cpl_propertylist_delete(header); header = NULL;
00427
00428
00429 if (cpl_frame_get_nextensions(frame) == 0 || bias != NULL) {
00430
00431
00432
00433
00434
00435
00436 variance = cpl_image_new(
00437 cpl_image_get_size_x(data),
00438 cpl_image_get_size_y(data),
00439 FORS_IMAGE_TYPE);
00440
00441 if (cpl_frame_get_group(frame) == CPL_FRAME_GROUP_CALIB) {
00442 cpl_msg_warning(cpl_func,
00443 "No error bars provided for calibration frame %s, "
00444 "assuming no errors. For complete error propagation, "
00445 "you may recreate this frame with this pipeline",
00446 filename);
00447 } else {
00448 cpl_image_add_scalar(variance, frame_setting->ron*frame_setting->ron);
00449 }
00450 }
00451 else {
00452
00453 extension = 1;
00454
00455
00456 variance = cpl_image_load(filename,
00457 FORS_IMAGE_TYPE, plane, extension);
00458
00459 assure( !cpl_error_get_code(), return image,
00460 "Could not load image from %s extension %d",
00461 filename, extension);
00462
00463 cpl_image_power(variance, 2);
00464
00465 assure( cpl_image_get_min(variance) >= 0, return image,
00466 "Illegal minimum variance: %g",
00467 cpl_image_get_min(variance));
00468
00469 cpl_image_delete(temp); temp = NULL;
00470 }
00471
00472 image = fors_image_new(data, variance);
00473
00474
00475 if (saturated != NULL) {
00476 double lo_cut = 0.5;
00477 double hi_cut = 65534.5;
00478 non_saturated = cpl_mask_threshold_image_create(
00479 data, lo_cut, hi_cut);
00480
00481 *saturated = (cpl_image_get_size_x(data) *
00482 cpl_image_get_size_y(data) -
00483 cpl_mask_count(non_saturated)) * 100.0 /
00484 (cpl_image_get_size_x(data) *
00485 cpl_image_get_size_y(data));
00486
00487 cpl_mask_delete(non_saturated); non_saturated = NULL;
00488
00489 cpl_msg_debug(cpl_func,
00490 "%f %% saturated pixels", *saturated);
00491 }
00492
00493 if (bias == NULL) {
00494
00495 }
00496 else {
00497 assure( cpl_frame_get_group(frame) == CPL_FRAME_GROUP_RAW,
00498 return image,
00499 "Refusing to subtract bias from non-raw (%s) input frame: %s",
00500 fors_frame_get_group_string(frame),
00501 filename);
00502
00503 cpl_msg_debug(cpl_func, "Subtracting bias from %s", filename);
00504
00505 fors_image_subtract(image, bias);
00506
00507 assure( !cpl_error_get_code(), return image,
00508 "Bias subtraction failed" );
00509
00510
00511
00512
00513
00514
00515
00516 double conad = 1.0 / frame_setting->average_gain;
00517
00518 temp = cpl_image_divide_scalar_create(image->data, conad);
00519 cpl_image_abs(temp);
00520
00521 cpl_image_add(image->variance, temp);
00522 }
00523
00524 cleanup;
00525 return image;
00526 }
00527
00528
00529 #undef cleanup
00530 #define cleanup \
00531 do { \
00532 cpl_image_delete(sigma); \
00533 } while(0)
00534
00543 void
00544 fors_image_save(const fors_image *image, const cpl_propertylist *header,
00545 const char *filename)
00546 {
00547 cpl_propertylist *extension_header = NULL;
00548 cpl_image *sigma = NULL;
00549
00550 assure( image != NULL, return, NULL );
00551
00552 assure( filename != NULL, return, NULL );
00553
00554 cpl_image_save(image->data, filename, CPL_BPP_IEEE_FLOAT, header,
00555 CPL_IO_DEFAULT);
00556 assure( !cpl_error_get_code(), return,
00557 "Cannot save product %s", filename);
00558
00559 sigma = cpl_image_power_create(image->variance, 0.5);
00560
00561
00562 cpl_image_save(sigma, filename, CPL_BPP_IEEE_FLOAT, extension_header,
00563 CPL_IO_EXTEND);
00564 assure( !cpl_error_get_code(), return,
00565 "Cannot save product %s", filename);
00566
00567 cleanup;
00568 return;
00569 }
00570
00571
00572 #undef cleanup
00573 #define cleanup \
00574 do { \
00575 cpl_image_delete(var_bkg); \
00576 cpl_image_delete(sigma_bkg); \
00577 } while(0)
00578
00589 void
00590 fors_image_save_sex(const fors_image *image, const cpl_propertylist *header,
00591 const char *filename_dat,
00592 const char *filename_var,
00593 int radius)
00594 {
00595 cpl_propertylist *extension_header = NULL;
00596 cpl_image *sigma_bkg = NULL;
00597 cpl_image *var_bkg = NULL;
00598
00599 assure( image != NULL, return, NULL );
00600
00601 assure( filename_dat != NULL, return, NULL );
00602 assure( filename_var != NULL, return, NULL );
00603
00604 cpl_image_save(image->data, filename_dat, CPL_BPP_IEEE_FLOAT, header,
00605 CPL_IO_DEFAULT);
00606 assure( !cpl_error_get_code(), return,
00607 "Cannot save product %s", filename_dat);
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617 cpl_msg_info(cpl_func, "Creating background error map");
00618
00619 bool filter_data = false;
00620 int xstep = radius/2;
00621
00622
00623
00624
00625
00626
00627 int ystep = radius/2;
00628 int xstart = 1;
00629 int ystart = 1;
00630 int xend = fors_image_get_size_x(image);
00631 int yend = fors_image_get_size_y(image);
00632
00633
00634 var_bkg = fors_image_filter_median_create(image,
00635 radius,
00636 radius,
00637 xstart, ystart,
00638 xend, yend,
00639 xstep, ystep,
00640 filter_data);
00641 assure( !cpl_error_get_code(), return,
00642 "Median filtering failed");
00643
00644 sigma_bkg = cpl_image_power_create(var_bkg, 0.5);
00645
00646 cpl_image_save(sigma_bkg, filename_var,
00647 CPL_BPP_IEEE_FLOAT, extension_header,
00648 CPL_IO_DEFAULT);
00649 assure( !cpl_error_get_code(), return,
00650 "Cannot save product %s", filename_var);
00651
00652 cleanup;
00653 return;
00654 }
00655
00656 #undef cleanup
00657 #define cleanup
00658
00663 int fors_image_get_size_x(const fors_image *image)
00664 {
00665 assure( image != NULL, return -1, NULL );
00666 return cpl_image_get_size_x(image->data);
00667 }
00668
00669 #undef cleanup
00670 #define cleanup
00671
00676 int fors_image_get_size_y(const fors_image *image)
00677 {
00678 assure( image != NULL, return -1, NULL );
00679 return cpl_image_get_size_y(image->data);
00680 }
00681
00682 #undef cleanup
00683 #define cleanup
00684
00688 const float *fors_image_get_data_const(const fors_image *image)
00689 {
00690 assure( image != NULL, return NULL, NULL );
00691
00692 assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return NULL, NULL );
00693
00694
00695
00696 return cpl_image_get_data_float(image->data);
00697 }
00698
00699 #undef cleanup
00700 #define cleanup
00701
00708 void
00709 fors_image_abs(fors_image *image)
00710 {
00711 assure( image != NULL, return, NULL );
00712
00713 cpl_image_abs(image->data);
00714
00715 return;
00716 }
00717
00718 #undef cleanup
00719 #define cleanup
00720
00727 void
00728 fors_image_square(fors_image *image)
00729 {
00730 assure( image != NULL, return, NULL );
00731
00732 cpl_image_multiply(image->data, image->data);
00733
00734
00735
00736
00737 cpl_image_multiply_scalar(image->variance, 2);
00738
00739 return;
00740 }
00741
00742
00743 #undef cleanup
00744 #define cleanup \
00745 do { \
00746 cpl_image_delete(temp); \
00747 } while(0)
00748
00756 void
00757 fors_image_subtract(fors_image *left, const fors_image *right)
00758 {
00759 cpl_image *temp = NULL;
00760 assure( left != NULL, return, NULL );
00761 assure( right != NULL, return, NULL );
00762
00763 cpl_image_subtract(left->data, right->data);
00764
00765
00766 cpl_image_add(left->variance, right->variance);
00767
00768 cleanup;
00769 return;
00770 }
00771
00772 #undef cleanup
00773 #define cleanup
00774
00784 void
00785 fors_image_multiply_noerr(fors_image *left, const cpl_image *right)
00786 {
00787 assure( left != NULL, return, NULL );
00788 assure( right != NULL, return, NULL );
00789 assure( cpl_image_get_size_x(left->data) == cpl_image_get_size_x(right) &&
00790 cpl_image_get_size_y(left->data) == cpl_image_get_size_y(right),
00791 return, "Incompatible image sizes: %dx%d and %dx%d",
00792 cpl_image_get_size_x(left->data),
00793 cpl_image_get_size_y(left->data),
00794 cpl_image_get_size_x(right),
00795 cpl_image_get_size_y(right));
00796
00797 cpl_image_multiply(left->data, right);
00798 cpl_image_multiply(left->variance, right);
00799 cpl_image_multiply(left->variance, right);
00800
00801 return;
00802 }
00803
00804 #undef cleanup
00805 #define cleanup
00806
00822 void
00823 fors_image_divide_noerr(fors_image *left, cpl_image *right)
00824 {
00825 assure( left != NULL, return, NULL );
00826 assure( right != NULL, return, NULL );
00827 assure( cpl_image_get_size_x(left->data) == cpl_image_get_size_x(right) &&
00828 cpl_image_get_size_y(left->data) == cpl_image_get_size_y(right),
00829 return, "Incompatible image sizes: %dx%d and %dx%d",
00830 cpl_image_get_size_x(left->data),
00831 cpl_image_get_size_y(left->data),
00832 cpl_image_get_size_x(right),
00833 cpl_image_get_size_y(right));
00834
00835 int x, y;
00836 int nx = cpl_image_get_size_x(right);
00837 int ny = cpl_image_get_size_y(right);
00838 float *datal = cpl_image_get_data_float(left->data);
00839 float *datav = cpl_image_get_data_float(left->variance);
00840 float *datar = cpl_image_get_data_float(right);
00841 for (y = 0; y < ny; y++) {
00842 for (x = 0; x < nx; x++) {
00843 if (datar[x + nx*y] == 0) {
00844 datar[x + nx*y] = 1;
00845 datal[x + nx*y] = 1;
00846
00847 datav[x + nx*y] = FORS_IMAGE_TYPE_MAX;
00848 }
00849 }
00850 }
00851
00852 cpl_image_divide(left->data, right);
00853 cpl_image_divide(left->variance, right);
00854 cpl_image_divide(left->variance, right);
00855
00856 return;
00857 }
00858
00859 #undef cleanup
00860 #define cleanup \
00861 do { \
00862 fors_image_delete(&dupl); \
00863 } while(0)
00864
00884 void
00885 fors_image_divide(fors_image *left, const fors_image *right)
00886 {
00887 fors_image *dupl = NULL;
00888
00889 assure( left != NULL, return, NULL );
00890 assure( right != NULL, return, NULL );
00891
00892 dupl = fors_image_duplicate(right);
00893
00894 cpl_image_divide(left->data, dupl->data);
00895
00896
00897 cpl_image_multiply(dupl->variance, left->data);
00898 cpl_image_multiply(dupl->variance, left->data);
00899
00900
00901
00902 cpl_image_add(left->variance, dupl->variance);
00903
00904
00905
00906 cpl_image_divide(left->variance, dupl->data);
00907 cpl_image_divide(left->variance, dupl->data);
00908
00909
00910
00911 int x, y;
00912 int nx = cpl_image_get_size_x(left->data);
00913 int ny = cpl_image_get_size_y(left->data);
00914 float *datal = cpl_image_get_data_float(left->data);
00915 float *datav = cpl_image_get_data_float(left->variance);
00916 float *datar = cpl_image_get_data_float(right->data);
00917 for (y = 0; y < ny; y++) {
00918 for (x = 0; x < nx; x++) {
00919 if (datar[x + nx*y] == 0) {
00920 datal[x + nx*y] = 1;
00921 datav[x + nx*y] = FORS_IMAGE_TYPE_MAX;
00922 }
00923 }
00924 }
00925
00926 cleanup;
00927 return;
00928 }
00929
00930 #undef cleanup
00931 #define cleanup \
00932 do { \
00933 cpl_image_delete(s22d12); \
00934 } while(0)
00935
00946 void
00947 fors_image_multiply(fors_image *left, const fors_image *right)
00948 {
00949 cpl_image *s22d12 = NULL;
00950
00951 assure( left != NULL, return, NULL );
00952 assure( right != NULL, return, NULL );
00953
00954 s22d12 = cpl_image_duplicate(right->variance);
00955 cpl_image_multiply(s22d12, left->data);
00956 cpl_image_multiply(s22d12, left->data);
00957
00958 cpl_image_multiply(left->variance, right->data);
00959 cpl_image_multiply(left->variance, right->data);
00960 cpl_image_add(left->variance, s22d12);
00961
00962 cpl_image_multiply(left->data, right->data);
00963
00964 cleanup;
00965 return;
00966 }
00967
00968
00969 #undef cleanup
00970 #define cleanup
00971
00983 void fors_image_subtract_scalar(fors_image *image, double s, double ds)
00984 {
00985 assure( image != NULL, return, NULL );
00986 assure( ds <= 0, return, "Unsupported");
00987
00988 cpl_image_subtract_scalar(image->data, s);
00989
00990 return;
00991 }
00992
00993
00994 #undef cleanup
00995 #define cleanup
00996
01008 void fors_image_divide_scalar(fors_image *image, double s, double ds)
01009 {
01010 assure( image != NULL, return, NULL );
01011 assure( s != 0, return, "Division by zero");
01012 assure( ds <= 0, return, "Unsupported");
01013
01014 cpl_image_divide_scalar(image->data, s);
01015 cpl_image_divide_scalar(image->variance, s*s);
01016
01017 return;
01018 }
01019
01020 #undef cleanup
01021 #define cleanup
01022
01034 void fors_image_multiply_scalar(fors_image *image, double s, double ds)
01035 {
01036 assure( image != NULL, return, NULL );
01037 assure( ds <= 0, return, "Unsupported");
01038
01039 cpl_image_multiply_scalar(image->data, s);
01040 cpl_image_multiply_scalar(image->variance, s*s);
01041
01042 return;
01043 }
01044
01045 #undef cleanup
01046 #define cleanup \
01047 do { \
01048 cpl_image_delete(temp); \
01049 } while(0)
01050
01063 void fors_image_exponential(fors_image *image, double b, double db)
01064 {
01065 cpl_image *temp = NULL;
01066
01067 assure( image != NULL, return, NULL );
01068 assure( b >= 0, return, "Negative base: %f", b);
01069 assure( db <= 0, return, "Unsupported");
01070
01071 cpl_image_exponential(image->data, b);
01072
01073 double lnb = log(b);
01074
01075 cpl_image_multiply_scalar(image->variance, lnb*lnb);
01076 cpl_image_multiply(image->variance, image->data);
01077 cpl_image_multiply(image->variance, image->data);
01078
01079 return;
01080 }
01081
01082
01083 #undef cleanup
01084 #define cleanup
01085
01090 double
01091 fors_image_get_min(const fors_image *image)
01092 {
01093 assure( image != NULL, return 0, NULL );
01094
01095 return cpl_image_get_min(image->data);
01096 }
01097
01098 #undef cleanup
01099 #define cleanup
01100
01105 double
01106 fors_image_get_max(const fors_image *image)
01107 {
01108 assure( image != NULL, return 0, NULL );
01109
01110 return cpl_image_get_max(image->data);
01111 }
01112
01113 #undef cleanup
01114 #define cleanup
01115
01121 double
01122 fors_image_get_mean(const fors_image *image, double *dmean)
01123 {
01124 assure( image != NULL, return 0, NULL );
01125 assure( dmean == NULL, return 0, "Unsupported");
01126
01127 return cpl_image_get_mean(image->data);
01128 }
01129
01130 #undef cleanup
01131 #define cleanup
01132
01138 double
01139 fors_image_get_median(const fors_image *image, double *dmedian)
01140 {
01141 assure( image != NULL, return 0, NULL );
01142 assure( dmedian == NULL, return 0, "Unsupported");
01143
01144 return cpl_image_get_median(image->data);
01145 }
01146
01147
01148 #undef cleanup
01149 #define cleanup
01150
01164 void fors_image_crop(fors_image *image,
01165 int xlo, int ylo,
01166 int xhi, int yhi)
01167 {
01168
01169
01170 assure( image != NULL, return, NULL );
01171 assure( 1 <= xlo && xlo <= xhi && xhi <= fors_image_get_size_x(image) &&
01172 1 <= ylo && ylo <= yhi && yhi <= fors_image_get_size_y(image),
01173 return, "Cannot extraction region (%d, %d) - (%d, %d) of "
01174 "%dx%d image",
01175 xlo, ylo, xhi, yhi,
01176 fors_image_get_size_x(image),
01177 fors_image_get_size_y(image));
01178
01179 cpl_image *new_data = cpl_image_extract(image->data,
01180 xlo, ylo,
01181 xhi, yhi);
01182 cpl_image_delete(image->data);
01183
01184 cpl_image* new_variance = cpl_image_extract(image->variance,
01185 xlo, ylo,
01186 xhi, yhi);
01187 cpl_image_delete(image->variance);
01188
01189 image->data = new_data;
01190 image->variance = new_variance;
01191
01192 return;
01193 }
01194
01220 cpl_image *
01221 fors_image_filter_median_create(const fors_image *image,
01222 int xradius,
01223 int yradius,
01224 int xstart,
01225 int ystart,
01226 int xend,
01227 int yend,
01228 int xstep,
01229 int ystep,
01230 bool use_data)
01231 {
01232 const cpl_image *input = NULL;
01233 cpl_image *smooth = NULL;
01234 int nx, ny;
01235
01236 assure( image != NULL, return smooth, NULL );
01237 passure( image->data != NULL, return smooth );
01238 passure( image->variance != NULL, return smooth );
01239
01240 input = (use_data) ? image->data : image->variance;
01241
01242 nx = cpl_image_get_size_x(input);
01243 ny = cpl_image_get_size_y(input);
01244
01245 if (xstep < 1) xstep = 1;
01246 if (ystep < 1) ystep = 1;
01247
01248 assure( 1 <= xstart && xstart <= xend && xend <= nx &&
01249 1 <= ystart && ystart <= yend && yend <= ny, return smooth,
01250 "Illegal region (%d, %d) - (%d, %d) of %dx%d image",
01251 xstart, ystart,
01252 xend, yend,
01253 nx, ny);
01254
01255 smooth = cpl_image_duplicate(input);
01256
01257
01258 assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return smooth, NULL );
01259
01260 const float *input_data = cpl_image_get_data_float_const(input);
01261 float *smooth_data = cpl_image_get_data_float(smooth);
01262 float *data = cpl_malloc((2*yradius + 1)*(2*xradius + 1)*sizeof(*data));
01263
01264 int y;
01265 for (y = ystart; y < yend; y++) {
01266
01267
01268
01269
01270
01271
01272 int ylo = y - (yradius/ystep) * ystep;
01273 int yhi = y + (yradius/ystep) * ystep;
01274
01275 while (ylo < ystart) ylo += ystep;
01276 while (yhi > yend ) yhi -= ystep;
01277
01278 int x;
01279 for (x = xstart; x < xend; x++) {
01280 int xlo = x - (xradius/xstep) * xstep;
01281 int xhi = x + (xradius/xstep) * xstep;
01282
01283 while (xlo < xstart) xlo += xstep;
01284 while (xhi > xend ) xhi -= xstep;
01285
01286
01287 int k = 0;
01288 int j, i;
01289 for (j = ylo; j <= yhi; j += ystep) {
01290 for (i = xlo; i <= xhi; i += xstep) {
01291 data[k++] = input_data[ (i-1) + (j-1)*nx ];
01292 }
01293 }
01294
01295
01296 smooth_data[ (x-1) + (y-1)*nx ] =
01297 fors_tools_get_median_float(data, k);
01298 }
01299 }
01300
01301 cpl_free(data);
01302 return smooth;
01303 }
01304
01305 #undef cleanup
01306 #define cleanup \
01307 do { \
01308 cpl_image_delete(input); \
01309 } while(0)
01310 cpl_image *
01311 fors_image_flat_fit_create(fors_image *image,
01312 int step,
01313 int degree,
01314 float level)
01315 {
01316 cpl_image *temp = NULL;
01317 cpl_image *input = NULL;
01318 cpl_image *smooth = NULL;
01319 int nx, ny;
01320
01321 assure( image != NULL, return smooth, NULL );
01322 passure( image->data != NULL, return smooth );
01323 assure( step > 0, return smooth, NULL );
01324 assure( degree >= 0, return smooth, NULL );
01325
01326
01327 temp = image->data;
01328
01329 nx = cpl_image_get_size_x(temp);
01330 ny = cpl_image_get_size_y(temp);
01331
01332
01333
01334
01335
01336 assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return smooth, NULL );
01337
01338
01339
01340
01341
01342 cpl_matrix *kernel = cpl_matrix_new(3, 3);
01343 cpl_matrix_fill(kernel, 1.0);
01344 input = cpl_image_filter_median(image->data, kernel);
01345 cpl_matrix_delete(kernel);
01346 const float *input_data = cpl_image_get_data_float_const(input);
01347
01348
01349
01350
01351
01352 int x, y, pos;
01353 int count = 0;
01354 for (y = 0; y < ny; y += step) {
01355 pos = y*nx;
01356 for (x = 0; x < nx; x += step, pos += step) {
01357 if (input_data[pos] > level) {
01358 count++;
01359 }
01360 }
01361 }
01362
01363 if (count < (degree+1)*(degree+2)) {
01364 step = sqrt((nx*nx)/((degree+1)*(degree+2))) / 2;
01365 if (step == 0)
01366 step = 1;
01367 cpl_msg_error(cpl_func, "Flat field image too small (%dx%d). "
01368 "Please provide a smaller resampling step (a good "
01369 "value would be %d)", nx, ny, step);
01370 cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT);
01371 cleanup;
01372 return NULL;
01373 }
01374
01375
01376
01377
01378
01379
01380 cpl_bivector *positions = cpl_bivector_new(count);
01381 double *xpos = cpl_bivector_get_x_data(positions);
01382 double *ypos = cpl_bivector_get_y_data(positions);
01383 cpl_vector *fluxes = cpl_vector_new(count);
01384 double *flux = cpl_vector_get_data(fluxes);
01385
01386 count = 0;
01387 for (y = 0; y < ny; y += step) {
01388 pos = y*nx;
01389 for (x = 0; x < nx; x += step, pos += step) {
01390 if (input_data[pos] > level) {
01391 xpos[count] = x;
01392 ypos[count] = y;
01393 flux[count] = input_data[pos];
01394 count++;
01395 }
01396 }
01397 }
01398
01399 cpl_image_delete(input); input = NULL;
01400
01401
01402
01403
01404
01405
01406 cpl_polynomial *model = cpl_polynomial_fit_2d_create(positions,
01407 fluxes,
01408 degree,
01409 NULL);
01410
01411 cpl_bivector_delete(positions);
01412 cpl_vector_delete(fluxes);
01413
01414 smooth = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
01415 float *smooth_data = cpl_image_get_data_float(smooth);
01416
01417 cpl_vector *point = cpl_vector_new(2);
01418 double *dpoint = cpl_vector_get_data(point);
01419
01420 for (y = 0; y < ny; y++) {
01421 pos = y*nx;
01422 dpoint[1] = y;
01423 for (x = 0; x < nx; x++, pos++) {
01424 dpoint[0] = x;
01425 smooth_data[pos] = cpl_polynomial_eval(model, point);
01426 }
01427 }
01428
01429 cpl_polynomial_delete(model);
01430 cpl_vector_delete(point);
01431
01432 cleanup;
01433 return smooth;
01434
01435 }
01436
01437 #undef cleanup
01438 #define cleanup
01439
01455 cpl_image *
01456 fors_image_filter_max_create(const fors_image *image,
01457 int xradius,
01458 int yradius,
01459 bool use_data)
01460 {
01461 const cpl_image *input = NULL;
01462 cpl_image *hmaxima = NULL;
01463 cpl_image *maxima = NULL;
01464 int nx, ny;
01465
01466 assure( image != NULL, return maxima, NULL );
01467 passure( image->data != NULL, return maxima );
01468 passure( image->variance != NULL, return maxima );
01469
01470 input = (use_data) ? image->data : image->variance;
01471
01472 nx = cpl_image_get_size_x(input);
01473 ny = cpl_image_get_size_y(input);
01474
01475
01476
01477
01478
01479 hmaxima = cpl_image_duplicate(input);
01480
01481
01482 assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return maxima, NULL );
01483
01484 float *input_data = (float *)cpl_image_get_data_float_const(input);
01485 float *maxima_data = cpl_image_get_data_float(hmaxima);
01486
01487 int y;
01488 for (y = 0; y < ny; y++) {
01489 const float *irow = input_data + y * nx;
01490 float *orow = maxima_data + y * nx;
01491 max_filter(irow, orow, nx, 2*xradius+1);
01492 }
01493
01494 cpl_image_turn(hmaxima, 1);
01495
01496
01497
01498
01499
01500 maxima = cpl_image_duplicate(hmaxima);
01501 input_data = cpl_image_get_data_float(hmaxima);
01502 maxima_data = cpl_image_get_data_float(maxima);
01503
01504
01505
01506
01507
01508 int x;
01509 for (x = 0; x < nx; x++) {
01510 const float *irow = input_data + x * ny;
01511 float *orow = maxima_data + x * ny;
01512 max_filter(irow, orow, ny, 2*yradius+1);
01513 }
01514
01515 cpl_image_delete(hmaxima);
01516
01517 cpl_image_turn(maxima, -1);
01518
01519 return maxima;
01520 }
01521
01522 #undef cleanup
01523 #define cleanup
01524
01530 double
01531 fors_image_get_stdev(const fors_image *image, double *dstdev)
01532 {
01533 assure( image != NULL, return 0, NULL );
01534 assure( dstdev == NULL, return 0, "Unsupported");
01535
01536 return cpl_image_get_stdev(image->data);
01537 }
01538 #undef cleanup
01539 #define cleanup \
01540 do { \
01541 cpl_mask_delete(rejected); \
01542 cpl_image_delete(im); \
01543 } while (0)
01544
01552 double fors_image_get_stdev_robust(const fors_image *image,
01553 double cut,
01554 double *dstdev)
01555 {
01556 cpl_mask *rejected = NULL;
01557 cpl_image *im = NULL;
01558
01559 assure( image != NULL, return 0, NULL );
01560 assure( cut > 0, return 0, "Illegal cut: %f", cut );
01561 assure( dstdev == NULL, return 0, "Unsupported");
01562
01563 double median = fors_image_get_median(image, NULL);
01564
01565 im = cpl_image_duplicate(image->data);
01566 cpl_image_subtract_scalar(im, median);
01567 cpl_image_power(im, 2);
01568
01569
01570 rejected = cpl_mask_threshold_image_create(image->data,
01571 median - cut,
01572 median + cut);
01573 cpl_mask_not(rejected);
01574 cpl_image_reject_from_mask(im, rejected);
01575
01576 double robust_stdev = sqrt(cpl_image_get_mean(im));
01577
01578 cleanup;
01579 return robust_stdev;
01580 }
01581
01582 #undef cleanup
01583 #define cleanup
01584
01594 double
01595 fors_image_get_error_mean(const fors_image *image, double *dmean)
01596 {
01597 double avg;
01598
01599 assure( image != NULL, return 0, NULL );
01600 assure( dmean == NULL, return 0, "Unsupported");
01601
01602 avg = cpl_image_get_mean(image->variance);
01603
01604
01605 assure( avg >= 0, return -1, "Average variance is %f", avg);
01606
01607 return sqrt(avg);
01608 }
01609
01610
01611 #undef cleanup
01612 #define cleanup \
01613 do { \
01614 cpl_imagelist_delete(datlist); \
01615 cpl_imagelist_delete(varlist); \
01616 } while (0)
01617
01626 fors_image *
01627 fors_image_collapse_create(const fors_image_list *images)
01628 {
01629 cpl_imagelist *datlist = NULL;
01630 cpl_imagelist *varlist = NULL;
01631 cpl_image *data = NULL;
01632 cpl_image *variance = NULL;
01633 const fors_image *i;
01634 int N = 0;
01635
01636 assure( images != NULL, return NULL, NULL );
01637 assure( fors_image_list_size(images) > 0, return NULL,
01638 "Cannot stack zero images");
01639
01640 i = fors_image_list_first_const(images);
01641
01642 datlist = cpl_imagelist_new();
01643 varlist = cpl_imagelist_new();
01644
01645 while(i != NULL) {
01646
01647
01648 cpl_imagelist_set(datlist,
01649 cpl_image_duplicate(i->data),
01650 cpl_imagelist_get_size(datlist));
01651 cpl_imagelist_set(varlist,
01652 cpl_image_duplicate(i->variance),
01653 cpl_imagelist_get_size(varlist));
01654 i = fors_image_list_next_const(images);
01655 N++;
01656 }
01657
01658 #ifdef CPL_IS_NOT_CRAP
01659 data = cpl_imagelist_collapse_create(datlist);
01660
01661 variance = cpl_imagelist_collapse_create(varlist);
01662 #else
01663 data = fors_imagelist_collapse_create(datlist);
01664
01665 variance = fors_imagelist_collapse_create(varlist);
01666 #endif
01667
01668 cpl_image_divide_scalar(variance, N);
01669
01670 cleanup;
01671 return fors_image_new(data, variance);
01672 }
01673
01674
01675 #undef cleanup
01676 #define cleanup \
01677 do { \
01678 cpl_imagelist_delete(datlist); \
01679 cpl_imagelist_delete(varlist); \
01680 } while (0)
01681
01692 fors_image *
01693 fors_image_collapse_minmax_create(const fors_image_list *images,
01694 int low, int high)
01695 {
01696 cpl_imagelist *datlist = NULL;
01697 cpl_imagelist *varlist = NULL;
01698 cpl_image *data = NULL;
01699 cpl_image *variance = NULL;
01700 const fors_image *i;
01701 int N = 0;
01702
01703 assure( images != NULL, return NULL, NULL );
01704 assure( fors_image_list_size(images) > low + high, return NULL,
01705 "Cannot reject more images than there are");
01706 assure( low*high >= 0 && low+high > 0, return NULL,
01707 "Invalid minmax rejection criteria");
01708
01709 i = fors_image_list_first_const(images);
01710
01711 datlist = cpl_imagelist_new();
01712 varlist = cpl_imagelist_new();
01713
01714 while(i != NULL) {
01715
01716
01717 cpl_imagelist_set(datlist,
01718 cpl_image_duplicate(i->data),
01719 cpl_imagelist_get_size(datlist));
01720 cpl_imagelist_set(varlist,
01721 cpl_image_duplicate(i->variance),
01722 cpl_imagelist_get_size(varlist));
01723 i = fors_image_list_next_const(images);
01724 N++;
01725 }
01726
01727 data = cpl_imagelist_collapse_minmax_create(datlist, low, high);
01728 variance = cpl_imagelist_collapse_minmax_create(varlist, low, high);
01729
01730 cpl_image_divide_scalar(variance, N);
01731
01732 cleanup;
01733 return fors_image_new(data, variance);
01734 }
01735
01747 fors_image *
01748 fors_image_collapse_median_create(const fors_image_list *images)
01749 {
01750 cpl_imagelist *datlist = NULL;
01751 cpl_imagelist *varlist = NULL;
01752 cpl_image *data = NULL;
01753 cpl_image *variance = NULL;
01754 const fors_image *i;
01755 int N = 0;
01756
01757 assure( images != NULL, return NULL, NULL );
01758 assure( fors_image_list_size(images) > 0, return NULL,
01759 "Cannot stack zero images");
01760
01761 i = fors_image_list_first_const(images);
01762
01763 datlist = cpl_imagelist_new();
01764 varlist = cpl_imagelist_new();
01765 while(i != NULL) {
01766
01767 cpl_imagelist_set(datlist,
01768 cpl_image_duplicate(i->data),
01769 cpl_imagelist_get_size(datlist));
01770 cpl_imagelist_set(varlist,
01771 cpl_image_duplicate(i->variance),
01772 cpl_imagelist_get_size(varlist));
01773
01774 i = fors_image_list_next_const(images);
01775 N++;
01776 }
01777
01778 #ifdef CPL_IS_NOT_CRAP
01779 data = cpl_imagelist_collapse_median_create(datlist);
01780
01781 variance = cpl_imagelist_collapse_create(varlist);
01782 #else
01783 data = fors_imagelist_collapse_median_create(datlist);
01784
01785 variance = fors_imagelist_collapse_create(varlist);
01786 #endif
01787
01788 cpl_image_divide_scalar(variance, N);
01789
01790 cpl_image_multiply_scalar(variance,
01791 fors_utils_median_corr(N) *
01792 fors_utils_median_corr(N));
01793
01794 cleanup;
01795 return fors_image_new(data, variance);
01796 }
01797
01798 #undef cleanup
01799 #define cleanup
01800
01817 void fors_image_draw(fors_image *image, int type,
01818 double x, double y,
01819 int radius, double color)
01820 {
01821 assure( image != NULL, return, NULL );
01822
01823 assure( type == 0 || type == 1 || type == 2,
01824 return , "Unsupported type %d", type);
01825
01826 assure( radius > 0, return, NULL );
01827
01828 if (type == 2) {
01829 int i;
01830 for (i = 0; i < 360; i++) {
01831
01832
01833 int px = x + radius*cos(i/(2*M_PI));
01834 int py = y + radius*sin(i/(2*M_PI));
01835
01836 if (1 <= px && px <= cpl_image_get_size_x(image->data) &&
01837 1 <= py && py <= cpl_image_get_size_y(image->data)) {
01838 cpl_image_set(image->data, px, py, color);
01839 cpl_image_set(image->variance, px, py, color > 0 ? color : 0);
01840 }
01841 }
01842 }
01843 else {
01844 int i;
01845
01846 for (i = -radius; i <= radius; i++) {
01847
01848 int px, py;
01849
01850 if (type == 0) {
01851 px = x + i;
01852 py = y;
01853 }
01854 else {
01855 px = x;
01856 py = y + i;
01857 }
01858
01859 if (1 <= px && px <= cpl_image_get_size_x(image->data) &&
01860 1 <= py && py <= cpl_image_get_size_y(image->data)) {
01861 cpl_image_set(image->data , px, py, color);
01862 cpl_image_set(image->variance, px, py, color > 0 ? color : 0);
01863 }
01864 }
01865 }
01866
01867 return;
01868 }
01869
01870 #define LIST_DEFINE
01871 #undef LIST_ELEM
01872 #define LIST_ELEM fors_image
01873 #include <list.h>
01874