#include "mrilib.h"
#include "extrema.h"
int main( int argc , char * argv[] )
{
THD_3dim_dataset * inset , * outset ;
int nx,ny,nz,nxyz,nval , ii,kk , nopt=1, nsum=0 ;
char * prefix = "edge3" , *insetname=NULL;
int datum=-1 , verb=0 , do_sd=0, do_sum=0 , do_sqr=0, firstds=0 ;
float ** sum , fsum;
float ** sd;
int border[3]={0,0,0};
int indims[3]={0,0,0};
int fscale=0 , gscale=0 , nscale=0 ;
float filterCoefs[3] = {1.0, 1.0, 1.0};
recursiveFilterType filterType = ALPHA_DERICHE;
/* recursiveFilterType filterType = GAUSSIAN_DERICHE; */
/*-- help? --*/
if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
printf("Usage: 3dedge3 [options] dset dset ...\n"
"Does 3D Edge detection using the library 3DEdge by;\n"
"by Gregoire Malandain (gregoire.malandain@sophia.inria.fr)\n"
"\n"
"Options :\n"
" -input iii = Input dataset\n"
" -verbose = Print out some information along the way.\n"
" -prefix ppp = Sets the prefix of the output dataset.\n"
" -datum ddd = Sets the datum of the output dataset.\n"
" -fscale = Force scaling of the output to the maximum integer range.\n"
" -gscale = Same as '-fscale', but also forces each output sub-brick to\n"
" to get the same scaling factor.\n"
" -nscale = Don't do any scaling on output to byte or short datasets.\n"
"\n"
"\n"
"References for the algorithms:\n"
" - Optimal edge detection using recursive filtering\n"
" R. Deriche, International Journal of Computer Vision,\n"
" pp 167-187, 1987.\n"
" - Recursive filtering and edge tracking: two primary tools\n"
" for 3-D edge detection, O. Monga, R. Deriche, G. Malandain\n"
" and J.-P. Cocquerez, Image and Vision Computing 4:9, \n"
" pp 203-214, August 1991.\n"
"\n"
) ;
exit(0) ;
}
/*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
mainENTRY("3dedge3 main"); machdep() ; PRINT_VERSION("3dedge3") ;
{ int new_argc ; char ** new_argv ;
addto_args( argc , argv , &new_argc , &new_argv ) ;
if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
}
AFNI_logger("3dedge3",argc,argv) ;
/*-- command line options --*/
while( nopt < argc ){
if( strcmp(argv[nopt],"-prefix") == 0 ){
if( ++nopt >= argc ){
fprintf(stderr,"** ERROR: need an argument after -prefix!\n"); exit(1);
}
prefix = argv[nopt] ;
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-input") == 0 ){
if( ++nopt >= argc ){
fprintf(stderr,"** ERROR: need an argument after -input!\n"); exit(1);
}
insetname = argv[nopt] ;
nopt++ ; continue ;
}
if( strncmp(argv[nopt],"-verbose",5) == 0 ){
verb++ ; nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-datum") == 0 ){
if( ++nopt >= argc ){
fprintf(stderr,"** ERROR: need an argument after -datum!\n"); exit(1);
}
if( strcmp(argv[nopt],"short") == 0 ){
datum = MRI_short ;
} else if( strcmp(argv[nopt],"float") == 0 ){
datum = MRI_float ;
} else if( strcmp(argv[nopt],"byte") == 0 ){
datum = MRI_byte ;
} else {
fprintf(stderr,"** ERROR -datum of type '%s' not supported in 3dedge3!\n",
argv[nopt] ) ;
exit(1) ;
}
nopt++ ; continue ; /* go to next arg */
}
if( strncmp(argv[nopt],"-nscale",6) == 0 ){
gscale = fscale = 0 ; nscale = 1 ;
nopt++ ; continue ;
}
if( strncmp(argv[nopt],"-fscale",6) == 0 ){
fscale = 1 ; nscale = 0 ;
nopt++ ; continue ;
}
if( strncmp(argv[nopt],"-gscale",6) == 0 ){
gscale = fscale = 1 ; nscale = 0 ;
nopt++ ; continue ;
}
fprintf(stderr,"** ERROR: unknown option %s\n",argv[nopt]) ;
exit(1) ;
}
if (!insetname) {
fprintf(stderr,"** ERROR: no input dset specified\n") ;
exit(1) ;
}
{
/*-- input dataset header --*/
inset = THD_open_dataset( insetname ) ;
if( !ISVALID_DSET(inset) ){
fprintf(stderr,"** ERROR: can't open dataset %s\n",argv[nopt]) ;
exit(1) ;
}
/*-- make workspace and empty output dataset --*/
if( nsum == 0 ){
nx = DSET_NX(inset) ;
ny = DSET_NY(inset) ;
nz = DSET_NZ(inset) ; nxyz= nx*ny*nz;
nval = DSET_NVALS(inset) ;
sum = (float **) malloc( sizeof(float *)*nval ) ; /* array of sub-bricks */
for( kk=0 ; kk < nval ; kk++ ){
sum[kk] = (float *) malloc(sizeof(float)*nxyz) ; /* kk-th sub-brick */
}
outset = EDIT_empty_copy( inset ) ;
if( datum < 0 ) datum = DSET_BRICK_TYPE(inset,0) ;
tross_Copy_History( inset , outset ) ;
tross_Make_History( "3dedge3" , argc,argv , outset ) ;
EDIT_dset_items( outset ,
ADN_prefix , prefix ,
ADN_datum_all , datum ,
ADN_none ) ;
if( THD_deathcon() && THD_is_file(outset->dblk->diskptr->header_name) ){
fprintf(stderr,
"*** Output file %s already exists -- cannot continue!\n",
outset->dblk->diskptr->header_name ) ;
exit(1) ;
}
}
/*-- read data from disk --*/
DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
if( verb ) fprintf(stderr," ++ read in dataset %s\n",insetname) ;
/*-- Edge detect each sub-brik --*/
indims[0] = DSET_NX(inset);
indims[1] = DSET_NY(inset);
indims[2] = DSET_NZ(inset);
border[0] = 50;
border[1] = 50;
border[2] = 50;
for( kk=0 ; kk < nval ; kk++ ){
if( verb )
fprintf(stderr," + sub-brick %d [%s]\n",
kk,MRI_TYPE_name[DSET_BRICK_TYPE(inset,kk)] ) ;
if (DSET_BRICK_TYPE(inset,kk) != DSET_BRICK_TYPE(inset,0)) {
fprintf(stderr,"ERROR: Sub-bricks of different types.\n"
"This is not splenda\n") ;
exit(1) ;
}
switch( DSET_BRICK_TYPE(inset,kk) ){
default:
fprintf(stderr,"ERROR: illegal input sub-brick datum\n") ;
exit(1) ;
case MRI_float:{
float *pp = (float *) DSET_ARRAY(inset,kk) ;
float fac = DSET_BRICK_FACTOR(inset,kk) ;
if( fac ) {
for( ii=0 ; ii < nxyz ; ii++ ) { pp[ii] *= fac; }
}
if ( Extract_Gradient_Maxima_3D( (void *)pp, FLOAT,
sum[kk], FLOAT,
indims,
border,
filterCoefs,
filterType ) == 0 ) {
fprintf( stderr, "ERROR: gradient extraction failed.\n" );
exit( 1 );
}
}
break ;
case MRI_short:{
short *pp = (short *) DSET_ARRAY(inset,kk) ;
float fac = DSET_BRICK_FACTOR(inset,kk) ;
if( fac ) {
for( ii=0 ; ii < nxyz ; ii++ ) { pp[ii] *= fac; }
}
if ( Extract_Gradient_Maxima_3D( (void *)pp, USHORT,
sum[kk], FLOAT,
indims,
border,
filterCoefs,
filterType ) == 0 ) {
fprintf( stderr, "ERROR: gradient extraction failed.\n" );
exit( 1 );
}
}
break ;
case MRI_byte:{
byte *pp = (byte *) DSET_ARRAY(inset,kk) ;
float fac = DSET_BRICK_FACTOR(inset,kk) ;
if( fac ) {
for( ii=0 ; ii < nxyz ; ii++ ) { pp[ii] *= fac; }
}
if ( Extract_Gradient_Maxima_3D( (void *)pp, UCHAR,
sum[kk], FLOAT,
indims,
border,
filterCoefs,
filterType ) == 0 ) {
fprintf( stderr, "ERROR: gradient extraction failed.\n" );
exit( 1 );
}
}
break ;
}
}
DSET_delete(inset) ;
}
switch( datum ){
default:
fprintf(stderr,
"*** Fatal Error ***\n"
"*** Somehow ended up with datum = %d\n",datum) ;
exit(1) ;
case MRI_float:{
for( kk=0 ; kk < nval ; kk++ ){
EDIT_substitute_brick(outset, kk, MRI_float, sum[kk]);
DSET_BRICK_FACTOR(outset, kk) = 0.0;
}
}
break ;
case MRI_byte:
case MRI_short:{
void ** dfim ;
float gtop , fimfac , gtemp ;
if( verb )
fprintf(stderr," ++ Scaling output to type %s brick(s)\n",
MRI_TYPE_name[datum] ) ;
dfim = (void **) malloc(sizeof(void *)*nval) ;
if( gscale ){ /* allow global scaling */
gtop = 0.0 ;
for( kk=0 ; kk < nval ; kk++ ){
gtemp = MCW_vol_amax( nxyz , 1 , 1 , MRI_float, sum[kk] ) ;
gtop = MAX( gtop , gtemp ) ;
if( gtemp == 0.0 )
fprintf(stderr," -- Warning: output sub-brick %d is all zeros!\n",kk) ;
}
}
for (kk = 0 ; kk < nval ; kk ++ ) {
if( ! gscale ){
gtop = MCW_vol_amax( nxyz , 1 , 1 , MRI_float, sum[kk] ) ;
if( gtop == 0.0 )
fprintf(stderr," -- Warning: output sub-brick %d is all zeros!\n",kk) ;
}
if( fscale ){
fimfac = (gtop > 0.0) ? MRI_TYPE_maxval[datum] / gtop : 0.0 ;
} else if( !nscale ){
fimfac = (gtop > MRI_TYPE_maxval[datum] || (gtop > 0.0 && gtop <= 1.0) )
? MRI_TYPE_maxval[datum]/ gtop : 0.0 ;
} else {
fimfac = 0.0 ;
}
if( verb ){
if( fimfac != 0.0 )
fprintf(stderr," ++ Sub-brick %d scale factor = %f\n",kk,fimfac) ;
else
fprintf(stderr," ++ Sub-brick %d: no scale factor\n" ,kk) ;
}
dfim[kk] = (void *) malloc( mri_datum_size(datum) * nxyz ) ;
if( dfim[kk] == NULL ){ fprintf(stderr,"*** malloc fails at output\n");exit(1); }
EDIT_coerce_scale_type( nxyz , fimfac ,
MRI_float, sum[kk] , datum,dfim[kk] ) ;
free( sum[kk] ) ;
EDIT_substitute_brick(outset, kk, datum, dfim[kk] );
DSET_BRICK_FACTOR(outset,kk) = (fimfac != 0.0) ? 1.0/fimfac : 0.0 ;
}
}
break ;
}
if( verb ) fprintf(stderr," ++ Computing output statistics\n") ;
THD_load_statistics( outset ) ;
THD_write_3dim_dataset( NULL,NULL , outset , True ) ;
if( verb ) fprintf(stderr," ++ Wrote output: %s\n",DSET_BRIKNAME(outset)) ;
exit(0) ;
}
syntax highlighted by Code2HTML, v. 0.9.1