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
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199 #ifdef HAVE_CONFIG_H
00200 # include <config.h>
00201 #endif
00202 #include "sinfo_vltPort.h"
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212 #include "sinfo_image_ops.h"
00213 #include "sinfo_resampling.h"
00214 #include "sinfo_local_types.h"
00223
00224
00225
00226
00235 double sinfo_new_my_median_image(cpl_image* im)
00236 {
00237 double m=0;
00238 register int i=0;
00239 int n=0;
00240 pixelvalue* pv=0;
00241 int ilx=0;
00242 int ily=0;
00243 float* pidata=NULL;
00244
00245
00246 if(im==NULL) sinfo_msg_error("Null Image");
00247 ilx=cpl_image_get_size_x(im);
00248 ily=cpl_image_get_size_y(im);
00249 pidata=cpl_image_get_data_float(im);
00250
00251 for ( i = 0 ; i < (int) ilx*ily ; i++ )
00252 {
00253 if ( isnan(pidata[i]) )
00254 {
00255
00256 } else {
00257 n++;
00258 }
00259 }
00260 pv = cpl_calloc(n,sizeof(pixelvalue));
00261 n=0;
00262 for ( i = 0 ; i < (int) ilx*ily ; i++ )
00263 {
00264 if ( isnan(pidata[i]) )
00265 {
00266
00267 } else {
00268 pv[n]=pidata[i];
00269 n++;
00270 }
00271 }
00272 if(pv == NULL || n == 0) {
00273 m=0;
00274 } else {
00275 m=sinfo_new_median(pv,n);
00276 }
00277 cpl_free(pv);
00278 if(isnan(m)){
00279 m=0;
00280 }
00281 return m;
00282 }
00283
00292 Vector * sinfo_new_mean_of_columns( cpl_image *im )
00293 {
00294 Vector * row=NULL ;
00295 int i=0;
00296 int j=0;
00297 int ilx=0;
00298 int ily=0;
00299 float* pidata=NULL;
00300
00301 if ( im == NULL )
00302 {
00303 sinfo_msg_error ("null image") ;
00304 return NullVector ;
00305 }
00306 ilx=cpl_image_get_size_x(im);
00307 ily=cpl_image_get_size_y(im);
00308 pidata=cpl_image_get_data_float(im);
00309
00310
00311 if ( NULL == (row = sinfo_new_vector (ilx)) )
00312 {
00313 sinfo_msg_error ("not able to allocate a sinfo_vector" ) ;
00314 return NullVector ;
00315 }
00316
00317 for ( i = 0 ; i < ilx ; i++ )
00318 {
00319 for ( j = 0 ; j < ily ; j++ )
00320 {
00321 if (!isnan(pidata[i+j*ilx]))
00322 {
00323 row->data[i] += pidata[i + j*(ilx)] ;
00324 }
00325 }
00326
00327 row->data[i] /= ily ;
00328 }
00329 return row ;
00330 }
00342 cpl_image * sinfo_new_clean_mean_of_columns( cpl_image *im,
00343 float lo_reject,
00344 float hi_reject)
00345 {
00346 cpl_image * row=NULL ;
00347 pixelvalue* buffer=NULL ;
00348 int i=0;
00349 int j=0;
00350 int k=0;
00351 int nv=0;
00352 int lo_n=0;
00353 int hi_n=0;
00354 int ilx=0;
00355 int ily=0;
00356 float* pidata=NULL;
00357 float* podata=NULL;
00358
00359 if ( im == NULL )
00360 {
00361 sinfo_msg_error ("null image") ;
00362 return NULL ;
00363 }
00364 ilx=cpl_image_get_size_x(im);
00365 ily=cpl_image_get_size_y(im);
00366 pidata=cpl_image_get_data_float(im);
00367
00368 if ((lo_reject + hi_reject) > 0.9)
00369 {
00370 sinfo_msg_error("illegal rejection thresholds: [%f] and [%f]",
00371 lo_reject, hi_reject) ;
00372 sinfo_msg_error("threshold sum should not be over "
00373 "0.90 aborting average") ;
00374 return NULL ;
00375 }
00376
00377 lo_n = (int) (ily * lo_reject + 0.5) ;
00378 hi_n = (int) (ily * hi_reject + 0.5) ;
00379 if (lo_n + hi_n >= ily)
00380 {
00381 sinfo_msg_error ("everything would be rejected") ;
00382 return NULL ;
00383 }
00384
00385
00386 if ( NULL == (row = cpl_image_new (ilx, 1,CPL_TYPE_FLOAT)) )
00387 {
00388 sinfo_msg_error ("cannot allocate new image") ;
00389 return NULL ;
00390 }
00391 podata=cpl_image_get_data_float(row);
00392
00393 buffer=(pixelvalue*) cpl_calloc(ily,sizeof(pixelvalue)) ;
00394
00395 for ( i = 0 ; i < ilx ; i++ )
00396 {
00397 for ( j = 0 ; j < ily ; j++ )
00398 {
00399 buffer[j] = pidata[i + j*(ilx)] ;
00400 }
00401 sinfo_pixel_qsort (buffer, ily) ;
00402
00403 nv = 0 ;
00404 for (k = lo_n ; k < ily - hi_n ; k ++)
00405 {
00406 if ( !isnan(buffer[k]) )
00407 {
00408 podata[i] += buffer[k] ;
00409 nv ++ ;
00410 }
00411 }
00412 podata[i] /= nv ;
00413
00414 }
00415 cpl_free(buffer);
00416 return row ;
00417 }
00418
00419
00429 cpl_image * sinfo_new_div_image_by_row( cpl_image *im, Vector *row )
00430 {
00431 cpl_image *image=NULL ;
00432 int i=0;
00433 int j=0;
00434 int ilx=0;
00435 int ily=0;
00436 float* pidata=NULL;
00437 float* podata=NULL;
00438
00439 if ( im == NULL || row == NULL )
00440 {
00441 sinfo_msg_error ("null image or null row") ;
00442 return NULL ;
00443 }
00444 ilx=cpl_image_get_size_x(im);
00445 ily=cpl_image_get_size_y(im);
00446 pidata=cpl_image_get_data_float(im);
00447
00448 if ( ilx != row -> n_elements )
00449 {
00450 sinfo_msg_error("image and row size not compatible") ;
00451 return NULL ;
00452 }
00453
00454 if ( NULL == (image = cpl_image_duplicate (im)) )
00455 {
00456 sinfo_msg_error ("cannot copy image") ;
00457 return NULL ;
00458 }
00459 podata=cpl_image_get_data_float(image);
00460
00461 for (i = 0 ; i < ilx ; i++ )
00462 {
00463 for (j = 0 ; j < ily ; j++)
00464 {
00465 if ( !isnan(pidata[i + j*(ilx)]) )
00466 {
00467 podata[i + j*(ilx)] = pidata[i + j*(ilx)] / row -> data[i] ;
00468 }
00469 }
00470 }
00471 return image ;
00472 }
00473
00474
00484 cpl_image * sinfo_new_mult_row_to_image( cpl_image *im, Vector *row )
00485 {
00486 cpl_image *image=NULL;
00487 int i=0;
00488 int j=0;
00489 int ilx=0;
00490 int ily=0;
00491 float* pidata=NULL;
00492 float* podata=NULL;
00493
00494
00495
00496
00497 if ( im == NULL || row == NULL )
00498 {
00499 sinfo_msg_error ("null image or null row") ;
00500 return NULL ;
00501 }
00502 ilx=cpl_image_get_size_x(im);
00503 ily=cpl_image_get_size_y(im);
00504 pidata=cpl_image_get_data_float(im);
00505
00506 if ( ilx != row -> n_elements )
00507 {
00508 sinfo_msg_error("image and row size not compatible") ;
00509 return NULL ;
00510 }
00511
00512 if ( NULL == (image = cpl_image_duplicate (im)) )
00513 {
00514 sinfo_msg_error ("cannot copy image") ;
00515 return NULL ;
00516 }
00517 podata=cpl_image_get_data_float(image);
00518
00519 for (i = 0 ; i < ilx ; i++ )
00520 {
00521 for (j = 0 ; j < ily ; j++)
00522 {
00523 if ( !isnan(pidata[i + j*(ilx)]) )
00524 {
00525 podata[i + j*(ilx)] = pidata[i + j*(ilx)] * row -> data[i] ;
00526 }
00527 }
00528 }
00529 return image ;
00530 }
00531
00532
00533
00534
00535
00536
00560 cpl_image * sinfo_new_col_tilt ( cpl_image * image, float sigmaFactor )
00561 {
00562 cpl_image * im=NULL;
00563 float * column=NULL ;
00564 double sum=0;
00565 double sum2=0;
00566 double mean=0;
00567 float sinfo_median=0;
00568 float noise=0 ;
00569 float * sig=NULL;
00570 float * dat=NULL;
00571 float a=0;
00572 float b=0;
00573 float siga=0;
00574 float sigb=0;
00575 float chi2=0;
00576 float q=0;
00577 int i=0;
00578 int j=0;
00579 int colnum=0;
00580 int npix=0;
00581 int mwt=0 ;
00582 int lx=0;
00583 int ly=0;
00584 float* p_in_data=NULL;
00585 float* p_ou_data=NULL;
00586
00587
00588 if ( image == NULL )
00589 {
00590 sinfo_msg_error ("no image given" ) ;
00591 return NULL ;
00592 }
00593
00594 if ( sigmaFactor <= 0. )
00595 {
00596 sinfo_msg_error ("no or negative sigma factor") ;
00597 return NULL ;
00598 }
00599 lx = cpl_image_get_size_x(image);
00600 ly = cpl_image_get_size_y(image);
00601
00602
00603
00604 if ( NULL == (im = cpl_image_new (lx,ly,CPL_TYPE_FLOAT )) )
00605 {
00606 sinfo_msg_error ("cannot allocate new image" ) ;
00607 return NULL ;
00608 }
00609
00610
00611 p_in_data = cpl_image_get_data_float(image);
00612 p_ou_data = cpl_image_get_data_float(im);
00613 for ( i = 0 ; i < lx ; i ++ )
00614 {
00615
00616 colnum = 0 ;
00617 column = (float *) cpl_calloc ( ly , sizeof (float *) ) ;
00618 sig = (float *) cpl_calloc ( ly , sizeof (float *) ) ;
00619 dat = (float *) cpl_calloc ( ly , sizeof (float *) ) ;
00620
00621
00622 for ( j = 0 ; j < ly ; j++ )
00623 {
00624 if ( !isnan(p_in_data[i + j*lx]) )
00625 {
00626 column[j] = p_in_data[i + j*lx] ;
00627 colnum ++ ;
00628 }
00629 }
00630 if ( colnum < 10 )
00631 {
00632
00633
00634 for ( j = 0 ; j < ly ; j++ )
00635 {
00636 p_ou_data[i + j*lx] = ZERO;
00637 }
00638
00639
00640
00641
00642
00643
00644 }
00645
00646
00647
00648
00649
00650
00651
00652 sinfo_pixel_qsort (column, colnum) ;
00653
00654 sum = 0. ;
00655 sum2 = 0. ;
00656 npix = 0 ;
00657
00658 for ( j = 0.1*colnum + 1 ; j <= 0.9*colnum ; j++ )
00659 {
00660 sum += column[j] ;
00661 sum2 += column[j] * column[j] ;
00662 npix ++ ;
00663 }
00664
00665 if (npix <= 1)
00666 {
00667 noise = sigmaFactor * 1000.;
00668 }
00669 else
00670 {
00671 mean = sum/(float)npix ;
00672 noise = sqrt( (sum2 - sum*mean)/(double)(npix -1) ) ;
00673 noise *= sigmaFactor ;
00674 }
00675
00676
00677
00678
00679
00680
00681
00682 if ( colnum % 2 == 1 )
00683 {
00684 sinfo_median = column[colnum/2] ;
00685 }
00686 else
00687 {
00688 sinfo_median = (column[colnum/2 - 1] + column[colnum/2])/2. ;
00689 }
00690
00691
00692
00693 colnum = 0 ;
00694 for ( j = 0; j < ly ; j ++ )
00695 {
00696 if ( !isnan(p_in_data[i+j*lx]) &&
00697 fabs ( (p_in_data[i+j*lx]) - sinfo_median) <= noise )
00698 {
00699 column[colnum] = p_in_data[i+j*lx] ;
00700 dat[colnum] = (float) j ;
00701 sig[colnum] = 1. ;
00702 colnum ++ ;
00703 }
00704 }
00705
00706 if ( colnum == 0 )
00707 {
00708
00709
00710
00711
00712
00713
00714
00715
00716 a=0./0.;
00717 b=0./0.;
00718 }
00719 else
00720 {
00721 mwt = 0 ;
00722 sinfo_my_fit ( dat, column, colnum, sig, mwt, &a,
00723 &b, &siga, &sigb, &chi2, &q ) ;
00724 }
00725 if ( fabs(b) >= SLOPE || fabs(a) >= SATURATION ||
00726 isnan(b) || isnan(a))
00727 {
00728 sinfo_msg_warning ("linear fit: slope is greater than limit: %f"
00729 " saturation level is reached: %f in column"
00730 " number %d ", b, a , i+1) ;
00731 }
00732
00733
00734 for ( j = 0; j < ly; j++ )
00735 {
00736 if ( !isnan(p_in_data[i+j*lx]) &&
00737 fabs(b) < SLOPE && fabs(a) < SATURATION )
00738 {
00739 p_ou_data[i+j*lx] = p_in_data[i+j*lx] - (a + b * (float)j) ;
00740 }
00741 else if ( isnan(p_in_data[i+j*lx]) )
00742 {
00743 p_ou_data[i+j*lx] = ZERO ;
00744 }
00745 else if ( (fabs(b) >= SLOPE ||
00746 fabs(a) >= SATURATION || isnan(a) || isnan(b)) &&
00747 !isnan(p_in_data[i+j*lx]) )
00748 {
00749 p_ou_data[i+j*lx] -= sinfo_median ;
00750 }
00751 else
00752 {
00753 sinfo_msg_error (" case is not possible! %f %f", b,a) ;
00754
00755
00756
00757
00758
00759 }
00760 }
00761 cpl_free (column) ;
00762 cpl_free (sig) ;
00763 cpl_free (dat) ;
00764 }
00765
00766 return im ;
00767 }
00768
00769
00770
00771
00772
00796 cpl_image * sinfo_new_median_image( cpl_image * im, float fmedian )
00797 {
00798 cpl_image * image=NULL ;
00799 pixelvalue * value=NULL ;
00800 pixelvalue sinfo_median=0 ;
00801 int * position=NULL ;
00802 int nposition=0 ;
00803 int n=0;
00804 int i=0;
00805 int j=0;
00806 int lx=0;
00807 int ly=0;
00808 float* p_in_data=NULL;
00809 float* p_ou_data=NULL;
00810
00811 if ( im == NULL )
00812 {
00813 sinfo_msg_error ("no image input") ;
00814 return NULL ;
00815 }
00816
00817 image = cpl_image_duplicate ( im ) ;
00818 lx=cpl_image_get_size_x(im);
00819 ly=cpl_image_get_size_y(im);
00820 p_in_data=cpl_image_get_data_float(im);
00821 p_ou_data=cpl_image_get_data_float(image);
00822
00823
00824
00825
00826
00827 for ( i = 0 ; i < (int) lx*ly ; i++ )
00828 {
00829
00830 if ( isnan(p_in_data[i]) )
00831 {
00832 continue ;
00833 }
00834
00835
00836 value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ;
00837 position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ;
00838
00839
00840
00841
00842
00843 position[0] = i + lx - 1 ;
00844 position[1] = i + lx ;
00845 position[2] = i + lx + 1 ;
00846 position[3] = i + 1 ;
00847 position[4] = i - lx + 1 ;
00848 position[5] = i - lx ;
00849 position[6] = i - lx - 1 ;
00850 position[7] = i - 1 ;
00851
00852
00853
00854
00855
00856
00857
00858 if ( i >= 0 && i < lx )
00859 {
00860 position[4] += 2 * lx ;
00861 position[5] += 2 * lx ;
00862 position[6] += 2 * lx ;
00863 }
00864 else if ( i >= ((int) lx*ly - lx ) && i < (int) lx*ly )
00865 {
00866 position[0] -= 2 * lx ;
00867 position[1] -= 2 * lx ;
00868 position[2] -= 2 * lx ;
00869 }
00870 else if ( i % lx == 0 )
00871 {
00872 position[0] += 2 ;
00873 position[6] += 2 ;
00874 position[7] += 2 ;
00875 }
00876 else if ( i % lx == lx - 1 )
00877 {
00878 position[2] -= 2 ;
00879 position[3] -= 2 ;
00880 position[4] -= 2 ;
00881 }
00882
00883
00884
00885
00886
00887
00888 nposition = 8 ;
00889 n = 0 ;
00890 for ( j = 0 ; j < nposition ; j ++ )
00891 {
00892 if ( !isnan(p_in_data[position[j]]) )
00893 {
00894 value[n] = p_in_data[position[j]] ;
00895 n ++ ;
00896 }
00897 }
00898 nposition = n ;
00899
00900 if ( nposition <= 1 )
00901 {
00902 p_ou_data[i] = ZERO ;
00903 cpl_free(value) ;
00904 cpl_free(position) ;
00905 continue ;
00906 }
00907
00908
00909
00910 sinfo_pixel_qsort ( value, nposition ) ;
00911 if ( nposition % 2 == 1 )
00912 {
00913 sinfo_median = value [ nposition/2 ] ;
00914 }
00915 else
00916 {
00917 sinfo_median = ( value [nposition/2 - 1] +
00918 value [nposition/2] ) / 2. ;
00919 }
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935 if ( fmedian == 0 )
00936 {
00937 p_ou_data[i] = sinfo_median ;
00938 }
00939 else if ( fmedian < 0 &&
00940 fabs ( sinfo_median - p_in_data[i] ) >= -fmedian )
00941 {
00942 p_ou_data[i] = sinfo_median ;
00943 }
00944 else if ( fmedian > 0 &&
00945 fabs ( sinfo_median - p_in_data[i] ) >= fmedian *
00946 sqrt(fabs(sinfo_median)) )
00947 {
00948 p_ou_data[i] = sinfo_median ;
00949 }
00950 else
00951 {
00952 cpl_free (value) ;
00953 cpl_free (position) ;
00954 continue ;
00955 }
00956
00957 cpl_free (value) ;
00958 cpl_free (position) ;
00959 }
00960 return image ;
00961 }
00962
00963
00964
00965
00976 cpl_image *
00977 sinfo_new_compare_images(cpl_image * im1,cpl_image * im2,cpl_image * origim )
00978 {
00979 cpl_image * image=NULL ;
00980 int i=0 ;
00981 int lx1=0;
00982 int ly1=0;
00983 int lx2=0;
00984 int ly2=0;
00985 float* p_in1_data=NULL;
00986 float* p_in2_data=NULL;
00987 float* p_ou_data=NULL;
00988 float* p_org_data=NULL;
00989
00990
00991 if ( im1 == NULL || im2 == NULL || origim == NULL )
00992 {
00993 sinfo_msg_error ("Null images as input" ) ;
00994 return NULL ;
00995 }
00996 lx1=cpl_image_get_size_x(im1);
00997 ly1=cpl_image_get_size_y(im1);
00998
00999 lx2=cpl_image_get_size_x(im2);
01000 ly2=cpl_image_get_size_y(im2);
01001
01002 p_in1_data=cpl_image_get_data_float(im1);
01003 p_in2_data=cpl_image_get_data_float(im2);
01004 p_org_data=cpl_image_get_data_float(origim);
01005 if ( lx1 != lx2 || ly1 != ly2 )
01006 {
01007 sinfo_msg_error ("incompatible image sizes" ) ;
01008 return NULL ;
01009 }
01010
01011
01012 if ( NULL == (image = cpl_image_new ( lx1, ly1, CPL_TYPE_FLOAT )) )
01013 {
01014 sinfo_msg_error ("cannot allocate new image" ) ;
01015 return NULL ;
01016 }
01017 p_ou_data=cpl_image_get_data_float(image);
01018 for ( i = 0 ; i < (int) lx1*ly1 ; i ++ )
01019 {
01020 if ( isnan(p_in1_data[i]) && isnan(p_in2_data[i]) )
01021 {
01022 p_ou_data[i] = ZERO ;
01023 }
01024 else
01025 {
01026 if ( p_in1_data[i] == p_in2_data[i] )
01027 {
01028 p_ou_data[i] = p_org_data[i] ;
01029 }
01030 else
01031 {
01032 p_ou_data[i] = ZERO ;
01033 }
01034 }
01035 }
01036 return image ;
01037 }
01038
01039
01040
01052 cpl_image *
01053 sinfo_new_promote_image_to_mask (cpl_image * im, int * n_badpixels )
01054 {
01055 cpl_image * reImage=NULL ;
01056 int i=0 ;
01057 int lx=0;
01058 int ly=0;
01059 float* p_in_data=NULL;
01060 float* p_ou_data=NULL;
01061
01062 if ( NULL == im )
01063 {
01064 sinfo_msg_error("no input image given!") ;
01065 return NULL ;
01066 }
01067 lx=cpl_image_get_size_x(im);
01068 ly=cpl_image_get_size_y(im);
01069 p_in_data=cpl_image_get_data_float(im);
01070
01071
01072 if ( NULL == (reImage = cpl_image_new (lx,ly,CPL_TYPE_FLOAT )) )
01073 {
01074 sinfo_msg_error ("cannot allocate new image!") ;
01075 return NULL ;
01076 }
01077 p_ou_data=cpl_image_get_data_float(reImage);
01078
01079 *n_badpixels = 0 ;
01080 for ( i = 0 ; i < (int) lx*ly ; i ++ )
01081 {
01082 if ( isnan(p_in_data[i]) )
01083 {
01084 p_ou_data[i] = 0. ;
01085 (*n_badpixels)++ ;
01086 }
01087 else
01088 {
01089 p_ou_data[i] = 1. ;
01090 }
01091 }
01092 return reImage ;
01093 }
01094
01095
01106 cpl_image * sinfo_new_mult_image_by_mask (cpl_image * im,cpl_image * mask )
01107 {
01108 cpl_image * reImage=NULL ;
01109 int i=0 ;
01110 int ix=0;
01111 int iy=0;
01112 int mx=0;
01113 int my=0;
01114
01115
01116 float* pmdata=NULL;
01117 float* podata=NULL;
01118
01119 if ( NULL == im )
01120 {
01121 sinfo_msg_error("no input image given!") ;
01122 return NULL ;
01123 }
01124 if ( NULL == mask )
01125 {
01126 sinfo_msg_error("no mask image given!") ;
01127 return NULL ;
01128 }
01129 ix=cpl_image_get_size_x(im);
01130 iy=cpl_image_get_size_y(im);
01131 mx=cpl_image_get_size_x(mask);
01132 my=cpl_image_get_size_y(mask);
01133
01134 if ( ix != mx || iy != my)
01135 {
01136 sinfo_msg_error("image sizes are not correspondent!") ;
01137 return NULL ;
01138 }
01139
01140 reImage = cpl_image_duplicate( im ) ;
01141 podata=cpl_image_get_data_float(reImage);
01142 pmdata=cpl_image_get_data_float(mask);
01143
01144 for ( i = 0 ; i < (int) ix*iy ; i ++ )
01145 {
01146 if ( pmdata[i] == 0. )
01147 {
01148 podata[i] = ZERO ;
01149 }
01150 }
01151
01152 return reImage ;
01153 }
01154
01155
01156
01166 cpl_image *
01167 sinfo_new_thresh_image (cpl_image * im, float lo_cut, float hi_cut )
01168 {
01169 cpl_image * image=NULL ;
01170 float* p_inp_data=NULL;
01171 float* p_out_data=NULL;
01172 int lx=0;
01173 int ly=0;
01174
01175 int i=0 ;
01176
01177 if (im == NULL)
01178 {
01179 sinfo_msg_error ("null image given") ;
01180 return NULL ;
01181 }
01182 lx=cpl_image_get_size_x(im);
01183 ly=cpl_image_get_size_y(im);
01184
01185 image = cpl_image_duplicate(im) ;
01186 p_inp_data=cpl_image_get_data(im);
01187 p_out_data=cpl_image_get_data(image);
01188 for ( i = 0 ; i < (int) lx*ly ; i ++ )
01189 {
01190 if ( p_inp_data[i] > (pixelvalue) hi_cut ||
01191 p_inp_data[i] < (pixelvalue) lo_cut )
01192 {
01193 p_out_data[i] = ZERO ;
01194 }
01195 }
01196 return image ;
01197 }
01198
01199
01200
01201
01226 cpl_image * sinfo_new_interpol_image ( cpl_image * im,
01227 cpl_image * mask,
01228 int max_radius,
01229 int n_pixels )
01230 {
01231 cpl_image * returnImage=NULL ;
01232 float* neighbors=NULL ;
01233 float sum=0;
01234 float mean=0;
01235 int i=0;
01236 int j=0;
01237 int k=0;
01238 int row=0;
01239 int col=0;
01240 int n_valid=0;
01241 int agreed=0;
01242
01243 int ilx=0;
01244 int ily=0;
01245 int mlx=0;
01246 int mly=0;
01247 float* pidata=NULL;
01248 float* podata=NULL;
01249 float* pmdata=NULL;
01250
01251 if ( NULL == im )
01252 {
01253 sinfo_msg_error("sorry, no input image given!") ;
01254 return NULL ;
01255 }
01256 ilx=cpl_image_get_size_x(im);
01257 ily=cpl_image_get_size_y(im);
01258 pidata=cpl_image_get_data_float(im);
01259
01260 if ( NULL == mask )
01261 {
01262 sinfo_msg_error("sorry, no mask image given!") ;
01263 return NULL ;
01264 }
01265
01266 mlx=cpl_image_get_size_x(mask);
01267 mly=cpl_image_get_size_y(mask);
01268 pmdata=cpl_image_get_data_float(mask);
01269
01270 if ( mlx != ilx || mly != mly )
01271 {
01272 sinfo_msg_error("images not compatible !") ;
01273 return NULL ;
01274 }
01275
01276 if ( max_radius <= 0 )
01277 {
01278 sinfo_msg_error("wrong number of pixels for maximal "
01279 "search radius given!") ;
01280 return NULL ;
01281 }
01282
01283 if ( n_pixels <= 2 )
01284 {
01285 sinfo_msg_error("wrong number of pixels used for interpolation given!") ;
01286 return NULL ;
01287 }
01288
01289 returnImage = cpl_image_duplicate ( im ) ;
01290 podata=cpl_image_get_data_float(returnImage);
01291
01292
01293 neighbors=cpl_calloc(4*max_radius*max_radius,sizeof(float)) ;
01294
01295 for ( col = 0 ; col < ilx ; col++ )
01296 {
01297 for ( row = 0 ; row < ily ; row++ )
01298 {
01299
01300 if ( isnan(pmdata[col+row*ilx]) || pmdata[col+row*ilx] == 0. )
01301 {
01302
01303 n_valid = 0 ;
01304 agreed = 0 ;
01305 for ( j = 1 ; j <= max_radius ; j++ )
01306 {
01307
01308
01309 for ( k = -j ; k < j ; k++ )
01310 {
01311 if ( col-j >= 0 && row+k < ily && row+k >= 0 )
01312 {
01313 if ( !isnan(pmdata[col-j+(row+k)*mlx]) ||
01314 pmdata[col-j+(row+k)*mlx] != 0 )
01315 {
01316 neighbors[n_valid]=pidata[col-j+(row+k)*ilx] ;
01317 n_valid++ ;
01318 }
01319 }
01320 }
01321
01322
01323 for ( k = -j ; k < j ; k++ )
01324 {
01325 if ( col+k < ilx && col+k >= 0 && row+j < ily )
01326 {
01327 if ( !isnan(pmdata[col+k+(row+j)*mlx]) ||
01328 pmdata[col+k+(row+j)*mlx] != 0. )
01329 {
01330 neighbors[n_valid]=pidata[col+k+(row+j)*ilx] ;
01331 n_valid++ ;
01332 }
01333 }
01334 }
01335
01336
01337 for ( k = -j ; k < j ; k++ )
01338 {
01339 if ( col+j < ilx && row-k >= 0 && row-k < ily )
01340 {
01341 if ( !isnan(pmdata[col+j+(row-k)*mlx]) ||
01342 pmdata[col+j+(row-k)*mlx] != 0. )
01343 {
01344 neighbors[n_valid]=pidata[col+j+(row-k)*ilx] ;
01345 n_valid++ ;
01346 }
01347 }
01348 }
01349
01350
01351 for ( k = -j ; k < j ; k++ )
01352 {
01353 if ( col-k >= 0 && col-k < ilx && row-j < ily )
01354 {
01355 if ( !isnan(pmdata[col-k+(row-j)*mlx]) ||
01356 pmdata[col-k+(row-j)*mlx] != 0. )
01357 {
01358 neighbors[n_valid]=pidata[col-k+(row-j)*ilx] ;
01359 n_valid++ ;
01360 }
01361 }
01362 }
01363
01364
01365 if ( n_valid >= n_pixels )
01366 {
01367 agreed = 1 ;
01368 break ;
01369 }
01370
01371 if ( j == 1 && n_valid >= 2 )
01372 {
01373 agreed = 1 ;
01374 break ;
01375 }
01376 }
01377 if ( n_valid < n_pixels && agreed == 0 )
01378 {
01379 sinfo_msg_error("not enough valid neighbors found to "
01380 "interpolate bad pixel in col: "
01381 "%d, row: %d", col, row ) ;
01382 return NULL ;
01383 }
01384 else
01385 {
01386
01387
01388
01389
01390
01391 if ( n_valid <= 8 )
01392 {
01393 sum = 0. ;
01394
01395 for ( i = 0 ; i < n_valid ; i++ )
01396 {
01397 sum += neighbors[i] ;
01398 }
01399 mean = sum / n_valid ;
01400
01401 podata[col+row*ilx] = mean ;
01402 }
01403 else
01404 {
01405 podata[col+row*ilx]=sinfo_new_median(neighbors,n_valid);
01406 }
01407 }
01408 }
01409 }
01410 }
01411 cpl_free(neighbors);
01412 return returnImage ;
01413 }
01414
01415
01438 cpl_image * sinfo_interpol_source_image ( cpl_image * im,
01439 cpl_image * mask,
01440 int max_rad,
01441 float ** slit_edges )
01442 {
01443 cpl_image * returnImage=NULL ;
01444 float validpixel[6] ;
01445 float sum=0 ;
01446 int n=0;
01447 int row=0;
01448 int col=0;
01449 int i=0;
01450 int k=0;
01451 int slitlet=0;
01452 int n_slitlets=0;
01453 int ilx=0;
01454 int ily=0;
01455 int mlx=0;
01456 int mly=0;
01457
01458 float* pidata=NULL;
01459 float* podata=NULL;
01460 float* pmdata=NULL;
01461
01462
01463 if ( NULL == im )
01464 {
01465 sinfo_msg_error("sorry, no input image given!") ;
01466 return NULL ;
01467 }
01468 ilx=cpl_image_get_size_x(im);
01469 ily=cpl_image_get_size_y(im);
01470 pidata=cpl_image_get_data_float(im);
01471
01472 if ( NULL == mask )
01473 {
01474 sinfo_msg_error("sorry, no bad pixel mask image given!") ;
01475 return NULL ;
01476 }
01477 mlx=cpl_image_get_size_x(mask);
01478 mly=cpl_image_get_size_y(mask);
01479 pmdata=cpl_image_get_data_float(mask);
01480
01481 if ( mlx != ilx || mly != ily )
01482 {
01483 sinfo_msg_error("images not compatible in size!") ;
01484 return NULL ;
01485 }
01486
01487 if ( max_rad < 1 )
01488 {
01489 sinfo_msg_error("sorry, wrong maximum distance given!") ;
01490 return NULL ;
01491 }
01492
01493 if ( slit_edges == NULL )
01494 {
01495 sinfo_msg_error("sorry, array slit_edges is empty!") ;
01496 return NULL ;
01497 }
01498
01499
01500 n_slitlets = N_SLITLETS ;
01501
01502
01503 returnImage = cpl_image_duplicate( im ) ;
01504 podata=cpl_image_get_data_float(returnImage);
01505
01506
01507
01508 for ( row = 0 ; row < ily ; row++ )
01509 {
01510 for ( col = 0 ; col < ilx ; col++ )
01511 {
01512 n = 0 ;
01513 if ( isnan(pmdata[col + row*mlx]) ||
01514 pmdata[col + row*mlx] == 0. ||
01515 isnan(pidata[col + row*mlx]) )
01516 {
01517
01518 slitlet = -1000 ;
01519 for ( k = 0 ; k < n_slitlets ; k++ )
01520 {
01521 if ( sinfo_new_nint(slit_edges[k][0]) <= col &&
01522 sinfo_new_nint(slit_edges[k][1]) >= col )
01523 {
01524 slitlet = k ;
01525 }
01526
01527
01528
01529
01530
01531
01532
01533 }
01534 for ( i = 0 ; i < 6 ; i++ )
01535 {
01536 validpixel[i] = 0. ;
01537 }
01538
01539
01540 for ( i = 1 ; i <= max_rad ; i++ )
01541 {
01542 if ( row + i < ily)
01543 {
01544 if ( !isnan(pmdata[col + (row+i) * mlx])
01545 && pmdata[col + (row+i) * mlx] != 0. &&
01546 !isnan(pidata[col + (row+i) * ilx]) )
01547 {
01548 validpixel[n] = pidata[col + (row+i) * ilx] ;
01549 n++ ;
01550 }
01551 }
01552 if ( row - i >= 0 )
01553 {
01554 if ( !isnan(pmdata[col + (row-i) * mlx])
01555 && pmdata[col + (row-i) * mlx] != 0. &&
01556 !isnan(pidata[col + (row-i) * ilx]) )
01557 {
01558 validpixel[n] = pidata[col + (row-i) * ilx] ;
01559 n++ ;
01560 }
01561 }
01562
01563
01564
01565 if ( col + i < ilx )
01566 {
01567 if ( slitlet != -1000 )
01568 {
01569 if (col+i <= sinfo_new_nint(slit_edges[slitlet][1]) &&
01570 !isnan(pmdata[col + i + row * mlx]) &&
01571 pmdata[col + i + row * mlx] != 0. &&
01572 !isnan(pidata[col + i + row * ilx]) )
01573 {
01574 validpixel[n] = pidata[col + i + row * ilx] ;
01575 n++ ;
01576 }
01577 }
01578 }
01579 if ( col - i >= 0 )
01580 {
01581 if ( slitlet != -1000 )
01582 {
01583 if (col-i >= sinfo_new_nint(slit_edges[slitlet][0]) &&
01584 !isnan(pmdata[col - i + row * mlx]) &&
01585 pmdata[col - i + row * mlx] != 0. &&
01586 !isnan(pidata[col - i + row * ilx]) )
01587 {
01588 validpixel[n] = pidata[col - i + row * ilx] ;
01589 n++ ;
01590 }
01591 }
01592 }
01593
01594 if ( i == 1 && n > 1 )
01595 {
01596 break ;
01597 }
01598 if ( n > 2 )
01599 {
01600 break ;
01601 }
01602 }
01603
01604 if ( n == 0 )
01605 {
01606 podata[col + row*ilx] = ZERO ;
01607
01608
01609
01610
01611 }
01612 else
01613 {
01614
01615
01616 sum = 0. ;
01617 for ( i = 0 ; i < n ; i++ )
01618 {
01619 sum += validpixel[i] ;
01620 }
01621 podata[col + row*ilx] = sum/n ;
01622 }
01623 }
01624 }
01625 }
01626
01627 return returnImage ;
01628 }
01629
01639 cpl_image * sinfo_new_stack_row_to_image ( Vector * row, int ly )
01640 {
01641 cpl_image * image=NULL;
01642 int col=0;
01643 int ro=0;
01644 float* podata=NULL;
01645
01646 if ( row == NullVector )
01647 {
01648 sinfo_msg_error ("Null sinfo_vector as input" ) ;
01649 return NULL ;
01650 }
01651 if ( ly <= 1 )
01652 {
01653 sinfo_msg_error ("wrong image length given" ) ;
01654 return NULL ;
01655 }
01656
01657
01658 if (NULL == (image = cpl_image_new(row->n_elements ,ly,CPL_TYPE_FLOAT )) )
01659 {
01660 sinfo_msg_error ("cannot allocate new image" ) ;
01661 return NULL ;
01662 }
01663 podata=cpl_image_get_data_float(image);
01664
01665 for ( col = 0 ; col < row -> n_elements ; col++ )
01666 {
01667 for ( ro = 0 ; ro < ly ; ro++ )
01668 {
01669 podata[col + ro*ly] = row -> data[col] ;
01670 }
01671 }
01672 return image ;
01673 }
01674
01690 Stats * sinfo_new_image_stats_on_rectangle ( cpl_image * im,
01691 float loReject,
01692 float hiReject,
01693 int llx,
01694 int lly,
01695 int urx,
01696 int ury )
01697 {
01698 Stats * retstats=NULL;
01699 int i=0 ;
01700 int row=0;
01701 int col=0;
01702 int n=0;
01703 int npix=0;
01704 int lo_n=0;
01705 int hi_n=0;
01706 double pix_sum=0;
01707 double sqr_sum=0;
01708 float * pix_array=NULL;
01709 int im_lx=0;
01710 int im_ly=0;
01711 float* pim=NULL;
01712
01713 if ( NULL == im )
01714 {
01715 sinfo_msg_error("sorry, no input image given!") ;
01716 return NULL ;
01717 }
01718 if ( loReject+hiReject >= 100. )
01719 {
01720 sinfo_msg_error("sorry, too much pixels rejected!") ;
01721 return NULL ;
01722 }
01723 if ( loReject < 0. || loReject >= 100. ||
01724 hiReject < 0. || hiReject >= 100. )
01725 {
01726 sinfo_msg_error("sorry, negative reject values!") ;
01727 return NULL ;
01728 }
01729
01730 im_lx=cpl_image_get_size_x(im);
01731 im_ly=cpl_image_get_size_x(im);
01732
01733 if ( llx < 0 || lly < 0 || urx < 0 || ury < 0 ||
01734 llx >= im_lx || lly >= im_ly || urx >= im_lx ||
01735 ury >= im_ly || ury <= lly || urx <= llx )
01736 {
01737 sinfo_msg_error("sorry, wrong pixel coordinates of rectangle!") ;
01738 return NULL ;
01739 }
01740
01741
01742 retstats = (Stats*) cpl_calloc(1, sizeof(Stats)) ;
01743 npix = (urx - llx + 1) * (ury - lly + 1) ;
01744 pix_array = (float*) cpl_calloc ( npix, sizeof(float) ) ;
01745
01746
01747
01748
01749 n = 0 ;
01750 pim = cpl_image_get_data_float(im);
01751 for ( row = lly ; row <= ury ; row++ )
01752 {
01753 for ( col = llx ; col <= urx ; col++ )
01754 {
01755 if ( !isnan(pim[col + row*im_lx]) )
01756 {
01757 pix_array[n] = pim[col + row*im_lx] ;
01758 n++ ;
01759 }
01760 }
01761 }
01762
01763 npix = n;
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774 if ( FLT_MAX == (retstats->cleanmean = sinfo_new_clean_mean(pix_array,
01775 npix, loReject, hiReject)) )
01776 {
01777 sinfo_msg_error("sinfo_new_clean_mean() did not work!") ;
01778 cpl_free(retstats) ;
01779 cpl_free(pix_array) ;
01780 return NULL ;
01781 }
01782
01783
01784
01785 lo_n = (int) (loReject / 100. * (float)npix) ;
01786 hi_n = (int) (hiReject / 100. * (float)npix) ;
01787 pix_sum = 0. ;
01788 sqr_sum = 0. ;
01789 n = 0 ;
01790 for ( i = lo_n ; i <= npix - hi_n ; i++ )
01791 {
01792 pix_sum += (double)pix_array[i] ;
01793 sqr_sum += ((double)pix_array[i] * (double)pix_array[i]) ;
01794 n++ ;
01795 }
01796
01797 if ( n == 0 )
01798 {
01799 sinfo_msg_error("number of clean pixels is zero!") ;
01800 cpl_free(retstats) ;
01801 cpl_free(pix_array) ;
01802 return NULL ;
01803 }
01804 retstats -> npix = n ;
01805 pix_sum /= (double) n ;
01806 sqr_sum /= (double) n ;
01807 retstats -> cleanstdev = (float)sqrt(sqr_sum - pix_sum * pix_sum) ;
01808 cpl_free (pix_array) ;
01809 return retstats ;
01810 }
01811
01812
01813
01822 cpl_image * sinfo_new_normalize_to_central_pixel ( cpl_image * image )
01823 {
01824 int col=0;
01825 int row=0;
01826 int i=0;
01827 int n=0;
01828 float* array=NULL ;
01829 float divisor=0;
01830 cpl_image * retImage=NULL;
01831 int ilx=0;
01832 int ily=0;
01833 float* pidata=NULL;
01834 float* podata=NULL;
01835
01836 if ( image == NULL )
01837 {
01838 sinfo_msg_error("no image given!") ;
01839 return NULL ;
01840 }
01841 ilx=cpl_image_get_size_x(image);
01842 ily=cpl_image_get_size_y(image);
01843 pidata=cpl_image_get_data_float(image);
01844
01845 retImage = cpl_image_duplicate(image) ;
01846 podata=cpl_image_get_data_float(retImage);
01847
01848 n = 0 ;
01849
01850
01851 array=cpl_calloc(2*ilx,sizeof(float)) ;
01852
01853 for ( row = ily/2 ; row < ily/2+1 ; row++ )
01854 {
01855 for ( col = 0 ; col < ilx ; col++ )
01856 {
01857 if ( !isnan(pidata[col+ilx*row]) )
01858 {
01859 array[n] = pidata[col+ilx*row] ;
01860 n++ ;
01861 }
01862 }
01863 }
01864
01865
01866 if ( isnan(divisor = sinfo_new_median(array, n) ) )
01867 {
01868 sinfo_msg_error("no sinfo_median possible!") ;
01869 return NULL ;
01870 }
01871 if ( 0 == divisor )
01872 {
01873 sinfo_msg_error("cannot divide by 0") ;
01874 return NULL ;
01875 }
01876
01877 for ( i = 0 ; i < (int) ilx*ily ; i++ )
01878 {
01879 if ( isnan(pidata[i]) )
01880 {
01881 podata[i] = ZERO ;
01882 }
01883 else
01884 {
01885 podata[i] = pidata[i]/divisor ;
01886 }
01887 }
01888 cpl_free(array);
01889 return retImage ;
01890 }
01891
01892
01921
01922
01923 cpl_image *
01924 sinfo_new_mpe_shift_image(
01925 cpl_image * image_in,
01926 double shift_x,
01927 double shift_y,
01928 double * interp_kernel)
01929 {
01930 cpl_image * shifted=NULL ;
01931 pixelvalue * first_pass=NULL ;
01932 pixelvalue * second_pass=NULL ;
01933 int samples = KERNEL_SAMPLES ;
01934 int i=0, j=0 ;
01935 double fx=0, fy=0 ;
01936 double rx=0, ry=0 ;
01937 int px=0, py=0 ;
01938 int tabx=0, taby=0 ;
01939 double value=0 ;
01940 size_t pos ;
01941 register pixelvalue * pix ;
01942 register pixelvalue * pixint ;
01943 int mid=0;
01944 double norm=0 ;
01945 double * ker=NULL ;
01946 int freeKernel = 1 ;
01947
01948 int ilx=0;
01949 int ily=0;
01950 float* pidata=NULL;
01951 float* psdata=NULL;
01952
01953
01954
01955 if (image_in==NULL) return NULL ;
01956
01957
01958 if ((fabs(shift_x)<1e-2) && (fabs(shift_y)<1e-2))
01959 return cpl_image_duplicate(image_in) ;
01960 ilx=cpl_image_get_size_x(image_in);
01961 ily=cpl_image_get_size_y(image_in);
01962 pidata=cpl_image_get_data_float(image_in);
01963
01964
01965
01966 if (interp_kernel == NULL) {
01967 ker = sinfo_generate_interpolation_kernel("default") ;
01968 if (ker == NULL) {
01969 sinfo_msg_error("kernel generation failure:aborting resampling") ;
01970 return NULL ;
01971 }
01972 } else {
01973 ker = interp_kernel ;
01974 freeKernel = 0 ;
01975 }
01976
01977 mid = (int)samples/(int)2 ;
01978 first_pass = cpl_calloc(ilx, ily*sizeof(pixelvalue)) ;
01979 shifted = cpl_image_new(ilx, ily,CPL_TYPE_FLOAT) ;
01980 psdata=cpl_image_get_data_float(shifted);
01981
01982 second_pass = psdata ;
01983
01984 pix = pidata ;
01985 if ( ilx != 1 )
01986 {
01987 for (j=0 ; j<ily ; j++)
01988 {
01989 for (i=0 ; i<ilx ; i++) {
01990 fx = (double)i-shift_x ;
01991 px = (int)fx ;
01992 rx = fx - (double)px ;
01993 pos = px + j * ilx ;
01994
01995 if ((px>1) && (px<(ilx-2)))
01996 {
01997 tabx = (int)(fabs((double)mid * rx)) ;
01998
01999 if (isnan(pix[pos]))
02000 {
02001 value = ZERO ;
02002 }
02003 else
02004 {
02005 if (isnan(pix[pos-1]))
02006 {
02007 pix[pos-1] = 0. ;
02008 }
02009 if (isnan(pix[pos+1]))
02010 {
02011 pix[pos+1] = 0. ;
02012 }
02013 if (isnan(pix[pos+2]))
02014 {
02015 pix[pos+2] = 0. ;
02016 }
02017
02018
02019
02020
02021
02022 value = (double)pix[pos-1] * ker[mid+tabx] +
02023 (double)pix[pos] * ker[tabx] +
02024 (double)pix[pos+1] * ker[mid-tabx] +
02025 (double)pix[pos+2] * ker[samples-tabx-1] ;
02026
02027
02028
02029
02030 norm = (double)ker[mid+tabx] +
02031 (double)ker[tabx] +
02032 (double)ker[mid-tabx] +
02033 (double)ker[samples-tabx-1] ;
02034 if (fabs(norm) > 1e-4) {
02035 value /= norm ;
02036 }
02037 }
02038 } else {
02039 value = ZERO ;
02040 }
02041
02042
02043
02044
02045 if ( isnan(value) )
02046 {
02047 first_pass[i+j*ilx] = ZERO ;
02048 }
02049 else
02050 {
02051 first_pass[i+j*ilx] = (pixelvalue)value ;
02052 }
02053 }
02054 }
02055 }
02056 else
02057 {
02058 memcpy(first_pass,pix,ily*sizeof(pixelvalue));
02059 }
02060
02061 pixint = first_pass ;
02062 for (i=0 ; i<ilx ; i++) {
02063 for (j=0 ; j<ily ; j++) {
02064 fy = (double)j - shift_y ;
02065 py = (int)fy ;
02066 ry = fy - (double)py ;
02067 pos = i + py * ilx ;
02068
02069 taby = (int)(fabs((double)mid * ry)) ;
02070
02071 if ((py>(int)1) && (py<(ily-2))) {
02072
02073 if (isnan(pixint[pos]) && ilx != 1 )
02074 {
02075 value = ZERO ;
02076 }
02077 else
02078 {
02079 if (isnan(pixint[pos-ilx]))
02080 {
02081 pixint[pos-ilx] = 0. ;
02082 }
02083 if (isnan(pixint[pos+ilx]))
02084 {
02085 pixint[pos+ilx] = 0. ;
02086 }
02087 if (isnan(pixint[pos+2*ilx]))
02088 {
02089 pixint[pos+2*ilx] = 0. ;
02090 }
02091
02092
02093
02094
02095 value = (double)pixint[pos-ilx] * ker[mid+taby] +
02096 (double)pixint[pos] * ker[taby] +
02097 (double)pixint[pos+ilx] * ker[mid-taby] +
02098 (double)pixint[pos+2*ilx]*ker[samples-taby-1];
02099
02100
02101
02102
02103 norm = (double)ker[mid+taby] +
02104 (double)ker[taby] +
02105 (double)ker[mid-taby] +
02106 (double)ker[samples-taby-1] ;
02107
02108 if (fabs(norm) > 1e-4) {
02109 value /= norm ;
02110 }
02111 }
02112 } else {
02113 value = ZERO ;
02114 }
02115 if (isnan(value))
02116 {
02117 second_pass[i+j*ilx] = ZERO ;
02118 }
02119 else
02120 {
02121 second_pass[i+j*ilx] = (pixelvalue)value ;
02122 }
02123 }
02124 }
02125
02126 cpl_free(first_pass) ;
02127 if (freeKernel)
02128 cpl_free(ker) ;
02129 return shifted ;
02130 }
02131
02132
02133
02164
02165
02166 void
02167 sinfo_new_shift_image_in_cube(
02168 cpl_image * image_in,
02169 double shift_x,
02170 double shift_y,
02171 double * interp_kernel,
02172 cpl_image * shifted,
02173 pixelvalue * first_pass)
02174 {
02175 pixelvalue * second_pass=NULL ;
02176 int samples = KERNEL_SAMPLES ;
02177 int i=0, j=0 ;
02178 double fx=0, fy=0 ;
02179 double rx=0, ry=0 ;
02180 int px=0, py=0 ;
02181 int tabx=0, taby=0 ;
02182 double value=0 ;
02183 size_t pos ;
02184 register pixelvalue * pix ;
02185 register pixelvalue * pixint ;
02186 int mid=0;
02187 double norm=0 ;
02188 double * ker=NULL ;
02189 int freeKernel = 1 ;
02190
02191 int ilx=0;
02192 int ily=0;
02193 int slx=0;
02194 int sly=0;
02195 float* pidata=NULL;
02196 float* psdata=NULL;
02197
02198
02199 if (image_in==NULL) shifted = NULL ;
02200 pidata=cpl_image_get_data_float(image_in);
02201 ilx=cpl_image_get_size_x(image_in);
02202 ily=cpl_image_get_size_y(image_in);
02203
02204 shifted=cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
02205 slx=ilx;
02206 sly=ily;
02207
02208 psdata=cpl_image_get_data_float(shifted);
02209
02210
02211 if ((fabs(shift_x)<1e-2) && (fabs(shift_y)<1e-2))
02212 memcpy(psdata,pidata, (size_t) slx*sly * sizeof(pixelvalue)) ;
02213
02214
02215 if (interp_kernel == NULL) {
02216 ker = sinfo_generate_interpolation_kernel("default") ;
02217 if (ker == NULL) {
02218 sinfo_msg_error("kernel generation failure:aborting resampling") ;
02219 shifted = NULL ;
02220 }
02221 } else {
02222 ker = interp_kernel ;
02223 freeKernel = 0 ;
02224 }
02225
02226 mid = (int)samples/(int)2 ;
02227 second_pass = psdata ;
02228
02229 pix = pidata ;
02230 for (j=0 ; j<ily ; j++) {
02231 for (i=1 ; i<ilx-2 ; i++) {
02232 fx = (double)i-shift_x ;
02233 px = (int)fx ;
02234 rx = fx - (double)px ;
02235
02236 pos = px + j * ilx ;
02237
02238 if ((px>1) && (px<(ilx-2))) {
02239 tabx = (int)(fabs((double)mid * rx)) ;
02240
02241 if (isnan(pix[pos]))
02242 {
02243 value = ZERO ;
02244 }
02245 else
02246 {
02247 if (isnan(pix[pos-1]))
02248 {
02249 pix[pos-1] = 0. ;
02250 }
02251 if (isnan(pix[pos+1]))
02252 {
02253 pix[pos+1] = 0. ;
02254 }
02255 if (isnan(pix[pos+2]))
02256 {
02257 pix[pos+2] = 0. ;
02258 }
02259
02260
02261
02262
02263
02264 value = (double)pix[pos-1] * ker[mid+tabx] +
02265 (double)pix[pos] * ker[tabx] +
02266 (double)pix[pos+1] * ker[mid-tabx] +
02267 (double)pix[pos+2] * ker[samples-tabx-1] ;
02268
02269
02270
02271
02272 norm = (double)ker[mid+tabx] +
02273 (double)ker[tabx] +
02274 (double)ker[mid-tabx] +
02275 (double)ker[samples-tabx-1] ;
02276 if (fabs(norm) > 1e-4) {
02277 value /= norm ;
02278 }
02279 }
02280 } else {
02281 value = 0.0 ;
02282 }
02283
02284
02285
02286
02287 if ( isnan(value) )
02288 {
02289 first_pass[i+j*ilx] = ZERO ;
02290 }
02291 else
02292 {
02293 first_pass[i+j*ilx] = (pixelvalue)value ;
02294 }
02295 }
02296 }
02297 pixint = first_pass ;
02298 for (i=0 ; i< ilx ; i++) {
02299 for (j=1 ; j< ily-2 ; j++) {
02300 fy = (double)j - shift_y ;
02301 py = (int)fy ;
02302 ry = fy - (double)py ;
02303 pos = i + py * ilx ;
02304
02305 taby = (int)(fabs((double)mid * ry)) ;
02306
02307 if ((py>(int)1) && (py<(ily-2))) {
02308
02309 if (isnan(pixint[pos]))
02310 {
02311 value = ZERO ;
02312 }
02313 else
02314 {
02315 if (isnan(pixint[pos-ilx]))
02316 {
02317 pixint[pos-ilx] = 0. ;
02318 }
02319 if (isnan(pixint[pos+ilx]))
02320 {
02321 pixint[pos+ilx] = 0. ;
02322 }
02323 if (isnan(pixint[pos+2*ilx]))
02324 {
02325 pixint[pos+2*ilx] = 0. ;
02326 }
02327
02328
02329
02330
02331 value = (double)pixint[pos-ilx] * ker[mid+taby] +
02332 (double)pixint[pos] * ker[taby] +
02333 (double)pixint[pos+ilx] * ker[mid-taby] +
02334 (double)pixint[pos+2*ilx]*ker[samples-taby-1];
02335
02336
02337
02338
02339 norm = (double)ker[mid+taby] +
02340 (double)ker[taby] +
02341 (double)ker[mid-taby] +
02342 (double)ker[samples-taby-1] ;
02343
02344 if (fabs(norm) > 1e-4) {
02345 value /= norm ;
02346 }
02347 }
02348 } else {
02349
02350 }
02351 if (isnan(value))
02352 {
02353 second_pass[i+j*ilx] = ZERO ;
02354 }
02355 else
02356 {
02357 second_pass[i+j*ilx] = (pixelvalue)value ;
02358 }
02359 }
02360 }
02361
02362 if (freeKernel)
02363 cpl_free(ker) ;
02364 }
02365
02366
02367 void sinfo_new_del_Stats( Stats * st)
02368 {
02369 cpl_free (st) ;
02370 }
02371
02378 cpl_image *
02379 sinfo_new_combine_masks ( cpl_image * firstMask, cpl_image * secondMask )
02380 {
02381 cpl_image * retMask=NULL ;
02382 int n=0 ;
02383 int olx=0;
02384 int oly=0;
02385 float* podata=NULL;
02386 float* pm1data=NULL;
02387 float* pm2data=NULL;
02388
02389 if ( firstMask == NULL || secondMask == NULL )
02390 {
02391 sinfo_msg_error ("no input mask image given!") ;
02392 return NULL ;
02393 }
02394 retMask = cpl_image_duplicate (firstMask) ;
02395 podata = cpl_image_get_data_float(retMask);
02396 pm1data = cpl_image_get_data_float(firstMask);
02397 pm2data = cpl_image_get_data_float(secondMask);
02398 olx=cpl_image_get_size_x(retMask);
02399 oly=cpl_image_get_size_y(retMask);
02400
02401 for ( n = 0 ; n < (int) olx*oly ; n++ )
02402 {
02403 if ( podata[n] == 0. || pm2data[n] == 0. )
02404 {
02405 podata[n] = 0. ;
02406 }
02407 else
02408 {
02409 podata[n] = 1. ;
02410 }
02411 }
02412 return retMask ;
02413 }
02414
02423 cpl_image * sinfo_new_slice_cube (cpl_imagelist * cube, int x, int y )
02424 {
02425 cpl_image * retImage=NULL ;
02426 int col=0, row=0, z=0 ;
02427 int inp=0;
02428 int ilx=0;
02429 int ily=0;
02430 cpl_image* img=NULL;
02431 float* podata=NULL;
02432 float* pidata=NULL;
02433
02434 if ( cube == NULL )
02435 {
02436 sinfo_msg_error("no cube given!") ;
02437 return NULL ;
02438 }
02439 if ( x > 31 || y > 31 )
02440 {
02441 sinfo_msg_warning("wrong x or y values!") ;
02442 }
02443
02444 img=cpl_imagelist_get(cube,0);
02445 ilx=cpl_image_get_size_x(img);
02446 ily=cpl_image_get_size_y(img);
02447 inp=cpl_imagelist_get_size(cube);
02448 if ( x < 0 )
02449 {
02450
02451 if ( NULL == (retImage = cpl_image_new(ilx, inp, CPL_TYPE_FLOAT)) )
02452 {
02453 sinfo_msg_error("could not allocate memory!") ;
02454 return NULL ;
02455 }
02456 podata=cpl_image_get_data_float(retImage);
02457 for ( z = 0 ; z < inp ; z++ )
02458 {
02459
02460 pidata=cpl_image_get_data_float(cpl_imagelist_get(cube,z));
02461 for ( col = 0 ; col < ilx ; col++ )
02462 {
02463 podata[col+z*ilx] = pidata[col+y*ilx] ;
02464 }
02465 }
02466 }
02467 else if ( y < 0 )
02468 {
02469
02470 if ( NULL == (retImage = cpl_image_new(ily, inp,CPL_TYPE_FLOAT)) )
02471 {
02472 sinfo_msg_error("could not allocate memory!") ;
02473 return NULL ;
02474 }
02475 podata=cpl_image_get_data_float(retImage);
02476
02477 for ( z = 0 ; z < inp ; z++ )
02478 {
02479 pidata=cpl_image_get_data_float(cpl_imagelist_get(cube,z));
02480 for ( row = 0 ; row < ily ; row++ )
02481 {
02482 podata[row+z*ily] = pidata[x+row*ily] ;
02483 }
02484 }
02485 }
02486 else
02487 {
02488 sinfo_msg_error("wrong input!") ;
02489 return NULL ;
02490 }
02491 return retImage ;
02492 }
02493
02505 cpl_image * sinfo_new_div_images_robust ( cpl_image * im1, cpl_image * im2 )
02506 {
02507 cpl_image * retIm=NULL ;
02508 float help=0 ;
02509 int i=0 ;
02510 int lx1=0;
02511 int ly1=0;
02512 int lx2=0;
02513 int ly2=0;
02514
02515 float* p1data=NULL;
02516 float* p2data=NULL;
02517 float* podata=NULL;
02518
02519 if ( im1 == NULL || im2 == NULL )
02520 {
02521 sinfo_msg_error("no input images given!") ;
02522 return NULL ;
02523 }
02524 lx1=cpl_image_get_size_x(im1);
02525 ly1=cpl_image_get_size_y(im1);
02526 lx2=cpl_image_get_size_x(im2);
02527 ly2=cpl_image_get_size_y(im2);
02528 p1data=cpl_image_get_data_float(im1);
02529 p2data=cpl_image_get_data_float(im2);
02530
02531 if ( lx1 != lx2 || ly1 != ly2 )
02532 {
02533 sinfo_msg_error("images not compatible!") ;
02534 return NULL ;
02535 }
02536 if ( NULL == (retIm = cpl_image_new(lx1, ly1, CPL_TYPE_FLOAT)) )
02537 {
02538 sinfo_msg_error("could not allocate memory!") ;
02539 return NULL ;
02540 }
02541 podata=cpl_image_get_data_float(retIm);
02542
02543 for ( i = 0 ; i < (int) lx2*ly2 ; i++ )
02544 {
02545 if ( !isnan(p2data[i]) )
02546 {
02547 help = 1./p2data[i] ;
02548 if (fabs( help )> THRESH )
02549 {
02550 help = 1. ;
02551 }
02552 }
02553 else
02554 {
02555 help = ZERO ;
02556 }
02557 if ( isnan(help) || isnan(p1data[i]) )
02558 {
02559 podata[i] = ZERO ;
02560 }
02561 else
02562 {
02563 podata[i] = p1data[i] * help ;
02564 }
02565 }
02566 return retIm ;
02567 }
02568
02569 cpl_image * sinfo_new_null_edges ( cpl_image * image)
02570 {
02571 cpl_image * new=NULL ;
02572 int i=0,j=0 ;
02573 int ilx=0;
02574 int ily=0;
02575 int olx=0;
02576 int oly=0;
02577
02578 float* pidata=NULL;
02579 float* podata=NULL;
02580
02581 if ( image == NULL )
02582 {
02583 sinfo_msg_error ("no input image given!\n") ;
02584 return NULL ;
02585 }
02586
02587
02588 new = cpl_image_duplicate (image) ;
02589 ilx=cpl_image_get_size_x(image);
02590 ily=cpl_image_get_size_y(image);
02591 olx=cpl_image_get_size_x(new);
02592 oly=cpl_image_get_size_y(new);
02593 pidata=cpl_image_get_data_float(image);
02594 podata=cpl_image_get_data_float(new);
02595
02596 for ( i = 0 ; i < olx ; i++ )
02597 {
02598 for ( j = 0 ; j < 4 ; j++)
02599 {
02600 podata[i+j*olx]=0;
02601 podata[i+(oly-j-1)*olx]=0;
02602 }
02603 }
02604 for ( i = 0 ; i < oly ; i++ )
02605 {
02606 for ( j = 0 ; j < 4 ; j++)
02607 {
02608 podata[j+i*olx]=0;
02609 podata[(olx-j-1)+i*olx]=0;
02610 }
02611 }
02612 return new ;
02613 }
02614
02615
02616 void sinfo_new_used_cor_map( cpl_image *im, cpl_image *map)
02617 {
02618 int i=0,j=0,index=0;
02619 float temp_array[2048];
02620 int lx=cpl_image_get_size_x(im);
02621 int ly=cpl_image_get_size_y(im);
02622 float* pidata=cpl_image_get_data_float(im);
02623 float* pmdata=cpl_image_get_data_float(map);
02624
02625 for( j=0; j<ly; j++)
02626 {
02627 for( i=0;i<lx;i++)
02628 {
02629 index = (int)pmdata[i+j*lx];
02630 temp_array[i] = pidata[index+j*lx];
02631 }
02632 for( i=0;i<lx;i++)
02633 {
02634 pidata[i+j*lx]= temp_array[i];
02635 }
02636 }
02637 }
02638
02639
02640
02641
02642
02665
02666
02667 cpl_image *
02668 sinfo_new_shift_image(
02669 cpl_image * image_in,
02670 double shift_x,
02671 double shift_y,
02672 double * interp_kernel)
02673 {
02674 cpl_image * shifted=NULL ;
02675 float * first_pass=NULL ;
02676 float * second_pass=NULL ;
02677 int samples = KERNEL_SAMPLES ;
02678 int i=0, j=0 ;
02679 double fx=0, fy=0 ;
02680 double rx=0, ry=0 ;
02681 int px=0, py=0 ;
02682 int tabx=0, taby=0 ;
02683 double value=0 ;
02684 size_t pos ;
02685 register float * pix=NULL ;
02686 register float * pixint=NULL ;
02687 int mid=0;
02688 double norm=0 ;
02689 double * ker=NULL ;
02690 int freeKernel = 1 ;
02691 int ilx=0;
02692 int ily=0;
02693
02694
02695 if (image_in==NULL) return NULL ;
02696
02697
02698 if ((fabs(shift_x)<1e-2) && (fabs(shift_y)<1e-2))
02699 return cpl_image_duplicate(image_in) ;
02700
02701
02702 if (interp_kernel == NULL) {
02703 ker = sinfo_generate_interpolation_kernel("default") ;
02704 if (ker == NULL) {
02705 sinfo_msg_error("kernel generation failure: aborting resampling") ;
02706 return NULL ;
02707 }
02708 } else {
02709 ker = interp_kernel ;
02710 freeKernel = 0 ;
02711 }
02712
02713 ilx=cpl_image_get_size_x(image_in);
02714 ily=cpl_image_get_size_y(image_in);
02715
02716 mid = (int)samples/(int)2 ;
02717 first_pass = cpl_calloc(ilx, ily*sizeof(float)) ;
02718 shifted = cpl_image_new(ilx, ily,CPL_TYPE_FLOAT) ;
02719 second_pass = cpl_image_get_data_float(shifted);
02720
02721 pix = cpl_image_get_data_float(image_in);
02722 for (j=0 ; j<ily ; j++) {
02723 for (i=1 ; i<ilx-2 ; i++) {
02724 fx = (double)i-shift_x ;
02725 px = (int)fx ;
02726 rx = fx - (double)px ;
02727
02728 pos = px + j * ilx ;
02729
02730 if ((px>1) && (px<(ilx-3))) {
02731 tabx = (int)(fabs((double)mid * rx)) ;
02732
02733
02734
02735
02736 value = (double)pix[pos-1] * ker[mid+tabx] +
02737 (double)pix[pos] * ker[tabx] +
02738 (double)pix[pos+1] * ker[mid-tabx] +
02739 (double)pix[pos+2] * ker[samples-tabx-1] ;
02740
02741
02742
02743
02744 norm = (double)ker[mid+tabx] +
02745 (double)ker[tabx] +
02746 (double)ker[mid-tabx] +
02747 (double)ker[samples-tabx-1] ;
02748 if (fabs(norm) > 1e-4) {
02749 value /= norm ;
02750 }
02751 } else {
02752 value = 0.0 ;
02753 }
02754
02755
02756
02757
02758 first_pass[i+j*ilx] = (float)value ;
02759 }
02760 }
02761 pixint = first_pass ;
02762 for (i=0 ; i<ilx ; i++) {
02763 for (j=1 ; j<ily-3 ; j++) {
02764 fy = (double)j - shift_y ;
02765 py = (int)fy ;
02766 ry = fy - (double)py ;
02767 pos = i + py * ilx ;
02768
02769 taby = (int)(fabs((double)mid * ry)) ;
02770
02771 if ((py>(int)1) && (py<(ily-2))) {
02772
02773
02774
02775
02776 value = (double)pixint[pos-ilx] * ker[mid+taby] +
02777 (double)pixint[pos] * ker[taby] +
02778 (double)pixint[pos+ilx] * ker[mid-taby] +
02779 (double)pixint[pos+2*ilx]*ker[samples-taby-1];
02780
02781
02782
02783
02784 norm = (double)ker[mid+taby] +
02785 (double)ker[taby] +
02786 (double)ker[mid-taby] +
02787 (double)ker[samples-taby-1] ;
02788
02789 if (fabs(norm) > 1e-4) {
02790 value /= norm ;
02791 }
02792 } else {
02793 value = 0.0 ;
02794 }
02795 second_pass[i+j*ilx] = (float)value ;
02796 }
02797 }
02798
02799 cpl_free(first_pass) ;
02800 if (freeKernel)
02801 cpl_free(ker) ;
02802 return shifted ;
02803 }
02804
02805