1 #include "mrilib.h"
2
3 /*-------------------------------------------------------------------*/
4 /* 3D sharpening filter, for program 3dSharpen3D.
5 * Only operates on positive values, and only on those
6 which have all 26 NN1, NN2, and NN3 points being positive
7 as well (to avoid weird stuff at edges).
8 * Uses a simple local difference approach, with NN1 & NN2.
9 * im must be in floating point, and is edited in place
10 * phi is between 0 and 1 (exclusive).
11 ***** RWCox -- 13 Feb 2017 *****
12 *//*-----------------------------------------------------------------*/
13
14 #define RFAC 0.5f
15
mri_sharpen3D_pos(MRI_IMAGE * im,float phi)16 void mri_sharpen3D_pos( MRI_IMAGE *im , float phi )
17 {
18 float *imar , *tar ;
19 float nphi , omphi , pval , psum , ppp , pbot,ptop , tbot,ttop ;
20 int nx,ny,nz, ii,jj,kk, joff,koff,jkoff,ijk, nx1,ny1,nz1,nxyz , ndon;
21
22 if( im == NULL || im->kind != MRI_float ) return ;
23 if( phi <= 0.0f || phi >= 1.0f ) return ;
24
25 nx = im->nx ; ny = im->ny ; nz = im->nz ;
26 if( nx < 3 || ny < 3 || nz < 3 ) return ;
27 nx1 = nx-1 ; ny1 = ny-1 ; nz1 = nz-1 ;
28 joff = nx ; koff = nx*ny ; nxyz = koff*nz ;
29
30 imar = MRI_FLOAT_PTR(im) ;
31 if( imar == NULL ) return ;
32
33 nphi = phi / 19.0f ;
34 omphi = 1.0/(1.0-phi) ;
35
36 tar = (float *)malloc(sizeof(float)*nxyz) ;
37 memcpy( tar , imar , sizeof(float)*nxyz) ;
38
39 tbot = 1.e+37 ; ttop = 0.0f ;
40 for( ijk=0 ; ijk < nxyz ; ijk++ ){
41 if( tar[ijk] > 0.0f ){
42 tbot = MIN(tbot,tar[ijk]) ;
43 ttop = MAX(ttop,tar[ijk]) ;
44 }
45 }
46 if( tbot >= ttop ){ free(tar) ; return ; }
47
48 #define DOFF(qq) \
49 { ppp = tar[ijk+(qq)] ; \
50 if( ppp <= 0.0f ) continue ; \
51 else { psum += ppp; pbot=MIN(pbot,ppp); ptop=MAX(ptop,ppp); } }
52
53 #define ZOFF(qq) \
54 { ppp = tar[ijk+(qq)] ; if( ppp <= 0.0f ) continue ; }
55
56 ndon = 0 ;
57 for( kk=1 ; kk < nz1 ; kk++ ){
58 for( jj=1 ; jj < ny1 ; jj++ ){
59 jkoff = jj*joff + kk*koff ;
60 for( ii=1 ; ii < nx1 ; ii++ ){
61 ijk = ii+jkoff ;
62 pbot = ptop = psum = pval = tar[ijk] ; if( pval <= 0.0f ) continue;
63 DOFF(-1) ; DOFF( 1) ; /* operate on the 18 */
64 DOFF(-joff) ; DOFF( joff) ; /* NN1 and NN2 points */
65 DOFF(-koff) ; DOFF( koff) ;
66 DOFF(-joff-koff) ; DOFF(-joff+koff) ;
67 DOFF( joff-koff) ; DOFF( joff-koff) ;
68 DOFF(-joff-1) ; DOFF(-joff+1) ;
69 DOFF( joff-1) ; DOFF( joff+1) ;
70 DOFF(-koff-1) ; DOFF(-koff+1) ;
71 DOFF( koff-1) ; DOFF( koff+1) ;
72 ZOFF( joff+koff+1); ZOFF( joff+koff-1); /* skip if any of the NN3 */
73 ZOFF( joff-koff+1); ZOFF( joff-koff-1); /* points are non-positive */
74 ZOFF(-joff+koff+1); ZOFF(-joff+koff-1);
75 ZOFF(-joff-koff+1); ZOFF(-joff-koff-1);
76
77 pbot -= RFAC*(pval-pbot) ; if( pbot < tbot ) pbot = tbot ;
78 ptop += RFAC*(ptop-pval) ; if( ptop > ttop ) ptop = ttop ;
79 pval = (pval - nphi*psum)*omphi ;
80 if( pval < pbot ) pval = pbot ; else if( pval > ptop ) pval = ptop ;
81 imar[ijk] = pval ; ndon++ ;
82 }}}
83
84 free(tar) ; return ;
85 }
86