1 #include "rhdd.c"
2
3 /*----------------------------------------------------------------------------*/
4
5 typedef struct {
6 float xc , yc , zc ; /* centroid of RHDD */
7 int num_ijk , num_in_mask ;
8 int *ijk ;
9 float *vijk ;
10 } RHDD_plist ;
11
12 #undef DESTROY_RHDD_plist
13 #define DESTROY_RHDD_plist(pl) do{ if( (pl)->ijk != NULL ) free((pl)->ijk ) ; \
14 if( (pl)->vijk != NULL ) free((pl)->vijk) ; \
15 free(pl) ; \
16 } while(0)
17
18 /*----------------------------------------------------------------------------*/
19
20 typedef struct {
21 float asiz ;
22 int num_rhdd ;
23 RHDD_plist **rhdd ;
24 } RHDD_array ;
25
26 #undef DESTROY_RHDD_array
27 #define DESTROY_RHDD_array(pa) do{ if( (pa)->rhdd != NULL ) free((pa)->rhdd) ; \
28 free(pa) ; \
29 } while(0)
30
31 /*----------------------------------------------------------------------------*/
32
33 typedef struct {
34 mat44 amat ; /* 12 params */
35 mat33 dmat ; /* 9 params */
36 int nrh ;
37 float *xrh , *yrh , *zrh ; /* 3*nrh params */
38 RHDD_array *rar ;
39 } RHDD_warp ;
40
41 /*----------------------------------------------------------------------------*/
42
RHDD_plist_create(int nx,int ny,int nz,float dx,float dy,float dz,float asiz,float xcen,float ycen,float zcen,byte * mask)43 RHDD_plist * RHDD_plist_create( int nx , int ny , int nz ,
44 float dx , float dy , float dz ,
45 float asiz, float xcen, float ycen, float zcen,
46 byte *mask )
47 {
48 int ii,jj,kk,ijk,num=0,nall=0,nin=0,ibot,itop,jbot,jtop,kbot,ktop,nxy=nx*ny ;
49 RHDD_plist *rpl ;
50 float xbot,xtop , ybot,ytop , zbot,ztop , xx,yy,zz , val ;
51
52 ENTRY("RHDD_plist_create") ;
53
54 #if 0
55 if( nx < 3 || ny < 3 || nz < 3 ) RETURN(NULL) ;
56 if( dx <= 0.0f ) dx = 1.0f ;
57 if( dy <= 0.0f ) dy = 1.0f ;
58 if( dz <= 0.0f ) dz = 1.0f ;
59 if( asiz <= dx || asiz <= dy || asiz <= dz ) RETURN(NULL) ;
60 #endif
61
62 rpl = (RHDD_plist *)calloc(1,sizeof(RHDD_plist)) ;
63 rpl->xc = xcen ; rpl->yc = ycen ; rpl->zc = zcen ;
64
65 xbot = xcen - 2.0f*asiz; ibot = (int)(xbot/dx) ; if( ibot < 0 ) ibot = 0 ;
66 xtop = xcen + 2.0f*asiz; itop = (int)(xtop/dx)+2; if( itop > nx ) itop = nx;
67
68 ybot = ycen - 2.0f*asiz; jbot = (int)(ybot/dy) ; if( jbot < 0 ) jbot = 0 ;
69 ytop = ycen + 2.0f*asiz; jtop = (int)(ytop/dy)+2; if( jtop > ny ) jtop = ny;
70
71 zbot = zcen - 2.0f*asiz; kbot = (int)(zbot/dz) ; if( kbot < 0 ) kbot = 0 ;
72 ztop = zcen + 2.0f*asiz; ktop = (int)(ztop/dz)+2; if( ktop > nz ) ktop = nz;
73
74 for( kk=kbot ; kk <= ktop ; kk++ ){
75 zz = (kk*dz - zcen) / asiz ;
76 for( jj=jbot ; jj <= jtop ; jj++ ){
77 yy = (jj*dy - ycen) / asiz ;
78 for( ii=ibot ; ii <= itop ; ii++ ){
79 xx = (ii*dx - xcen) / asiz ;
80 val = rhddc2(xx,yy,zz) ;
81 if( fabsf(val) > 0.01f ){
82 if( num == nall ){
83 nall = 2*num + 16 ;
84 rpl->ijk = (int * )realloc(rpl->ijk ,sizeof(int )*nall) ;
85 rpl->vijk = (float *)realloc(rpl->vijk,sizeof(float)*nall) ;
86 }
87 rpl->ijk [num] = ijk = ii + jj*nx + kk*nxy ;
88 rpl->vijk[num] = val ; num++ ;
89 if( mask == NULL || mask[ijk] != 0 ) nin++ ;
90 }
91 }
92 }
93 }
94
95 rpl->num_ijk = num ; rpl->num_in_mask = nin ;
96
97 if( num == 0 ){ free(rpl) ; rpl = NULL ; }
98 else if( num < nall ){
99 rpl->ijk = (int * )realloc(rpl->ijk ,sizeof(int )*num) ;
100 rpl->vijk = (float *)realloc(rpl->vijk,sizeof(float)*num) ;
101 }
102
103 RETURN(rpl) ;
104 }
105
106 /*----------------------------------------------------------------------------*/
107
RHDD_array_create(int nx,int ny,int nz,float dx,float dy,float dz,float asiz,byte * mask)108 RHDD_array * RHDD_array_create( int nx , int ny , int nz ,
109 float dx , float dy , float dz ,
110 float asiz , byte *mask )
111 {
112 int pb,pt , qb,qt , rb,rt , pp,qq,rr , np,nq,nr,npq,npqr , nmax,ngood ;
113 float ah=0.5f*asiz , xt,yt,zt ;
114 THD_mat33 latmat , invlatmat ; THD_fvec3 pqr , xyz ;
115 RHDD_plist *rpl ; RHDD_array *rar ;
116
117 ENTRY("RHDD_array_create") ;
118
119 if( nx < 9 || ny < 9 || nz < 9 ) RETURN(NULL) ;
120 if( dx <= 0.0f ) dx = 1.0f ;
121 if( dy <= 0.0f ) dy = 1.0f ;
122 if( dz <= 0.0f ) dz = 1.0f ;
123 xt = MAX(dx,dy) ; xt = MAX(xt,dz) ;
124 if( asiz < 3.0f*xt ) RETURN(NULL) ;
125
126 /* find range of (p,q,r) indexes needed to cover volume,
127 by checking out all 7 corners besides (0,0,0) (where p=q=r=0) */
128
129 LOAD_MAT( latmat , -ah,ah,ah , ah,-ah,ah , ah,ah,-ah ) ;
130 invlatmat = MAT_INV(latmat) ;
131
132 xt = (nx-1)*dx ; yt = (ny-1)*dy ; zt = (nz-1)*dz ;
133 pb = pt = qb = qt = rb = rt = 0 ; /* initialize (p,q,r) bot, top values */
134
135 LOAD_FVEC3(xyz , xt,0.0f,0.0f ); pqr = MATVEC( invlatmat , xyz ) ;
136 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
137 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
138 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
139
140 LOAD_FVEC3(xyz , xt,yt,0.0f ) ; pqr = MATVEC( invlatmat , xyz ) ;
141 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
142 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
143 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
144
145 LOAD_FVEC3(xyz , xt,0.0f,zt ) ; pqr = MATVEC( invlatmat , xyz ) ;
146 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
147 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
148 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
149
150 LOAD_FVEC3(xyz , xt,yt,zt ) ; pqr = MATVEC( invlatmat , xyz ) ;
151 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
152 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
153 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
154
155 LOAD_FVEC3(xyz , 0.0f,yt,0.0f ); pqr = MATVEC( invlatmat , xyz ) ;
156 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
157 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
158 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
159
160 LOAD_FVEC3(xyz , 0.0f,0.0f,zt ); pqr = MATVEC( invlatmat , xyz ) ;
161 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
162 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
163 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
164
165 LOAD_FVEC3(xyz , 0.0f,yt,zt ) ; pqr = MATVEC( invlatmat , xyz ) ;
166 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
167 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
168 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
169
170 pb-- ; qb-- ; rb-- ; /* push outwards, for luck */
171 pt++ ; qt++ ; rt++ ;
172
173 /* Lattice index range is (p,q,r) = (pb..pt,qb..qt,rb..rt) inclusive */
174
175 np = pt-pb+1 ; /* number of p values to consider */
176 nq = qt-qb+1 ; npq = np*nq ;
177 nr = rt-rb+1 ; npqr = npq*nr ;
178
179 rar = (RHDD_array *)calloc(1,sizeof(RHDD_array)) ;
180 rar->asiz = asiz ;
181 rar->num_rhdd = npqr ;
182 rar->rhdd = (RHDD_plist **)calloc(npqr,sizeof(RHDD_plist *)) ;
183
184 for( nmax=0,rr=rb ; rr <= rt ; rr++ ){
185 for( qq=qb ; qq <= qt ; qq++ ){
186 for( pp=pb ; pp <= pt ; pp++ ){
187 LOAD_FVEC3(pqr,pp,qq,rr) ; xyz = MATVEC(latmat,pqr) ;
188 rar->rhdd[(pp-pb)+(qq-qb)*np+(rr-rb)*npq] = rpl =
189 RHDD_plist_create( nx,ny,nz , dx,dy,dz , asiz ,
190 xyz.xyz[0] , xyz.xyz[1] , xyz.xyz[2] , mask ) ;
191 if( rpl != NULL && rpl->num_in_mask > nmax ) nmax = rpl->num_in_mask ;
192 }}}
193
194 if( nmax <= 9 ){
195 for( pp=0 ; pp < npqr ; pp++ ){
196 if( rar->rhdd[pp] != NULL ) DESTROY_RHDD_plist(rar->rhdd[pp]) ;
197 }
198 DESTROY_RHDD_array(rar) ; RETURN(NULL) ;
199 }
200
201 nmax = (int)(0.3456789f*nmax) ;
202 for( ngood=pp=0 ; pp < npqr ; pp++ ){
203 if( rar->rhdd[pp] != NULL ){
204 if( rar->rhdd[pp]->num_in_mask < nmax ){
205 DESTROY_RHDD_plist(rar->rhdd[pp]) ; rar->rhdd[pp] = NULL ;
206 } else {
207 ngood++ ;
208 }
209 }
210 }
211
212 if( ngood < npqr ){
213 RHDD_array *qar ;
214 qar = (RHDD_array *)calloc(1,sizeof(RHDD_array)) ;
215 qar->asiz = asiz ;
216 qar->num_rhdd = ngood ;
217 qar->rhdd = (RHDD_plist **)calloc(ngood,sizeof(RHDD_plist *)) ;
218 for( pp=qq=0 ; pp < npqr ; pp++ ){
219 if( rar->rhdd[pp] != NULL ) qar->rhdd[qq++] = rar->rhdd[pp] ;
220 }
221 DESTROY_RHDD_array(rar) ; rar = qar ;
222 }
223
224 RETURN(rar) ;
225 }
226