1 #include "mrilib.h"
2
3 #undef MANUAL_ORIENT
4
5 /*----------------------------------------------------------------*/
6 /* returns a new string, which can be free()-ed later if you like */
7 /*----------------------------------------------------------------*/
8
EDIT_convert_to_geometry_string(char * str)9 char * EDIT_convert_to_geometry_string( char *str )
10 {
11 char *lstr=NULL , *cpt ;
12 float a11,a12,a13,a14 ;
13 float a21,a22,a23,a24 ;
14 float a31,a32,a33,a34 ;
15 float dx,dy,dz , xorg,yorg,zorg ;
16 int nx,ny,nz , ii, dicomorigin = 0;
17
18 ENTRY("EDIT_convert_to_geometry_string") ;
19
20 if( str == NULL || strlen(str) < 4 ) RETURN(NULL) ;
21
22 lstr = strdup(str) ; /* copy input */
23
24 if( ISVALID_GEOMETRY_STRING(lstr) ) RETURN(lstr) ; /* that was easy */
25
26 /*--- convert lstr to MATRIX() type of string ---*/
27
28 if( strncasecmp(lstr,"TLRC",4) == 0 ){ /* for very old types of data */
29
30 free(lstr) ;
31 lstr = strdup("MATRIX(1,0,0,-80 , 0,1,0,-80 , 0,0,1,-65):161,191,151") ;
32
33 } else if ( strncasecmp(lstr,"MNI2009",7) == 0 ){
34
35 free(lstr) ;
36 lstr = strdup("MATRIX(-1,0,0,96,0,-1,0,132,0,0,1,-78):193,229,193") ;
37
38 } else if( lstr[3] == ':' ){ /* "XYZ:" where XYZ are some coordinate order */
39
40 THD_coorder cord ;
41
42 THD_coorder_fill( lstr , &cord ) ;
43
44 /* convert commas to blanks */
45
46 for( cpt=lstr ; *cpt != '\0' ; cpt++ ) if( *cpt == ',' ) *cpt = ' ' ;
47 nx = ny = nz = -1; dx = dy = dz = -1.0f; xorg = yorg = zorg = 0.0f ;
48
49 if (!strncmp(lstr+4,"D:",2)) { /* scan for 9 values (DICOM) */
50 dicomorigin = 1;
51 ii = sscanf(lstr+6,"%d%f%f%d%f%f%d%f%f",
52 &nx,&xorg,&dx , &ny,&yorg,&dy , &nz,&zorg,&dz ) ;
53 if( ii < 9 ||
54 nx <= 0 || ny <= 0 || nz <= 0 ||
55 dx <= 0 || dy <= 0 || dz <= 0 ){
56 ERROR_message("Negative or 0 voxel counts or voxel sizes");
57 free(lstr); RETURN(NULL);
58 }
59
60 } else { /* scan for 9 values (not DICOM) */
61
62 dicomorigin = 0;
63 ii = sscanf(lstr+4,"%d%f%f%d%f%f%d%f%f",
64 &nx,&xorg,&dx , &ny,&yorg,&dy , &nz,&zorg,&dz ) ;
65 if( ii < 9 ||
66 nx <= 0 || ny <= 0 || nz <= 0 ||
67 dx <= 0 || dy <= 0 || dz <= 0 ){ free(lstr); RETURN(NULL); }
68 }
69
70 a11 = dx ; a21 = 0.0f ; a31 = 0.0f ;
71 THD_coorder_to_dicom( &cord , &a11,&a21,&a31 ) ;
72
73 a12 = 0.0f ; a22 = dy ; a32 = 0.0f ;
74 THD_coorder_to_dicom( &cord , &a12,&a22,&a32 ) ;
75
76 a13 = 0.0f ; a23 = 0.0f ; a33 = dz ;
77 THD_coorder_to_dicom( &cord , &a13,&a23,&a33 ) ;
78
79 a14 = xorg ; a24 = yorg ; a34 = zorg ;
80 if (!dicomorigin) THD_coorder_to_dicom( &cord , &a14,&a24,&a34 ) ;
81
82 cpt = (char *)malloc(sizeof(char)*666) ;
83 sprintf(cpt,
84 "MATRIX(%.4f,%.4f,%.4f,%.4f , %.4f,%.4f,%.4f,%.4f , %.4f,%.4f,%.4f,%.4f):%d,%d,%d" ,
85 a11,a12,a13,a14 , a21,a22,a23,a24, a31,a32,a33,a34 , nx,ny,nz ) ;
86
87 free(lstr) ; lstr = cpt ;
88
89 } else {
90
91 free(lstr) ; lstr = NULL ;
92
93 }
94
95 RETURN(lstr) ;
96 }
97
98
99 /*-------------------------------------------------------------------------*/
100 /*! Create an empty dataset with geometry given by a string. Examples:
101 - "tlrc"
102 - "MNI2009"
103 - "RAI:nx,xorg,dx,ny,yorg,dy,nz,zorg,dz"
104 - "RAI:D:nx,xorg,dx,ny,yorg,dy,nz,zorg,dz" if xorg, yorg, zorg are DICOM
105 - "MATRIX(a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34):nx,ny,nz"
106 *//*-----------------------------------------------------------------------*/
107
EDIT_geometry_constructor(char * gstr,char * prefix)108 THD_3dim_dataset * EDIT_geometry_constructor( char *gstr , char *prefix )
109 {
110 THD_3dim_dataset *dset = NULL ;
111 THD_ivec3 orixyz , nxyz ;
112 THD_fvec3 dxyz , orgxyz ;
113 mat44 ijk_to_dicom44 ; THD_mat33 R ;
114 int view=VIEW_ORIGINAL_TYPE ;
115 float dx,dy,dz , xorg,yorg,zorg ;
116 int nx,ny,nz , ii, dicomorigin = 0;
117 char *lstr , *cpt ;
118 float a11,a12,a13,a14 ;
119 float a21,a22,a23,a24 ;
120 float a31,a32,a33,a34 ;
121 float orgx,orgy,orgz ;
122
123 ENTRY("EDIT_geometry_constructor") ;
124
125 lstr = EDIT_convert_to_geometry_string(gstr) ;
126
127 if( !ISVALID_GEOMETRY_STRING(lstr) ){
128 if( lstr != NULL ) free(lstr) ;
129 RETURN(NULL) ;
130 }
131
132 if( strncasecmp(gstr,"TLRC",4) == 0 ||
133 strncasecmp(gstr,"MNI" ,3) == 0 ) view = VIEW_TALAIRACH_TYPE ;
134
135 /*--- decode MATRIX geometry string ---*/
136
137 #if 0
138 INFO_message("EDIT_geometry_constructor: string = %s",lstr) ;
139 #endif
140
141 /* convert commas to blanks */
142
143 for( cpt=lstr ; *cpt != '\0' ; cpt++ ) if( *cpt == ',' ) *cpt = ' ' ;
144
145 /* scan after the 'MATRIX:' part to get the geometry values */
146
147 ii = sscanf(lstr+7,"%f%f%f%f%f%f%f%f%f%f%f%f):%d%d%d",
148 &a11,&a12,&a13,&a14 ,
149 &a21,&a22,&a23,&a24 , &a31,&a32,&a33,&a34 , &nx,&ny,&nz ) ;
150 free(lstr) ;
151 if( ii < 15 ) RETURN(NULL) ; /* didn't get 15 values from string :( */
152
153 /*--- create dataset and put stuff into it ---*/
154
155 dset = EDIT_empty_copy(NULL) ;
156
157 /*-- orientation stuff --*/
158
159 LOAD_MAT44( ijk_to_dicom44 ,
160 a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34 ) ;
161
162 #if 0
163 DUMP_MAT44("EDIT_geometry_constructor ijk_to_dicom",ijk_to_dicom44) ;
164 #endif
165
166 dset->daxes->ijk_to_dicom_real = ijk_to_dicom44 ;
167 dset->daxes->ijk_to_dicom = ijk_to_dicom44 ;
168
169 nxyz.ijk[0] = nx ; nxyz.ijk[1] = ny ; nxyz.ijk[2] = nz ;
170 EDIT_dset_items( dset , ADN_nxyz,nxyz , ADN_none ) ;
171
172 #ifdef MANUAL_ORIENT /* This one does not work.
173 Changed #ifndef to #ifdef
174 One below seems to work.
175 ZSS, RKR, Jan 22 10 */
176 THD_daxes_from_mat44( dset->daxes ) ;
177 #else
178 LOAD_MAT( R , a11,a12,a13 , a21,a22,a23 , a31,a32,a33 ) ;
179
180 orixyz = THD_matrix_to_orientation( R ) ; /* compute orientation codes */
181
182 #if 0
183 INFO_message("EDIT_geometry_constructor: orientation codes = %d %d %d",orixyz.ijk[0],orixyz.ijk[1],orixyz.ijk[2]) ;
184 #endif
185
186 orgx = ijk_to_dicom44.m[ORIENT_xyzint[orixyz.ijk[0]]-1][3] ; /* somewhat */
187 orgy = ijk_to_dicom44.m[ORIENT_xyzint[orixyz.ijk[1]]-1][3] ; /* confusing */
188 orgz = ijk_to_dicom44.m[ORIENT_xyzint[orixyz.ijk[2]]-1][3] ; /* I admit! */
189 LOAD_FVEC3( orgxyz , orgx,orgy,orgz ) ;
190
191 dx = sqrtf( a11*a11 + a21*a21 + a31*a31 ) ;
192 dy = sqrtf( a12*a12 + a22*a22 + a32*a32 ) ;
193 dz = sqrtf( a13*a13 + a23*a23 + a33*a33 ) ;
194 if( ORIENT_sign[orixyz.ijk[0]] == '-' ) dx = -dx ;
195 if( ORIENT_sign[orixyz.ijk[1]] == '-' ) dy = -dy ;
196 if( ORIENT_sign[orixyz.ijk[2]] == '-' ) dz = -dz ;
197 LOAD_FVEC3( dxyz , dx,dy,dz ) ;
198 EDIT_dset_items( dset ,
199 ADN_xyzdel , dxyz ,
200 ADN_xyzorg , orgxyz ,
201 ADN_xyzorient , orixyz ,
202 ADN_none ) ;
203 #endif
204
205 dset->idcode.str[0] = 'G' ; /* this is just for fun */
206 dset->idcode.str[1] = 'E' ;
207 dset->idcode.str[2] = 'O' ;
208
209 if( !THD_filename_ok(prefix) ) prefix = "gggeom" ; /* loser user */
210
211 EDIT_dset_items( dset ,
212 ADN_prefix , prefix ,
213 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
214 ADN_view_type , view ,
215 ADN_type , HEAD_FUNC_TYPE ,
216 ADN_func_type , FUNC_BUCK_TYPE ,
217 ADN_none ) ;
218
219 RETURN(dset) ;
220 }
221
222 /*-------------------------------------------------------------------------*/
223 /* create a random dataset [16 Mar 2016] */
224
jRandomDataset(int nx,int ny,int nz,int nt)225 THD_3dim_dataset * jRandomDataset( int nx, int ny, int nz, int nt )
226 {
227 THD_3dim_dataset *dset ;
228 char gstr[128] ;
229 int iv,jj,nvox; float *far , zz=0.0f ;
230
231 if( nx < 2 ) return NULL ;
232 if( ny < 1 ) ny = nx ;
233 if( nz < 1 ) nz = 1 ;
234 if( nt < 1 ) nt = 1 ;
235
236 sprintf(gstr,"RAI:%d,0,1.0,%d,0,1.0,%d,0,1.0",nx,ny,nz) ;
237 dset = EDIT_geometry_constructor(gstr,"jRandomDataset") ;
238
239 EDIT_dset_items( dset ,
240 ADN_nvals , nt ,
241 ADN_none ) ;
242 if( nt > 1 ){
243 EDIT_dset_items( dset ,
244 ADN_ntt , nt ,
245 ADN_ttdel , 1.0f ,
246 ADN_none ) ;
247 }
248
249 nvox = nx*ny*nz ;
250 for( iv=0 ; iv < nt ; iv++ ){
251 EDIT_substitute_brick( dset , iv , MRI_float , NULL ) ;
252 far = DSET_ARRAY( dset , iv ) ;
253 for( jj=0 ; jj < nvox ; jj++ ) far[jj] = 2.0f*(float)(drand48())-1.0f ;
254 if( nvox%32 == 0 && iv < nt-1 ){
255 for( jj=0 ; jj < 17 ; jj++ ) zz += drand48() ; /* avoid artifacts */
256 }
257 }
258
259 return dset ;
260 }
261
262 /*-------------------------------------------------------------------------*/
263 /* create a random 1D file [17 Mar 2016] */
264
jRandom1D(int nx,int ny)265 MRI_IMAGE * jRandom1D( int nx , int ny )
266 {
267 MRI_IMAGE *im ; int ii,jj,kk,nxy ; float *far , zz=0.0f ;
268
269 if( nx < 1 ) return NULL ;
270 if( ny < 1 ) ny = 1 ;
271
272 im = mri_new( nx , ny , MRI_float ) ;
273 far = MRI_FLOAT_PTR(im) ;
274 nxy = nx*ny ;
275 for( kk=jj=0 ; jj < ny ; jj++ ){
276 for( ii=0 ; ii < nx ; ii++,kk++ )
277 far[kk] = 2.0f*(float)(drand48())-1.0f ;
278 if( nx%8 == 0 && jj < ny-1 ){
279 int qq ; for( qq=0 ; qq < 11 ; qq++ ) zz += drand48() ;
280 }
281 }
282
283 return im ;
284 }
285
286 /*-------------------------------------------------------------------------*/
287
EDIT_get_geometry_string(THD_3dim_dataset * dset)288 char * EDIT_get_geometry_string( THD_3dim_dataset *dset )
289 {
290 char *ggg ;
291
292 if( !ISVALID_MAT44(dset->daxes->ijk_to_dicom_real) ){
293 THD_daxes_to_mat44(dset->daxes) ;
294 dset->daxes->ijk_to_dicom_real = dset->daxes->ijk_to_dicom ;
295 }
296
297 #if 0
298 DUMP_MAT44("EDIT_get_geometry_string: using ijk_to_dicom_real",dset->daxes->ijk_to_dicom_real) ;
299 #endif
300
301 ggg = EDIT_imat_to_geometry_string( dset->daxes->ijk_to_dicom_real ,
302 DSET_NX(dset),DSET_NY(dset),DSET_NZ(dset) ) ;
303 return ggg ;
304 }
305
306 /*-------------------------------------------------------------------------*/
307
EDIT_imat_to_geometry_string(mat44 imat,int nx,int ny,int nz)308 char * EDIT_imat_to_geometry_string( mat44 imat , int nx,int ny,int nz )
309 {
310 static char gstr[666] ;
311 float a11,a12,a13,a14 ;
312 float a21,a22,a23,a24 ;
313 float a31,a32,a33,a34 ;
314 char b11[32],b12[32],b13[32],b14[32] ;
315 char b21[32],b22[32],b23[32],b24[32] ;
316 char b31[32],b32[32],b33[32],b34[32] ;
317
318 UNLOAD_MAT44( imat ,
319 a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34 ) ;
320
321 MV_fval_to_char(a11,b11) ; MV_fval_to_char(a12,b12) ;
322 MV_fval_to_char(a13,b13) ; MV_fval_to_char(a14,b14) ;
323 MV_fval_to_char(a21,b21) ; MV_fval_to_char(a22,b22) ;
324 MV_fval_to_char(a23,b23) ; MV_fval_to_char(a24,b24) ;
325 MV_fval_to_char(a31,b31) ; MV_fval_to_char(a32,b32) ;
326 MV_fval_to_char(a33,b33) ; MV_fval_to_char(a34,b34) ;
327
328 sprintf( gstr ,
329 "MATRIX(%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s):%d,%d,%d" ,
330 b11,b12,b13,b14,b21,b22,b23,b24,b31,b32,b33,b34 , nx,ny,nz ) ;
331
332 return strdup(gstr) ;
333 }
334
335 /*-------------------------------------------------------------------------*/
336
EDIT_geometry_string_to_mat44(char * gstr)337 mat44_nxyz EDIT_geometry_string_to_mat44( char *gstr )
338 {
339 mat44_nxyz omat ;
340 char *cpt, *lstr ;
341 float a11,a12,a13,a14 ;
342 float a21,a22,a23,a24 ;
343 float a31,a32,a33,a34 ;
344 int nx,ny,nz , ii ;
345
346 LOAD_IDENT_MAT44(omat.mat) ; omat.nx = omat.ny = omat.nz = 0 ;
347 if( !ISVALID_GEOMETRY_STRING(gstr) ) return omat ;
348
349 lstr = strdup(gstr) ;
350
351 for( cpt=lstr ; *cpt != '\0' ; cpt++ ) if( *cpt == ',' ) *cpt = ' ' ;
352 ii = sscanf(lstr+7,"%f%f%f%f%f%f%f%f%f%f%f%f):%d%d%d",
353 &a11,&a12,&a13,&a14 ,
354 &a21,&a22,&a23,&a24 , &a31,&a32,&a33,&a34 , &nx,&ny,&nz ) ;
355 free(lstr) ;
356 if( ii < 15 ) return omat ;
357
358 LOAD_MAT44( omat.mat ,
359 a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34 ) ;
360 omat.nx = nx ; omat.ny = ny ; omat.nz = nz ;
361 return omat ;
362 }
363
364 /*-------------------------------------------------------------------------*/
365
EDIT_geometry_string_to_delxyz(char * gstr)366 float_triple EDIT_geometry_string_to_delxyz( char *gstr )
367 {
368 mat44_nxyz omat ; float_triple dxyz={0.0f,0.0f,0.0f} ;
369
370 omat = EDIT_geometry_string_to_mat44(gstr) ;
371 if( omat.nx == 0 ) return dxyz ;
372
373 dxyz.a = MAT44_COLNORM(omat.mat,0) ;
374 dxyz.b = MAT44_COLNORM(omat.mat,1) ;
375 dxyz.c = MAT44_COLNORM(omat.mat,2) ; return dxyz ;
376 }
377
378 /*-------------------------------------------------------------------------*/
379
EDIT_geometry_string_diff(char * astr,char * bstr)380 float EDIT_geometry_string_diff( char *astr , char *bstr )
381 {
382 mat44_nxyz amat , bmat ;
383 float err=0.0f ;
384
385 if( astr == NULL || bstr == NULL ) return 666.0f ; /* bad inputs */
386
387 if( strcasecmp(astr,bstr) == 0 ) return err ; /* identical strings! */
388
389 amat = EDIT_geometry_string_to_mat44(astr) ;
390 bmat = EDIT_geometry_string_to_mat44(bstr) ;
391
392 if( amat.nx != bmat.nx ) err += 1000.0f ; /* grid size mismatch! */
393 if( amat.ny != bmat.ny ) err += 1000.0f ;
394 if( amat.nz != bmat.nz ) err += 1000.0f ;
395
396 err += fabsf(amat.mat.m[0][0]-bmat.mat.m[0][0]) /* check for */
397 +fabsf(amat.mat.m[0][1]-bmat.mat.m[0][1]) /* matrix mismatch */
398 +fabsf(amat.mat.m[0][2]-bmat.mat.m[0][2])
399 +fabsf(amat.mat.m[0][3]-bmat.mat.m[0][3])
400 +fabsf(amat.mat.m[1][0]-bmat.mat.m[1][0])
401 +fabsf(amat.mat.m[1][1]-bmat.mat.m[1][1])
402 +fabsf(amat.mat.m[1][2]-bmat.mat.m[1][2])
403 +fabsf(amat.mat.m[1][3]-bmat.mat.m[1][3])
404 +fabsf(amat.mat.m[2][0]-bmat.mat.m[2][0])
405 +fabsf(amat.mat.m[2][1]-bmat.mat.m[2][1])
406 +fabsf(amat.mat.m[2][2]-bmat.mat.m[2][2])
407 +fabsf(amat.mat.m[2][3]-bmat.mat.m[2][3]) ;
408
409 return err ;
410 }
411
412 /*-------------------------------------------------------------------------*/
413 /* pad a geometry string (same in all directions) */
414
EDIT_geometry_string_pad(char * gsin,int npad)415 char * EDIT_geometry_string_pad( char *gsin , int npad )
416 {
417 char *gsout ; mat44 cmat ; mat44_nxyz cmat_nxyz ;
418 float dx,dy,dz , xorg,yorg,zorg ; int nx,ny,nz ;
419
420 if( npad <= 0 ) return NULL ;
421
422 cmat_nxyz = EDIT_geometry_string_to_mat44( gsin ) ;
423 if( cmat_nxyz.nx <= 0 ) return NULL ;
424 cmat = cmat_nxyz.mat; nx = cmat_nxyz.nx; ny = cmat_nxyz.ny; nz = cmat_nxyz.nz;
425 dx = MAT44_CLEN(cmat,0); dy = MAT44_CLEN(cmat,1); dz = MAT44_CLEN(cmat,2);
426
427 MAT44_VEC(cmat,-npad,-npad,-npad,xorg,yorg,zorg) ;
428 LOAD_MAT44_VEC(cmat,xorg,yorg,zorg) ;
429
430 gsout = EDIT_imat_to_geometry_string( cmat, nx+2*npad,ny+2*npad,nz+2*npad ) ;
431 return gsout ;
432 }
433
434 /*-------------------------------------------------------------------------*/
435 /* Return a geomstring that includes the entire volume of all of them. */
436
EDIT_geomstring_from_collection(int nstr,char ** gsin)437 char * EDIT_geomstring_from_collection( int nstr , char **gsin )
438 {
439 char *gs , *hs ;
440 int ii , ndif=0 ; float val ;
441 float xxbot,xxtop, yybot,yytop, zzbot,zztop, dxyzbot ;
442 float xxmin,xxmax, yymin,yymax, zzmin,zzmax ;
443 THD_3dim_dataset *qset ; mat44 cmat ; int nxnew,nynew,nznew ;
444
445 ENTRY("EDIT_geomstring_from_collection") ;
446
447 if( nstr <= 0 || gsin == NULL ) RETURN(NULL) ; /* ztoopid caller */
448 gs = gsin[0] ;
449 if( nstr == 1 ) RETURN(strdup(gs)) ; /* the easy case */
450
451 for( ii=1 ; ii < nstr ; ii++ ){
452 hs = gsin[ii] ; val = EDIT_geometry_string_diff(gs,hs) ;
453 if( val > 0.01f ){
454 #if 0
455 INFO_message("Warp geometries '%s' and '%s' differ by %f",gs,hs,val) ;
456 #endif
457 ndif++ ;
458 }
459 }
460 if( ndif == 0 ){ /* the easy case */
461 #if 0
462 INFO_message("No warp geometry differences") ;
463 #endif
464 RETURN(strdup(gs)) ;
465 }
466
467 /* find a box big enough to hold everybody */
468
469 xxbot = yybot = zzbot = WAY_BIG ; dxyzbot = WAY_BIG ;
470 xxtop = yytop = zztop = -WAY_BIG ;
471 for( ii=0 ; ii < nstr ; ii++ ){
472 gs = gsin[ii] ; qset = EDIT_geometry_constructor(gs,"Junk") ;
473 THD_set_dicom_box(qset->daxes) ;
474 /* get the DICOM min/max values of coords */
475 xxmin = qset->daxes->dicom_xxmin ; yymin = qset->daxes->dicom_yymin ; zzmin = qset->daxes->dicom_zzmin ;
476 xxmax = qset->daxes->dicom_xxmax ; yymax = qset->daxes->dicom_yymax ; zzmax = qset->daxes->dicom_zzmax ;
477 /* get the bot/top values of DICOM coords seen thus far */
478 if( xxmin < xxbot ) xxbot = xxmin; if( yymin < yybot ) yybot = yymin; if( zzmin < zzbot ) zzbot = zzmin;
479 if( xxmax > xxtop ) xxtop = xxmax; if( yymax > yytop ) yytop = yymax; if( zzmax > zztop ) zztop = zzmax;
480 /* get the grid spacings */
481 xxmin = fabsf(DSET_DX(qset)); yymin = fabsf(DSET_DY(qset)); zzmin = fabsf(DSET_DY(qset));
482 /* get the smallest grid spacing seen thus far */
483 if( xxmin < dxyzbot ) dxyzbot = xxmin; if( yymin < dxyzbot ) dxyzbot = yymin; if( zzmin < dxyzbot ) dxyzbot = zzmin;
484 }
485 nxnew = 1 + (int)((xxtop-xxbot)/dxyzbot) ; /* we use the smallest grid spacing */
486 nynew = 1 + (int)((yytop-yybot)/dxyzbot) ;
487 nznew = 1 + (int)((zztop-zzbot)/dxyzbot) ;
488 LOAD_MAT44(cmat , dxyzbot , 0.0f , 0.0f , xxbot ,
489 0.0f , dxyzbot , 0.0f , yybot ,
490 0.0f , 0.0f , dxyzbot , zzbot ) ;
491 gs = EDIT_imat_to_geometry_string( cmat , nxnew,nynew,nznew ) ;
492
493 RETURN(gs) ;
494 }
495
496 /*-------------------------------------------------------------------------*/
497 /* Return a geomstring that includes the DICOM corner points given. */
498
EDIT_geomstring_from_corners(float xxbot,float xxtop,float yybot,float yytop,float zzbot,float zztop,float dx,float dy,float dz)499 char * EDIT_geomstring_from_corners( float xxbot, float xxtop ,
500 float yybot, float yytop ,
501 float zzbot, float zztop ,
502 float dx, float dy, float dz )
503 {
504 char *gs ; mat44 cmat ; int nxnew,nynew,nznew ; float ttt ;
505
506 ENTRY("EDIT_geomstring_from_corners") ;
507
508 if( dx == 0.0f || dy == 0.0f || dz == 0.0f ) RETURN(NULL ) ;
509
510 if( yybot > yytop ){ ttt = yybot ; yybot = yytop ; yytop = ttt ; }
511 if( zzbot > zztop ){ ttt = zzbot ; zzbot = zztop ; zztop = ttt ; }
512
513 if( xxbot > xxtop ){ ttt = xxbot ; xxbot = xxtop ; xxtop = ttt ; }
514 nxnew = (dx > 0.0f) ? 1 + (int)((xxtop-xxbot)/dx)
515 : 1 + (int)((xxbot-xxtop)/dx) ;
516
517 if( yybot > yytop ){ ttt = yybot ; yybot = yytop ; yytop = ttt ; }
518 nynew = (dy > 0.0f) ? 1 + (int)((yytop-yybot)/dy)
519 : 1 + (int)((yybot-yytop)/dy) ;
520
521 if( zzbot > zztop ){ ttt = zzbot ; zzbot = zztop ; zztop = ttt ; }
522 nznew = (dz > 0.0f) ? 1 + (int)((zztop-zzbot)/dz)
523 : 1 + (int)((zzbot-zztop)/dz) ;
524
525 LOAD_MAT44(cmat , dx , 0.0f , 0.0f , xxbot ,
526 0.0f , dy , 0.0f , yybot ,
527 0.0f , 0.0f , dz , zzbot ) ;
528 gs = EDIT_imat_to_geometry_string( cmat , nxnew,nynew,nznew ) ;
529
530 RETURN(gs) ;
531 }
532