1 #ifndef _MRI_NWARP_HEADER_
2 #define _MRI_NWARP_HEADER_
3 
4 /********** This file defines stuff for nonlinear warping
5             in AFNIland, and is intended to be included into mrilib.h *********/
6 
7 #ifdef  __cplusplus
8 extern "C" {
9 #endif
10 
11 /*----------------------------------------------------------------------------*/
12 /* for mri_nwarp.c */
13 
14  /* The basic struct for a nonlinear warp in index space,
15     with hooks to allow it to be converted to a coordinate space 3D dataset */
16 
17  /*----------
18     April 2021 - Significant changes:
19     Modify warp struct so that the grid displacements are now stored in
20     two parts:
21       little wiggly stuff (xd,yd,zd),
22       affine transformation for long-range stuff/points outside grid (amat)
23     The mapping of a point at index (i,j,k) is to index
24            [ i + xd(i,j,k) ]
25        [M] [ j + yd(i,j,k) ] + [c]
26            [ k + zd(i,j,k) ]
27     where [amat] is the catenation [M c] of 3x3 matrix M and 3-vector c.
28 
29     The intention is that the "pure" index warp of [i,j,k] + [xd,yd,zd]
30     is bounded to the 3D grid, and for (i,j,k) outside of the grid, the
31     warp transformation is purely affine [M][i,j,k]+[c]. In the past,
32     [M] was implicitly the identity and [c] was zero -- that is, all the
33     displacements were carried by [xd,yd,zd], which lead to various
34     practical problems brought out by DRG and PAT (those devils). ----------*/
35 
36 typedef struct {
37   int    nx ,  ny ,  nz ;  /* Grid dimensions */
38   float *xd , *yd , *zd ;  /* Displacments (in index space) */
39   float *hv , *je , *se ;  /* Various auxiliary volumes */
40   int   use_amat ;         /* Whether to use amat */
41   mat44 amat , amati ;     /* Affine component of warp (and its inverse) */
42    /*- stuff below here is for conversion to/from 3D dataset format -*/
43   mat44 cmat , imat ;      /* cmat: i->x ; imat: x->i */
44   char *geomstring ;       /* grid geometry of 3D dataset */
45   int view ;               /* view/space of the dataset */
46 } IndexWarp3D ;
47 
48  /* an array of nonlinear warps [Unused at this time] */
49 
50 typedef struct {
51   int nwarp ;          /* number of them */
52   IndexWarp3D **warp ;
53 } IndexWarp3DArray ;
54 
55  /* an image plus a warp [Used in 3dQwarp = warped image plus the warp] */
56 
57 typedef struct {
58   MRI_IMAGE *im ;
59   IndexWarp3D *warp ;
60 } Image_plus_Warp ;
61 
62  /* a warp and its inverse, packaged together for fun */
63 
64 typedef struct {
65   IndexWarp3D *fwarp ;   /* forward warp */
66   IndexWarp3D *iwarp ;   /* inverse warp */
67 } IndexWarp3D_pair ;
68 
69  /* a 4x4 matrix and its inverse [Unused] */
70 
71 typedef struct {
72   mat44 fwarp ;
73   mat44 iwarp ;
74 } mat44_pair ;
75 
76  /* a collection of 4x4 matrices (e.g., time dependent rotations) */
77 
78 typedef struct { /* 17 Oct 2014 */
79   int   nmar ;
80   char  fname[128] ;
81   mat44 *mar ;
82 } mat44_vec ;
83 
84  /* get the iii-th mat44 element from a mat44_vec struct */
85 
86 #define M44V_mat(mmm,iii) ( ((iii) < (mmm)->nmar) ? (mmm)->mar[iii]             \
87                                                   : (mmm)->mar[(mmm)->nmar-1] )
88 
89  /* kill a mat44_vec struct (all of it, Frank) */
90 
91 #define DESTROY_mat44_vec(mv)                  \
92  do{ if( (mv)->mar != NULL ) free((mv)->mar) ; \
93      free(mv) ;                                \
94  } while(0) ;
95 
96  /* list of warps to catenate, some nonlinear (nwarp) and some affine (awarp) */
97 
98 typedef struct { /* 17 Oct 2014 */
99   int ncat , nvar , flags ;
100   THD_3dim_dataset **nwarp ;
101   mat44_vec        **awarp ;
102   char              *actual_geomstring ;
103   char              *master_geomstring ;
104   mat44              actual_cmat , actual_imat ;
105   int              xpad  ,ypad  ,zpad ;
106   float            xshift,yshift,zshift ;
107 } Nwarp_catlist ;
108 
109 #define NWC_INVERT_MASK 1  /* for flags field above, when catenating warps */
110 
111  /* macros to get components from a Nwarp_catlist struct */
112 
113 #define NWC_nwarp(nnn,iii) ( ((nnn)->nwarp != NULL) ? (nnn)->nwarp[iii] : NULL )
114 #define NWC_awarp(nnn,iii) ( ((nnn)->awarp != NULL) ? (nnn)->awarp[iii] : NULL )
115 #define NWC_null(nnn,iii)  ( NWC_nwarp(nnn,iii)==NULL && NWC_awarp(nnn,iii)==NULL )
116 
117  /* prototypes for warping utilities and functions */
118 
119 extern THD_3dim_dataset * IW3D_from_nwarp_catlist( Nwarp_catlist * , int ) ;
120 extern void IW3D_destroy_nwarp_catlist( Nwarp_catlist * ) ;
121 extern int IW3D_reduce_nwarp_catlist( Nwarp_catlist * ) ;
122 extern Nwarp_catlist * IW3D_read_nwarp_catlist( char * ) ;
123 extern void THD_set_nwarp_apply_prefix( char *ppp ) ; /* 15 Mar 2021 */
124 
125 extern IndexWarp3D * IW3D_create( int nx , int ny , int nz ) ;
126 extern void IW3D_destroy( IndexWarp3D *AA ) ;
127 extern float IW3D_normL1  ( IndexWarp3D *AA , IndexWarp3D *BB ) ;
128 extern float IW3D_normL2  ( IndexWarp3D *AA , IndexWarp3D *BB ) ;
129 extern float IW3D_normLinf( IndexWarp3D *AA , IndexWarp3D *BB ) ;
130 extern int_sextet IW3D_warpbox( IndexWarp3D *AA , float fac , float dthr ) ; /* 18 Mar 2021 */
131 extern IndexWarp3D * IW3D_empty_copy( IndexWarp3D *AA ) ;
132 extern IndexWarp3D * IW3D_copy( IndexWarp3D *AA , float fac ) ;
133 extern IndexWarp3D * IW3D_sum( IndexWarp3D *AA, float Afac, IndexWarp3D *BB, float Bfac ) ;
134 extern void IW3D_scale( IndexWarp3D *AA , float fac ) ;
135 extern IndexWarp3D * IW3D_from_dataset( THD_3dim_dataset *dset , int empty , int ivs ) ;
136 extern THD_3dim_dataset * IW3D_to_dataset( IndexWarp3D *AA , char *prefix ) ;
137 extern float IW3D_load_hexvol( IndexWarp3D *AA , float *hv ) ;
138 extern float IW3D_load_energy( IndexWarp3D *AA ) ;
139 extern void IW3D_load_bsv( IndexWarp3D *AA , float,float,float, float *bb , float *ss , float *vv ) ;
140 extern IndexWarp3D * IW3D_compose( IndexWarp3D *AA , IndexWarp3D *BB     , int icode ) ;
141 extern IndexWarp3D * IW3D_invert ( IndexWarp3D *AA , IndexWarp3D *BBinit , int icode ) ;
142 extern IndexWarp3D * IW3D_sqrtinv( IndexWarp3D *AA , int icode ) ;
143 extern IndexWarp3D * IW3D_sqrt   ( IndexWarp3D *AA , int icode ) ;
144 extern IndexWarp3D * IW3D_from_poly( int npar, float *par, IndexWarp3D *WW ) ;
145 extern THD_3dim_dataset * NwarpCalcRPN( char *expr, char *prefix, int icode, int acode ) ;
146 extern void NwarpCalcRPN_verb(int i) ;
147 
148 extern void IW3D_bounding_box_clear( IndexWarp3D *AA , float thresh ) ; /* 29 Mar 2021 */
149 extern int_sextet IW3D_bounding_box( IndexWarp3D *AA , float thresh ) ;
150 
151 extern void THD_interp_floatim( MRI_IMAGE *fim ,
152                                 int np , float *ip , float *jp , float *kp ,
153                                 int code, float *outar ) ;
154 extern void THD_interp_complexim( MRI_IMAGE *fim ,
155                                   int np , float *ip , float *jp , float *kp ,
156                                   int code, complex *outar ) ; /* 27 Mar 2018 */
157 extern MRI_IMARR * THD_setup_nwarp( MRI_IMARR *bimar,
158                                     int use_amat    , mat44 amat ,
159                                     mat44 cmat_bim  ,
160                                     int incode      , float wfac ,
161                                     mat44 cmat_src  ,
162                                     mat44 cmat_out  ,
163                                     int nx_out      , int ny_out , int nz_out  ) ;
164 extern THD_3dim_dataset * THD_nwarp_dataset( THD_3dim_dataset *dset_nwarp ,
165                                              THD_3dim_dataset *dset_src  ,
166                                              THD_3dim_dataset *dset_mast ,
167                                              char *prefix , int wincode , int dincode ,
168                                              float dxyz_mast , float wfac , int nvlim ,
169                                              MRI_IMAGE *amatim ) ;
170 
171 extern THD_3dim_dataset * THD_nwarp_dataset_NEW( Nwarp_catlist    *nwc       ,
172                                                  THD_3dim_dataset *dset_src  ,
173                                                  THD_3dim_dataset *dset_mast ,
174                                                  char *prefix, int wincode, int dincode,
175                                                  float dxyz_mast, float wfac, int nvlim ) ;
176 
177 extern int THD_nwarp_forward_xyz( THD_3dim_dataset *dset_nwarp ,
178                                   float dfac , int npt ,
179                                   float *xin , float *yin , float *zin ,
180                                   float *xut , float *yut , float *zut  ) ;
181 
182 extern int THD_nwarp_inverse_xyz( THD_3dim_dataset *dset_nwarp ,
183                                   float dfac , int npt ,
184                                   float *xin , float *yin , float *zin ,
185                                   float *xut , float *yut , float *zut  ) ;
186 /*----------------------------------------------------------------------------*/
187 
188 #ifdef  __cplusplus
189 }
190 #endif
191 
192 #endif /* _MCW_MRILIB_HEADER_ */
193