00001
00002
00003
00004
00005
00006
00007 #include <stdlib.h>
00008
00009 #include "util.h"
00010 #include "imcore.h"
00011
00012
00013
00014 static void sortm (float ia[], int ib[], int n);
00015 static void quicksort (float x[], int point[], int l, int nfilt);
00016
00017
00018
00019 extern void bfilt (float **xbuf, int nx, int ny) {
00020 float *ybuf, *save;
00021 int mfilt = 5, j, k;
00022
00023
00024 ybuf = (float *) malloc(MAX(nx,ny) * sizeof(float));
00025 save = (float *) malloc((nx+1) * ny * sizeof(float));
00026
00027
00028
00029 for(k = 0; k < ny; k++) {
00030 for(j = 0; j < nx; j++) {
00031 save[(nx+1)*k+j] = xbuf[k][j];
00032 ybuf[j] = xbuf[k][j];
00033 }
00034 filt1d(ybuf, nx, mfilt);
00035 for(j = 0; j < nx; j++) xbuf[k][j] = ybuf[j];
00036 }
00037
00038
00039 for(k = 0; k < nx; k++) {
00040 for(j = 0; j < ny; j++) ybuf[j] = xbuf[j][k];
00041 filt1d(ybuf, ny, mfilt);
00042 for(j = 0; j < ny; j++)
00043
00044 if(save[(nx+1)*j+k] > -1000.0)
00045 xbuf[j][k] = MIN(save[(nx+1)*j+k], ybuf[j]);
00046 }
00047
00048
00049 for(k = 0; k < ny; k++) {
00050 for(j = 0; j < nx; j++) ybuf[j] = xbuf[k][j];
00051 hanning(ybuf, nx);
00052 for(j = 0; j < nx; j++) xbuf[k][j] = ybuf[j];
00053 }
00054
00055
00056 for(k = 0; k < nx; k++) {
00057 for(j = 0; j < ny; j++) ybuf[j] = xbuf[j][k];
00058 hanning(ybuf, ny);
00059 for(j = 0; j < ny; j++) xbuf[j][k] = ybuf[j];
00060 }
00061
00062
00063 free((void *) ybuf);
00064 free((void *) save);
00065 }
00066
00067
00068
00069 void filt1d (float ybuf[], int mpt, int mfilt) {
00070 float *wbuf;
00071 int i, irc;
00072
00073 wbuf = (float *) malloc(mpt * sizeof(float));
00074
00075
00076 irc = 0;
00077 for(i = 0; i < mpt; i++){
00078 if(ybuf[i] > -1000.0){
00079 wbuf[irc] = ybuf[i];
00080 irc++;
00081 }
00082 }
00083 if(irc == 0) {
00084 free((void *) wbuf);
00085 return;
00086 }
00087 median(wbuf, irc, mfilt);
00088 irc = 0;
00089 for(i = 0; i < mpt; i++){
00090 if(ybuf[i] > -1000.0){
00091 ybuf[i] = wbuf[irc];
00092 irc++;
00093 }
00094 }
00095 padext(ybuf, mpt);
00096
00097 free((void *) wbuf);
00098 }
00099
00100
00101
00102
00103 void padext (float x[], int n) {
00104 int i, j, ilow, ihih=0, ic;
00105 float xlow, xhih, slope, t1 ,t2;
00106
00107 i = 0;
00108 while(x[i] <= -1000.0) i++;
00109 ilow = i;
00110 for(i = ilow+1; i < n; i++){
00111 if(x[i] <= -1000.0) {
00112 ic = 1;
00113 if (i < n - 1) {
00114 while(x[i+ic] <= -1000.0) {
00115 ic++;
00116 if (i+ic >= n-1)
00117 break;
00118 }
00119 }
00120 if(i+ic < n-1){
00121 xlow = x[i-1];
00122 xhih = x[i+ic];
00123 for(j = 0; j < ic; j++){
00124 t2 = ((float) j+1)/((float) ic+1);
00125 t1 = 1.0 - t2;
00126 x[i+j] = t1*xlow+t2*xhih;
00127 }
00128 }
00129 } else {
00130 ihih = i;
00131 }
00132 }
00133
00134 if(ilow > 0){
00135 slope = x[ilow+1]-x[ilow];
00136 for(i = 0; i < ilow; i++) x[i] = x[ilow]-slope*(ilow-i);
00137 }
00138 if(ihih < n-1) {
00139 slope = x[ihih]-x[ihih-1];
00140 for(i = ihih+1; i < n; i++) x[i] = x[ihih]+slope*(i-ihih);
00141 }
00142 }
00143
00144
00145
00146 void hanning (float xbuf[], int npt) {
00147 float *ybuf;
00148 float sum = 0.0, xmns, xmnf;
00149 int nfilt = 3, i, il, ilow, nelem;
00150
00151 if(npt <= nfilt)
00152 return;
00153
00154
00155 il = nfilt/2;
00156 ilow = MAX(3,nfilt/4);
00157 ilow = (ilow/2)*2 + 1;
00158
00159 for(i = 0; i < ilow; i++)
00160 sum += xbuf[i];
00161
00162 xmns = sum/((float) ilow);
00163
00164 sum=0.0;
00165 for(i = 0; i < ilow; i++)
00166 sum += xbuf[npt-1-i];
00167
00168 xmnf = sum/((float) ilow);
00169
00170
00171 nelem = npt + nfilt;
00172
00173 ybuf = (float *) malloc(nelem * sizeof(float));
00174
00175
00176
00177
00178 for(i = 0; i < il; i++) {
00179 ybuf[i] = 2.0 * xmns - xbuf[il+ilow-1-i];
00180 ybuf[npt+i+il] = 2.0 * xmnf - xbuf[npt-i-ilow-1];
00181 }
00182
00183 for(i = 0; i < npt; i++)
00184 ybuf[i+il] = xbuf[i];
00185
00186
00187 for(i = 0; i < npt; i++)
00188 xbuf[i] = 0.25 * (ybuf[i] + 2.0 * ybuf[i+1] + ybuf[i+2]);
00189
00190 free((void *) ybuf);
00191 }
00192
00193
00194
00195 void median (float xbuf[], int npt, int nfilt) {
00196 float *ybuf, *array;
00197 float xmns, xmnf;
00198 int *point;
00199 int nfo2p1, i, il, ilow, j, jl, jh, nelem, l=0;
00200
00201 if((nfilt/2)*2 == nfilt) nfilt++;
00202 if(npt <= nfilt) return;
00203 nfo2p1 = nfilt/2;
00204
00205
00206 nelem = npt + nfilt;
00207 ybuf = (float *) malloc(nelem * sizeof(float));
00208
00209
00210 array = (float *) malloc(nfilt * sizeof(float));
00211 point = (int *) malloc(nfilt * sizeof(int));
00212
00213
00214
00215
00216 il = nfilt/2;
00217 ilow = MAX(3, nfilt/4);
00218 ilow = (ilow/2)*2 + 1;
00219
00220 for(i = 0; i < ilow; i++) array[i] = xbuf[i];
00221 sortm(array, point, ilow);
00222 xmns = array[ilow/2];
00223
00224 for(i = 0; i < ilow; i++) array[i] = xbuf[npt-1-i];
00225 sortm(array, point, ilow);
00226 xmnf = array[ilow/2];
00227
00228
00229 for(i = 0; i < il; i++) {
00230 ybuf[i] = 2.0 * xmns - xbuf[il+ilow-1-i];
00231 ybuf[npt+i+il] = 2.0 * xmnf - xbuf[npt-i-ilow-1];
00232 }
00233 for(i = 0; i < npt; i++) ybuf[i+il] = xbuf[i];
00234
00235
00236 for(i = 0; i < nfilt; i++) {
00237 array[i] = ybuf[i];
00238 point[i] = i+1;
00239 }
00240
00241 sortm(array, point, nfilt);
00242
00243 xbuf[0] = array[nfo2p1];
00244 jl = nfilt;
00245 jh = nfilt+npt-1;
00246 for(j = jl; j < jh; j++) {
00247
00248 for(i = 0; i < nfilt; i++) {
00249 if(point[i] != 1) {
00250 point[i]--;
00251 continue;
00252 }
00253 point[i] = nfilt;
00254 array[i] = ybuf[j];
00255 l = i;
00256 }
00257 quicksort(array, point, l, nfilt);
00258 xbuf[j-jl+1] = array[nfo2p1];
00259 }
00260
00261
00262 free((void *) point);
00263 free((void *) array);
00264 free((void *) ybuf);
00265 }
00266
00267 static void sortm (float ia[], int ib[], int n) {
00268 int i,j, ii, jj, ifin, iu;
00269 float it;
00270
00271 jj = 2;
00272 while(jj < n) jj = 2 * jj;
00273 jj = MIN(n,(3 * jj)/4 - 1);
00274 while(jj > 1) {
00275 jj = jj/2;
00276 ifin = n - jj;
00277 for(ii = 0; ii < ifin; ii++) {
00278 i = ii;
00279 j = i + jj;
00280 if(ia[i] <= ia[j]) continue;
00281 it = ia[j];
00282 iu = ib[j];
00283 do {
00284 ia[j] = ia[i];
00285 ib[j] = ib[i];
00286 j = i;
00287 i = i - jj;
00288 if (i < 0) break;
00289 } while(ia[i] > it);
00290 ia[j] = it;
00291 ib[j] = iu;
00292 }
00293 }
00294 }
00295
00296 static void quicksort (float x[], int point[], int l, int nfilt) {
00297 float test, temp;
00298 int i, it, j, npt, ii;
00299
00300 test = x[l];
00301 j = nfilt;
00302 for(i = 0; i < nfilt; i++) {
00303 if(i != l && test <= x[i]) {
00304 j = i;
00305 break;
00306 }
00307 }
00308 if(j - 1 == l) return;
00309
00310 if(j < l) {
00311 temp = x[l];
00312 it = point[l];
00313 npt = l - j;
00314 for(i = 0; i < npt; i++) {
00315 ii = l - i - 1;
00316 x[ii+1] = x[ii];
00317 point[ii+1] = point[ii];
00318 }
00319 x[j] = temp;
00320 point[j] = it;
00321 }
00322 else if(j > l) {
00323 temp = x[l];
00324 it = point[l];
00325 j--;
00326 npt = j - l;
00327 if(npt != 0) {
00328 for(i = 0; i < npt; i++) {
00329 ii = l + i + 1;
00330 x[ii-1] = x[ii];
00331 point[ii-1] = point[ii];
00332 }
00333 }
00334 x[j] = temp;
00335 point[j] = it;
00336 }
00337 }
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360