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 substitute one brick array for another in a 3D dataset.
11 - dset = dataset to be edited
12 - ival = index of sub-brick to be replaced
13 - ftype = datum type of substitute array
14 - fim = array of substitute data -- it may be NULL. If it is not NULL,
15 then should contain mri_datum_size(ftype) * nxx*nyy*nzz bytes.
16 * If it is NULL, then space will be calloc-ed.
17 * Do NOT free this memory after this function - the
18 data in the dataset IS this memory, not a copy of it!
19
20 Notes:
21 - The original brick (an MRI_IMAGE within an MRI_IMARR) is deleted
22 from memory. A new brick is put in its place.
23 - This can only be done on a brick that is malloc-ed, not mmap-ed!
24 - The brick_bytes and total_bytes field of the datablock are
25 patched, but the brick_fac field is NOT changed here.
26 -----------------------------------------------------------------------------*/
27
EDIT_substitute_brick(THD_3dim_dataset * dset,int ival,int ftype,void * fim)28 void EDIT_substitute_brick(THD_3dim_dataset *dset, int ival, int ftype,void *fim)
29 {
30 THD_datablock *dblk ;
31 MRI_IMAGE *newim , *oldim ;
32 int nullfim = (fim==NULL) ;
33 int64_t nbytes ;
34
35 ENTRY("EDIT_substitute_brick") ;
36
37 /**-- Sanity Checks --**/
38
39 if( ! ISVALID_3DIM_DATASET(dset) ) EXRETURN; /* error! */
40 if( dset->dblk->brick == NULL ) EXRETURN; /* ditto! */
41 if( dset->dblk->malloc_type != DATABLOCK_MEM_MALLOC )EXRETURN; /* ditto! */
42 if( ival >= dset->dblk->nvals || ival < 0 ) EXRETURN; /* ditto! */
43 if( ftype < 0 || ftype > LAST_MRI_TYPE ) EXRETURN; /* ditto! */
44
45 oldim = DSET_BRICK(dset,ival) ; if( oldim == NULL ) EXRETURN; /* ditto! */
46
47 newim = mri_empty_conforming( oldim , ftype ) ; /* new sub-brick */
48 nbytes = (int64_t)newim->nvox
49 * (int64_t)newim->pixel_size ; /* how big it is */
50 mri_free( oldim ) ; /* kill old one */
51
52 if( nullfim ){ /* if needed, */
53 fim = calloc( 1,nbytes ) ; /* make array */
54 if( fim == NULL )
55 ERROR_exit("malloc (out of memory) error for dataset sub-brick #%d",ival) ;
56 }
57 mri_fix_data_pointer( fim , newim ) ; /* attach new data */
58 DSET_BRICK(dset,ival) = newim ; /* put in dataset */
59
60 /** change the byte count for this sub-brick and the total dataset **/
61
62 dset->dblk->total_bytes += (nbytes - dset->dblk->brick_bytes[ival]) ;
63 dset->dblk->brick_bytes[ival] = nbytes ;
64
65 DSET_CRUSH_BSTAT(dset,ival) ;
66
67 THD_patch_dxyz_one(dset,ival) ; /* 05 Jun 2007 */
68
69 EXRETURN ;
70 }
71
72 /*---------------------------------------------------------------------------*/
73 /*! Similar to EDIT_substitute_brick(), but allows the data type to be
74 converted/scaled. The brick_fac field will be set. [18 Sep 2006]
75
76 Free fim after the function returns UNLESS ftype == stype
77 -----------------------------------------------------------------------------*/
78
EDIT_substscale_brick(THD_3dim_dataset * dset,int ival,int ftype,void * fim,int stype,float fac)79 void EDIT_substscale_brick(THD_3dim_dataset *dset, int ival,
80 int ftype,void *fim , int stype,float fac )
81 {
82 float *far ;
83 double *dar;
84 int ii,nvox ;
85
86 ENTRY("EDIT_substscale_brick") ;
87
88 /**-- Sanity Checks --**/
89
90 if( ! ISVALID_3DIM_DATASET(dset) ) EXRETURN; /* error! */
91 if( dset->dblk->brick == NULL ) EXRETURN; /* ditto! */
92 if( dset->dblk->malloc_type != DATABLOCK_MEM_MALLOC )EXRETURN; /* ditto! */
93 if( ival >= dset->dblk->nvals || ival < 0 ) EXRETURN; /* ditto! */
94 if( ftype < 0 || ftype > LAST_MRI_TYPE ) EXRETURN; /* ditto! */
95
96 /** the trivial operation? **/
97
98 if( stype < 0 || stype > LAST_MRI_TYPE || stype == ftype ){
99 EDIT_substitute_brick( dset,ival,ftype,fim ) ;
100 EDIT_BRICK_FACTOR( dset,ival,0.0f ) ;
101 EXRETURN ;
102 }
103 if( fim == NULL ){
104 EDIT_substitute_brick( dset,ival,stype,NULL ) ;
105 EDIT_BRICK_FACTOR( dset,ival,0.0f ) ;
106 EXRETURN ;
107 }
108
109 /** at this time, can only scale float inputs to shorts or bytes **/
110
111 if( ftype != MRI_float && ftype != MRI_double ){ /* ZSS Dec 2010 */
112 ERROR_message("EDIT_substscale_brick: non-float and non-double input! "
113 "(%d %d %d)", ftype, stype, MRI_short);
114 EXRETURN;
115 }
116 if( stype != MRI_short && stype != MRI_byte ){
117 ERROR_message("EDIT_substscale_brick: non-short/byte output!"); EXRETURN;
118 }
119
120 if (ftype == MRI_float) {
121 far = (float *)fim ; nvox = DSET_NVOX(dset) ;
122
123 /** compute factor, if not supplied by user **/
124
125 if( fac <= 0.0f ){
126 float bot,top , abot,atop,mmm ; int isin ;
127 bot = top = far[0] ;
128 for( ii=1 ; ii < nvox ; ii++ ){
129 if( far[ii] < bot ) bot = far[ii] ;
130 else if( far[ii] > top ) top = far[ii] ;
131 }
132 abot = fabsf(bot); atop = fabsf(top); mmm = MAX(abot,atop);
133 if( mmm == 0.0f ){ /** data values are all zero! **/
134 fac = 1.0f ;
135 } else if( stype == MRI_short ){
136 isin = is_integral_data( nvox , MRI_float , far ) ;
137 fac = (isin && mmm <= 32767.0f) ? 1.0f : 32767.0f / mmm ;
138 } else if( stype == MRI_byte ){
139 if( bot < 0.0f ){
140 for( ii=0 ; ii < nvox ; ii++ ) if( far[ii] < 0.0f ) far[ii] = 0.0f ;
141 }
142 if( top > 0.0f ){
143 isin = is_integral_data( nvox , MRI_float , far ) ;
144 fac = (isin && top <= 255.0f) ? 1.0f : 255.0f / top ;
145 } else {
146 WARNING_message("EDIT_substscale_brick: "
147 "no positive data for -> byte");
148 fac = 1.0f ;
149 }
150 }
151 }
152
153 /*-- now do the scaling and substitution --*/
154
155 if( stype == MRI_short ){
156
157 short *sar = (short *)malloc(sizeof(short)*nvox) ;
158 float val ;
159 if( fac == 1.0f ){
160 STATUS("storing to shorts") ;
161 for( ii=0 ; ii < nvox ; ii++ ){
162 sar[ii] = SHORTIZE(far[ii]) ;
163 }
164 } else {
165 STATUS("scaling to shorts") ;
166 for( ii=0 ; ii < nvox ; ii++ ){
167 val = fac*far[ii] ; sar[ii] = SHORTIZE(val) ;
168 }
169 }
170 STATUS("putting into dataset") ;
171 EDIT_substitute_brick( dset,ival,MRI_short,sar ) ;
172 fac = (fac==1.0f) ? 0.0f : 1.0f/fac ;
173 STATUS("setting scale factor") ;
174 EDIT_BRICK_FACTOR( dset,ival,fac ) ;
175
176 EDIT_misfit_report( DSET_FILECODE(dset), ival,
177 nvox , fac , sar , far ) ;
178
179 } else if( stype == MRI_byte ){
180
181 byte *bar = (byte *)malloc(sizeof(byte)*nvox) ;
182 float val ;
183 if( fac == 1.0f ){
184 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = BYTEIZE(far[ii]) ;
185 } else {
186 for( ii=0 ; ii < nvox ; ii++ ){
187 val = fac*far[ii] ; bar[ii] = BYTEIZE(val) ;
188 }
189 }
190 EDIT_substitute_brick( dset,ival,MRI_byte,bar ) ;
191 fac = (fac==1.0f) ? 0.0f : 1.0f/fac ;
192 EDIT_BRICK_FACTOR( dset,ival,fac ) ;
193 }
194 } else if (ftype == MRI_double) { /* ZSS Dec 2010 */
195 dar = (double *)fim ; nvox = DSET_NVOX(dset) ;
196
197 /** compute factor, if not supplied by user **/
198
199 if( fac <= 0.0f ){
200 double bot,top , abot,atop,mmm ; int isin ;
201 bot = top = dar[0] ;
202 for( ii=1 ; ii < nvox ; ii++ ){
203 if( dar[ii] < bot ) bot = dar[ii] ;
204 else if( dar[ii] > top ) top = dar[ii] ;
205 }
206 abot = fabs(bot); atop = fabs(top); mmm = MAX(abot,atop);
207 if( mmm == 0.0f ){ /** data values are all zero! **/
208 fac = 1.0f ;
209 } else if( stype == MRI_short ){
210 isin = is_integral_data( nvox , MRI_double , dar ) ;
211 fac = (isin && mmm <= 32767.0f) ? 1.0f : 32767.0f / mmm ;
212 } else if( stype == MRI_byte ){
213 if( bot < 0.0f ){
214 for( ii=0 ; ii < nvox ; ii++ ) if( dar[ii] < 0.0f ) dar[ii] = 0.0f ;
215 }
216 if( top > 0.0f ){
217 isin = is_integral_data( nvox , MRI_double , dar ) ;
218 fac = (isin && top <= 255.0f) ? 1.0f : 255.0f / top ;
219 } else {
220 WARNING_message("EDIT_substscale_brick: "
221 "no positive data for -> byte");
222 fac = 1.0f ;
223 }
224 }
225 }
226
227 /*-- now do the scaling and substitution --*/
228
229 if( stype == MRI_short ){
230
231 short *sar = (short *)malloc(sizeof(short)*nvox) ;
232 double val ;
233 if( fac == 1.0f ){
234 STATUS("storing to shorts") ;
235 for( ii=0 ; ii < nvox ; ii++ ){
236 sar[ii] = SHORTIZE(dar[ii]) ;
237 }
238 } else {
239 STATUS("scaling to shorts") ;
240 for( ii=0 ; ii < nvox ; ii++ ){
241 val = fac*dar[ii] ; sar[ii] = SHORTIZE(val) ;
242 }
243 }
244 STATUS("putting into dataset") ;
245 EDIT_substitute_brick( dset,ival,MRI_short,sar ) ;
246 fac = (fac==1.0f) ? 0.0f : 1.0f/fac ;
247 STATUS("setting scale factor") ;
248 EDIT_BRICK_FACTOR( dset,ival,fac ) ;
249
250 /* Skip the misfit report. None exists for double
251 input at the moment.
252 EDIT_misfit_report( ) ;*/
253
254 } else if( stype == MRI_byte ){
255
256 byte *bar = (byte *)malloc(sizeof(byte)*nvox) ;
257 double val ;
258 if( fac == 1.0f ){
259 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = BYTEIZE(dar[ii]) ;
260 } else {
261 for( ii=0 ; ii < nvox ; ii++ ){
262 val = fac*dar[ii] ; bar[ii] = BYTEIZE(val) ;
263 }
264 }
265 EDIT_substitute_brick( dset,ival,MRI_byte,bar ) ;
266 fac = (fac==1.0f) ? 0.0f : 1.0f/fac ;
267 EDIT_BRICK_FACTOR( dset,ival,fac ) ;
268 }
269 }
270
271 EXRETURN ;
272 }
273
274 /*---------------------------------------------------------------------------*/
275 static char EDIT_index_prefix = '#' ;
EDIT_set_index_prefix(char c)276 void EDIT_set_index_prefix(char c){ EDIT_index_prefix = c ;
277 #if 0
278 INFO_message("EDIT_set_index_prefix('%c') --> %c",c,EDIT_get_index_prefix() ) ;
279 #endif
280 } /* 11 Jul 2019 */
281
EDIT_get_index_prefix(void)282 char EDIT_get_index_prefix(void) { return EDIT_index_prefix ; } /* 11 Jul 2019 */
283 /*---------------------------------------------------------------------------*/
284