1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8 #include "thd.h"
9
10 /*-----------------------------------------------------------------------
11 Create FD_bricks for viewing purposes.
12
13 FD_bricks are a relic of the very earliest code that evolved into
14 AFNI. By the time I invented THD_3dim_datasets, I had enough code
15 in place to keep this data structure alive. Basically, it exists
16 only to provide a structure that enables fast copying of data out
17 of a 3D array into 2D arrays.
18 -------------------------------------------------------------------------*/
19
THD_setup_bricks(THD_3dim_dataset * dset)20 FD_brick ** THD_setup_bricks( THD_3dim_dataset * dset )
21 {
22 int r2l=0 , a2p=0 , i2s=0 ;
23 THD_dataxes *daxes ;
24 FD_brick **br ;
25
26 ENTRY("THD_setup_bricks") ;
27
28 if( ! ISVALID_3DIM_DATASET(dset) ) RETURN(NULL) ;
29
30 daxes = CURRENT_DAXES(dset) ;
31 if( ! ISVALID_DATAXES(daxes) ) RETURN(NULL) ;
32
33 /*----- create FD_bricks for viewing purposes -----*/
34
35 switch( daxes->xxorient ){
36 case ORI_R2L_TYPE: r2l = 1 ; break ;
37 case ORI_L2R_TYPE: r2l = -1 ; break ;
38 case ORI_P2A_TYPE: a2p = -1 ; break ;
39 case ORI_A2P_TYPE: a2p = 1 ; break ;
40 case ORI_I2S_TYPE: i2s = 1 ; break ;
41 case ORI_S2I_TYPE: i2s = -1 ; break ;
42 }
43
44 switch( daxes->yyorient ){
45 case ORI_R2L_TYPE: r2l = 2 ; break ;
46 case ORI_L2R_TYPE: r2l = -2 ; break ;
47 case ORI_P2A_TYPE: a2p = -2 ; break ;
48 case ORI_A2P_TYPE: a2p = 2 ; break ;
49 case ORI_I2S_TYPE: i2s = 2 ; break ;
50 case ORI_S2I_TYPE: i2s = -2 ; break ;
51 }
52
53 switch( daxes->zzorient ){
54 case ORI_R2L_TYPE: r2l = 3 ; break ;
55 case ORI_L2R_TYPE: r2l = -3 ; break ;
56 case ORI_P2A_TYPE: a2p = -3 ; break ;
57 case ORI_A2P_TYPE: a2p = 3 ; break ;
58 case ORI_I2S_TYPE: i2s = 3 ; break ;
59 case ORI_S2I_TYPE: i2s = -3 ; break ;
60 }
61
62 if( r2l==0 || a2p==0 || i2s==0 ){
63 char buf[256] ;
64 sprintf(buf,"Illegal orientation codes: %d %d %d",
65 daxes->xxorient,daxes->yyorient,daxes->zzorient ) ;
66 THD_FATAL_ERROR(buf) ;
67 }
68
69 /* now we can set up views: axial, sagittal, coronal;
70 the top option is the way I think it ought to be, so that
71 in the axial and coronal views, left is left and right is right;
72 the bottom options are the radiologists conventions, so sue me! */
73
74 br = (FD_brick **) RwcMalloc( sizeof(FD_brick *) * 3 ) ;
75
76 #undef LEFT_IS_LEFT
77 #ifdef LEFT_IS_LEFT
78 br[0] = THD_3dim_dataset_to_brick(dset,-r2l, a2p,-i2s); /* axi */
79 br[1] = THD_3dim_dataset_to_brick(dset, a2p,-i2s,-r2l); /* sag */
80 br[2] = THD_3dim_dataset_to_brick(dset,-r2l,-i2s,-a2p); /* cor */
81 #else
82 br[0] = THD_3dim_dataset_to_brick(dset, r2l, a2p, i2s); /* axi */
83 br[1] = THD_3dim_dataset_to_brick(dset, a2p,-i2s,-r2l); /* sag */
84 br[2] = THD_3dim_dataset_to_brick(dset, r2l,-i2s, a2p); /* cor */
85 #endif
86
87 MCW_strncpy( br[0]->namecode , "Axial" , 32 ) ;
88 MCW_strncpy( br[1]->namecode , "Sagittal" , 32 ) ;
89 MCW_strncpy( br[2]->namecode , "Coronal" , 32 ) ;
90
91 RETURN(br) ;
92 }
93
94 /*----------------------------------------------------------------------
95 07 Dec 2001 - orient an FD brick in any legal way
96 (e.g., orients = "RAI" for standard axial)
97 ------------------------------------------------------------------------*/
98
THD_oriented_brick(THD_3dim_dataset * dset,char * orients)99 FD_brick * THD_oriented_brick( THD_3dim_dataset *dset , char *orients )
100 {
101 int r2l=0 , a2p=0 , i2s=0 , xx,yy,zz , pp=0,qq=0,rr=0 ;
102 THD_dataxes *daxes ;
103 FD_brick *br ;
104
105 ENTRY("THD_oriented_brick") ;
106
107 if( !ISVALID_DSET(dset) ||
108 orients == NULL ||
109 strlen(orients) < 3 ) RETURN(NULL) ;
110
111 daxes = CURRENT_DAXES(dset) ;
112 if( !ISVALID_DATAXES(daxes) ) RETURN(NULL) ;
113
114 xx = ORCODE( orients[0] ) ;
115 yy = ORCODE( orients[1] ) ;
116 zz = ORCODE( orients[2] ) ;
117 if( !OR3OK(xx,yy,zz) ) RETURN(NULL) ;
118
119 /*----- create FD_bricks for viewing purposes -----*/
120
121 switch( daxes->xxorient ){
122 case ORI_R2L_TYPE: r2l = 1 ; break ;
123 case ORI_L2R_TYPE: r2l = -1 ; break ;
124 case ORI_P2A_TYPE: a2p = -1 ; break ;
125 case ORI_A2P_TYPE: a2p = 1 ; break ;
126 case ORI_I2S_TYPE: i2s = 1 ; break ;
127 case ORI_S2I_TYPE: i2s = -1 ; break ;
128 }
129
130 switch( daxes->yyorient ){
131 case ORI_R2L_TYPE: r2l = 2 ; break ;
132 case ORI_L2R_TYPE: r2l = -2 ; break ;
133 case ORI_P2A_TYPE: a2p = -2 ; break ;
134 case ORI_A2P_TYPE: a2p = 2 ; break ;
135 case ORI_I2S_TYPE: i2s = 2 ; break ;
136 case ORI_S2I_TYPE: i2s = -2 ; break ;
137 }
138
139 switch( daxes->zzorient ){
140 case ORI_R2L_TYPE: r2l = 3 ; break ;
141 case ORI_L2R_TYPE: r2l = -3 ; break ;
142 case ORI_P2A_TYPE: a2p = -3 ; break ;
143 case ORI_A2P_TYPE: a2p = 3 ; break ;
144 case ORI_I2S_TYPE: i2s = 3 ; break ;
145 case ORI_S2I_TYPE: i2s = -3 ; break ;
146 }
147
148 if( r2l==0 || a2p==0 || i2s==0 ) RETURN(NULL) ;
149
150 switch( xx ){
151 case ORI_R2L_TYPE: pp = r2l ; break ;
152 case ORI_L2R_TYPE: pp = -r2l ; break ;
153 case ORI_P2A_TYPE: pp = -a2p ; break ;
154 case ORI_A2P_TYPE: pp = a2p ; break ;
155 case ORI_I2S_TYPE: pp = i2s ; break ;
156 case ORI_S2I_TYPE: pp = -i2s ; break ;
157 }
158
159 switch( yy ){
160 case ORI_R2L_TYPE: qq = r2l ; break ;
161 case ORI_L2R_TYPE: qq = -r2l ; break ;
162 case ORI_P2A_TYPE: qq = -a2p ; break ;
163 case ORI_A2P_TYPE: qq = a2p ; break ;
164 case ORI_I2S_TYPE: qq = i2s ; break ;
165 case ORI_S2I_TYPE: qq = -i2s ; break ;
166 }
167
168 switch( zz ){
169 case ORI_R2L_TYPE: rr = r2l ; break ;
170 case ORI_L2R_TYPE: rr = -r2l ; break ;
171 case ORI_P2A_TYPE: rr = -a2p ; break ;
172 case ORI_A2P_TYPE: rr = a2p ; break ;
173 case ORI_I2S_TYPE: rr = i2s ; break ;
174 case ORI_S2I_TYPE: rr = -i2s ; break ;
175 }
176
177 if( pp==0 || qq==0 || rr==0 ) RETURN(NULL) ;
178
179 br = THD_3dim_dataset_to_brick(dset,pp,qq,rr) ;
180 RETURN(br) ;
181 }
182
183 /*======================================================================
184 routines adapted from original afni code for display of a 3D dataset
185 ========================================================================*/
186
187 /*----------------------------------------------------------------------
188 prepare a "display brick" for a 3D dataset;
189 ax_n = 1,2,3 for displaying the x,y,z axis along the n-th display
190 direction (n=1 --> across, n=2 --> down, n=3 --> slices )
191 negative value means reverse
192
193 a return of NULL means something bad happened
194 ------------------------------------------------------------------------*/
195
THD_3dim_dataset_to_brick(THD_3dim_dataset * dset,int ax_1,int ax_2,int ax_3)196 FD_brick * THD_3dim_dataset_to_brick( THD_3dim_dataset *dset ,
197 int ax_1, int ax_2, int ax_3 )
198 {
199 FD_brick *br ; /* will be output */
200 THD_dataxes *daxes ; /* connection to actual axes */
201
202 int xyz_dim[4] , xyz_stp[4] , xyz_dir[4] ;
203 float xyz_del[4] , xyz_org[4] ;
204
205 int x_dir,y_dir,z_dir , sx,sy,sz , aax_1,aax_2,aax_3 , nx,ny,nz ;
206
207 /*-- sanity check --*/
208
209 ENTRY("THD_3dim_dataset_to_brick") ;
210
211 if( ! ISVALID_3DIM_DATASET(dset) ) RETURN(NULL) ;
212
213 daxes = CURRENT_DAXES(dset) ;
214
215 aax_1 = abs(ax_1) ; /* which axes, regardless of + or - */
216 aax_2 = abs(ax_2) ;
217 aax_3 = abs(ax_3) ;
218
219 if( aax_1 < 1 || aax_1 > 3 || /* range checks */
220 aax_2 < 1 || aax_2 > 3 ||
221 aax_3 < 1 || aax_3 > 3 ) RETURN(NULL) ;
222
223 xyz_dir[1] = xyz_dir[2] = xyz_dir[3] = 0 ;
224
225 xyz_dir[aax_1] = ax_1 ; /* assign to original directions */
226 xyz_dir[aax_2] = ax_2 ;
227 xyz_dir[aax_3] = ax_3 ;
228
229 x_dir = xyz_dir[1] ; /* if any |ax_n| is duplicated */
230 y_dir = xyz_dir[2] ; /* then one of these will end */
231 z_dir = xyz_dir[3] ; /* up as zero --> bad inputs! */
232
233 if( x_dir == 0 || y_dir == 0 || z_dir == 0 ) RETURN(NULL) ;
234
235 /*-- the inputs are good, so create a brick: --*/
236
237 br = myRwcNew(FD_brick) ; /* new brick */
238 br->dset = dset ; /* dataset */
239 br->resam_code = RESAM_NN_TYPE ; /* crudest type */
240 br->parent = NULL ;
241 br->brother = NULL ;
242 br->deltival = 0 ; /* 23 Feb 2011 */
243
244 br->thr_resam_code = RESAM_NN_TYPE ; /* 09 Dec 1997 */
245
246 /*-- at this point, x_dir is +1 or -1, y_dir is +2 or -2, etc. --*/
247
248 nx = daxes->nxx ; ny = daxes->nyy ; nz = daxes->nzz ;
249
250 sx = (x_dir > 0) ? (0) : (nx-1) ; /* starting voxel indices */
251 sy = (y_dir > 0) ? (0) : (ny-1) ; /* for each original dimension */
252 sz = (z_dir > 0) ? (0) : (nz-1) ;
253
254 br->start = sx + sy*nx + sz*nx*ny ; /* overall starting voxel index */
255
256 /*-- assign original dimensions to arrays,
257 then pick out the permuted dimensions --*/
258
259 xyz_dim[1] = nx ; /* dimensions */
260 xyz_dim[2] = ny ;
261 xyz_dim[3] = nz ;
262
263 LOAD_IVEC3( br->nxyz , nx,ny,nz ) ; /* save stuff in br */
264 LOAD_IVEC3( br->sxyz , sx,sy,sz ) ;
265 LOAD_IVEC3( br->a123 , ax_1,ax_2,ax_3 ) ;
266
267 xyz_stp[1] = 1 ; /* index step sizes */
268 xyz_stp[2] = nx ;
269 xyz_stp[3] = nx * ny ;
270
271 xyz_del[1] = daxes->xxdel ; /* voxel physical step sizes (mm) */
272 xyz_del[2] = daxes->yydel ;
273 xyz_del[3] = daxes->zzdel ;
274
275 xyz_org[1] = daxes->xxorg ; /* voxel origins (mm) */
276 xyz_org[2] = daxes->yyorg ;
277 xyz_org[3] = daxes->zzorg ;
278
279 br->n1 = xyz_dim[aax_1] ; /* permute dimensions, etc. */
280 br->n2 = xyz_dim[aax_2] ;
281 br->n3 = xyz_dim[aax_3] ;
282
283 br->d1 = (ax_1 > 0) ? (xyz_stp[aax_1]) : (-xyz_stp[aax_1]) ;
284 br->d2 = (ax_2 > 0) ? (xyz_stp[aax_2]) : (-xyz_stp[aax_2]) ;
285 br->d3 = (ax_3 > 0) ? (xyz_stp[aax_3]) : (-xyz_stp[aax_3]) ;
286
287 br->e1 = br->n1 * br->d1 ; /* last indices for readout */
288 br->e2 = br->n2 * br->d2 ;
289
290 br->del1 = fabs(xyz_del[aax_1]) ; /* dimensions */
291 br->del2 = fabs(xyz_del[aax_2]) ;
292 br->del3 = fabs(xyz_del[aax_3]) ;
293
294 br->namecode[0] = '\0' ;
295
296 br->tmask = NULL ;
297 br->ntmask = -666 ;
298
299 RETURN(br) ;
300 }
301