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