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
9 /*---------------------------------------------------------------------------*/
10 /*! Routine to attach sub-bricks to the end of an already existing 3D dataset.
11 - nbr = number of extra bricks
12 - tbr = array of data types in each brick (defaults to MRI_short)
13 - fbr = array of new brick_fac factors in each brick
14 if NULL, new factors will be set to zero
15 - sbr = array of pointers to data in each brick
16 if sbr[i]==NULL, then the i'th extra brick will not have data
17 if sbr==NULL, then all new extra bricks will not have data
18 (in which case you will have to use mri_fix_data_pointer to
19 attach data at a later time, if you plan to use the data)
20
21 Note that this can only be done on a brick that is malloc-ed, not mmap-ed!
22 -----------------------------------------------------------------------------*/
23
EDIT_add_bricklist(THD_3dim_dataset * dset,int nbr,int * tbr,float * fbr,void * sbr[])24 void EDIT_add_bricklist( THD_3dim_dataset *dset ,
25 int nbr, int *tbr, float *fbr , void *sbr[] )
26 {
27 int ibr , typ , nx,ny,nz , nvals,new_nvals ;
28 THD_datablock *dblk ;
29 MRI_IMAGE *qim ;
30 char str[32] ;
31
32 ENTRY("EDIT_add_bricklist") ;
33
34 /**-- Sanity Checks --**/
35
36 if( ! ISVALID_3DIM_DATASET(dset) || nbr <= 0 ) EXRETURN; /* error! */
37 if( dset->dblk->brick == NULL ) EXRETURN; /* error! */
38 if( dset->dblk->malloc_type != DATABLOCK_MEM_MALLOC )EXRETURN; /* error! */
39
40 dblk = dset->dblk ;
41 nvals = dblk->nvals ;
42 nx = dblk->diskptr->dimsizes[0] ;
43 ny = dblk->diskptr->dimsizes[1] ;
44 nz = dblk->diskptr->dimsizes[2] ;
45
46 /**-- reallocate the brick control information --**/
47
48 new_nvals = nvals + nbr ;
49 dblk->brick_bytes = (int64_t *) RwcRealloc( (char *) dblk->brick_bytes ,
50 sizeof(int64_t) * new_nvals ) ;
51
52 dblk->brick_fac = (float *) RwcRealloc( (char *) dblk->brick_fac ,
53 sizeof(float) * new_nvals ) ;
54
55 dblk->nvals = dblk->diskptr->nvals = new_nvals ;
56
57 /** allocate new sub-brick images **/
58
59 for( ibr=0 ; ibr < nbr ; ibr++ ){
60 typ = (tbr != NULL ) ? tbr[ibr] : MRI_short ;
61 qim = mri_new_vol_empty( nx,ny,nz , typ ) ; /* image with no data */
62
63 if( sbr != NULL && sbr[ibr] != NULL ) /* attach data to image */
64 mri_fix_data_pointer( sbr[ibr] , qim ) ;
65
66 ADDTO_IMARR( dblk->brick , qim ) ; /* attach image to dset */
67
68 dblk->brick_fac[nvals+ibr] = (fbr != NULL) ? fbr[ibr] : 0.0 ;
69 dblk->brick_bytes[nvals+ibr] = (int64_t)qim->pixel_size * (int64_t)qim->nvox ;
70 dblk->total_bytes += dblk->brick_bytes[ibr] ;
71 }
72
73 /** allocate new sub-brick auxiliary data: labels **/
74
75 if( dblk->brick_lab == NULL )
76 THD_init_datablock_labels( dblk ) ;
77 else
78 dblk->brick_lab = (char **) RwcRealloc( (char *) dblk->brick_lab ,
79 sizeof(char *) * new_nvals ) ;
80 for( ibr=0 ; ibr < nbr ; ibr++ ){
81 sprintf( str , "#%d" , nvals+ibr ) ;
82 dblk->brick_lab[nvals+ibr] = NULL ;
83 THD_store_datablock_label( dblk , nvals+ibr , str ) ;
84 }
85
86 /** keywords **/
87
88 if( dblk->brick_keywords == NULL )
89 THD_init_datablock_keywords( dblk ) ;
90 else
91 dblk->brick_keywords = (char **) RwcRealloc( (char *) dblk->brick_keywords ,
92 sizeof(char *) * new_nvals ) ;
93 for( ibr=0 ; ibr < nbr ; ibr++ ){
94 dblk->brick_keywords[nvals+ibr] = NULL ;
95 THD_store_datablock_keywords( dblk , nvals+ibr , NULL ) ;
96 }
97
98 /** stataux **/
99
100 if( dblk->brick_statcode != NULL ){
101 dblk->brick_statcode = (int *) RwcRealloc( (char *) dblk->brick_statcode ,
102 sizeof(int) * new_nvals ) ;
103 dblk->brick_stataux = (float **) RwcRealloc( (char *) dblk->brick_stataux ,
104 sizeof(float *) * new_nvals ) ;
105
106 for( ibr=0 ; ibr < nbr ; ibr++ ){
107 dblk->brick_statcode[nvals+ibr] = 0 ;
108 dblk->brick_stataux[nvals+ibr] = NULL ;
109 }
110 }
111
112 /** fdrcurve **/
113
114 if( dblk->brick_fdrcurve != NULL ){
115 dblk->brick_fdrcurve = (floatvec **) realloc( (void *)dblk->brick_fdrcurve ,
116 sizeof(floatvec *) * new_nvals ) ;
117 for( ibr=0 ; ibr < nbr ; ibr++ ) dblk->brick_fdrcurve[nvals+ibr] = NULL ;
118 }
119
120 if( dblk->brick_mdfcurve != NULL ){ /* 22 Oct 2008 */
121 dblk->brick_mdfcurve = (floatvec **) realloc( (void *)dblk->brick_mdfcurve ,
122 sizeof(floatvec *) * new_nvals ) ;
123 for( ibr=0 ; ibr < nbr ; ibr++ ) dblk->brick_mdfcurve[nvals+ibr] = NULL ;
124 }
125
126 EXRETURN;
127 }
128
129 /*---------------------------------------------------------------------------*/
130 /*! Add one brick to the end of a dataset. */
131
EDIT_add_brick(THD_3dim_dataset * dset,int typ,float fac,void * br)132 void EDIT_add_brick( THD_3dim_dataset *dset, int typ , float fac , void *br )
133 {
134 int ttt = typ ;
135 float fff = fac ;
136 void *bbb = br ;
137
138 EDIT_add_bricklist( dset , 1 , &ttt , &fff , &bbb ) ;
139 return ;
140 }
141
142 /*---------------------------------------------------------------------------*/
143 /*!
144 Turn float arrays into sub-bricks of a preset type
145 (based on code in 3dMean)
146
147 dset (THD_3dim_dataset *) new dset to which arrays will be added
148 far (float **) each far[i] is to become one sub-brick in dset
149 nval (int) the number of arrays in far
150 otype (int) the sub-brick type. Supported options are:
151 MRI_float (so far
152 scaleopt (char) scaling options:
153 'A' scale if needed
154 'F' do scale each sub-brick
155 'G' scale all sub-bricks with the same factor
156 'N' Do not scale
157 verb (int) loquaciousness
158 returns 1 if all is well
159 0 all hell broke loose
160
161 */
EDIT_add_bricks_from_far(THD_3dim_dataset * dset,float ** far,int nval,int otype,char scaleopt,int verb)162 int EDIT_add_bricks_from_far(THD_3dim_dataset *dset,
163 float **far, int nval,
164 int otype, char scaleopt,
165 int verb)
166 {
167 int ii=0, kk=0, nxyz;
168
169 ENTRY("EDIT_add_bricks_from_far");
170
171 if (scaleopt != 'A' && scaleopt != 'F' && scaleopt != 'G' && scaleopt != 'N'){
172 ERROR_message("Bad scaleopt value of %c", scaleopt);
173 RETURN(0);
174 }
175
176 if (!dset) {
177 ERROR_message("NULL input");
178 RETURN(0);
179 }
180
181 nxyz = DSET_NVOX(dset);
182
183 switch( otype ){
184
185 default:
186 ERROR_message("Somehow ended up with otype = %d\n",otype) ;
187 RETURN(0) ;
188
189 case MRI_float:{
190 for( kk=0 ; kk < nval ; kk++ ){
191 EDIT_substitute_brick(dset, kk, MRI_float, far[kk]);
192 DSET_BRICK_FACTOR(dset, kk) = 0.0;
193 far[kk] = NULL;
194 }
195 }
196 break ;
197
198 case MRI_byte:
199 case MRI_short:{
200 void ** dfim ;
201 float gtop=0.0 , fimfac , gtemp ;
202
203 if( verb )
204 fprintf(stderr," ++ Scaling output to type %s brick(s)\n",
205 MRI_TYPE_name[otype] ) ;
206
207 dfim = (void **) malloc(sizeof(void *)*nval) ;
208
209 if( scaleopt == 'G' ){ /* allow global scaling */
210 gtop = 0.0 ;
211 for( kk=0 ; kk < nval ; kk++ ){
212 gtemp = MCW_vol_amax( nxyz , 1 , 1 , MRI_float, far[kk] ) ;
213 gtop = MAX( gtop , gtemp ) ;
214 if( gtemp == 0.0 )
215 WARNING_message("output sub-brick %d is all zeros!\n",kk) ;
216 }
217 }
218
219 for (kk = 0 ; kk < nval ; kk ++ ) {
220
221 if( scaleopt != 'G' && scaleopt != 'N'){
222 /* compute max value in this sub-brick */
223 gtop = MCW_vol_amax( nxyz , 1 , 1 , MRI_float, far[kk] ) ;
224 if( gtop == 0.0 )
225 WARNING_message("output sub-brick %d is all zeros!\n",kk) ;
226
227 }
228
229 if( scaleopt == 'F' || scaleopt == 'G'){ /* scaling needed */
230
231 fimfac = (gtop > 0.0) ? MRI_TYPE_maxval[otype] / gtop : 0.0 ;
232
233 } else if( scaleopt == 'A' ){ /* only if needed */
234
235 fimfac = ( gtop > MRI_TYPE_maxval[otype] ||
236 (gtop > 0.0 && gtop < 1.0) )
237 ? MRI_TYPE_maxval[otype]/ gtop : 0.0 ;
238
239 if( fimfac == 0.0 && gtop > 0.0 ){ /* 14 May 2010 */
240 float fv,iv ; /* force scaling if */
241 for( ii=0 ; ii < nxyz ; ii++ ){ /* non-integers are inside */
242 fv = far[kk][ii] ; iv = rint(fv) ;
243 if( fabsf(fv-iv) >= 0.01 ){
244 fimfac = MRI_TYPE_maxval[otype] / gtop ; break ;
245 }
246 }
247 }
248
249 } else if( scaleopt == 'N') { /* no scaling allowed */
250 fimfac = 0.0 ;
251 } else {
252 ERROR_message("Should not see this one");
253 RETURN(0);
254 }
255
256
257 if( verb ){
258 if( fimfac != 0.0 )
259 INFO_message("Sub-brick %d scale factor = %f\n",kk,fimfac) ;
260 else
261 INFO_message("Sub-brick %d: no scale factor\n" ,kk) ;
262 }
263
264 dfim[kk] = (void *) malloc( mri_datum_size(otype) * nxyz ) ;
265 if( dfim[kk] == NULL ){
266 ERROR_message("malloc fails at output\n");
267 exit(1);
268 }
269
270 EDIT_coerce_scale_type( nxyz , fimfac ,
271 MRI_float, far[kk] , otype,dfim[kk] ) ;
272 if( otype == MRI_short )
273 EDIT_misfit_report( DSET_FILECODE(dset) , kk ,
274 nxyz , (fimfac != 0.0f) ? 1.0f/fimfac : 0.0f ,
275 dfim[kk] , far[kk] ) ;
276 free( far[kk] ) ; far[kk] = NULL;
277 EDIT_substitute_brick(dset, kk, otype, dfim[kk] );
278
279 DSET_BRICK_FACTOR(dset,kk) = (fimfac != 0.0) ? 1.0/fimfac : 0.0 ;
280
281 dfim[kk]=NULL;
282 }
283 free(dfim); dfim = NULL;
284 }
285 break ;
286 }
287
288
289 RETURN(1);
290 }
291