1 #include "mrilib.h"
2
3 /* prototypes */
4
5 int THD_conformant_dataxes( THD_dataxes *ax, THD_dataxes *bx ) ;
6 THD_dataxes * THD_superset_dataxes( THD_dataxes *ax, THD_dataxes *bx ) ;
7 static int_pair zpadax_pm( int nx_super, float xorg_super,
8 int nx_input, float xorg_input, float dx ) ;
9 int THD_conformist( int ndset, THD_3dim_dataset **dset, int flags, int *ijkpad ) ;
10
11 #define CONFORM_REWRITE 1
12 #define CONFORM_NOREFIT 2
13
14 /*----------------------------------------------------------------------------*/
15 /* For finding and extending a collection of datasets to have the same grid. */
16
THD_conformist(int ndset,THD_3dim_dataset ** dset,int flags,int * ijkpad)17 int THD_conformist( int ndset, THD_3dim_dataset **dset, int flags , int *ijkpad )
18 {
19 int iset , xpad_m,xpad_p , ypad_m,ypad_p , zpad_m,zpad_p , nwrit ;
20 THD_dataxes *cx , *dx ;
21 THD_3dim_dataset *qset ;
22 int_pair pm ;
23 int do_norefit = (flags & CONFORM_NOREFIT) ;
24 int do_rewrite = (flags & CONFORM_REWRITE) && !do_norefit ;
25 int do_ijkpad = (ijkpad != NULL) ;
26
27 ENTRY("THD_conformist") ;
28
29 /* check for good inputs */
30
31 if( ndset <= 0 || dset == NULL ) RETURN(-1) ;
32 for( iset=0 ; iset < ndset ; iset++ )
33 if( !ISVALID_DSET(dset[iset]) ) RETURN(-2) ;
34
35 /* check if all inputs are on the same grid (i.e., no work to do) */
36
37 for( iset=1 ; iset < ndset ; iset++ ){
38 if( ! EQUIV_DATAXES(dset[0]->daxes,dset[iset]->daxes) ) break ;
39 }
40 if( iset == ndset ){
41 if( do_ijkpad ){
42 for( iset=0 ; iset < 6*ndset ; iset++ ) ijkpad[iset] = 0 ;
43 }
44 RETURN(0) ;
45 }
46
47 /* construct the dataxes that encloses all the input datasets */
48
49 cx = THD_superset_dataxes( dset[0]->daxes , dset[1]->daxes ) ;
50 if( cx == NULL ){
51 ERROR_message("3dConformist: '%s' and '%s' aren't compatible",
52 DSET_HEADNAME(dset[0]) , DSET_HEADNAME(dset[1]) ) ;
53 RETURN(-3) ;
54 }
55
56 for( iset=2 ; iset < ndset ; iset++ ){
57 dx = THD_superset_dataxes( cx , dset[iset]->daxes ) ;
58 if( dx == NULL ){
59 ERROR_message("3dConformist: '%s' is not compatible with others",
60 DSET_HEADNAME(dset[iset]) ) ;
61 myRwcFree(cx) ; RETURN(-3) ;
62 }
63 myRwcFree(cx) ; cx = dx ; dx = NULL ;
64 }
65
66 /* now, re-create and re-write all datasets */
67
68 if( do_rewrite )
69 fprintf(stderr," + thd_conformist re-write loop: ") ;
70
71 for( nwrit=iset=0 ; iset < ndset ; iset++ ){
72 if( EQUIV_DATAXES(cx,dset[iset]->daxes) ){ /* already OK */
73 if( do_rewrite ) fprintf(stderr,"-") ;
74 if( do_ijkpad )
75 ijkpad[6*iset+0] = ijkpad[6*iset+1] = ijkpad[6*iset+2]
76 = ijkpad[6*iset+3] = ijkpad[6*iset+4] = ijkpad[6*iset+5] = 0 ;
77 continue ;
78 }
79 pm = zpadax_pm( cx->nxx , cx->xxorg ,
80 dset[iset]->daxes->nxx , dset[iset]->daxes->xxorg ,
81 cx->xxdel ) ;
82 xpad_m = pm.i ; xpad_p = pm.j ;
83 pm = zpadax_pm( cx->nyy , cx->yyorg ,
84 dset[iset]->daxes->nyy , dset[iset]->daxes->yyorg ,
85 cx->yydel ) ;
86 ypad_m = pm.i ; ypad_p = pm.j ;
87 pm = zpadax_pm( cx->nzz , cx->zzorg ,
88 dset[iset]->daxes->nzz , dset[iset]->daxes->zzorg ,
89 cx->zzdel ) ;
90 zpad_m = pm.i ; zpad_p = pm.j ;
91 if( do_ijkpad ){
92 ijkpad[6*iset+0] = xpad_m; ijkpad[6*iset+1] = xpad_p;
93 ijkpad[6*iset+2] = ypad_m; ijkpad[6*iset+3] = ypad_p;
94 ijkpad[6*iset+4] = zpad_m; ijkpad[6*iset+5] = zpad_p;
95 }
96 if( do_norefit ) continue ; /* just wanted ijkpad, I guess */
97 qset = THD_zeropad( dset[iset] ,
98 xpad_m,xpad_p , ypad_m,ypad_p , zpad_m,zpad_p ,
99 "BertieWooster" , ZPAD_PURGE | ZPAD_IJK ) ;
100 if( qset == NULL ){ /* this should never happen */
101 if( do_rewrite ) fprintf(stderr,"\n") ;
102 ERROR_message("thd_conformist: skipping dataset %s",
103 DSET_HEADNAME(dset[iset])) ;
104 continue ;
105 }
106 qset->idcode = dset[iset]->idcode ;
107 EDIT_dset_items( qset , ADN_prefix , DSET_PREFIX(dset[iset]) , ADN_none ) ;
108 if( do_rewrite ){
109 THD_delete_3dim_dataset( dset[iset] , True ) ;
110 DSET_overwrite(qset) ; DSET_unload(qset) ;
111 fprintf(stderr,"+") ;
112 } else {
113 THD_delete_3dim_dataset( dset[iset] , False ) ;
114 }
115 dset[iset] = qset ; nwrit++ ;
116 }
117 if( do_rewrite ) fprintf(stderr,"\n") ;
118
119 RETURN(nwrit) ;
120 }
121
122 /*----------------------------------------------------------------------------*/
123
THD_conformant_dataxes(THD_dataxes * ax,THD_dataxes * bx)124 int THD_conformant_dataxes( THD_dataxes *ax , THD_dataxes *bx )
125 {
126 double xo,yo,zo ;
127
128 if( ax->xxorient != bx->xxorient ||
129 ax->yyorient != bx->yyorient ||
130 ax->zzorient != bx->zzorient ) return 0 ;
131
132 if( fabsf(ax->xxdel-bx->xxdel) > 0.001f ) return 0 ;
133 if( fabsf(ax->yydel-bx->yydel) > 0.001f ) return 0 ;
134 if( fabsf(ax->zzdel-bx->zzdel) > 0.001f ) return 0 ;
135
136 xo = ((double)ax->xxorg - (double)bx->xxorg) / (double)ax->xxdel ;
137 yo = ((double)ax->yyorg - (double)bx->yyorg) / (double)ax->yydel ;
138 zo = ((double)ax->zzorg - (double)bx->zzorg) / (double)ax->zzdel ;
139
140 if( fabs(xo-rint(xo)) > 0.01 ||
141 fabs(yo-rint(yo)) > 0.01 ||
142 fabs(zo-rint(zo)) > 0.01 ) return 0 ;
143
144 return 1 ;
145 }
146
147 /*----------------------------------------------------------------------------*/
148
zpadax_pm(int nx_super,float xorg_super,int nx_input,float xorg_input,float dx)149 static int_pair zpadax_pm( int nx_super , float xorg_super ,
150 int nx_input , float xorg_input , float dx )
151 {
152 int_pair pm ; float ts , ti ;
153
154 ts = xorg_super + nx_super * dx ;
155 ti = xorg_input + nx_input * dx ;
156 pm.i = (int)rintf((xorg_input-xorg_super)/dx) ;
157 pm.j = (int)rintf((ts-ti)/dx) ;
158
159 return pm ;
160 }
161
162 /*----------------------------------------------------------------------------*/
163
THD_superset_dataxes(THD_dataxes * ax,THD_dataxes * bx)164 THD_dataxes * THD_superset_dataxes( THD_dataxes *ax , THD_dataxes *bx )
165 {
166 THD_dataxes *cx ;
167 float dx,dy,dz , axo,ayo,azo , bxo,byo,bzo , ae,be ;
168 int nxa,nya,nza , nxb,nyb,nzb , ndif ;
169 float cxo,cyo,czo ;
170 int nxc,nyc,nzc ;
171
172 if( !THD_conformant_dataxes(ax,bx) ) return NULL ;
173
174 /* create new dataxes as copy of first one */
175
176 cx = myRwcNew(THD_dataxes) ; *cx = *ax ;
177 cx->parent = NULL ;
178 /* if( EQUIV_DATAXES(ax,bx) ) return cx ; */
179
180 /* load some variables from the input structs */
181
182 dx = ax->xxdel ; dy = ax->yydel ; dz = ax->zzdel ; /* same for ax & bx */
183 axo = ax->xxorg ; ayo = ax->yyorg ; azo = ax->zzorg ; /* ax origins */
184 bxo = bx->xxorg ; byo = bx->yyorg ; bzo = bx->zzorg ; /* bx origins */
185 nxa = ax->nxx ; nya = ax->nyy ; nza = ax->nzz ; /* ax grid lengths */
186 nxb = bx->nxx ; nyb = bx->nyy ; nzb = bx->nzz ; /* bx grid lengths */
187
188 /* extend origins to the outermost */
189
190 ndif = (int)rintf((axo-bxo)/dx) ; cxo = (ndif <= 0) ? axo : bxo ;
191 ndif = (int)rintf((ayo-byo)/dy) ; cyo = (ndif <= 0) ? ayo : byo ;
192 ndif = (int)rintf((azo-bzo)/dz) ; czo = (ndif <= 0) ? azo : bzo ;
193
194 cx->xxorg = cxo ; cx->yyorg = cyo ; cx->zzorg = czo ;
195
196 /* extend grid lengths to the outermost */
197
198 ae = axo + nxa*dx ; be = bxo + nxb*dx ;
199 ndif = (int)rintf((ae-be)/dx) ; if( ndif < 0 ) ae = be ;
200 nxc = (int)rintf((ae-cxo)/dx) ;
201
202 ae = ayo + nya*dy ; be = byo + nyb*dy ;
203 ndif = (int)rintf((ae-be)/dy) ; if( ndif < 0 ) ae = be ;
204 nyc = (int)rintf((ae-cyo)/dy) ;
205
206 ae = azo + nza*dz ; be = bzo + nzb*dz ;
207 ndif = (int)rintf((ae-be)/dz) ; if( ndif < 0 ) ae = be ;
208 nzc = (int)rintf((ae-czo)/dz) ;
209
210 cx->nxx = nxc ; cx->nyy = nyc ; cx->nzz = nzc ;
211
212 #if 0
213 fprintf(stderr,"\n") ;
214 INFO_message("ax: nxyz=%d %d %d org=%g %g %g del=%g %g %g",ax->nxx,ax->nyy,ax->nzz,ax->xxorg,ax->yyorg,ax->zzorg,ax->xxdel,ax->yydel,ax->zzdel) ;
215 INFO_message("bx: nxyz=%d %d %d org=%g %g %g del=%g %g %g",bx->nxx,bx->nyy,bx->nzz,bx->xxorg,bx->yyorg,bx->zzorg,bx->xxdel,bx->yydel,bx->zzdel) ;
216 INFO_message("cx: nxyz=%d %d %d org=%g %g %g del=%g %g %g",cx->nxx,cx->nyy,cx->nzz,cx->xxorg,cx->yyorg,cx->zzorg,cx->xxdel,cx->yydel,cx->zzdel) ;
217 #endif
218
219 /* fix the matrices etc */
220
221 LOAD_ZERO_MAT(cx->to_dicomm) ;
222 THD_daxes_to_mat44(cx) ;
223 THD_set_daxes_bbox(cx) ;
224 cx->ijk_to_dicom_real = cx->ijk_to_dicom ;
225
226 return cx ;
227 }
228