1 #include "afni.h"
2
3 /*--------------- Sample 0D function: log10 of each input point ------------*/
4
log10_func(int num,float * vec)5 void log10_func( int num , float *vec )
6 {
7 int ii ;
8 float vmax , vmin ;
9
10 if( num <= 0 || vec == NULL ) return ;
11
12 /* find largest element */
13
14 vmax = vec[0] ;
15 for( ii=1 ; ii < num ; ii++ ) vmax = MAX( vmax , vec[ii] ) ;
16
17 /* if all are nonpositive, return all zeros */
18
19 if( vmax <= 0.0 ){
20 for( ii=0 ; ii < num ; ii++ ) vec[ii] = 0.0 ;
21 return ;
22 }
23
24 /* find smallest positive element */
25
26 vmin = vmax ;
27 for( ii=0 ; ii < num ; ii++ )
28 if( vec[ii] > 0.0 ) vmin = MIN( vmin , vec[ii] ) ;
29
30 /* take log10 of each positive element;
31 nonpositive elements get the log10 of vmin instead */
32
33 vmin = log10(vmin) ;
34 for( ii=0 ; ii < num ; ii++ )
35 vec[ii] = (vec[ii] > 0.0) ? log10(vec[ii]) : vmin ;
36
37 return ;
38 }
39
40 /*------------ Sample 0D function: Signed sqrt of each input point ---------*/
41
ssqrt_func(int num,float * vec)42 void ssqrt_func( int num , float *vec )
43 {
44 int ii ;
45 double val ;
46
47 if( num <= 0 || vec == NULL ) return ;
48
49 for( ii=0 ; ii < num ; ii++ ){
50 val = sqrt(fabs(vec[ii])) ; /* will be positive */
51 vec[ii] = (vec[ii] >= 0.0) ? val : -val ; /* output sign = input sign */
52 }
53
54 return ;
55 }
56
57 /*------------ Sample 0D function: Absolute value of each input point --------*/
58
absval_func(int num,float * vec)59 void absval_func( int num , float *vec ) /* 20 Oct 2020 */
60 {
61 int ii ;
62 double val ;
63
64 if( num <= 0 || vec == NULL ) return ;
65
66 for( ii=0 ; ii < num ; ii++ )
67 vec[ii] = fabsf(vec[ii]) ;
68
69 return ;
70 }
71
72 /*--------------------------------------------------------------------------*/
73 /* adaptive = downweight things a long ways from median */
74
adaptive_weighted_mean(int num,float * x)75 float adaptive_weighted_mean( int num , float *x )
76 {
77 float med,mad, wt,wsum, xsum ; int ii ;
78
79 if( num <= 0 || x == NULL ) return (0.0f) ;
80 else if( num == 1 ) return (x[0]) ;
81 else if( num == 2 ) return (0.5f*(x[0]+x[1])) ;
82
83 qmedmad_float( num , x , &med , &mad ) ;
84 if( mad <= 0.0f ) return (med) ;
85
86 wsum = xsum = 0.0f ; mad = 0.56789f / mad ;
87 for( ii=0 ; ii < num ; ii++ ){
88 wt = mad*fabsf(x[ii]-med); wt = 1.0f / (1.0f+wt*wt*wt); wsum += wt;
89 xsum += wt * x[ii] ;
90 }
91 return (xsum/wsum) ;
92 }
93
94 /*--------------------------------------------------------------------------*/
95 /* 1D adaptive filter 9 points wide */
96
adpt_wt_mn9(int num,double to,double dt,float * vec)97 void adpt_wt_mn9( int num , double to,double dt, float *vec )
98 {
99 float x[9] , *nv ; int ii,jj,kk , n1=num-1 ;
100
101 nv = (float *)malloc(sizeof(float)*num) ;
102
103 for( ii=0 ; ii < num ; ii++ ){
104
105 for( jj=-4 ; jj <= 4 ; jj++ ){
106 kk = ii+jj ; if( kk < 0 ) kk = 0 ; else if( kk > n1 ) kk = n1 ;
107 x[jj+4] = vec[kk] ;
108 }
109
110 nv[ii] = adaptive_weighted_mean( 9 , x ) ;
111 }
112
113 memcpy(vec,nv,sizeof(float)*num) ; free(nv) ; return ;
114 }
115
116 /*--------------------------------------------------------------------------*/
117 /* 1D adaptive filter 19 points wide */
118
adpt_wt_mn19(int num,double to,double dt,float * vec)119 void adpt_wt_mn19( int num , double to,double dt, float *vec )
120 {
121 float x[19] , *nv ; int ii,jj,kk , n1=num-1 ;
122
123 nv = (float *)malloc(sizeof(float)*num) ;
124
125 for( ii=0 ; ii < num ; ii++ ){
126
127 for( jj=-9 ; jj <= 9 ; jj++ ){
128 kk = ii+jj ; if( kk < 0 ) kk = 0 ; else if( kk > n1 ) kk = n1 ;
129 x[jj+9] = vec[kk] ;
130 }
131
132 nv[ii] = adaptive_weighted_mean( 19 , x ) ;
133 }
134
135 memcpy(vec,nv,sizeof(float)*num) ; free(nv) ; return ;
136 }
137
138 /*--------------------------------------------------------------------------*/
139 /* 1D adaptive filter user-initialized points wide (must be odd) */
140
adpt_wt_mnXX(int num,double to,double dt,float * vec)141 void adpt_wt_mnXX( int num , double to,double dt, float *vec )
142 {
143 static int nXX=0,nHH=0 ; static float *XX=NULL ;
144 float *nv ; int ii,jj,kk , n1=num-1 ;
145
146 if( vec == NULL ){ /* setup call */
147 if( num > 0 && num != nXX ){
148 nXX = num ; nHH = nXX/2 ;
149 XX = (float *)realloc((void *)XX,sizeof(float)*(2*nHH+1)) ;
150 }
151 return ;
152 }
153
154 nv = (float *)malloc(sizeof(float)*num) ;
155
156 for( ii=0 ; ii < num ; ii++ ){
157
158 for( jj=-nHH ; jj <= nHH ; jj++ ){
159 kk = ii+jj ; if( kk < 0 ) kk = 0 ; else if( kk > n1 ) kk = n1 ;
160 XX[jj+nHH] = vec[kk] ;
161 }
162
163 nv[ii] = adaptive_weighted_mean( nXX , XX ) ;
164 }
165
166 memcpy(vec,nv,sizeof(float)*num) ; free(nv) ; return ;
167 }
168
169 /*--------------- Sample 1D function: Order Statistics Filter -------------*/
170
osfilt3_func(int num,double to,double dt,float * vec)171 void osfilt3_func( int num , double to,double dt, float *vec )
172 {
173 int ii ;
174 float aa,bb,cc ;
175
176 bb = vec[0] ; cc = vec[1] ;
177 for( ii=1 ; ii < num-1 ; ii++ ){
178 aa = bb ; bb = cc ; cc = vec[ii+1] ;
179 vec[ii] = OSFILT(aa,bb,cc) ; /* see mrilib.h */
180 }
181
182 return ;
183 }
184
185 /*--------------- Sample 1D function: Median of 3 Filter ----------------*/
186
median3_func(int num,double to,double dt,float * vec)187 void median3_func( int num , double to,double dt, float *vec )
188 {
189 int ii ;
190 float aa,bb,cc ;
191
192 bb = vec[0] ; cc = vec[1] ;
193 for( ii=1 ; ii < num-1 ; ii++ ){
194 aa = bb ; bb = cc ; cc = vec[ii+1] ;
195 vec[ii] = MEDIAN(aa,bb,cc) ; /* see mrilib.h */
196 }
197
198 return ;
199 }
200
201 /*-------- Sample 1D function: Despike7 Filter [07 Oct 2010] ----------*/
202
203 #define MTHR 6.789f /* threshold parameter */
204
205 #if 0
206 #undef mmm7
207 #define mmm7(j) \
208 { float qqq[7] ; int jj = (j)-3 ; \
209 if( jj < 0 ) jj = 0; else if( jj+7 > num ) jj = num-7; \
210 memcpy(qqq,vec+jj,sizeof(float)*7) ; \
211 med = qmed_float(7,qqq); qqq[0] = fabsf(qqq[0]-med); \
212 qqq[1] = fabsf(qqq[1]-med); qqq[2] = fabsf(qqq[2]-med); \
213 qqq[3] = fabsf(qqq[3]-med); qqq[4] = fabsf(qqq[4]-med); \
214 qqq[5] = fabsf(qqq[5]-med); qqq[6] = fabsf(qqq[6]-med); \
215 mad = qmed_float(7,qqq); }
216
217 void despike7_func( int num , double to,double dt , float *vec )
218 {
219 int ii ; float *zma,*zme , med,mad,val ;
220
221 if( num < 7 ) return ;
222 zme = (float *)malloc(sizeof(float)*num) ;
223 zma = (float *)malloc(sizeof(float)*num) ;
224
225 for( ii=0 ; ii < num ; ii++ ){
226 mmm7(ii) ; zme[ii] = med ; zma[ii] = mad ;
227 }
228 mad = qmed_float(num,zma) ; free(zma) ;
229 if( mad <= 0.0f ){ free(zme) ; return ; } /* should not happen */
230 mad *= MTHR ; /* threshold value */
231
232 for( ii=0 ; ii < num ; ii++ )
233 if( fabsf(vec[ii]-zme[ii]) > mad ) vec[ii] = zme[ii] ;
234
235 free(zme) ; return ;
236 }
237 #undef mmm7
238 #endif
239
240 /** this next func doesn't work well -- don't include it in AFNI **/
241 #if 0
242 void despike7pp_func( int num , double to,double dt , float *vec )
243 {
244 int ii ; float *dvv ;
245
246 despike7_func(num,to,dt,vec) ;
247 if( num < 9 ) return ;
248 dvv = malloc(sizeof(float)*num) ;
249 for( ii=1 ; ii < num ; ii++ ) dvv[ii-1] = vec[ii] -vec[ii-1] ;
250 despike7_func( num-1 , to,dt , dvv ) ;
251 for( ii=1 ; ii < num ; ii++ ) vec[ii] = vec[ii-1]+dvv[ii-1] ;
252 free(dvv) ; return ;
253 }
254 #endif
255
256 /*-------- Sample 1D function: Despike9 Filter [08 Oct 2010] ----------*/
257
258 #undef mmm9
259 #define mmm9(j) \
260 { float qqq[9] ; int jj = (j)-4 ; \
261 if( jj < 0 ) jj = 0; else if( jj+9 > num ) jj = num-9; \
262 memcpy(qqq,vec+jj,sizeof(float)*9) ; \
263 med = qmed_float(9,qqq); qqq[0] = fabsf(qqq[0]-med); \
264 qqq[1] = fabsf(qqq[1]-med); qqq[2] = fabsf(qqq[2]-med); \
265 qqq[3] = fabsf(qqq[3]-med); qqq[4] = fabsf(qqq[4]-med); \
266 qqq[5] = fabsf(qqq[5]-med); qqq[6] = fabsf(qqq[6]-med); \
267 qqq[7] = fabsf(qqq[7]-med); qqq[8] = fabsf(qqq[8]-med); \
268 mad = qmed_float(9,qqq); }
269
despike9_func(int num,double to,double dt,float * vec)270 void despike9_func( int num , double to,double dt , float *vec )
271 {
272 int ii ; float *zma,*zme , med,mad,val ;
273
274 if( num < 9 ) return ;
275 zme = (float *)malloc(sizeof(float)*num) ;
276 zma = (float *)malloc(sizeof(float)*num) ;
277
278 for( ii=0 ; ii < num ; ii++ ){
279 mmm9(ii) ; zme[ii] = med ; zma[ii] = mad ;
280 }
281 mad = qmed_float(num,zma) ; free(zma) ;
282 if( mad <= 0.0f ){ free(zme) ; return ; } /* should not happen */
283 mad *= MTHR ; /* threshold value */
284
285 for( ii=0 ; ii < num ; ii++ )
286 if( fabsf(vec[ii]-zme[ii]) > mad ) vec[ii] = zme[ii] ;
287
288 free(zme) ; return ;
289 }
290 #undef mmm9
291
292 /*----------------------------------------------------------------------------*/
293 /* Despiking-25 filter */
294
295 #undef SWAP
296 #define SWAP(x,y) (temp=x,x=y,y=temp)
297 #undef SORT2
298 #define SORT2(a,b) if(a>b) SWAP(a,b)
299
median25f(float * p)300 static INLINE float median25f(float *p)
301 {
302 register float temp ;
303 SORT2(p[0], p[1]) ; SORT2(p[3], p[4]) ; SORT2(p[2], p[4]) ;
304 SORT2(p[2], p[3]) ; SORT2(p[6], p[7]) ; SORT2(p[5], p[7]) ;
305 SORT2(p[5], p[6]) ; SORT2(p[9], p[10]) ; SORT2(p[8], p[10]) ;
306 SORT2(p[8], p[9]) ; SORT2(p[12], p[13]) ; SORT2(p[11], p[13]) ;
307 SORT2(p[11], p[12]) ; SORT2(p[15], p[16]) ; SORT2(p[14], p[16]) ;
308 SORT2(p[14], p[15]) ; SORT2(p[18], p[19]) ; SORT2(p[17], p[19]) ;
309 SORT2(p[17], p[18]) ; SORT2(p[21], p[22]) ; SORT2(p[20], p[22]) ;
310 SORT2(p[20], p[21]) ; SORT2(p[23], p[24]) ; SORT2(p[2], p[5]) ;
311 SORT2(p[3], p[6]) ; SORT2(p[0], p[6]) ; SORT2(p[0], p[3]) ;
312 SORT2(p[4], p[7]) ; SORT2(p[1], p[7]) ; SORT2(p[1], p[4]) ;
313 SORT2(p[11], p[14]) ; SORT2(p[8], p[14]) ; SORT2(p[8], p[11]) ;
314 SORT2(p[12], p[15]) ; SORT2(p[9], p[15]) ; SORT2(p[9], p[12]) ;
315 SORT2(p[13], p[16]) ; SORT2(p[10], p[16]) ; SORT2(p[10], p[13]) ;
316 SORT2(p[20], p[23]) ; SORT2(p[17], p[23]) ; SORT2(p[17], p[20]) ;
317 SORT2(p[21], p[24]) ; SORT2(p[18], p[24]) ; SORT2(p[18], p[21]) ;
318 SORT2(p[19], p[22]) ; SORT2(p[8], p[17]) ; SORT2(p[9], p[18]) ;
319 SORT2(p[0], p[18]) ; SORT2(p[0], p[9]) ; SORT2(p[10], p[19]) ;
320 SORT2(p[1], p[19]) ; SORT2(p[1], p[10]) ; SORT2(p[11], p[20]) ;
321 SORT2(p[2], p[20]) ; SORT2(p[2], p[11]) ; SORT2(p[12], p[21]) ;
322 SORT2(p[3], p[21]) ; SORT2(p[3], p[12]) ; SORT2(p[13], p[22]) ;
323 SORT2(p[4], p[22]) ; SORT2(p[4], p[13]) ; SORT2(p[14], p[23]) ;
324 SORT2(p[5], p[23]) ; SORT2(p[5], p[14]) ; SORT2(p[15], p[24]) ;
325 SORT2(p[6], p[24]) ; SORT2(p[6], p[15]) ; SORT2(p[7], p[16]) ;
326 SORT2(p[7], p[19]) ; SORT2(p[13], p[21]) ; SORT2(p[15], p[23]) ;
327 SORT2(p[7], p[13]) ; SORT2(p[7], p[15]) ; SORT2(p[1], p[9]) ;
328 SORT2(p[3], p[11]) ; SORT2(p[5], p[17]) ; SORT2(p[11], p[17]) ;
329 SORT2(p[9], p[17]) ; SORT2(p[4], p[10]) ; SORT2(p[6], p[12]) ;
330 SORT2(p[7], p[14]) ; SORT2(p[4], p[6]) ; SORT2(p[4], p[7]) ;
331 SORT2(p[12], p[14]) ; SORT2(p[10], p[14]) ; SORT2(p[6], p[7]) ;
332 SORT2(p[10], p[12]) ; SORT2(p[6], p[10]) ; SORT2(p[6], p[17]) ;
333 SORT2(p[12], p[17]) ; SORT2(p[7], p[17]) ; SORT2(p[7], p[10]) ;
334 SORT2(p[12], p[18]) ; SORT2(p[7], p[12]) ; SORT2(p[10], p[18]) ;
335 SORT2(p[12], p[20]) ; SORT2(p[10], p[20]) ; SORT2(p[10], p[12]) ;
336 return (p[12]);
337 }
338
339 /*--- get the local median and MAD of values vec[j-12 .. j+12] ---*/
340
341 #undef mead25
342 #define mead25(j) \
343 { float qqq[25] ; int jj=(j)-12 ; register int pp; \
344 if( jj < 0 ) jj = 0; else if( jj+25 > num ) jj = num-25; \
345 for( pp=0 ; pp < 25 ; pp++ ) qqq[pp] = vec[jj+pp] ; \
346 med = median25f(qqq) ; \
347 for( pp=0 ; pp < 25 ; pp++ ) qqq[pp] = fabsf(qqq[pp]-med) ; \
348 mad = median25f(qqq); }
349
DES_despike25(int num,double to,double dt,float * vec)350 void DES_despike25( int num , double to,double dt , float *vec )
351 {
352 int ii , nsp ; float *zma,*zme , med,mad,val ;
353 static float *deswks = NULL ;
354 static int ndeswks = 0 ;
355
356 if( vec == NULL ) return ;
357 if( num < 25 ) { despike9_func( num , to,dt , vec ) ; return ; }
358
359 if( ndeswks < num ){
360 deswks = (float *)realloc(deswks,sizeof(float)*(4*num)) ;
361 ndeswks = num ;
362 }
363
364 zme = deswks ; zma = zme + num ;
365
366 for( ii=0 ; ii < num ; ii++ ){
367 mead25(ii) ; zme[ii] = med ; zma[ii] = mad ;
368 }
369 mad = qmed_float(num,zma) ;
370 if( mad <= 0.0f ) return ;
371 mad *= 6.1234f ; /* threshold value */
372
373 for( nsp=ii=0 ; ii < num ; ii++ )
374 if( fabsf(vec[ii]-zme[ii]) > mad ){ vec[ii] = zme[ii]; nsp++; }
375
376 return ;
377 }
378 #undef mead25
379 #undef SORT2
380 #undef SWAP
381
382 /*-------------- Sample 1D function: HRF Decon [29 Oct 2010] ---------------*/
383
lhs_legendre(float x,float bot,float top,float n)384 static float lhs_legendre( float x, float bot, float top, float n )
385 {
386 double xx ;
387 xx = 2.0*(x-bot)/(top-bot) - 1.0 ; /* now in range -1..1 */
388 return (float)Plegendre(xx,n) ;
389 }
390
hrfdecon_func(int num,double to,double dt,float * vec)391 void hrfdecon_func( int num , double to,double dt , float *vec )
392 {
393 int nkern , nbase , ii,pp ;
394 float *kern , **base , ksum=0.0f ;
395 floatvec *bfit ;
396
397 if( num < 49 || dt >= 4.0f ) return ;
398 if( dt <= 0.0f ) dt = 1.0f ;
399
400 nkern = 1+(int)(16.0f/dt) ; if( nkern >= num/2 ) return ;
401 kern = (float *)calloc(sizeof(float),nkern) ;
402 for( ii=1 ; ii < nkern ; ii++ ){
403 kern[ii] = powf(ii*dt,8.6f)*exp(-ii*dt/0.547) ; ksum += kern[ii] ;
404 }
405 for( ii=1 ; ii < nkern ; ii++ ) kern[ii] /= ksum ;
406 for( ii=nkern-1 ; ii > 2 ; ii-- ) if( kern[ii] >= 0.0444f ) break ;
407 nkern = ii ;
408
409 nbase = 1+(int)(num*dt/100.0f) ; if( nbase < 3 ) nbase = 3 ;
410 base = (float **)malloc(sizeof(float *)*nbase) ;
411 for( pp=0 ; pp < nbase ; pp++ ){
412 base[pp] = (float *)malloc(sizeof(float)*num) ;
413 for( ii=0 ; ii < num ; ii++ )
414 base[pp][ii] = lhs_legendre( (float)ii, 0.0f, num-1.0f, pp ) ;
415 }
416
417 despike9_func( num , to,dt , vec ) ;
418
419 { float *vlam = (float *)malloc(sizeof(float)*(nbase+num)) ;
420 for( ii=0 ; ii < num ; ii++ ) vlam[ii] = 9.999f ;
421 for( ii=0 ; ii < nbase ; ii++ ) vlam[ii+num] = 0.0f ;
422 THD_lasso_fixlam(9.999f) ;
423 THD_lasso_setlamvec( nbase+num , vlam ) ; free(vlam) ;
424 }
425
426 bfit = THD_deconvolve( num , vec ,
427 0 , nkern-1 , kern ,
428 nbase , base , -1 , NULL , 0 , 7 , -3.333f ) ;
429
430 if( bfit != NULL ){
431 memcpy(vec,bfit->ar,sizeof(float)*num) ;
432 KILL_floatvec(bfit) ;
433 }
434
435 for( pp=0 ; pp < nbase ; pp++ ) free(base[pp]) ;
436 free(base) ; free(kern) ; return ;
437 }
438
439 /*---------------- Sample 1D function: abs(FFT) [30 Jun 2000] --------------*/
440
absfft_func(int num,double to,double dt,float * vec)441 void absfft_func( int num , double to,double dt, float *vec )
442 {
443 static complex *cx=NULL ;
444 static int ncx=0 , numold=0 ;
445 float f0,f1,f2 ;
446 int ii ;
447
448 if( num < 2 ) return ;
449 if( num != numold ){
450 numold = num ;
451 ncx = csfft_nextup_even(numold) ;
452 cx = (complex *)realloc(cx,sizeof(complex)*ncx) ;
453 INFO_message("1D FFT: ndata=%d nfft=%d",num,ncx) ;
454 }
455
456 get_quadratic_trend( num , vec , &f0,&f1,&f2 ) ; /* thd_detrend.c */
457
458 for( ii=0 ; ii < num ; ii++ ){ cx[ii].r = vec[ii]-(f0+f1*ii+f2*ii*ii); cx[ii].i = 0.0; }
459 for( ; ii < ncx ; ii++ ){ cx[ii].r = cx[ii].i = 0.0 ; }
460
461 csfft_cox( -1 , ncx , cx ) ; /* csfft.c */
462
463 vec[0] = 0.0 ;
464 for( ii=1 ; ii < num ; ii++ ) vec[ii] = CABS(cx[ii]) ;
465
466 return ;
467 }
468
469 /*---------------- Sample 1D function: 0..1 scaling [02 Sep 2009] ----------*/
470
ztone_func(int num,double to,double dt,float * vec)471 void ztone_func( int num , double to,double dt, float *vec )
472 {
473 int ii ; float vbot,vtop ;
474
475 if( num < 2 || vec == NULL ) return ;
476 vbot = vtop = vec[0] ;
477 for( ii=1 ; ii < num ; ii++ ){
478 if( vec[ii] < vbot ) vbot = vec[ii] ;
479 else if( vec[ii] > vtop ) vtop = vec[ii] ;
480 }
481 if( vbot == vtop ){
482 for( ii=0 ; ii < num ; ii++ ) vec[ii] = 0.5f ;
483 } else {
484 float fac = 1.0f / (vtop-vbot) ;
485 for( ii=0 ; ii < num ; ii++ ) vec[ii] = (vec[ii]-vbot)*fac ;
486 }
487
488 return ;
489 }
490
491 /*---------------- Sample 1D function: L1 normalize [03 Sep 2009] ----------*/
492
L1normalize_func(int num,double to,double dt,float * vec)493 void L1normalize_func( int num , double to,double dt, float *vec )
494 {
495 int ii ; float vsum ;
496
497 if( num < 2 || vec == NULL ) return ;
498
499 vsum = 0.0f ;
500 for( ii=0 ; ii < num ; ii++ ) vsum += fabsf(vec[ii]) ;
501 if( vsum == 0.0f ) return ;
502 vsum = 1.0f / vsum ;
503 for( ii=0 ; ii < num ; ii++ ) vec[ii] *= vsum ;
504 return ;
505 }
506
507 /*---------------- Sample 1D function: L2 normalize [03 Sep 2009] ----------*/
508
L2normalize_func(int num,double to,double dt,float * vec)509 void L2normalize_func( int num , double to,double dt, float *vec )
510 {
511 int ii ; float vsum ;
512
513 if( num < 2 || vec == NULL ) return ;
514
515 vsum = 0.0f ;
516 for( ii=0 ; ii < num ; ii++ ) vsum += vec[ii]*vec[ii] ;
517 if( vsum <= 0.0f ) return ;
518 vsum = 1.0f / sqrtf(vsum) ;
519 for( ii=0 ; ii < num ; ii++ ) vec[ii] *= vsum ;
520 return ;
521 }
522
523 /*----------- Sample slice projection functions [31 Jan 2002] -----------*/
524
min_proj(int n,float * ar)525 float min_proj( int n , float *ar )
526 {
527 float v = ar[0] ;
528 int ii ;
529 for( ii=1 ; ii < n ; ii++ ) if( v > ar[ii] ) v = ar[ii] ;
530 return v ;
531 }
532
max_proj(int n,float * ar)533 float max_proj( int n , float *ar )
534 {
535 float v = ar[0] ;
536 int ii ;
537 for( ii=1 ; ii < n ; ii++ ) if( v < ar[ii] ) v = ar[ii] ;
538 return v ;
539 }
540
mean_proj(int n,float * ar)541 float mean_proj( int n , float *ar )
542 {
543 float v=0.0 ;
544 int ii ;
545 for( ii=0 ; ii < n ; ii++ ) v += ar[ii] ;
546 return (v/n) ;
547 }
548
extreme_proj(int n,float * ar)549 float extreme_proj( int n , float *ar ) /* 02 Feb 2002 */
550 {
551 float vv,ww , med=qmed_float(n,ar) ;
552 int ii , jj ;
553
554 jj = 0 ; vv = fabs(ar[0]-med) ; /* Find the value */
555 for( ii=1 ; ii < n ; ii++ ){ /* furthest from */
556 ww = fabs(ar[ii]-med) ; /* the median. */
557 if( ww > vv ){ vv=ww; jj=ii; }
558 }
559 return ar[jj] ;
560 }
561
osfilt_proj(int n,float * ar)562 float osfilt_proj( int n , float *ar ) /* 07 Dec 2007 */
563 {
564 int ii , n2 ; float v=0.0f , d=0.0f ;
565 qsort_float( n , ar ) ;
566 n2 = n/2 ;
567 for( ii=0 ; ii < n2 ; ii++ ){
568 v += (ii+1.0f)*(ar[ii]+ar[n-1-ii]) ;
569 d += 2.0f*(ii+1.0f) ;
570 }
571 v += (n2+1.0f)*ar[n2] ; d += (n2+1.0f) ;
572 return (v/d) ;
573 }
574
mad_proj(int n,float * ar)575 float mad_proj( int n , float *ar ) /* 07 Dec 2007 */
576 {
577 float v ;
578 qmedmad_float( n , ar , NULL , &v ) ; return v ;
579 }
580
581 /*======================================================================*/
582 /*----------------------- Sample 2D transformations --------------------*/
583 /*======================================================================*/
584
585 static float *atemp = NULL ;
586 static int natemp = -666 ;
587
588 #define MAKE_ATEMP(nvox) \
589 do{ if( natemp < (nvox) ){ \
590 if( atemp != NULL ) free(atemp) ; \
591 natemp = (nvox) ; \
592 atemp = (float *) malloc( sizeof(float) * natemp ) ; } } while(0)
593
594 #define AT(i,j) atemp[(i)+(j)*nx]
595
596 /*------------------------------------------------------------------------*/
597
median9_box_func(int nx,int ny,double dx,double dy,float * ar)598 void median9_box_func( int nx , int ny , double dx, double dy, float *ar )
599 {
600 int ii , jj , nxy , joff ;
601 float aa[9] ;
602 float *ajj , *ajm , *ajp ;
603
604 if( nx < 3 || ny < 3 ) return ;
605
606 /** make space and copy input into it **/
607
608 nxy = nx * ny ;
609 MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
610 memcpy(atemp,ar,sizeof(float)*nxy) ;
611
612 /** process copy of input back into the input array **/
613
614 for( jj=0 ; jj < ny ; jj++ ){
615
616 joff = jj * nx ; /* offset into this row */
617 ajj = atemp + joff ; /* pointer to this row */
618
619 ajm = (jj==0 ) ? ajj : ajj-nx ; /* pointer to last row */
620 ajp = (jj==ny-1) ? ajj : ajj+nx ; /* pointer to next row */
621
622 /* do interior points of this row */
623
624 for( ii=1 ; ii < nx-1 ; ii++ ){
625 aa[0] = ajm[ii-1] ; aa[1] = ajm[ii] ; aa[2] = ajm[ii+1] ;
626 aa[3] = ajj[ii-1] ; aa[4] = ajj[ii] ; aa[5] = ajj[ii+1] ;
627 aa[6] = ajp[ii-1] ; aa[7] = ajp[ii] ; aa[8] = ajp[ii+1] ;
628 #if 0
629 isort_float( 9 , aa ) ; ar[ii+joff] = aa[4] ;
630 #else
631 ar[ii+joff] = qmed_float(9,aa) ; /* 25 Oct 2000 */
632 #endif
633 }
634
635 /* do leading edge point (ii=0) */
636
637 aa[0] = ajm[0] ; aa[1] = ajm[0] ; aa[2] = ajm[1] ;
638 aa[3] = ajj[0] ; aa[4] = ajj[0] ; aa[5] = ajj[1] ;
639 aa[6] = ajp[0] ; aa[7] = ajp[0] ; aa[8] = ajp[1] ;
640 #if 0
641 isort_float( 9 , aa ) ; ar[joff] = aa[4] ;
642 #else
643 ar[joff] = qmed_float(9,aa) ; /* 25 Oct 2000 */
644 #endif
645
646 /* do trailing edge point (ii=nx-1) */
647
648 aa[0] = ajm[nx-2] ; aa[1] = ajm[nx-1] ; aa[2] = ajm[nx-1] ;
649 aa[3] = ajj[nx-2] ; aa[4] = ajj[nx-1] ; aa[5] = ajj[nx-1] ;
650 aa[6] = ajp[nx-2] ; aa[7] = ajp[nx-1] ; aa[8] = ajp[nx-1] ;
651 #if 0
652 isort_float( 9 , aa ) ; ar[nx-1+joff] = aa[4] ;
653 #else
654 ar[nx-1+joff] = qmed_float(9,aa) ; /* 25 Oct 2000 */
655 #endif
656 }
657 return ;
658 }
659
660 /*------------------------------------------------------------------------*/
661
winsor9_box_func(int nx,int ny,double dx,double dy,float * ar)662 void winsor9_box_func( int nx , int ny , double dx, double dy, float *ar )
663 {
664 int ii , jj , nxy , joff ;
665 float aa[9] ;
666 float *ajj , *ajm , *ajp ;
667
668 if( nx < 3 || ny < 3 ) return ;
669
670 /** make space and copy input into it **/
671
672 nxy = nx * ny ;
673 MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
674 memcpy(atemp,ar,sizeof(float)*nxy) ;
675
676 /** process copy of input back into the input array **/
677
678 for( jj=0 ; jj < ny ; jj++ ){
679
680 joff = jj * nx ; /* offset into this row */
681 ajj = atemp + joff ; /* pointer to this row */
682
683 ajm = (jj==0 ) ? ajj : ajj-nx ; /* pointer to last row */
684 ajp = (jj==ny-1) ? ajj : ajj+nx ; /* pointer to next row */
685
686 /* do interior points of this row */
687
688 for( ii=1 ; ii < nx-1 ; ii++ ){
689 aa[0] = ajm[ii-1] ; aa[1] = ajm[ii] ; aa[2] = ajm[ii+1] ;
690 aa[3] = ajj[ii-1] ; aa[4] = ajj[ii] ; aa[5] = ajj[ii+1] ;
691 aa[6] = ajp[ii-1] ; aa[7] = ajp[ii] ; aa[8] = ajp[ii+1] ;
692 isort_float( 9 , aa ) ;
693 if( ar[ii+joff] < aa[2] ) ar[ii+joff] = aa[2] ;
694 else if( ar[ii+joff] > aa[6] ) ar[ii+joff] = aa[6] ;
695 }
696
697 /* do leading edge point (ii=0) */
698
699 aa[0] = ajm[0] ; aa[1] = ajm[0] ; aa[2] = ajm[1] ;
700 aa[3] = ajj[0] ; aa[4] = ajj[0] ; aa[5] = ajj[1] ;
701 aa[6] = ajp[0] ; aa[7] = ajp[0] ; aa[8] = ajp[1] ;
702 isort_float( 9 , aa ) ;
703 if( ar[joff] < aa[2] ) ar[joff] = aa[2] ;
704 else if( ar[joff] > aa[6] ) ar[joff] = aa[6] ;
705
706 /* do trailing edge point (ii=nx-1) */
707
708 aa[0] = ajm[nx-2] ; aa[1] = ajm[nx-1] ; aa[2] = ajm[nx-1] ;
709 aa[3] = ajj[nx-2] ; aa[4] = ajj[nx-1] ; aa[5] = ajj[nx-1] ;
710 aa[6] = ajp[nx-2] ; aa[7] = ajp[nx-1] ; aa[8] = ajp[nx-1] ;
711 isort_float( 9 , aa ) ;
712 if( ar[nx-1+joff] < aa[2] ) ar[nx-1+joff] = aa[2] ;
713 else if( ar[nx-1+joff] > aa[6] ) ar[nx-1+joff] = aa[6] ;
714 }
715 return ;
716 }
717
718 /*------------------------------------------------------------------------*/
719
osfilt9_box_func(int nx,int ny,double dx,double dy,float * ar)720 void osfilt9_box_func( int nx , int ny , double dx, double dy, float *ar )
721 {
722 int ii , jj , nxy , joff ;
723 float aa[9] ;
724 float *ajj , *ajm , *ajp ;
725
726 if( nx < 3 || ny < 3 ) return ;
727
728 /** make space and copy input into it **/
729
730 nxy = nx * ny ;
731 MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
732 memcpy(atemp,ar,sizeof(float)*nxy) ;
733
734 /** process copy of input back into the input array **/
735
736 for( jj=0 ; jj < ny ; jj++ ){
737
738 joff = jj * nx ; /* offset into this row */
739 ajj = atemp + joff ; /* pointer to this row */
740
741 ajm = (jj==0 ) ? ajj : ajj-nx ; /* pointer to last row */
742 ajp = (jj==ny-1) ? ajj : ajj+nx ; /* pointer to next row */
743
744 /* do interior points of this row */
745
746 #undef OSUM
747 #define OSUM(a,b,c,d,e) ( 0.1*((a)+(e)) + 0.2*((b)+(d)) + 0.4*(c) )
748
749 for( ii=1 ; ii < nx-1 ; ii++ ){
750 aa[0] = ajm[ii-1] ; aa[1] = ajm[ii] ; aa[2] = ajm[ii+1] ;
751 aa[3] = ajj[ii-1] ; aa[4] = ajj[ii] ; aa[5] = ajj[ii+1] ;
752 aa[6] = ajp[ii-1] ; aa[7] = ajp[ii] ; aa[8] = ajp[ii+1] ;
753 isort_float( 9 , aa ) ;
754 ar[ii+joff] = OSUM( aa[2],aa[3],aa[4],aa[5],aa[6] ) ;
755 }
756
757 /* do leading edge point (ii=0) */
758
759 aa[0] = ajm[0] ; aa[1] = ajm[0] ; aa[2] = ajm[1] ;
760 aa[3] = ajj[0] ; aa[4] = ajj[0] ; aa[5] = ajj[1] ;
761 aa[6] = ajp[0] ; aa[7] = ajp[0] ; aa[8] = ajp[1] ;
762 isort_float( 9 , aa ) ;
763 ar[joff] = OSUM( aa[2],aa[3],aa[4],aa[5],aa[6] ) ;
764
765 /* do trailing edge point (ii=nx-1) */
766
767 aa[0] = ajm[nx-2] ; aa[1] = ajm[nx-1] ; aa[2] = ajm[nx-1] ;
768 aa[3] = ajj[nx-2] ; aa[4] = ajj[nx-1] ; aa[5] = ajj[nx-1] ;
769 aa[6] = ajp[nx-2] ; aa[7] = ajp[nx-1] ; aa[8] = ajp[nx-1] ;
770 isort_float( 9 , aa ) ;
771 ar[nx-1+joff] = OSUM( aa[2],aa[3],aa[4],aa[5],aa[6] ) ;
772 }
773 return ;
774 }
775
776 /*------------------------------------------------------------------------*/
777
median21_box_func(int nx,int ny,double dx,double dy,float * ar)778 void median21_box_func( int nx , int ny , double dx, double dy, float *ar )
779 {
780 int ii , jj , nxy , joff ;
781 float aa[21] ;
782 float *ajj , *ajm , *ajp , *ajmm , *ajpp ;
783
784 if( nx < 5 || ny < 5 ) return ;
785
786 /** make space and copy input into it **/
787
788 nxy = nx * ny ;
789 MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
790 memcpy(atemp,ar,sizeof(float)*nxy) ;
791
792 /** process copy of input back into the input array **/
793
794 for( jj=1 ; jj < ny-1 ; jj++ ){
795
796 joff = jj * nx ; /* offset into this row */
797 ajj = atemp + joff ; /* pointer to this row */
798
799 ajm = ajj-nx ; /* pointer to last row */
800 ajp = ajj+nx ; /* pointer to next row */
801
802 ajmm = (jj == 1 ) ? ajm : ajm-nx ; /* to last last row */
803 ajpp = (jj ==ny-2) ? ajp : ajp+nx ; /* to next next row */
804
805 /* do interior points of this row */
806
807 for( ii=2 ; ii < nx-2 ; ii++ ){
808 aa[0]=ajmm[ii-1]; aa[1]=ajmm[ii]; aa[2]=ajmm[ii+1];
809
810 aa[ 3]=ajm[ii-2]; aa[ 4]=ajm[ii-1]; aa[ 5]=ajm[ii]; aa[ 6]=ajm[ii+1]; aa[ 7]=ajm[ii+2];
811 aa[ 8]=ajj[ii-2]; aa[ 9]=ajj[ii-1]; aa[10]=ajj[ii]; aa[11]=ajj[ii+1]; aa[12]=ajj[ii+2];
812 aa[13]=ajp[ii-2]; aa[14]=ajp[ii-1]; aa[15]=ajp[ii]; aa[16]=ajp[ii+1]; aa[17]=ajp[ii+2];
813
814 aa[18]=ajpp[ii-1]; aa[19]=ajpp[ii]; aa[20]=ajpp[ii+1];
815
816 #if 0
817 isort_float( 21 , aa ) ; ar[ii+joff] = aa[10] ;
818 #else
819 ar[ii+joff] = qmed_float(21,aa) ; /* 25 Oct 2000 */
820 #endif
821 }
822
823 }
824 return ;
825 }
826
827 /*------------------------------------------------------------------------*/
828
winsor21_box_func(int nx,int ny,double dx,double dy,float * ar)829 void winsor21_box_func( int nx , int ny , double dx, double dy, float *ar )
830 {
831 int ii , jj , nxy , joff ;
832 float aa[21] ;
833 float *ajj , *ajm , *ajp , *ajmm , *ajpp ;
834
835 static int kbot=-1 , ktop ;
836
837 if( nx < 5 || ny < 5 ) return ;
838
839 /** initialize cutoffs [07 Dec 1999] **/
840
841 if( kbot < 0 ){
842 char *ee = my_getenv("AFNI_WINSOR21_CUTOFF") ;
843 kbot = 6 ; /* default */
844 if( ee != NULL ){
845 ii = strtol( ee , NULL , 10 ) ;
846 if( ii > 0 && ii < 10 ) kbot = ii ;
847 }
848 ktop = 20 - kbot ;
849 }
850
851 /** make space and copy input into it **/
852
853 nxy = nx * ny ;
854 MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
855 memcpy(atemp,ar,sizeof(float)*nxy) ;
856
857 /** process copy of input back into the input array **/
858
859 for( jj=1 ; jj < ny-1 ; jj++ ){
860
861 joff = jj * nx ; /* offset into this row */
862 ajj = atemp + joff ; /* pointer to this row */
863
864 ajm = ajj-nx ; /* pointer to last row */
865 ajp = ajj+nx ; /* pointer to next row */
866
867 ajmm = (jj == 1 ) ? ajm : ajm-nx ; /* to last last row */
868 ajpp = (jj ==ny-2) ? ajp : ajp+nx ; /* to next next row */
869
870 /* do interior points of this row */
871
872 for( ii=2 ; ii < nx-2 ; ii++ ){
873 aa[0]=ajmm[ii-1]; aa[1]=ajmm[ii]; aa[2]=ajmm[ii+1];
874
875 aa[ 3]=ajm[ii-2]; aa[ 4]=ajm[ii-1]; aa[ 5]=ajm[ii]; aa[ 6]=ajm[ii+1]; aa[ 7]=ajm[ii+2];
876 aa[ 8]=ajj[ii-2]; aa[ 9]=ajj[ii-1]; aa[10]=ajj[ii]; aa[11]=ajj[ii+1]; aa[12]=ajj[ii+2];
877 aa[13]=ajp[ii-2]; aa[14]=ajp[ii-1]; aa[15]=ajp[ii]; aa[16]=ajp[ii+1]; aa[17]=ajp[ii+2];
878
879 aa[18]=ajpp[ii-1]; aa[19]=ajpp[ii]; aa[20]=ajpp[ii+1];
880
881 isort_float( 21 , aa ) ;
882
883 if( ar[ii+joff] < aa[kbot] ) ar[ii+joff] = aa[kbot] ;
884 else if( ar[ii+joff] > aa[ktop] ) ar[ii+joff] = aa[ktop] ;
885 }
886
887 }
888 return ;
889 }
890
891 /*------------------------------------------------------------------------*/
892
adapt_mean_21_box_func(int nx,int ny,double dx,double dy,float * ar)893 void adapt_mean_21_box_func( int nx, int ny, double dx, double dy, float *ar )
894 {
895 int ii , jj , nxy , joff ;
896 float aa[21] ;
897 float *ajj , *ajm , *ajp , *ajmm , *ajpp ;
898
899 if( nx < 5 || ny < 5 ) return ;
900
901 /** make space and copy input into it **/
902
903 nxy = nx * ny ;
904 MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
905 memcpy(atemp,ar,sizeof(float)*nxy) ;
906
907 /** process copy of input back into the input array **/
908
909 for( jj=1 ; jj < ny-1 ; jj++ ){
910
911 joff = jj * nx ; /* offset into this row */
912 ajj = atemp + joff ; /* pointer to this row */
913
914 ajm = ajj-nx ; /* pointer to last row */
915 ajp = ajj+nx ; /* pointer to next row */
916
917 ajmm = (jj == 1 ) ? ajm : ajm-nx ; /* to last last row */
918 ajpp = (jj ==ny-2) ? ajp : ajp+nx ; /* to next next row */
919
920 /* do interior points of this row */
921
922 for( ii=2 ; ii < nx-2 ; ii++ ){
923 aa[0]=ajmm[ii-1]; aa[1]=ajmm[ii]; aa[2]=ajmm[ii+1];
924
925 aa[ 3]=ajm[ii-2]; aa[ 4]=ajm[ii-1]; aa[ 5]=ajm[ii]; aa[ 6]=ajm[ii+1]; aa[ 7]=ajm[ii+2];
926 aa[ 8]=ajj[ii-2]; aa[ 9]=ajj[ii-1]; aa[10]=ajj[ii]; aa[11]=ajj[ii+1]; aa[12]=ajj[ii+2];
927 aa[13]=ajp[ii-2]; aa[14]=ajp[ii-1]; aa[15]=ajp[ii]; aa[16]=ajp[ii+1]; aa[17]=ajp[ii+2];
928
929 aa[18]=ajpp[ii-1]; aa[19]=ajpp[ii]; aa[20]=ajpp[ii+1];
930
931 ar[ii+joff] = adaptive_weighted_mean( 21 , aa ) ;
932 }
933 }
934 return ;
935 }
936
937 /*------- [30 Jun 2000: abs(2D FFT) function] ----------------------------*/
938
fft2D_absfunc(int nx,int ny,double dx,double dy,float * ar)939 void fft2D_absfunc( int nx , int ny , double dx, double dy, float *ar )
940 {
941 complex *cxar , *cpt ;
942 int nxup,nyup , ii,jj ;
943 float fi,fj , *fpt ;
944
945 if( nx < 5 || ny < 5 ) return ;
946
947 nxup = csfft_nextup_even(nx) ; /* get FFT size */
948 nyup = csfft_nextup_even(ny) ;
949
950 cxar = (complex *) malloc(sizeof(complex)*nxup*nyup) ;
951
952 /* copy input to output, sign-alternating and zero-padding along the way */
953
954 cpt = cxar ;
955 fpt = ar ;
956 fj = 1.0 ;
957 for( jj=0 ; jj < ny ; jj++ ){
958 fi = fj ; fj = -fj ;
959 for(ii=0; ii<nx ; ii++){cpt->r=*fpt*fi; cpt->i=0.0; cpt++;fpt++;fi=-fi;}
960 for( ; ii<nxup; ii++){cpt->r=cpt->i=0.0; cpt++;}
961 }
962 for( ; jj < nyup ; jj++ ){cpt->r=cpt->i=0.0; cpt++;}
963
964 /* row FFTs */
965
966 for( jj=0 ; jj < ny ; jj++ )
967 csfft_cox( -1 , nxup , cxar+jj*nxup ) ;
968
969 /* column FFTs */
970
971 cpt = (complex *) malloc(sizeof(complex)*nyup) ;
972
973 for( ii=0 ; ii < nxup ; ii++ ){
974 for( jj=0 ; jj < nyup ; jj++ ) cpt[jj] = cxar[ii+jj*nxup] ;
975 csfft_cox( -1 , nyup , cpt ) ;
976 for( jj=0 ; jj < nyup ; jj++ ) cxar[ii+jj*nxup] = cpt[jj] ;
977 }
978
979 /* copy to output */
980
981 for( jj=0 ; jj < ny ; jj++ )
982 for( ii=0 ; ii < nx ; ii++ )
983 ar[ii+jj*nx] = CABS(cxar[ii+jj*nxup]) ;
984
985 free(cxar) ; free(cpt) ; return ;
986 }
987
988 /*------- [26 Mar 2008: arg(2D FFT) function] ----------------------------*/
989
fft2D_phasefunc(int nx,int ny,double dx,double dy,float * ar)990 void fft2D_phasefunc( int nx , int ny , double dx, double dy, float *ar )
991 {
992 complex *cxar , *cpt ;
993 int nxup,nyup , ii,jj ;
994 float fi,fj , *fpt ;
995
996 if( nx < 5 || ny < 5 ) return ;
997
998 nxup = csfft_nextup_even(nx) ; /* get FFT size */
999 nyup = csfft_nextup_even(ny) ;
1000
1001 cxar = (complex *) malloc(sizeof(complex)*nxup*nyup) ;
1002
1003 /* copy input to output, sign-alternating and zero-padding along the way */
1004
1005 cpt = cxar ;
1006 fpt = ar ;
1007 fj = 1.0 ;
1008 for( jj=0 ; jj < ny ; jj++ ){
1009 fi = fj ; fj = -fj ;
1010 for(ii=0; ii<nx ; ii++){cpt->r=*fpt*fi; cpt->i=0.0; cpt++;fpt++;fi=-fi;}
1011 for( ; ii<nxup; ii++){cpt->r=cpt->i=0.0; cpt++;}
1012 }
1013 for( ; jj < nyup ; jj++ ){cpt->r=cpt->i=0.0; cpt++;}
1014
1015 /* row FFTs */
1016
1017 for( jj=0 ; jj < ny ; jj++ )
1018 csfft_cox( -1 , nxup , cxar+jj*nxup ) ;
1019
1020 /* column FFTs */
1021
1022 cpt = (complex *) malloc(sizeof(complex)*nyup) ;
1023
1024 for( ii=0 ; ii < nxup ; ii++ ){
1025 for( jj=0 ; jj < nyup ; jj++ ) cpt[jj] = cxar[ii+jj*nxup] ;
1026 csfft_cox( -1 , nyup , cpt ) ;
1027 for( jj=0 ; jj < nyup ; jj++ ) cxar[ii+jj*nxup] = cpt[jj] ;
1028 }
1029
1030 /* copy to output */
1031
1032 for( jj=0 ; jj < ny ; jj++ )
1033 for( ii=0 ; ii < nx ; ii++ )
1034 ar[ii+jj*nx] = CARG(cxar[ii+jj*nxup]) ;
1035
1036 free(cxar) ; free(cpt) ; return ;
1037 }
1038
1039 /*----- 28 Oct 2014: sharpness estimate for an image -------------------------*/
1040
sharpness2D_func(int nx,int ny,double dx,double dy,float * ar)1041 void sharpness2D_func( int nx , int ny , double dx, double dy, float *ar )
1042 {
1043 MRI_IMAGE *inim , *outim ;
1044
1045 if( ar == NULL || nx < 5 || ny < 5 ) return ;
1046
1047 inim = mri_new_vol_empty( nx,ny,1 , MRI_float ) ;
1048 mri_set_data_pointer(inim,ar) ;
1049 outim = mri_sharpness(inim) ;
1050 if( outim != NULL ){
1051 memcpy( ar , MRI_FLOAT_PTR(outim) , sizeof(float)*nx*ny ) ;
1052 mri_free(outim) ; mri_clear_and_free(inim) ;
1053 }
1054 return ;
1055 }
1056