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 "afni.h"
8 #include "maxima.h"
9
10 char grMessage[ R_MESSAGE_L ];
11
12 /*----------------------------------------------------------------------
13 **
14 ** Find local maxima, storing the extrema as a mask.
15 **
16 **----------------------------------------------------------------------
17 */
process_data(maxima_s * M)18 int process_data( maxima_s * M )
19 {
20 ENTRY("process_data");
21 ( void )find_local_maxima( M );
22
23 if ( ! create_point_list( M ) )
24 RETURN(0);
25
26 gr_orig_data = M->sdata; /* global needed for sorting */
27 if ( M->negatives )
28 qsort( M->P.plist, M->P.used, sizeof( int ), point_comp_neg );
29 else
30 qsort( M->P.plist, M->P.used, sizeof( int ), point_comp_pos );
31
32 if ( M->debug > 1 )
33 show_point_list_s( "+d point list sorted: ", &M->P, M->debug );
34
35 if ( ( M->min_dist > 1.0 ) && ! apply_min_dist( M ) )
36 RETURN(0);
37
38 if ( M->debug > 1 )
39 show_point_list_s( "+d point list cleaned: ", &M->P, M->debug );
40
41 if ( M->outfile )
42 apply_fill_radius( M );
43
44 RETURN(1);
45 }
46
47
48 /*----------------------------------------------------------------------
49 ** Display the contents of the point_list_s struct.
50 **----------------------------------------------------------------------
51 */
show_point_list_s(char * mesg,point_list_s * p,int debug)52 void show_point_list_s( char * mesg, point_list_s * p, int debug )
53 {
54 int c;
55
56 ENTRY("show_point_list_s");
57
58 if ( mesg ) fputs( mesg, stderr );
59
60 fprintf(stderr, "point_list_s @ %p, used = %d, M = %d\n",
61 (void *)p, p->used, p->M);
62
63 if ( debug <= 0 ) EXRETURN; /* we're done */
64
65 fprintf(stderr," plist starting @ %p:", (void *)p->plist );
66
67 for ( c = 0; c < p->used; c++ )
68 fprintf(stderr," %d", p->plist[c] );
69 fprintf(stderr,"\n");
70
71 EXRETURN;
72 }
73
74
75 /*----------------------------------------------------------------------
76 **
77 ** Display the contents of the maxima structure.
78 **
79 **----------------------------------------------------------------------
80 */
show_maxima_s(char * mesg,maxima_s * M)81 void show_maxima_s( char * mesg, maxima_s * M )
82 {
83 ENTRY("show_maxima_s");
84
85 if ( mesg ) fputs( mesg, stderr );
86
87 fprintf( stderr,
88 "------------------------------\n"
89 "dset * : %p\n"
90 "sdata * : %p\n"
91 "result * : %p\n"
92 "nx : %d\n"
93 "ny : %d\n"
94 "nz : %d\n"
95 "nxy : %d\n"
96 "nvox : %d\n"
97
98 "P.plist : %p\n"
99 "P.used : %d\n"
100 "P.M : %d\n"
101
102 "extrema count : %d\n"
103
104 "data_type : %d\n"
105 "adn_type : %d\n"
106 "func_type : %d\n"
107
108 "outfile : %s\n"
109 "sval_style : %d\n"
110
111 "cutoff : %f\n"
112 "min_dist : %f\n"
113 "out_rad : %f\n"
114
115 "negatives : %d\n"
116 "ngbr_style : %d\n"
117 "overwrite : %d\n"
118 "quiet : %d\n"
119 "coords_only : %d\n"
120 "true_max : %d\n"
121 "dicom_coords : %d\n"
122 "debug : %d\n"
123 "------------------------------\n",
124
125 (void *)M->dset, (void *)M->sdata, (void *)M->result,
126 M->nx, M->ny, M->nz, M->nxy, M->nvox,
127 (void *)M->P.plist, M->P.used, M->P.M,
128 M->extrema_count,
129 M->data_type, M->adn_type, M->func_type,
130 M->outfile, M->sval_style,
131 M->cutoff, M->min_dist, M->out_rad,
132 M->negatives, M->ngbr_style, M->overwrite,
133 M->quiet, M->coords_only, M->true_max, M->dicom_coords, M->debug
134 );
135
136 EXRETURN;
137 }
138
139
140 /*----------------------------------------------------------------------
141 **
142 ** Remove mask points around lower intensity extrema.
143 **
144 **----------------------------------------------------------------------
145 */
146 static int
apply_min_dist(maxima_s * M)147 apply_min_dist( maxima_s * M )
148 {
149 int * iptr, count;
150 point_list_s newP = { NULL, 0, 0 };
151
152
153 for ( count = 0, iptr = M->P.plist; count < M->P.used; count++, iptr++ )
154 clear_around_point( *iptr, M, &newP );
155
156 free( M->P.plist ); /* replace old point list with new */
157 M->P = newP;
158
159 return 1;
160 }
161
162
163 /*----------------------------------------------------------------------
164 **
165 ** Simply clear to a radius of min_dist around the given point.
166 **
167 ** Clear all memory in the given radius and reset the point.
168 **
169 **----------------------------------------------------------------------
170 */
171 static int
clear_around_point(int p,maxima_s * M,point_list_s * newP)172 clear_around_point( int p, maxima_s * M, point_list_s * newP )
173 {
174 int X, Y, Z;
175 int xmin, xmax, ymin, ymax, zmin, zmax;
176 int yc, zc, xrad, yrad, yrad2;
177 int xbase, ybase, zbase;
178
179 short * optr;
180 float radius = M->min_dist;
181
182 static point_list_s P = { NULL, 0, 0 }; /* for allocation speed */
183 P.used = 0;
184
185
186 X = p % M->nx;
187 Y = (p % M->nxy) / M->nx;
188 Z = p / M->nxy;
189
190 if ( ! *(M->result + ( Z*M->ny + Y ) * M->nx + X ) )
191 return 1;
192
193 zmin = ( Z < radius ) ? Z : radius;
194 zmax = ( Z + radius >= M->nz ) ? ( M->nz-Z-1 ) : radius;
195
196 if ( M->debug > 1 )
197 fprintf(stderr,"+d index %d, val %f\n", p, M->sdata[p]*gr_fac);
198
199 for ( zc = -zmin; zc <= zmax; zc++ )
200 {
201 zbase = ( Z + zc ) * M->nx * M->ny;
202 yrad2 = radius * radius - zc * zc;
203 yrad = (int)sqrt( yrad2 );
204
205 ymin = ( Y < yrad ) ? Y : yrad;
206 ymax = ( Y + yrad >= M->ny ) ? ( M->ny - Y - 1 ) : yrad;
207
208 for ( yc = -ymin; yc <= ymax; yc++ )
209 {
210 ybase = ( Y + yc ) * M->nx;
211 xrad = (int)sqrt( yrad2 - yc * yc );
212
213 xmin = ( X < xrad ) ? X : xrad;
214 xmax = ( X + xrad >= M->nx ) ? ( M->nx - X - 1 ) : xrad;
215
216 optr = M->result + ybase + zbase;
217
218 for ( xbase = X-xmin; xbase <= X+xmax; xbase++ )
219 if ( *( optr + xbase ) )
220 {
221 *(optr + xbase) = 0;
222 /* maybe a switch later */
223 if ( M->ngbr_style == MAX_WEIGHTED_AVE_STYLE )
224 {
225 if ( ! add_point_to_list( &P, xbase+ybase+zbase ) )
226 return 0;
227 if ( M->debug > 2 )
228 fprintf(stderr," coords %d %d %d [%d], value %f\n",
229 xbase, Y+yc, Z+zc, xbase+ybase+zbase,
230 M->sdata[xbase+ybase+zbase]*gr_fac);
231 }
232 }
233 }
234 }
235
236 switch( M->ngbr_style )
237 {
238 case MAX_SORT_N_REMOVE_STYLE :
239 if ( ! add_point_to_list( newP, p ) )
240 return 0;
241
242 break;
243
244 case MAX_WEIGHTED_AVE_STYLE :
245 if ( ( p = weighted_index( &P, M ) ) < 0 )
246 return 0;
247
248 if ( ! add_point_to_list( newP, p ) )
249 return 0;
250
251 break;
252
253 default :
254 sprintf( grMessage, "Error: cap_00\nInvalid ngbr_style %d.",
255 M->ngbr_style );
256 rERROR( grMessage );
257 return 0;
258
259 break;
260 }
261
262 return 1;
263 }
264
265
266 /*----------------------------------------------------------------------
267 **
268 ** Given a list of points and the data in M, calculate the index
269 ** of the weighted average of those points. Values are from M->sdata.
270 **
271 **----------------------------------------------------------------------
272 */
273 static int
weighted_index(point_list_s * P,maxima_s * M)274 weighted_index( point_list_s * P, maxima_s * M )
275 {
276 double total_x = 0.0, total_y = 0.0, total_z = 0.0; /* weight*position */
277 double weight = 0.0;
278 double value;
279 int x, y, z;
280 int count, index;
281 int * iptr;
282
283 if ( ( P->plist == NULL ) || ( P->used <= 0 ) )
284 {
285 rERROR( "Error wi_00\nEmpty point list." );
286 return( -1 );
287 }
288
289 if ( P->used == 1 ) /* no weighting necessary */
290 return( P->plist[0] );
291
292 for ( count = 0, iptr = P->plist; count < P->used; count++, iptr++ )
293 {
294 index = *iptr;
295
296 x = index % M->nx;
297 y = ( index % M->nxy ) / M->nx;
298 z = index / M->nxy;
299
300 value = M->sdata[ index ];
301
302 weight += value;
303 total_x += value * x;
304 total_y += value * y;
305 total_z += value * z;
306 }
307
308 if ( M->debug > 1 )
309 fprintf(stderr, "-d nvals, weight, ave = %d, %f, %f\n",
310 P->used, weight*gr_fac, weight*gr_fac/P->used);
311
312 if ( weight <= 0.0 )
313 {
314 sprintf( grMessage, "Error: wi_10\nunexpected weight of %f", weight );
315 rERROR( grMessage );
316 }
317
318 x = ( int )( total_x / weight + 0.4 ); /* ~rounded average */
319 y = ( int )( total_y / weight + 0.4 );
320 z = ( int )( total_z / weight + 0.4 );
321
322 index = ( z * M->ny + y ) * M->nx + x;
323
324 if ( M->debug > 1 )
325 fprintf(stderr, "-d weighted i,j,k, ind, val = %d, %d, %d, %d, %f\n",
326 x, y, z, index, M->sdata[index]*gr_fac);
327
328 return index;
329 }
330
331
332 /*----------------------------------------------------------------------
333 **
334 ** Create a point list of all set extremas.
335 **
336 ** Walk through the result data, and for any set point, insert into list.
337 **
338 **----------------------------------------------------------------------
339 */
340 static int
create_point_list(maxima_s * M)341 create_point_list( maxima_s * M )
342 {
343 short * mptr;
344 int count;
345 point_list_s * P = &M->P;
346
347 ENTRY("create_pint_list");
348
349 mptr = M->result;
350 for ( count = 0; count < M->nvox; count++ )
351 if ( *mptr++ )
352 if ( ! add_point_to_list( P, count ) )
353 RETURN(0);
354
355 if ( M->debug > 0 )
356 show_point_list_s( "+d point list created: ", &M->P, M->debug );
357
358 RETURN(1);
359 }
360
361
362 /*----------------------------------------------------------------------
363 **
364 ** Add a point to the list.
365 **
366 **----------------------------------------------------------------------
367 */
368 static int
add_point_to_list(point_list_s * P,int offset)369 add_point_to_list( point_list_s * P, int offset )
370 {
371 ENTRY("add_point_to_list");
372 if ( ! P->plist )
373 {
374 P->M = 100;
375 P->used = 0;
376
377 if ( ( P->plist = (int *)malloc( P->M * sizeof(int) ) ) == NULL )
378 {
379 rERROR( "Error: aptl_10\n"
380 "Failed to allocate memory for initial point list.\n" );
381 RETURN(0);
382 }
383 }
384 else if ( P->used == P->M )
385 {
386 P->M = P->M + 100;
387
388 if ( ( P->plist = (int *)realloc( P->plist, P->M*sizeof(int))) == NULL )
389 {
390 sprintf( grMessage, "Error: aptl_20\n"
391 "Failed to reallocate %d ints for point list", P->M );
392 RETURN(0);
393 }
394 }
395
396 P->plist[P->used] = offset;
397 P->used++;
398
399 RETURN(1);
400 }
401
402
403 /*----------------------------------------------------------------------
404 **
405 ** Fill a sphere around any set point.
406 **
407 **----------------------------------------------------------------------
408 */
409 static int
apply_fill_radius(maxima_s * M)410 apply_fill_radius( maxima_s * M )
411 {
412 int count, outval;
413 int x, y, z;
414 int * iptr;
415
416 ENTRY("apply_fill_radius");
417
418 for ( count = 0, iptr = M->P.plist; count < M->P.used; count++, iptr++ )
419 {
420 outval = (M->sval_style == 0) ? MAX_MASK_FILL_VAL :
421 (M->sval_style == 1) ? (count+1) :
422 (M->P.used - count);
423
424 if ( M->out_rad < 1.0 )
425 {
426 M->result[ *iptr ] = outval;
427 continue;
428 }
429
430 x = *iptr % M->nx;
431 y = ( *iptr % M->nxy ) / M->nx;
432 z = *iptr / M->nxy;
433
434 radial_fill( x, y, z, M, outval );
435 }
436
437 RETURN(1);
438 }
439
440
441 /*----------------------------------------------------------------------
442 **
443 ** Fill a radius around the current point.
444 **
445 **----------------------------------------------------------------------
446 */
447 static int
radial_fill(int X,int Y,int Z,maxima_s * M,int val)448 radial_fill( int X, int Y, int Z, maxima_s * M, int val )
449 {
450 int xmin, xmax, ymin, ymax, zmin, zmax;
451 int yc, zc, xrad, yrad, yrad2;
452 int xbase, ybase, zbase;
453
454 short * sptr, * optr;
455 float radius = M->out_rad;
456
457 ENTRY("radial_fill");
458
459 zmin = ( Z < radius ) ? Z : radius;
460 zmax = ( Z + radius >= M->nz ) ? ( M->nz-Z-1 ) : radius;
461
462 for ( zc = -zmin; zc <= zmax; zc++ )
463 {
464 zbase = ( Z + zc ) * M->nx * M->ny;
465 yrad2 = radius * radius - zc * zc;
466 yrad = (int)sqrt( yrad2 );
467
468 ymin = ( Y < yrad ) ? Y : yrad;
469 ymax = ( Y + yrad >= M->ny ) ? ( M->ny - Y - 1 ) : yrad;
470
471 for ( yc = -ymin; yc <= ymax; yc++ )
472 {
473 ybase = ( Y + yc ) * M->nx;
474 xrad = (int)sqrt( yrad2 - yc * yc );
475
476 xmin = ( X < xrad ) ? X : xrad;
477 xmax = ( X + xrad >= M->nx ) ? ( M->nx - X - 1 ) : xrad;
478
479 optr = M->result + ybase + zbase;
480
481 for ( xbase = X-xmin; xbase <= X+xmax; xbase++ )
482 {
483 sptr = optr + xbase;
484
485 if ( ! *sptr )
486 *sptr = val;
487 }
488 }
489 }
490
491 RETURN(1);
492 }
493
494
495 /*----------------------------------------------------------------------
496 **
497 ** Find local maxima in short brick.
498 **
499 **----------------------------------------------------------------------
500 */
501 static int
find_local_maxima(maxima_s * M)502 find_local_maxima( maxima_s * M )
503 {
504 short * sourcep, * destp;
505 int c, cx, cy, cz;
506 int maxx = M->nx - 1;
507 int maxy = M->ny - 1;
508 int nx = M->nx, nxy = M->nx * M->ny;
509 int offset[ 26 ]; /* for speed */
510
511 ENTRY("find_local_maxima");
512
513 offset[ 0] = +1;
514 offset[ 1] = -1;
515 offset[ 2] = nx;
516 offset[ 3] = nx+1;
517 offset[ 4] = nx-1;
518 offset[ 5] = -nx;
519 offset[ 6] = -nx+1;
520 offset[ 7] = -nx-1;
521 offset[ 8] = nxy;
522 offset[ 9] = nxy+1;
523 offset[10] = nxy-1;
524 offset[11] = nxy+nx;
525 offset[12] = nxy+nx+1;
526 offset[13] = nxy+nx-1;
527 offset[14] = nxy-nx;
528 offset[15] = nxy-nx+1;
529 offset[16] = nxy-nx-1;
530 offset[17] = -nxy;
531 offset[18] = -nxy+1;
532 offset[19] = -nxy-1;
533 offset[20] = -nxy+nx;
534 offset[21] = -nxy+nx+1;
535 offset[22] = -nxy+nx-1;
536 offset[23] = -nxy-nx;
537 offset[24] = -nxy-nx+1;
538 offset[25] = -nxy-nx-1;
539
540 sourcep = M->sdata + nxy; /* skip first plane */
541 destp = M->result + nxy;
542
543 for ( cz = 0; cz < M->nz-2; cz++ ) /* and skip last plane (1->2) [rickr] */
544 {
545 for ( cy = 0; cy < M->ny; cy++ )
546 {
547 if ( ( cy == 0 ) || ( cy == maxy ) )
548 {
549 sourcep += nx;
550 destp += nx;
551
552 continue;
553 }
554
555 for ( cx = 0; cx < M->nx; cx++ )
556 {
557 if ( ( cx == 0 ) || ( cx == maxx ) )
558 {
559 sourcep++;
560 destp++;
561
562 continue;
563 }
564
565 if ( ! M->negatives && ( *sourcep < M->cutoff ) )
566 {
567 sourcep++;
568 destp++;
569
570 continue;
571 }
572
573 if ( M->negatives && ( *sourcep > M->cutoff ) )
574 {
575 sourcep++;
576 destp++;
577
578 continue;
579 }
580
581 *destp = MAX_MASK_FILL_VAL;
582 M->extrema_count++;
583
584 if ( M->true_max )
585 {
586 if ( ! M->negatives )
587 {
588 for ( c = 0; c < 26; c++ )
589 if ( *sourcep <= sourcep[offset[c]] )
590 {
591 *destp = 0;
592 M->extrema_count--;
593
594 break;
595 }
596 }
597 else
598 {
599 for ( c = 0; c < 26; c++ )
600 if ( *sourcep >= sourcep[offset[c]] )
601 {
602 *destp = 0;
603 M->extrema_count--;
604
605 break;
606 }
607 }
608 }
609 else
610 {
611 if ( ! M->negatives )
612 {
613 for ( c = 0; c < 26; c++ )
614 if ( *sourcep < sourcep[offset[c]] )
615 {
616 *destp = 0;
617 M->extrema_count--;
618
619 break;
620 }
621 }
622 else
623 {
624 for ( c = 0; c < 26; c++ )
625 if ( *sourcep > sourcep[offset[c]] )
626 {
627 *destp = 0;
628 M->extrema_count--;
629
630 break;
631 }
632 }
633 }
634
635 sourcep++;
636 destp++;
637 }
638 }
639 }
640
641 if ( M->debug > 3 ) show_maxima_s( "post find local maxima: ", M );
642
643 RETURN(1);
644 }
645
646
647 /*----------------------------------------------------------------------
648 **
649 ** Display extrema coordinates and write output BRIK.
650 **
651 ** Put any new BRIK into memory.
652 **
653 **----------------------------------------------------------------------
654 */
set_results(r_afni_s * A,maxima_s * M,THD_3dim_dataset * dset)655 int set_results( r_afni_s * A, maxima_s * M, THD_3dim_dataset * dset )
656 {
657 ENTRY("write_results");
658
659 if ( ! dset ) RETURN(1);
660
661 EDIT_dset_items( dset,
662 ADN_prefix, M->outfile,
663 ADN_label1, M->outfile,
664 ADN_nvals, 1,
665 ADN_ntt, 0,
666 ADN_type, HEAD_FUNC_TYPE,
667 ADN_func_type, FUNC_FIM_TYPE,
668 ADN_none
669 );
670
671 EDIT_substitute_brick( dset, 0, M->data_type, M->result );
672 EDIT_BRICK_FACTOR ( dset, 0, 0.0 );
673
674 RETURN(1);
675 }
676
677
678 /*----------------------------------------------------------------------
679 * print a list of strings
680 **----------------------------------------------------------------------
681 */
disp_str_list(char * list[],int len)682 int disp_str_list( char * list[], int len )
683 {
684 int c;
685 for( c = 0; c < len; c++ ) puts(list[c]);
686 fflush(stdout);
687 return 0;
688 }
689
690
691 /*----------------------------------------------------------------------
692 **
693 ** For every point in list, display the talairach coordinates.
694 **
695 **----------------------------------------------------------------------
696 */
display_coords(r_afni_s * A,maxima_s * M)697 int display_coords( r_afni_s * A, maxima_s * M )
698 {
699 THD_fvec3 f3, t3;
700 THD_ivec3 i3;
701 float prod, factor = A->factor[0];
702 short * optr;
703 short * mptr;
704 int * iptr;
705 int X, Y, Z, count;
706
707 THD_coorder cord; /* 16 Apr 2013 */
708 char orcode[4];
709
710 ENTRY("display_coords");
711
712 /* For dset based coordinates (including signs), use THD_coorder struct
713 and THD_dicom_to_coorder(). Previously (i.e. for past 15 years), the
714 dataset coord order just permuted the axes, but did not apply signs.
715 This issue was reported by G Pagnoni. 16 Apr 2013 [rickr]*/
716 orcode[0] = ORIENT_first[M->dset->daxes->xxorient];
717 orcode[1] = ORIENT_first[M->dset->daxes->yyorient];
718 orcode[2] = ORIENT_first[M->dset->daxes->zzorient];
719 orcode[3] = '\0';
720 THD_coorder_fill(orcode, &cord);
721
722 point_list_s * P = &M->P;
723
724 if( !M->coords_only ) /* 18 Aug 2006 [rickr] */
725 {
726 printf( "---------------------------------------------\n" );
727 if ( M->dicom_coords )
728 printf( "RAI mm coordinates:\n\n" );
729 else
730 printf( "%3s mm coordinates:\n\n", orcode);
731 }
732
733 for ( count = 0, iptr = P->plist; count < P->used; count++, iptr++ )
734 {
735 X = *iptr % M->nx;
736 Y = (*iptr % M->nxy) / M->nx;
737 Z = *iptr / M->nxy;
738 i3.ijk[0] = X; i3.ijk[1] = Y; i3.ijk[2] = Z;
739 f3 = THD_3dind_to_dicomm_no_wod(M->dset, i3);
740 if ( ! M->dicom_coords ) /* first dicom, then invert 16 Apr 2013 */
741 THD_dicom_to_coorder(&cord, f3.xyz, f3.xyz+1, f3.xyz+2);
742
743 optr = M->sdata + *iptr;
744 mptr = M->result + *iptr;
745
746 if ( M->coords_only )
747 printf( "%7.2f %7.2f %7.2f\n", f3.xyz[0], f3.xyz[1], f3.xyz[2] );
748 else
749 {
750 prod = *optr * factor;
751 printf( "(%7.2f %7.2f %7.2f) : val = %f\n",
752 f3.xyz[0], f3.xyz[1], f3.xyz[2], prod );
753 }
754 }
755
756 if( !M->coords_only )
757 {
758 if ( P->used )
759 printf( "\nnumber of extrema = %d\n", P->used );
760 else
761 printf( "No extrema found.\n" );
762 printf( "---------------------------------------------\n" );
763 }
764
765 RETURN(1);
766 }
767
768
769 /*----------------------------------------------------------------------
770 **
771 ** Initialize the afni structure.
772 **
773 ** These are the fields that need to be set before reading a dataset.
774 **
775 **----------------------------------------------------------------------
776 */
init_afni_s(r_afni_s * A)777 int init_afni_s( r_afni_s * A )
778 {
779 ENTRY("init_afni_s");
780
781 memset( A, 0, sizeof( r_afni_s ) );
782
783 A->must_be_short = 1;
784 A->want_floats = 1;
785 A->subs_must_equal = 1;
786 A->max_subs = 1;
787 A->sub_brick = 0;
788
789 RETURN(1);
790 }
791
792
793 /*----------------------------------------------------------------------
794 **
795 ** Initialize the maxima structure.
796 **
797 **----------------------------------------------------------------------
798 */
init_maxima_s(maxima_s * M,r_afni_s * A,char * outprefix)799 int init_maxima_s( maxima_s * M, r_afni_s * A, char * outprefix )
800 {
801 ENTRY("init_maxima_s");
802
803 M->dset = A->dset[0];
804
805 M->sdata = A->simage[0];
806
807 if ( ( M->result = (short *)calloc( A->nvox, sizeof( short ) ) ) == NULL )
808 {
809 sprintf( grMessage, "Error: ims_05\n"
810 "Failed to allocate M for %d shorts.", A->nvox );
811 rERROR( grMessage );
812 RETURN(0);
813 }
814
815 M->nx = A->nx;
816 M->ny = A->ny;
817 M->nz = A->nz;
818 M->nxy = A->nx * A->ny;
819 M->nvox = A->nvox;
820
821 M->P.plist = NULL;
822 M->P.used = 0;
823 M->P.M = 0;
824
825 M->extrema_count = 0;
826
827 M->data_type = MRI_short; /* output will be short */
828 M->adn_type = HEAD_FUNC_TYPE;
829 M->func_type = FUNC_FIM_TYPE;
830
831 if ( outprefix && strlen( outprefix ) > R_FILE_L )
832 {
833 sprintf( grMessage, "Error: ims_10\n"
834 "Outfile prefix exceeds %d characters.", R_FILE_L );
835 rERROR( grMessage );
836 RETURN(0);
837 }
838
839 if ( outprefix )
840 strcpy( M->outfile, outprefix );
841 else
842 *M->outfile = 0;
843
844 M->cutoff = 0.0;
845 M->min_dist = 0.0;
846 M->out_rad = 0.0;
847
848 M->negatives = 0;
849 M->ngbr_style = MAX_SORT_N_REMOVE_STYLE;
850 M->overwrite = 0;
851 M->quiet = 0;
852 M->coords_only = 0;
853 M->true_max = 0;
854 M->debug = 0;
855
856 RETURN(1);
857 }
858
859
860 /*----------------------------------------------------------------------
861 **
862 ** Comparasin function of positives for qsort.
863 **
864 **----------------------------------------------------------------------
865 */
866 int
point_comp_pos(const void * p1,const void * p2)867 point_comp_pos( const void * p1, const void * p2 )
868 {
869 short v1 = *( gr_orig_data + *(int *)p1 );
870 short v2 = *( gr_orig_data + *(int *)p2 );
871
872 if ( v1 < v2 ) /* reverse order to get large numbers first */
873 return 1;
874 else if ( v1 > v2 )
875 return -1;
876
877 return 0;
878 }
879
880
881 /*----------------------------------------------------------------------
882 **
883 ** Comparasin function of negatives for qsort.
884 **
885 **----------------------------------------------------------------------
886 */
887 int
point_comp_neg(const void * p1,const void * p2)888 point_comp_neg( const void * p1, const void * p2 )
889 {
890 short v1 = *( gr_orig_data + *(int *)p1 );
891 short v2 = *( gr_orig_data + *(int *)p2 );
892
893 /* more negative is greater AND reverse the order */
894 if ( v1 > v2 )
895 return 1;
896 else if ( v1 < v2 )
897 return -1;
898
899
900 return 0;
901 }
902
903
904 /*----------------------------------------------------------------------
905 **
906 ** Fill the afni structure from an existing dataset.
907 **
908 **----------------------------------------------------------------------
909 */
910 int
r_set_afni_s_from_dset(r_afni_s * A,THD_3dim_dataset * dset,int debug)911 r_set_afni_s_from_dset( r_afni_s * A, THD_3dim_dataset * dset, int debug )
912 {
913 ENTRY("r_set_afni_s_from_dset");
914
915 if( debug > 3 ) disp_r_afni_s( "-d initial struct", A);
916
917 if ( A->num_dsets >= R_MAX_AFNI_DSETS )
918 {
919 sprintf( grMessage, "Error: rsasfd_00\n"
920 "We only have memory to hold %d datasets. exiting...\n",
921 R_MAX_AFNI_DSETS );
922 rERROR( grMessage );
923
924 RETURN(0);
925 }
926
927 A->dset[ 0 ] = dset; /* rickr - use sub-brick */
928 A->simage[ 0 ] = ( short * )DSET_ARRAY( dset, A->sub_brick );
929
930 if ( !A->simage[0] )
931 {
932 sprintf(grMessage,
933 "** data not available, is this in warp-on-demand mode?\n");
934 rERROR(grMessage);
935 RETURN(0);
936 }
937
938 if ((A->factor[0] = DSET_BRICK_FACTOR(dset, A->sub_brick)) == 0.0 )
939 A->factor[0] = 1.0;
940
941 A->subs [ 0 ] = DSET_NVALS( dset );
942
943 A->nx = dset->daxes->nxx;
944 A->ny = dset->daxes->nyy;
945 A->nz = dset->daxes->nzz;
946 A->nvox = A->nx * A->ny * A->nz;
947
948 if ( A->want_floats )
949 {
950 int count;
951 short * sptr;
952 float * fptr;
953 float factor = A->factor[ 0 ]; /* just for speed */
954
955 if ( ( A->fimage[ 0 ] =
956 ( float * )calloc( A->nvox, sizeof( float ) ) ) == NULL )
957 {
958 sprintf( grMessage, "Error: rsasfd_10\n"
959 "Failed to allocate memory for %d floats.\n", A->nvox );
960 rERROR( grMessage );
961
962 RETURN(0);
963 }
964
965 fptr = A->fimage[ 0 ];
966 sptr = A->simage[ 0 ];
967 for ( count = 0; count < A->nvox; count++ )
968 *fptr++ = *sptr++ * factor;
969 }
970
971 A->max_u_short = r_get_max_u_short( (u_short *)A->simage[0], A->nvox );
972
973 /* A->num_dsets++; not using more than one */
974
975 if( debug > 1 ) disp_r_afni_s( "-d initial struct", A);
976
977 RETURN(1);
978 }
979
980
981 /*----------------------------------------------------------------------
982 **
983 ** Compute the maximum unsigned short in an array.
984 **
985 **----------------------------------------------------------------------
986 */
987 u_short
r_get_max_u_short(u_short * S,int size)988 r_get_max_u_short( u_short * S, int size )
989 {
990 u_short * usptr, max = *S;
991 int c = 0;
992
993 for ( c = 0, usptr = S; c < size; c++, usptr++ )
994 {
995 if ( *usptr > max )
996 max = *usptr;
997 }
998
999 return max;
1000 }
1001
1002
1003 /*----------------------------------------------------------------------
1004 **
1005 ** Free all local memory.
1006 **
1007 **----------------------------------------------------------------------
1008 */
free_memory(r_afni_s * A,maxima_s * M)1009 void free_memory( r_afni_s * A, maxima_s * M )
1010 {
1011 ENTRY("free_memory");
1012 if ( A->want_floats && A->fimage[0] )
1013 free( A->fimage[0] );
1014
1015 if ( M->result && !M->outfile[0] )
1016 free( M->result );
1017
1018 if ( M->P.plist )
1019 free( M->P.plist );
1020
1021 EXRETURN;
1022 }
1023
1024
1025 /*----------------------------------------------------------------------
1026 * display contents of r_afni_s struct
1027 *----------------------------------------------------------------------
1028 */
disp_r_afni_s(char * mesg,r_afni_s * A)1029 int disp_r_afni_s( char * mesg, r_afni_s * A )
1030 {
1031 ENTRY("disp_r_afni_s");
1032
1033 if( mesg ) puts(mesg);
1034
1035 printf("-d r_afni_s @ %p\n"
1036 " must_be_short, want_floats = %d, %d\n"
1037 " subs_must_equal, max_subs = %d, %d\n"
1038 " dset, simage = %p, %p\n"
1039 " factor = %f\n"
1040 " subs, sub_brick = %d, %d\n"
1041 " nx, ny, nz, nvox = %d, %d, %d, %d\n"
1042 " fimage = %p\n"
1043 " max_u_short, num_dsets = %d, %d\n",
1044 A,
1045 A->must_be_short, A->want_floats,
1046 A->subs_must_equal, A->max_subs,
1047 (void *)A->dset[0], (void *)A->simage[0],
1048 A->factor[0], A->subs[0], A->sub_brick,
1049 A->nx, A->ny, A->nz, A->nvox, (void *)A->fimage,
1050 A->max_u_short, A->num_dsets);
1051
1052 return 0;
1053 }
1054
1055