1 #include "mrilib.h"
2 #include "extrema.h"
3
4 #define ALLOW_INEDGE
5
6 #ifdef ALLOW_INEDGE
7 #include "mri_intedge.c" /* Jul 2018 */
8 #endif
9
main(int argc,char * argv[])10 int main( int argc , char * argv[] )
11 {
12 THD_3dim_dataset * inset , * outset ;
13 int nx,ny,nz,nxyz,nval , ii,kk , nopt=1, nsum=0 ;
14 char * prefix = "edge3" , *insetname=NULL;
15 int datum=-1 , verb=0 , do_sd=0, do_sum=0 , do_sqr=0, firstds=0 ;
16 float ** sum , fsum, ffac = 0.0;
17 float ** sd;
18 int border[3]={0,0,0};
19 int indims[3]={0,0,0};
20 int fscale=0 , gscale=0 , nscale=0 ;
21 float filterCoefs[3] = {1.0, 1.0, 1.0};
22 recursiveFilterType filterType = ALPHA_DERICHE;
23 /* recursiveFilterType filterType = GAUSSIAN_DERICHE; */
24
25 int do_inedge=0 ; /* Jul 2018 */
26
27 /*-- help? --*/
28
29 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
30 printf(
31 "Usage: 3dedge3 [options] dset dset ...\n"
32 "Does 3D Edge detection using the library 3DEdge by;\n"
33 "by Gregoire Malandain (gregoire.malandain@sophia.inria.fr)\n"
34 "\n"
35 "Options :\n"
36 " -input iii = Input dataset\n"
37 " -verbose = Print out some information along the way.\n"
38 " -prefix ppp = Sets the prefix of the output dataset.\n"
39 " -datum ddd = Sets the datum of the output dataset.\n"
40 " -fscale = Force scaling of the output to the maximum integer range.\n"
41 " -gscale = Same as '-fscale', but also forces each output sub-brick to\n"
42 " to get the same scaling factor.\n"
43 " -nscale = Don't do any scaling on output to byte or short datasets.\n"
44 " -scale_floats VAL = Multiply input by VAL, but only if the input datum is\n"
45 " float. This is needed when the input dataset\n"
46 " has a small range, like 0 to 2.0 for instance.\n"
47 " With such a range, very few edges are detected due to\n"
48 " what I suspect to be truncation problems.\n"
49 " Multiplying such a dataset by 10000 fixes the problem\n"
50 " and the scaling is undone at the output.\n"
51 #if 0
52 "\n"
53 " -inedge = Process the #0 sub-brick like the '-inedge' option of\n"
54 " 3dQwarp and write the result out. Mostly for checking\n"
55 " things out, not for any other practical purpose. [Jul 2018]\n"
56 "\n"
57 #endif
58 "\n"
59 "References for the algorithms:\n"
60 " - Optimal edge detection using recursive filtering\n"
61 " R. Deriche, International Journal of Computer Vision,\n"
62 " pp 167-187, 1987.\n"
63 " - Recursive filtering and edge tracking: two primary tools\n"
64 " for 3-D edge detection, O. Monga, R. Deriche, G. Malandain\n"
65 " and J.-P. Cocquerez, Image and Vision Computing 4:9, \n"
66 " pp 203-214, August 1991.\n"
67 "\n"
68 ) ;
69 PRINT_COMPILE_DATE ; exit(0) ;
70 }
71
72 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
73
74 mainENTRY("3dedge3 main"); machdep() ; PRINT_VERSION("3dedge3") ;
75
76 { int new_argc ; char ** new_argv ;
77 addto_args( argc , argv , &new_argc , &new_argv ) ;
78 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
79 }
80
81 AFNI_logger("3dedge3",argc,argv) ;
82
83 /*-- command line options --*/
84
85 while( nopt < argc ){
86
87
88 if( strcmp(argv[nopt],"-prefix") == 0 ){
89 if( ++nopt >= argc ){
90 fprintf(stderr,
91 "** ERROR: need an argument after -prefix!\n"); exit(1);
92 }
93 prefix = argv[nopt] ;
94 nopt++ ; continue ;
95 }
96 if( strcmp(argv[nopt],"-input") == 0 ){
97 if( ++nopt >= argc ){
98 fprintf(stderr,"** ERROR: need an argument after -input!\n");
99 exit(1);
100 }
101 insetname = argv[nopt] ;
102 nopt++ ; continue ;
103 }
104
105 if( strcmp(argv[nopt],"-inedge") == 0 ){ /* Jul 2018 */
106 #ifndef ALLOW_INEDGE
107 WARNING_message("-inedge is disabled - ignoring this option :(") ;
108 #endif
109 do_inedge = 1 ; nopt++ ; continue ;
110 }
111
112 if( strcmp(argv[nopt],"-scale_floats") == 0 ){
113 if( ++nopt >= argc ){
114 fprintf(stderr,"** ERROR: a factor after -scale_floats!\n");
115 exit(1);
116 }
117 ffac = (float)strtod(argv[nopt],NULL) ;
118 nopt++ ; continue ;
119 }
120
121 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
122 verb++ ; nopt++ ; continue ;
123 }
124
125 if( strcmp(argv[nopt],"-datum") == 0 ){
126 if( ++nopt >= argc ){
127 fprintf(stderr,"** ERROR: need an argument after -datum!\n");
128 exit(1);
129 }
130 if( strcmp(argv[nopt],"short") == 0 ){
131 datum = MRI_short ;
132 } else if( strcmp(argv[nopt],"float") == 0 ){
133 datum = MRI_float ;
134 } else if( strcmp(argv[nopt],"byte") == 0 ){
135 datum = MRI_byte ;
136 } else {
137 fprintf(stderr,
138 "** ERROR -datum of type '%s' not supported in 3dedge3!\n",
139 argv[nopt] ) ;
140 exit(1) ;
141 }
142 nopt++ ; continue ; /* go to next arg */
143 }
144
145 if( strncmp(argv[nopt],"-nscale",6) == 0 ){
146 gscale = fscale = 0 ; nscale = 1 ;
147 nopt++ ; continue ;
148 }
149
150 if( strncmp(argv[nopt],"-fscale",6) == 0 ){
151 fscale = 1 ; nscale = 0 ;
152 nopt++ ; continue ;
153 }
154
155 if( strncmp(argv[nopt],"-gscale",6) == 0 ){
156 gscale = fscale = 1 ; nscale = 0 ;
157 nopt++ ; continue ;
158 }
159
160 fprintf(stderr,"** ERROR: unknown option %s\n",argv[nopt]) ;
161 exit(1) ;
162 }
163
164 if (!insetname) {
165 fprintf(stderr,"** ERROR: no input dset specified\n") ;
166 exit(1) ;
167 }
168
169 {
170
171 /*-- input dataset header --*/
172
173 inset = THD_open_dataset( insetname ) ;
174 if( !ISVALID_DSET(inset) ){
175 fprintf(stderr,"** ERROR: can't open dataset %s\n",argv[nopt]) ;
176 exit(1) ;
177 }
178
179 /*-- make workspace and empty output dataset --*/
180
181 if( nsum == 0 ){
182
183 nx = DSET_NX(inset) ;
184 ny = DSET_NY(inset) ;
185 nz = DSET_NZ(inset) ; nxyz= nx*ny*nz;
186 nval = DSET_NVALS(inset) ;
187
188 sum = (float **) malloc( sizeof(float *)*nval ) ;
189 /* array of sub-bricks */
190 for( kk=0 ; kk < nval ; kk++ ){
191 sum[kk] = (float *) malloc(sizeof(float)*nxyz) ;
192 /* kk-th sub-brick */
193 }
194
195 outset = EDIT_empty_copy( inset ) ;
196
197 if( datum < 0 ) datum = DSET_BRICK_TYPE(inset,0) ;
198
199 tross_Copy_History( inset , outset ) ;
200 tross_Make_History( "3dedge3" , argc,argv , outset ) ;
201
202 EDIT_dset_items( outset ,
203 ADN_prefix , prefix ,
204 ADN_datum_all , datum ,
205 ADN_none ) ;
206
207 if( THD_deathcon() && THD_is_file(outset->dblk->diskptr->header_name) ){
208 fprintf(stderr,
209 "*** Output file %s already exists -- cannot continue!\n",
210 outset->dblk->diskptr->header_name ) ;
211 exit(1) ;
212 }
213
214 }
215
216
217 /*-- read data from disk --*/
218
219 DSET_mallocize(inset);
220 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
221
222 if( verb ) fprintf(stderr," ++ read in dataset %s\n",insetname) ;
223
224 /*-- Do the -inedge option now [Jul 2018] --*/
225
226
227 #ifdef ALLOW_INEDGE
228 if( do_inedge ){
229 MRI_IMAGE *bim = THD_extract_float_brick(0,inset) ;
230 if( bim == NULL ) ERROR_exit("-inedge: Cannot get sub-brick #0") ;
231 mri_interior_edgeize(bim,4,0.222f) ;
232
233 EDIT_dset_items( outset, ADN_nvals,1 , ADN_datum_all,MRI_float , ADN_none ) ;
234 EDIT_substitute_brick(outset, 0, MRI_float, MRI_FLOAT_PTR(bim) ) ;
235 DSET_write(outset) ;
236 INFO_message("wrote out -inedge dataset: %s",DSET_BRIKNAME(outset)) ;
237 exit(0) ;
238 }
239 #endif
240
241 /*-- Edge detect each sub-brik --*/
242
243 indims[0] = DSET_NX(inset);
244 indims[1] = DSET_NY(inset);
245 indims[2] = DSET_NZ(inset);
246 border[0] = 50;
247 border[1] = 50;
248 border[2] = 50;
249 for( kk=0 ; kk < nval ; kk++ ){
250
251 if( verb )
252 fprintf(stderr," + sub-brick %d [%s]\n",
253 kk,MRI_TYPE_name[DSET_BRICK_TYPE(inset,kk)] ) ;
254 if (DSET_BRICK_TYPE(inset,kk) != DSET_BRICK_TYPE(inset,0)) {
255 fprintf(stderr,"ERROR: Sub-bricks of different types.\n"
256 "This is not splenda\n") ;
257 exit(1) ;
258 }
259 switch( DSET_BRICK_TYPE(inset,kk) ){
260 default:
261 fprintf(stderr,"ERROR: illegal input sub-brick datum\n") ;
262 exit(1) ;
263
264 case MRI_float:{
265 float *pp = (float *) DSET_ARRAY(inset,kk) ;
266 float fac = DSET_BRICK_FACTOR(inset,kk) ;
267 if(ffac) {
268 if (fac) fac = fac*ffac;
269 else fac = ffac;
270 }
271 if( fac ) {
272 for( ii=0 ; ii < nxyz ; ii++ ) { pp[ii] *= fac; }
273 }
274 if ( Extract_Gradient_Maxima_3D( (void *)pp, FLOAT,
275 sum[kk], FLOAT,
276 indims,
277 border,
278 filterCoefs,
279 filterType ) == 0 ) {
280 fprintf( stderr, "ERROR: gradient extraction failed.\n" );
281 exit( 1 );
282 }
283
284 }
285 break ;
286
287 case MRI_short:{
288 short *pp = (short *) DSET_ARRAY(inset,kk) ;
289 float fac = DSET_BRICK_FACTOR(inset,kk) ;
290
291 if( fac && verb) {
292 INFO_message(
293 "Ignoring brick factor of %f for Gradient Extraction",
294 fac);
295 }
296 if (verb) INFO_message("Going to Extract_Gradient_Maxima_3D\n");
297 if ( Extract_Gradient_Maxima_3D( (void *)pp, SSHORT,
298 sum[kk], FLOAT,
299 indims,
300 border,
301 filterCoefs,
302 filterType ) == 0 ) {
303 fprintf( stderr, "ERROR: gradient extraction failed.\n" );
304 exit( 1 );
305 }
306 }
307 break ;
308
309 case MRI_byte:{
310 byte *pp = (byte *) DSET_ARRAY(inset,kk) ;
311 float fac = DSET_BRICK_FACTOR(inset,kk) ;
312
313 if( fac && verb) {
314 INFO_message(
315 "Ignoring brick factor of %f for Gradient Extraction",
316 fac);
317 }
318 if ( Extract_Gradient_Maxima_3D( (void *)pp, UCHAR,
319 sum[kk], FLOAT,
320 indims,
321 border,
322 filterCoefs,
323 filterType ) == 0 ) {
324 fprintf( stderr, "ERROR: gradient extraction failed.\n" );
325 exit( 1 );
326 }
327 }
328 break ;
329 }
330 }
331
332 DSET_delete(inset) ;
333 }
334
335
336
337 switch( datum ){
338
339 default:
340 fprintf(stderr,
341 "*** Fatal Error ***\n"
342 "*** Somehow ended up with datum = %d\n",datum) ;
343 exit(1) ;
344
345 case MRI_float:{
346 for( kk=0 ; kk < nval ; kk++ ){
347 EDIT_substitute_brick(outset, kk, MRI_float, sum[kk]);
348 /* Undo ffac scaling.*/
349 if (ffac != 0.0) {
350 DSET_BRICK_FACTOR(outset, kk) = 1.0/ffac;
351 } else {
352 DSET_BRICK_FACTOR(outset, kk) = 0.0;
353 }
354 }
355 }
356 break ;
357
358 case MRI_byte:
359 case MRI_short:{
360 void ** dfim ;
361 float gtop=0.0 , fimfac , gtemp ;
362
363 if( verb )
364 fprintf(stderr," ++ Scaling output to type %s brick(s)\n",
365 MRI_TYPE_name[datum] ) ;
366
367 dfim = (void **) malloc(sizeof(void *)*nval) ;
368
369 if( gscale ){ /* allow global scaling */
370 gtop = 0.0 ;
371 for( kk=0 ; kk < nval ; kk++ ){
372 gtemp = MCW_vol_amax( nxyz , 1 , 1 , MRI_float, sum[kk] ) ;
373 gtop = MAX( gtop , gtemp ) ;
374 if( gtemp == 0.0 )
375 fprintf(stderr," -- Warning: output sub-brick %d is all zeros!\n",kk) ;
376 }
377 }
378
379 for (kk = 0 ; kk < nval ; kk ++ ) {
380
381 if( ! gscale ){
382 gtop = MCW_vol_amax( nxyz , 1 , 1 , MRI_float, sum[kk] ) ;
383 if( gtop == 0.0 )
384 fprintf(stderr," -- Warning: output sub-brick %d is all zeros!\n",kk) ;
385 }
386
387 if( fscale ){
388 fimfac = (gtop > 0.0) ? MRI_TYPE_maxval[datum] / gtop : 0.0 ;
389 } else if( !nscale ){
390 fimfac = (gtop > MRI_TYPE_maxval[datum] || (gtop > 0.0 && gtop <= 1.0) )
391 ? MRI_TYPE_maxval[datum]/ gtop : 0.0 ;
392 } else {
393 fimfac = 0.0 ;
394 }
395
396 if( verb ){
397 if( fimfac != 0.0 )
398 fprintf(stderr," ++ Sub-brick %d scale factor = %f\n",kk,fimfac) ;
399 else
400 fprintf(stderr," ++ Sub-brick %d: no scale factor\n" ,kk) ;
401 }
402
403 dfim[kk] = (void *) malloc( mri_datum_size(datum) * nxyz ) ;
404 if( dfim[kk] == NULL ){ fprintf(stderr,"*** malloc fails at output\n");exit(1); }
405
406 EDIT_coerce_scale_type( nxyz , fimfac ,
407 MRI_float, sum[kk] , datum,dfim[kk] ) ;
408 free( sum[kk] ) ;
409 EDIT_substitute_brick(outset, kk, datum, dfim[kk] );
410
411 DSET_BRICK_FACTOR(outset,kk) = (fimfac != 0.0) ? 1.0/fimfac : 0.0 ;
412 }
413 }
414 break ;
415 }
416
417 if( verb ) fprintf(stderr," ++ Computing output statistics\n") ;
418 THD_load_statistics( outset ) ;
419
420 THD_write_3dim_dataset( NULL,NULL , outset , True ) ;
421 if( verb ) fprintf(stderr," ++ Wrote output: %s\n",DSET_BRIKNAME(outset)) ;
422
423 exit(0) ;
424 }
425