1 #include "mrilib.h"
2
3 static int verb = 0 ;
mri_polyfit_verb(int ii)4 void mri_polyfit_verb( int ii ){ verb = ii ; }
5
6 /*----------------------------------------------------------------------------*/
7 /* define some polynomials over the domain [-1..1] */
8
9 #undef LP1
10 #undef LP2
11 #undef LP3
12 #undef LP4
13 #undef LP5
14 #undef LP6
15 #undef LP7
16 #undef LP8
17 #undef LP9
18
19 #define LP1(x) (x)
20 #define LP2(x) ((x)*(x)-0.3333333f)
21 #define LP3(x) (((x)*(x)-0.6f)*(x))*1.5f
22 #define LP4(x) ((x)*(x)*((x)*(x)-0.857143f)+0.0857143f)*2.0f
23 #define LP5(x) (((x)*(x)*((x)*(x)-1.11111f)+0.238095f)*(x))*3.0f
24 #define LP6(x) ((x)*(x)*((x)*(x)*((x)*(x)-1.36364f)+0.454545f)-0.021645f)*6.0f
25 #define LP7(x) (((x)*(x)*((x)*(x)*((x)*(x)-1.61538f)+0.734266f)-0.081585f)*(x))*10.0f
26
27 #define LP8(x) ( (x)*(x) * \
28 ( (x)*(x) * \
29 ( (x)*(x) * \
30 ( (x)*(x)-1.86667f )+1.07692f )-0.195804f )+0.0054390f ) * 18.0f
31
32 #define LP9(x) ( ( (x)*(x) * \
33 ( (x)*(x) * \
34 ( (x)*(x) * \
35 ( (x)*(x)-2.11765f )+1.48235f )-0.380090f )+0.0259153f ) \
36 * (x) )*32.0f
37
38 /*----------------------------------------------------------------------------*/
39 /* Evaluate product polynomial px(x)*py(y)*pz(z) at a bunch of points. */
40
poly3D(int px,int py,int pz,int nxyz,float * x,float * y,float * z,float * val)41 static void poly3D( int px, int py, int pz,
42 int nxyz, float *x, float *y, float *z, float *val )
43 {
44 register int ii ;
45
46 ENTRY("poly3D") ;
47
48 switch( px ){
49 default: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = 1.0f ; break ;
50 case 1: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = LP1(x[ii]) ; break ;
51 case 2: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = LP2(x[ii]) ; break ;
52 case 3: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = LP3(x[ii]) ; break ;
53 case 4: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = LP4(x[ii]) ; break ;
54 case 5: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = LP5(x[ii]) ; break ;
55 case 6: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = LP6(x[ii]) ; break ;
56 case 7: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = LP7(x[ii]) ; break ;
57 case 8: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = LP8(x[ii]) ; break ;
58 case 9: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = LP9(x[ii]) ; break ;
59 }
60
61 switch( py ){
62 case 1: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP1(y[ii]) ; break ;
63 case 2: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP2(y[ii]) ; break ;
64 case 3: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP3(y[ii]) ; break ;
65 case 4: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP4(y[ii]) ; break ;
66 case 5: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP5(y[ii]) ; break ;
67 case 6: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP6(y[ii]) ; break ;
68 case 7: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP7(y[ii]) ; break ;
69 case 8: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP8(y[ii]) ; break ;
70 case 9: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP9(y[ii]) ; break ;
71 }
72
73 switch( pz ){
74 case 1: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP1(z[ii]) ; break ;
75 case 2: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP2(z[ii]) ; break ;
76 case 3: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP3(z[ii]) ; break ;
77 case 4: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP4(z[ii]) ; break ;
78 case 5: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP5(z[ii]) ; break ;
79 case 6: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP6(z[ii]) ; break ;
80 case 7: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP7(z[ii]) ; break ;
81 case 8: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP8(z[ii]) ; break ;
82 case 9: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= LP9(z[ii]) ; break ;
83 }
84
85 EXRETURN ;
86 }
87
88 /*----------------------------------------------------------------------------*/
89 /* Evaluate product Hermite functions. */
90
91 #define HP1(x) (x)
92 #define HP2(x) ((x)*(x)-2.0f)
93 #define HP3(x) ((x)*(2.0f*(x)*(x)-3.0f))
94 #define HP4(x) ((4.0f*(x)*(x)-12.0f)*(x)*(x)+3.0f)
95 #define HP5(x) ((x)*((4.0f*(x)*(x)-20.0f)*(x)*(x)+15.0f))
96 #define HP6(x) (((2.0f*(x)*(x)-15.0f)*(x)*(x)+22.5f)*(x)*(x)-3.75f)
97 #define HP7(x) (0.5f*(x)*(((8.0f*(x)*(x)-84.0f)*(x)*(x)+210.0f)*(x)*(x)-105.0f))
98 #define HP8(x) (((((x)*(x)-14.0f)*(x)*(x)+52.5f)*(x)*(x)-52.5f)*(x)*(x)+6.5625f)
99 #define HP9(x) (0.02f*(x)*((((16.0f*(x)*(x)-288.0f)*(x)*(x)+1512.0f)*(x)*(x)-2520.0f)*(x)*(x)+945.0f))
100
101 #define HH1(x) HP1(3.0f*(x))*exp(-(x)*(x))
102 #define HH2(x) HP2(3.0f*(x))*exp(-2.0f*(x)*(x))
103 #define HH3(x) HP3(3.0f*(x))*exp(-3.0f*(x)*(x))
104 #define HH4(x) HP4(3.0f*(x))*exp(-5.0f*(x)*(x))
105 #define HH5(x) HP5(3.0f*(x))*exp(-6.0f*(x)*(x))
106 #define HH6(x) HP6(3.0f*(x))*exp(-6.0f*(x)*(x))
107 #define HH7(x) HP7(3.0f*(x))*exp(-7.0f*(x)*(x))
108 #define HH8(x) HP8(3.0f*(x))*exp(-7.0f*(x)*(x))
109 #define HH9(x) HP9(3.0f*(x))*exp(-7.0f*(x)*(x))
110
herm3D(int px,int py,int pz,int nxyz,float * x,float * y,float * z,float * val)111 static void herm3D( int px, int py, int pz,
112 int nxyz, float *x, float *y, float *z, float *val )
113 {
114 register int ii ;
115
116 ENTRY("herm3D") ;
117
118 switch( px ){
119 default: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = 1.0f ; break ;
120 case 1: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = HH1(x[ii]) ; break ;
121 case 2: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = HH2(x[ii]) ; break ;
122 case 3: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = HH3(x[ii]) ; break ;
123 case 4: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = HH4(x[ii]) ; break ;
124 case 5: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = HH5(x[ii]) ; break ;
125 case 6: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = HH6(x[ii]) ; break ;
126 case 7: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = HH7(x[ii]) ; break ;
127 case 8: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = HH8(x[ii]) ; break ;
128 case 9: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] = HH9(x[ii]) ; break ;
129 }
130
131 switch( py ){
132 case 1: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH1(y[ii]) ; break ;
133 case 2: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH2(y[ii]) ; break ;
134 case 3: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH3(y[ii]) ; break ;
135 case 4: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH4(y[ii]) ; break ;
136 case 5: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH5(y[ii]) ; break ;
137 case 6: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH6(y[ii]) ; break ;
138 case 7: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH7(y[ii]) ; break ;
139 case 8: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH8(y[ii]) ; break ;
140 case 9: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH9(y[ii]) ; break ;
141 }
142
143 switch( pz ){
144 case 1: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH1(z[ii]) ; break ;
145 case 2: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH2(z[ii]) ; break ;
146 case 3: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH3(z[ii]) ; break ;
147 case 4: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH4(z[ii]) ; break ;
148 case 5: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH5(z[ii]) ; break ;
149 case 6: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH6(z[ii]) ; break ;
150 case 7: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH7(z[ii]) ; break ;
151 case 8: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH8(z[ii]) ; break ;
152 case 9: for( ii=0 ; ii < nxyz ; ii++ ) val[ii] *= HH9(z[ii]) ; break ;
153 }
154
155 EXRETURN ;
156 }
157
158 /*----------------------------------------------------------------------------*/
159
160 static void (*BFUN3D)(int,int,int,int,float *,float *,float*,float *) = poly3D ;
161
mri_polyfit_set_basis(char * str)162 void mri_polyfit_set_basis( char *str )
163 {
164 if( str != NULL && strcasecmp(str,"Hermite") == 0 )
165 BFUN3D = herm3D ;
166 else
167 BFUN3D = poly3D ; /* default */
168 return ;
169 }
170
171 /*----------------------------------------------------------------------------*/
172 /* for getting the fit coefficients [26 Feb 2019] */
173
174 static floatvec *pfit_vec = NULL ;
mri_polyfit_get_fitvec(void)175 floatvec * mri_polyfit_get_fitvec(void){ return pfit_vec ; }
176
177 /*----------------------------------------------------------------------------*/
178 /* Fit a polynomial to a 3D image and return the fitted image:
179 nord = maximum order of polynomial (0..9)
180 mask = optional mask of voxels to actually use
181 mrad = optional preliminary median filter radius (voxels)
182 meth = 2 for least squares fit, 1 for L1 fit
183 *//*--------------------------------------------------------------------------*/
184
mri_polyfit(MRI_IMAGE * imin,int nord,MRI_IMARR * exar,byte * mask,float mrad,int meth)185 MRI_IMAGE * mri_polyfit( MRI_IMAGE *imin, int nord, MRI_IMARR *exar, byte *mask, float mrad, int meth )
186 {
187 MRI_IMAGE *fim ;
188 float *far,*fdat , **ref=NULL , *xx,*yy,*zz,*vv ; floatvec *fvit ;
189 int *pir=NULL , *pjr=NULL , *pkr=NULL ;
190 int nmask , ii,jj,kk,pp,qq , nx,ny,nz,nxyz , ibot,itop,jbot,jtop,kbot,ktop ;
191 int pi,pj,pk , pitop=0,pjtop=0,pktop=0 , nref=0 , npref=0 ;
192 float imid=0.0,jmid=0.0,kmid=0.0 , ifac=0.0f,jfac=0.0f,kfac=0.0f , rfac ;
193 int nex=0 ;
194
195 ENTRY("mri_polyfit") ;
196
197 KILL_floatvec(pfit_vec) ;
198
199 if( imin == NULL || (nord < 0 && exar == NULL) || nord > 9 ) RETURN(NULL) ;
200
201 if( nord < -1 ) nord = -1 ;
202
203 ii = mri_dimensionality(imin) ;
204 if( ii < 2 || ii > 3 ){
205 if( verb ) ERROR_message("mri_polyfit: dimensionality = %d",ii) ;
206 RETURN(NULL) ;
207 }
208
209 if( exar != NULL ) nex = IMARR_COUNT(exar) ;
210
211 if( mask != NULL ) nmask = THD_countmask( imin->nvox , mask ) ;
212 else nmask = imin->nvox ;
213
214 jj = 8*( nex + (int)rint( pow( (double)(nord+1) , (double)ii ) ) ) ;
215 if( nmask < jj ){
216 if( verb ) ERROR_message("mri_polyfit: #points=%d < %d",nmask,jj) ;
217 RETURN(NULL) ;
218 }
219
220 #undef GOOD
221 #define GOOD(i) (mask == NULL || mask[i])
222
223 nx = imin->nx ; ny = imin->ny ; nz = imin->nz ; nxyz = nx*ny*nz ;
224
225 /* find min and max values of indexes in the mask */
226
227 if( mask == NULL ){
228 ibot = jbot = kbot = 0 ; itop = nx-1 ; jtop = ny-1 ; ktop = nz-1 ;
229 } else {
230 ibot = jbot = kbot = 9999999 ; itop = jtop = ktop = -9999999 ;
231 for( qq=kk=0 ; kk < nz ; kk++ ){
232 for( jj=0 ; jj < ny ; jj++ ){
233 for( ii=0 ; ii < nx ; ii++,qq++ ){
234 if( mask[qq] ){
235 if( ii < ibot ) ibot = ii ;
236 if( jj < jbot ) jbot = jj ;
237 if( kk < kbot ) kbot = kk ;
238 if( ii > itop ) itop = ii ;
239 if( jj > jtop ) jtop = jj ;
240 if( kk > ktop ) ktop = kk ;
241 }
242 }}}
243 if( ((ibot>=itop) + (jbot>=jtop) + (kbot>=ktop)) > 1 ) RETURN(NULL) ;
244 }
245
246 if( nord >= 0 ){
247 imid = 0.5f*(ibot+itop) ; jmid = 0.5f*(jbot+jtop) ; kmid = 0.5f*(kbot+ktop) ;
248
249 if( ibot < itop ){ kk = (itop-ibot)/3; pitop = MIN(nord,kk); ifac = 2.0f/(itop-ibot); }
250 if( jbot < jtop ){ kk = (jtop-jbot)/3; pjtop = MIN(nord,kk); jfac = 2.0f/(jtop-jbot); }
251 if( kbot < ktop ){ kk = (ktop-kbot)/3; pktop = MIN(nord,kk); kfac = 2.0f/(ktop-kbot); }
252
253 xx = (float *)malloc(sizeof(float)*nmask) ;
254 yy = (float *)malloc(sizeof(float)*nmask) ;
255 zz = (float *)malloc(sizeof(float)*nmask) ;
256 for( pp=qq=kk=0 ; kk < nz ; kk++ ){
257 for( jj=0 ; jj < ny ; jj++ ){
258 for( ii=0 ; ii < nx ; ii++,qq++ ){
259 if( GOOD(qq) ){
260 xx[pp] = (ii-imid)*ifac ;
261 yy[pp] = (jj-jmid)*jfac ;
262 zz[pp] = (kk-kmid)*kfac ; pp++ ;
263 }
264 }}}
265
266 pir = (int *)malloc(sizeof(int)*(pktop+1)*(pjtop+1)*(pitop+1)) ;
267 pjr = (int *)malloc(sizeof(int)*(pktop+1)*(pjtop+1)*(pitop+1)) ;
268 pkr = (int *)malloc(sizeof(int)*(pktop+1)*(pjtop+1)*(pitop+1)) ;
269 for( nref=pk=0 ; pk <= pktop ; pk++ ){
270 for( pj=0 ; pj <= pjtop ; pj++ ){
271 for( pi=0 ; pi <= pitop ; pi++ ){
272 if( pi+pj+pk <= nord ){
273 pir[nref] = pi ; pjr[nref] = pj ; pkr[nref] = pk ; nref++ ;
274 }
275 }}}
276 npref = nref ;
277
278 if( verb )
279 ININFO_message("create %d polynomial regressors with %d points each",nref,nmask) ;
280
281 ref = (float **)malloc(sizeof(float *)*nref) ;
282 for( kk=0 ; kk < nref ; kk++ ){
283 ref[kk] = (float *)malloc(sizeof(float)*nmask) ;
284 BFUN3D( pir[kk],pjr[kk],pkr[kk] , nmask , xx,yy,zz , ref[kk] ) ;
285 }
286 free(zz) ; free(yy) ; free(xx) ;
287 }
288
289 if( nex > 0 ){ /* 27 Dec 2012 */
290 float *ex, *rr ;
291 if( verb ) ININFO_message("add %d image regressors",nex) ;
292 ref = (float **)realloc(ref,sizeof(float *)*(npref+nex)) ;
293 for( kk=0 ; kk < nex ; kk++ ){
294 ex = MRI_FLOAT_PTR(IMARR_SUBIM(exar,kk)) ;
295 rr = ref[npref+kk] = (float *)malloc(sizeof(float)*nmask) ;
296 for( qq=ii=0 ; ii < nxyz ; ii++ ){ if( GOOD(ii) ) rr[qq++] = ex[ii]; }
297 }
298 nref = npref + nex ;
299 }
300
301 fim = mri_to_float(imin) ; /* convert input to float */
302 if( mrad > 0.0f ){
303 MRI_IMAGE *qim ;
304 if( verb ) ININFO_message("median filter input image") ;
305 qim = mri_medianfilter(fim,mrad,mask,verb) ;
306 if( qim != NULL ){ mri_free(fim) ; fim = qim ; }
307 }
308 far = MRI_FLOAT_PTR(fim) ; /* will later be output image */
309
310 if( nmask < nxyz ){
311 fdat = (float *)malloc(sizeof(float)*nmask) ;
312 for( qq=ii=0 ; ii < nxyz ; ii++ ){ if( GOOD(ii) ) fdat[qq++] = far[ii]; }
313 } else {
314 fdat = far ;
315 }
316
317 if( meth < 1 ) meth = 1 ;
318 else if( meth > 2 ) meth = 2 ;
319 if( verb ) ININFO_message("L%d fit polynomial field to data",meth) ;
320
321 fvit = THD_fitter( nmask , fdat , nref , ref , meth , NULL ) ;
322
323 if( fdat != far ) free(fdat) ;
324 for( kk=0 ; kk < nref ; kk++ ) free(ref[kk]) ;
325 free(ref) ;
326
327 if( fvit == NULL ){
328 if( pir != NULL ){ free(pir); free(pjr); free(pkr); }
329 ERROR_message("Can't calculate fit in mri_polyfit?! :-(") ;
330 RETURN(NULL) ;
331 }
332
333 COPY_floatvec(pfit_vec,fvit) ; /* save coefficients */
334
335 memset( far , 0 , sizeof(float)*nxyz ) ;
336
337 if( nord >= 0 ){
338 if( verb ) ININFO_message("compute final fit polynomial field") ;
339 xx = (float *)malloc(sizeof(float)*nxyz) ;
340 yy = (float *)malloc(sizeof(float)*nxyz) ;
341 zz = (float *)malloc(sizeof(float)*nxyz) ;
342 vv = (float *)malloc(sizeof(float)*nxyz) ;
343 for( qq=kk=0 ; kk < nz ; kk++ ){
344 for( jj=0 ; jj < ny ; jj++ ){
345 for( ii=0 ; ii < nx ; ii++,qq++ ){
346 xx[qq] = (ii-imid)*ifac ;
347 yy[qq] = (jj-jmid)*jfac ;
348 zz[qq] = (kk-kmid)*kfac ;
349 }}}
350
351 for( kk=0 ; kk < npref ; kk++ ){
352 BFUN3D( pir[kk],pjr[kk],pkr[kk] , nxyz , xx,yy,zz , vv ) ;
353 rfac = fvit->ar[kk] ;
354 for( ii=0 ; ii < nxyz ; ii++ ) far[ii] += rfac * vv[ii] ;
355 }
356 free(pir); free(pjr); free(pkr); free(vv); free(zz); free(yy); free(xx);
357 }
358
359 if( nex > 0 ){
360 float *ex ;
361 for( kk=0 ; kk < nex ; kk++ ){
362 ex = MRI_FLOAT_PTR(IMARR_SUBIM(exar,kk)) ;
363 rfac = fvit->ar[kk+npref] ;
364 for( ii=0 ; ii < nxyz ; ii++ ) far[ii] += rfac * ex[ii] ;
365 }
366 }
367
368 KILL_floatvec(fvit) ;
369
370 RETURN(fim) ;
371 }
372
373 /*----------------------------------------------------------------------------*/
374
mri_polyfit_byslice(MRI_IMAGE * imin,int nord,MRI_IMARR * exar,byte * mask,float mrad,int meth)375 MRI_IMAGE * mri_polyfit_byslice( MRI_IMAGE *imin, int nord, MRI_IMARR *exar,
376 byte *mask, float mrad, int meth )
377 {
378 MRI_IMAGE *sl_imin , *sl_exim , *sl_out , *out_im ;
379 MRI_IMARR *sl_exar=NULL , *out_imar=NULL ;
380 byte *sl_mask=NULL ;
381 int nx=imin->nx , ny=imin->ny , nz=imin->nz , nex=0 ;
382 int kk ;
383
384 ENTRY("mri_polyfit_byslice") ;
385
386 if( nz == 1 ){
387 out_im = mri_polyfit( imin,nord,exar , mask,mrad,meth ) ;
388 RETURN(out_im) ;
389 }
390
391 INIT_IMARR(out_imar) ;
392
393 for( kk=0 ; kk < nz ; kk++ ){
394 sl_imin = mri_cut_3D( imin , 0,nx-1 , 0,ny-1 , kk,kk ) ;
395 if( mask != NULL ) sl_mask = mask + kk*nx*ny ;
396 if( exar != NULL ){
397 int nex=IMARR_COUNT(exar) , ee ; MRI_IMAGE *eim ;
398 INIT_IMARR(sl_exar) ;
399 for( ee=0 ; ee < nex ; ee++ ){
400 eim = mri_cut_3D( IMARR_SUBIM(exar,ee) , 0,nx-1 , 0,ny-1 , kk,kk ) ;
401 ADDTO_IMARR(sl_exar,eim) ;
402 }
403 }
404 sl_out = mri_polyfit( sl_imin,nord,sl_exar , sl_mask,mrad,meth ) ;
405 ADDTO_IMARR(out_imar,sl_out) ;
406 if( sl_exar != NULL ) DESTROY_IMARR(sl_exar) ;
407 mri_free(sl_imin) ;
408 }
409
410 out_im = mri_catvol_1D( out_imar , 3 ) ;
411 DESTROY_IMARR(out_imar) ;
412 RETURN(out_im) ;
413 }
414