1 #include "mrilib.h"
2
3 #include "thd_dset_to_grayplot.c"
4 #ifdef USE_OMP
5 # include "cs_pv.c"
6 #endif
7
8 void show_help(void) ;
9
10 #define IS_NUMERIC(sss) \
11 ( ( isdigit((sss)[0]) ) || \
12 ( (sss)[0] == '.' && isdigit((sss)[1]) ) || \
13 ( (sss)[0] == '-' && isdigit((sss)[1]) ) || \
14 ( (sss)[0] == '-' && (sss)[1] == '.' && isdigit((sss)[2]) ) )
15
16 /*===========================================================================*/
17
main(int argc,char * argv[])18 int main( int argc , char *argv[] )
19 {
20 THD_3dim_dataset *dset=NULL ;
21 MRI_IMAGE *imout ;
22 char *prefix = "Grayplot.png" ;
23 byte *mask=NULL ; int mask_nvox=0 ;
24 int iarg = 1 ;
25 int polort=2 , nxout=0,nyout=0 ; float fwhm=0.0f ;
26 float bnd_bot = 1000000.0, bnd_top = -1000000.0;
27 int DO_PERCENT = 0, DO_RANGE = 0, DO_RAW = 0;
28
29 /*----------------------------------------------------------------*/
30
31 if( argc < 2 ){
32 show_help() ;
33 exit(0) ;
34 }
35
36 /*----------------------------------------------------------------*/
37
38 mainENTRY("3dGrayplot"); machdep(); (void)COX_clock_time();
39 if( strcasecmp(argv[1],"-help") != 0 ){
40 PRINT_VERSION("3dGrayplot") ;
41 }
42
43 while( iarg < argc && argv[iarg][0] == '-' ){
44
45 /*-----*/
46
47 if( strcasecmp(argv[iarg],"-help") == 0 ){
48 show_help() ;
49 exit(0) ;
50 }
51
52 /*-----*/
53
54 if( strcasecmp(argv[iarg],"-dimen") == 0 ){
55 if( ++iarg >= argc-1 )
56 ERROR_exit("'-dimen' needs 2 arguments") ;
57 nxout = (int)strtod(argv[iarg++],NULL) ;
58 nyout = (int)strtod(argv[iarg++],NULL) ;
59 continue ;
60 }
61
62 /*-----*/
63
64 if( strcasecmp(argv[iarg],"-polort") == 0 ){
65 if( ++iarg >= argc )
66 ERROR_exit("'-polort' needs an argument") ;
67 polort = (int)strtod(argv[iarg++],NULL) ;
68 continue ;
69 }
70
71 /*-----*/
72
73 if( strcasecmp(argv[iarg],"-pvorder") == 0 ){ /* 04 May 2018 */
74 grayplot_order_by_pvmap(1) ; iarg++ ; continue ;
75 }
76 if( strcasecmp(argv[iarg],"-peelorder") == 0 ){ /* 08 May 2018 */
77 grayplot_order_by_peels(1) ; iarg++ ; continue ;
78 }
79 if( strcasecmp(argv[iarg],"-ijkorder") == 0 ){ /* 09 May 2018 */
80 grayplot_order_by_ijk(1) ; iarg++ ; continue ;
81 }
82 if( strcasecmp(argv[iarg],"-LJorder") == 0 ){ /* 05 Feb 2020 */
83 grayplot_order_by_lj(1) ; iarg++ ; continue ;
84 }
85
86 /*-----*/
87
88 if( strncasecmp(argv[iarg],"-oldresam",7) == 0 ||
89 strncasecmp(argv[iarg],"-resamold",7) == 0 ){ /* 14 Aug 2018 */
90 ININFO_message("Using old (deprecated) resampling method") ;
91 grayplot_set_use_old_resam() ; iarg++ ; continue ;
92 }
93
94 /*-----*/
95
96 if( strcasecmp(argv[iarg],"-fwhm") == 0 ){
97 if( ++iarg >= argc )
98 ERROR_exit("'-fwhm' needs an argument") ;
99 fwhm = (float)strtod(argv[iarg++],NULL) ;
100 continue ;
101 }
102
103 /*-----*/
104
105 if( strcasecmp(argv[iarg],"-prefix") == 0 ){
106 if( ++iarg >= argc )
107 ERROR_exit("'-prefix' needs an argument") ;
108 if( STRING_HAS_SUFFIX(argv[iarg],".png") ||
109 STRING_HAS_SUFFIX(argv[iarg],".pgm") ||
110 STRING_HAS_SUFFIX(argv[iarg],".jpg") ){
111 prefix = strdup(argv[iarg]) ;
112 } else {
113 prefix = malloc(sizeof(char)*(strlen(argv[iarg])+16)) ;
114 sprintf(prefix,"%s.png",argv[iarg]) ;
115 }
116 if( !THD_filename_ok(prefix) )
117 ERROR_exit("prefix '%s' is illegal :(",prefix) ;
118 iarg++ ; continue ;
119 }
120
121 /*-----*/
122
123 if( strcasecmp(argv[iarg],"-mask") == 0 ){
124 THD_3dim_dataset *mset ; int mmm ;
125 if( ++iarg >= argc ) ERROR_exit("Need argument after '-mask'") ;
126 if( mask != NULL ) ERROR_exit("Can't have two mask inputs") ;
127 mset = THD_open_dataset( argv[iarg] ) ;
128 CHECK_OPEN_ERROR(mset,argv[iarg]) ;
129 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
130 mask_nvox = DSET_NVOX(mset) ;
131 if( DSET_BRICK_TYPE(mset,0) != MRI_byte &&
132 DSET_BRICK_TYPE(mset,0) != MRI_short )
133 ERROR_exit("-mask dataset is NOT byte- or short-valued :(") ;
134
135 if( DSET_BRICK_TYPE(mset,0) == MRI_byte ){
136 mask = DSET_BRICK_ARRAY(mset,0) ;
137 } else {
138 short *smm = DSET_BRICK_ARRAY(mset,0) ;
139 mask = (byte *)calloc(sizeof(byte),mask_nvox) ;
140 EDIT_coerce_type( mask_nvox , MRI_short,smm , MRI_byte,mask ) ;
141 DSET_unload(mset) ;
142 }
143
144 mmm = THD_countmask( mask_nvox , mask ) ;
145 ININFO_message("Number of voxels in mask = %d",mmm) ;
146 if( mmm < 19 ) ERROR_exit("Mask is too small to process") ;
147 iarg++ ; continue ;
148 }
149
150 /*-----*/
151
152 if( strcasecmp(argv[iarg],"-input") == 0 ){
153 if( ++iarg >= argc ) ERROR_exit("Need argument after '-input'") ;
154 if( dset != NULL ) ERROR_exit("Can't have two -input options") ;
155 dset = THD_open_dataset( argv[iarg] ) ;
156 CHECK_OPEN_ERROR(dset,argv[iarg]) ;
157 if( DSET_NVALS(dset) < 9 )
158 ERROR_exit("Dataset %s has too few time points",DSET_HEADNAME(dset)) ;
159 ININFO_message("Loading dataset %s",DSET_HEADNAME(dset)) ;
160 DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
161 iarg++ ; continue ;
162 }
163
164 /*-----*/
165
166 if( strcasecmp(argv[iarg],"-percent") == 0 ){ /* 30 Jul 2018 */
167 grayplot_set_percent() ;
168 grayplot_norming_none() ;
169 DO_PERCENT = 1;
170 iarg++ ; continue ;
171 }
172
173 /*-----*/
174
175 if( strcasecmp(argv[iarg],"-range") == 0 ){
176 if( ++iarg < argc && IS_NUMERIC(argv[iarg]) ){
177 float val = (float)strtod(argv[iarg],NULL) ;
178 if( val > 0.0f ) grayplot_set_range(val) ;
179 DO_RANGE = 1;
180 iarg++ ;
181 }
182 continue ;
183 }
184
185 /*-----*/
186
187 // [PT: Sep 20, 2021] all display of raw data, with user-selected
188 // interval
189 if( strcasecmp(argv[iarg],"-raw_with_bounds") == 0 ){
190 if( ++iarg < argc && IS_NUMERIC(argv[iarg]) ){
191 bnd_bot = (float)strtod(argv[iarg], NULL) ;
192 }
193 else
194 ERROR_exit("Need first arg after '-raw_with_bounds' to be numeric\n"
195 "\t-> this doesn't look numeric: %s", argv[iarg]);
196 if( ++iarg < argc && IS_NUMERIC(argv[iarg]) ){
197 bnd_top = (float)strtod(argv[iarg], NULL) ;
198 INFO_message("OK: %s", argv[iarg]);
199 }
200 else
201 ERROR_exit("Need second arg after '-raw_with_bounds' to be numeric\n"
202 "\t-> this doesn't look numeric: %s", argv[iarg]);
203
204 if( bnd_top > bnd_bot )
205 grayplot_set_raw_range(bnd_bot, bnd_top) ;
206 else
207 ERROR_exit("Prob with '-raw_with_bounds BOT TOP': TOP <= BOT\n"
208 "\t-> Please check your input again; current range is:\n"
209 "\t %.6f %.6f", bnd_bot, bnd_top) ;
210 grayplot_norming_none() ; // don't scale/transform data
211 DO_RAW = 1;
212
213 iarg++ ;
214 continue ;
215 }
216
217 /*-----*/
218
219 ERROR_exit("Unknown option '%s'",argv[iarg]) ;
220 }
221
222 /*----------------------------------------------------------------*/
223
224 if ( (DO_RAW && DO_PERCENT) || (DO_RAW && DO_RANGE) ){
225 ERROR_exit("Do not combine '-raw_with_bounds ..' with either "
226 "'-percent' or '-range ..'");
227 }
228
229 /*----- no dataset given yet? -----*/
230
231 if( dset == NULL ){
232 if( iarg >= argc ) ERROR_exit("No input dataset?") ;
233 dset = THD_open_dataset( argv[iarg] ) ;
234 CHECK_OPEN_ERROR(dset,argv[iarg]) ;
235 if( DSET_NVALS(dset) < 9 )
236 ERROR_exit("Dataset %s has too few time points",DSET_HEADNAME(dset)) ;
237 ININFO_message("Loading dataset %s",DSET_HEADNAME(dset)) ;
238 DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
239 }
240
241 /*----- create default mask -----*/
242
243 if( mask == NULL ){
244 int ii,jj,kk , nx,ny,nz,nxy , nx1,ny1,nz1 , vv ;
245 ININFO_message("No mask input ==> using all interior voxels") ;
246 mask_nvox = DSET_NVOX(dset) ;
247 mask = (byte *)malloc(sizeof(byte)*mask_nvox) ;
248 nx = DSET_NX(dset); ny = DSET_NY(dset); nz = DSET_NZ(dset); nxy = nx*ny;
249 nx1 = nx-1 ; ny1 = ny-1 ; nz1 = nz-1 ;
250 /* allow surface datasets, if not peel ordering [18 Jun 2019 rickr] */
251 if( grayplot_is_ordered_by_peels() && (nx1 == 0 || ny1 == 0 || nz1 == 0) )
252 ERROR_exit("Cannot peel automask on dataset with only one slice :(") ;
253 for( vv=kk=0 ; kk < nz ; kk++ ){
254 for( jj=0 ; jj < ny ; jj++ ){
255 for( ii=0 ; ii < nx ; ii++,vv++){
256 mask[vv] = (ii==0 || jj==0 || kk==0 || ii==nx1 || jj==ny1 || kk==nz1)
257 ? 0 : 1 ;
258 }}}
259 vv = THD_countmask( mask_nvox , mask ) ;
260 ININFO_message("Number of voxels in mask = %d",vv) ;
261 } else if( mask_nvox != DSET_NVOX(dset) ){
262 ERROR_exit("mask and dataset voxel counts don't match :(") ;
263 }
264
265 /*------- do all the work -------*/
266
267 imout = THD_dset_to_grayplot( dset,mask , nxout,nyout , polort,fwhm ) ;
268
269 mri_write_png( prefix , imout ) ;
270
271 ININFO_message("3dGrayplot: Elapsed = %.1f s\n", COX_clock_time() ) ;
272
273 exit(0) ;
274 }
275
276 /*===========================================================================*/
277
show_help(void)278 void show_help(void)
279 {
280 printf("\n"
281 "Make a grayplot from a 3D+time dataset, sort of like Jonathan Power:\n"
282 " https://www.ncbi.nlm.nih.gov/pubmed/27510328\n"
283 " https://www.jonathanpower.net/2017-ni-the-plot.html\n"
284 "Result is saved to a PNG image for your viewing delight.\n"
285 "\n"
286 " * This style of plot is also called a carpet plot,\n"
287 " but REAL carpets are much more attractive, IMHO.\n"
288 "\n"
289 " * The horizontal axis of the grayplot is time, and the\n"
290 " vertical axis is all 3 spatial dimensions collapsed into 1.\n"
291 "\n"
292 " * Also see AFNI script @grayplot, as well as the QC output\n"
293 " generated by afni_proc.py.\n"
294 "\n"
295 "Usage:\n"
296 "\n"
297 " 3dGrayplot [options] inputdataset\n"
298 "\n"
299 "OPTIONS: [lots of them]\n"
300 "--------\n"
301 "\n"
302 " -mask mset = Name of mask dataset\n"
303 " * Voxels that are 0 in mset will not be processed.\n"
304 " * Dataset must be byte-valued (8 bits: 0..255);\n"
305 " shorts (16 bits) are also acceptable, but only\n"
306 " values from 1.255 will be processed.\n"
307 " * Each distinct value from 1..255 will be processed\n"
308 " separately, and the output image will be ordered\n"
309 " with the mask=1 voxels on top, mask=2 voxels next,\n"
310 " and so on down the image.\n"
311 " * A partition (e.g., mask=3) with fewer than 9 voxels\n"
312 " will not be processed at all.\n"
313 " * If there is more than one partition, horizontal dashed\n"
314 " lines will drawn between them.\n"
315 " * If '-mask' is not given, then all voxels will be used,\n"
316 " except those at the very edge of a volume. Doing this is\n"
317 " usually not a good idea, as the non-brain tissue will\n"
318 " take up a lot of useless space in the output grayplot.\n"
319 "\n"
320 " -input dataset = Alternative way to input the dataset to process.\n"
321 "\n"
322 " -prefix ppp.png = Name for output file.\n"
323 " * Default is Grayplot.png (if you don't use this option)\n"
324 " * If the filename ends in '.jpg', a JPEG file is output.\n"
325 " * If the filename ends in '.pgm', a PGM file is output.\n"
326 " [PGM files can be manipulated with the NETPBM package.]\n"
327 " * If the filename does not end in '.jpg' OR in '.png'\n"
328 " OR in '.pgm', then '.png' will be added at the end.\n"
329 "\n"
330 " -dimen X Y = Output size of image in pixels.\n"
331 " * X = width = time axis direction\n"
332 " * Y = height = voxel/space dimensions\n"
333 " * Defaults are X=1024 Y=512 -- suitable for screen display.\n"
334 " * For publication, you might want more pixels, as in\n"
335 " -dimen 1800 1200\n"
336 " which would be 6 inches wide by 4 inches high, at the usual\n"
337 " 300 dots-per-inch (dpi) of high resolution image printing.\n"
338 " ** Note that there are usually many more voxels in the Y direction\n"
339 " than there are pixels in the output image. This fact requires\n"
340 " coarsening the Y output grid and resampling the data to match.\n"
341 " See the next option for a little more information about\n"
342 " how this resampling is implemented.\n"
343 "\n"
344 " -oldresam = The method for resampling the processed dataset to the final\n"
345 " grayscale image size was changed/improved in a major way.\n"
346 " If you want to use the original method, then give this option.\n"
347 " * The only reason for using this option is for\n"
348 " comparison with the new method.\n"
349 " * The new resampling method uses minimum-sidelobe local averaging\n"
350 " when coarsening the grid (vertical direction Y = voxels/space)\n"
351 " -- whose purpose is to reduce aliasing artifacts\n"
352 " * And uses cubic interpolation when refining the grid\n"
353 " (usually horizontal direction = time) -- whose purpose\n"
354 " is purely beauty -- compared to the older linear interpolation.\n"
355 " * Note that the collapsing of multiple voxels into one pixel in\n"
356 " the Y direction will tend to cancel out signals that change sign\n"
357 " within neighbors in the voxel ordering method you choose.\n"
358 " (See the 'order' options below.)\n"
359 "\n"
360 " -polort p = Order of polynomials for detrending.\n"
361 " * Default value is 2 (mean, slope, quadratic curve).\n"
362 " * Use '-1' if data is already detrended and de-meaned.\n"
363 " (e.g., is an AFNI errts.* file or other residual dataset)\n"
364 "\n"
365 " -fwhm f = FWHM of blurring radius to use in the dataset before\n"
366 " making the image.\n"
367 " * Each partition (i.e., mask=1, mask=2, ...) is blurred\n"
368 " independently, as in program 3dBlurInMask.\n"
369 " * Default value is 0 mm = no blurring.\n"
370 " [In the past, the default value was 6.]\n"
371 " * If the dataset was NOT previously blurred, a little\n"
372 " spatial blurring here will help bring out larger scale\n"
373 " features in the times series, which might otherwise\n"
374 " look very noisy.\n"
375 "\n"
376 " ** The following four options control the ordering of **\n"
377 " ** voxels in the grayplot, in the vertical direction. **\n"
378 "\n"
379 " -pvorder = Within each mask partition, order the voxels (top to\n"
380 " bottom) by how well they match the two leading principal\n"
381 " components of that partition. The result is to make the\n"
382 " top part of each partition be made up of voxels with\n"
383 " similar time series, and the bottom part will be more\n"
384 " 'random looking'.\n"
385 " ++ The presence of a lot of temporal structure in a\n"
386 " grayplot of a 'errts' residual dataset indicates\n"
387 " that the 'removal' of unwanted time series components\n"
388 " did not work well.\n"
389 " ++ Using '-pvorder' to put all the structured time series\n"
390 " close together will make such problems more visible.\n"
391 " ++ IMHO, this is the most useful ordering.\n"
392 "\n"
393 " -LJorder = Within each mask partition, order the voxels (top to\n"
394 " bottom) by their Ljung-Box statistics, which is a measure\n"
395 " of temporal correlation.\n"
396 " ++ Experimental; probably not useful.\n"
397 "\n"
398 " -peelorder = Within each mask partition, order the voxels by how\n"
399 " many 'peel' steps are needed to get from the partition\n"
400 " boundary to a given voxel.\n"
401 " ++ This ordering puts voxels in 'similar' geometrical\n"
402 " positions sort-of close together in the image.\n"
403 " And is usually not very interesting, IMHO.\n"
404 "\n"
405 " -ijkorder = Set the intra-partition ordering to the default, by\n"
406 " dataset 3D index ('ijk').\n"
407 " ++ In AFNI's +tlrc ordering, this ordering primarily will\n"
408 " be from Inferior to Superior in the brain (from top to\n"
409 " bottom in the grayplot image).\n"
410 " ++ This is the default ordering method, but not the best.\n"
411 "\n"
412 " ** These options control the scaling from voxel value to gray level **\n"
413 "\n"
414 " -range X = Set the range of the data to be plotted to be 'X'.\n"
415 " Each time series is first normalized by its values to:\n"
416 " Z[i] = (t[i] - mean_t)/stdev_t.\n"
417 " When this option is used, then:\n"
418 " * a value of 0 will be plotted as middle-gray\n"
419 " * a value of +X (or above) will be plotted as white\n"
420 " * a value of -X (or below) will be plotted as black\n"
421 " Thus, this option should be used with data that is centered\n"
422 " around zero -- or will be so after '-polort' detrending.\n"
423 " * For example, if you are applying this option to an\n"
424 " afni_proc.py 'errts' (residuals) dataset, a good value\n"
425 " of X to use is 3 or 4, since those values are in percents.\n"
426 " * The @grayplot script uses '-range 3.89' since that is the\n"
427 " value at which a standard normal N(0,1) deviate has a 1e-4\n"
428 " two-sided tail probability. (If nothing else, this sounds cool.)\n"
429 " If you do NOT use '-range', then the data will be automatically\n"
430 " normalized so each voxel time series has RMS value 1, and then\n"
431 " the grayscale plot will be black-to-white being the min-to-max,\n"
432 " where the min and max computed over the entire detrended\n"
433 " and normalized dataset.\n"
434 " * This default automatic normalizing and scaling makes it\n"
435 " almost impossible to directly compare grayplots from\n"
436 " different datasets. This difficulty is why the '-range'\n"
437 " and '-percent' options were added.\n"
438 "\n"
439 " -percent = Use this option on 'raw' time series datasets, to compute\n"
440 " the mean of each voxel timeseries and then use that value\n"
441 " to scale the values to percent differences from the mean.\n"
442 " * NOT suitable for use with a residual 'errts' dataset!\n"
443 " * Should be combined with '-range'.\n"
444 " * Detrending will be applied while calculating the mean.\n"
445 " By default, that will be quadratic detrending of each\n"
446 " voxel time series, but that can be changed with the\n"
447 " '-polort' option.\n"
448 "\n"
449 " -raw_with_bounds A B\n"
450 " = Use this option on 'raw' time series datasets, map values\n"
451 " <= A to black, those >= B to white, and intermediate values\n"
452 " to grays.\n"
453 " * Can be used with any kind of dataset, but probably makes\n"
454 " most sense to use with scaled ones (errts, fitts or\n"
455 " all_runs).\n"
456 " * Should NOT be combined with '-range' or '-percent'.\n"
457 "\n"
458 "** Quick hack for Cesar Caballero-Gaudes, April 2018, by @AFNIman.\n"
459 " As such, this program may be modified in the future to be more useful,\n"
460 " or at least more beautifully gorgeous.\n"
461 "\n"
462 "** Applied to 'raw' EPI data, the results may not be very informative.\n"
463 " It seems to be more useful to look at the grayplot calculated from\n"
464 " pre-processed data (e.g., time series registered, filtered, etc.).\n"
465 "\n"
466 "** See also the script @grayplot, which can process the results from\n"
467 " afni_proc.py and produce an image with the grayplot combined with\n"
468 " a graph of the motion magnitude, and with the GM, WM, and CSF in\n"
469 " different partitions.\n"
470 "\n"
471 "** afni_proc.py uses this program to create grayplots of the residuals\n"
472 " from regression analysis, as part of its Quality Control (QC) output.\n"
473 "\n"
474 "--------\n"
475 "EXAMPLE:\n"
476 "--------\n"
477 "The following commands first generate a time series dataset,\n"
478 "then create grayplots using each of the ordering methods\n"
479 "(so you can compare them). No mask is given.\n"
480 "\n"
481 " 3dcalc -a jRandomDataset:64:64:30:256 -datum float \\\n"
482 " -prefix Qsc.nii -expr 'abs(.3+cos(0.1*i))*sin(0.1*t+0.1*i)+gran(0,3)'\n"
483 " 3dGrayplot -pvorder -prefix QscPV.png -input Qsc.nii -fwhm 8\n"
484 " 3dGrayplot -ijkorder -prefix QscIJK.png -input Qsc.nii -fwhm 8\n"
485 " 3dGrayplot -peelorder -prefix QscPEEL.png -input Qsc.nii -fwhm 8\n"
486 "\n"
487 "\n"
488 ) ;
489 }
490