diff --git a/CGAL_ImageIO/src/CGAL_ImageIO/mincio.cpp b/CGAL_ImageIO/src/CGAL_ImageIO/mincio.cpp index 19a22154c93..8bf8b3351d7 100644 --- a/CGAL_ImageIO/src/CGAL_ImageIO/mincio.cpp +++ b/CGAL_ImageIO/src/CGAL_ImageIO/mincio.cpp @@ -19,479 +19,9 @@ // // Author(s) : ASCLEPIOS Project (INRIA Sophia-Antipolis), Laurent Rineau -#ifdef _MSC_VER -#pragma warning ( disable : 4068 4786 4081 ) -#endif - -#ifdef MINC_FILES +#ifndef CGAL_HEADER_ONLY #include "mincio.h" -#include -#include -#include -#include +#include "mincio_impl.h" - - -/** Magic header for MINC (MNI NetCDF) file format */ -#ifdef MINC_FILES -# define MINC_MAGIC "CDF" -#endif - - - -int readMincHeader(_image *im, const char* filename, - double *startx, double *starty, double *startz, - double *stepx, double *stepy, double *stepz, - double *Xcosines, double *Ycosines, double *Zcosines) { - int fin, varid, minc_icv, id, ndims, i, dim[256], length, - minfound = 0, maxfound = 0; - long start[3] = {0, 0, 0}, count[3]; - unsigned int plane1, plane2, line, shift[3], ct[3]; - nc_type type; - char sign_type[MI_MAX_ATTSTR_LEN], name[MAX_NC_NAME]; - void *data; - double range[2], swap; - - /* open minc file */ - fin = miopen((char *) filename, NC_NOWRITE); - if(fin == MI_ERROR) return 0; - - /* get number of dimensions and type */ - id = ncvarid(fin, MIimage); - (void) ncvarinq(fin, id, NULL, &type, &ndims, dim, NULL); - if(ndims != 3) { - fprintf(stderr, "unsupported %i dimensional minc file\n", ndims); - return 0; - } - - /* get sign */ - if ((miattgetstr(fin, id, MIsigntype, MI_MAX_ATTSTR_LEN, sign_type) - == NULL) || ((strcmp(sign_type, MI_UNSIGNED)!=0) && - (strcmp(sign_type, MI_SIGNED)!=0))) { - if (type == NC_BYTE) - (void) strcpy(sign_type, MI_UNSIGNED); - else - (void) strcpy(sign_type, MI_SIGNED); - } - im->sign = (!strcmp(sign_type, MI_UNSIGNED)) ? SGN_UNSIGNED : SGN_SIGNED; - - switch(type) { - case NC_CHAR: - case NC_BYTE: - im->wdim = 1; - im->wordKind = WK_FIXED; - im->sign = SGN_UNSIGNED; - break; - case NC_SHORT: - im->wdim = 2; - im->wordKind = WK_FIXED; - break; - case NC_LONG: - im->wdim = 4; - im->wordKind = WK_FIXED; - break; - case NC_FLOAT: - im->wdim = 4; - im->wordKind = WK_FLOAT; - im->sign = SGN_UNKNOWN; - break; - case NC_DOUBLE: - im->wdim = 8; - im->wordKind = WK_FLOAT; - im->sign = SGN_UNKNOWN; - break; - default: - fprintf(stderr, "unsupported minc type\n"); - return 0; - break; - } - - - /* get dimensions length and voxel size */ - ncopts = 0; - for(i = 0; i < ndims; i++) { - (void) ncdiminq(fin, dim[i], name, &count[i]); - varid = ncvarid(fin, name); - if(!strcmp(name, MIxspace)) { - im->xdim = count[i]; - if(miattget1(fin, varid, MIstep, NC_DOUBLE, stepx) == MI_ERROR) - *stepx = 1.0; - im->vx = fabs(*stepx); - if(miattget1(fin, varid, MIstart, NC_DOUBLE, startx) == MI_ERROR) - *startx = 0.0; - if(miattget(fin, varid, MIdirection_cosines, NC_DOUBLE, 3, - (void *) Xcosines, (int *) 0) == MI_ERROR ) { - Xcosines[0] = 0.0; - Xcosines[1] = 0.0; - Xcosines[2] = 0.0; - } - } - else if(!strcmp(name, MIyspace)) { - im->ydim = count[i]; - if(miattget1(fin, varid, MIstep, NC_DOUBLE, stepy) == MI_ERROR) - *stepy = 1.0; - im->vy = fabs(*stepy); - if(miattget1(fin, varid, MIstart, NC_DOUBLE, starty) == MI_ERROR) - *starty = 0.0; - if(miattget(fin, varid, MIdirection_cosines, NC_DOUBLE, 3, - (void *) Ycosines, (int *) 0) == MI_ERROR ) { - Ycosines[0] = 0.0; - Ycosines[1] = 0.0; - Ycosines[2] = 0.0; - } - } - else if(!strcmp(name, MIzspace)) { - im->zdim = count[i]; - if(miattget1(fin, varid, MIstep, NC_DOUBLE, stepz) == MI_ERROR) - *stepz = 1.0; - im->vz = fabs(*stepz); - if(miattget1(fin, varid, MIstart, NC_DOUBLE, startz) == MI_ERROR) - *startz = 0.0; - if(miattget(fin, varid, MIdirection_cosines, NC_DOUBLE, 3, - (void *) Zcosines, (int *) 0) == MI_ERROR ) { - Zcosines[0] = 0.0; - Zcosines[1] = 0.0; - Zcosines[2] = 0.0; - } - } else { - fprintf(stderr, "unsupported dimension %s\n", name); - return 0; - } - } - - - if( miattget( fin, id, MIvalid_range, NC_DOUBLE, 2, - (void *) range, &length) == MI_ERROR || length != 2 ) { - if( miattget1( fin, id, MIvalid_min, NC_DOUBLE, - (void *) &range[0] ) != MI_ERROR ) - minfound = 1; - if( miattget1( fin, id, MIvalid_max, NC_DOUBLE, - (void *) &range[1] ) != MI_ERROR ) - maxfound = 1; - } - else { - if( range[0] > range[1] ) { - swap = range[0]; - range[0] = range[1]; - range[1] = swap; - } - minfound = 1; - maxfound = 1; - } - - if(!minfound) { - switch(type) { - case NC_CHAR: - case NC_BYTE: - range[0] = 0; - break; - case NC_SHORT: - if(im->sign == SGN_SIGNED) range[0] = -32768.0; - else range[0] = 0; - break; - case NC_LONG: - if(im->sign == SGN_SIGNED) range[0] = -2147483648.0; - else range[0] = 0; - break; - default: - range[0] = 0.0; - break; - } - } - if(!maxfound) { - switch(type) { - case NC_CHAR: - case NC_BYTE: - range[1] = 255.0; - break; - case NC_SHORT: - if(im->sign == SGN_SIGNED) range[1] = 32767.0; - else range[1] = 65535.0; - break; - case NC_LONG: - if(im->sign == SGN_SIGNED) range[1] = 2147483647.0; - else range[1] = 4294967295.0; - break; - default: - range[1] = 1.0; - break; - } - } - - ncopts = NC_VERBOSE | NC_FATAL; - im->vdim = 1; - - /* read data bytes */ - im->data = ImageIO_alloc(im->xdim * im->ydim * im->zdim * im->vdim * - im->wdim); - - minc_icv = miicv_create(); - (void) miicv_setint( minc_icv, MI_ICV_TYPE, type ); - (void) miicv_setstr( minc_icv, MI_ICV_SIGN, sign_type ); - (void) miicv_setint( minc_icv, MI_ICV_DO_NORM, 1 ); - (void) miicv_setint( minc_icv, MI_ICV_DO_FILLVALUE, 1 ); - (void) miicv_setdbl( minc_icv, MI_ICV_VALID_MIN, range[0] ); - (void) miicv_setdbl( minc_icv, MI_ICV_VALID_MAX, range[1] ); - - (void) miicv_attach( minc_icv, fin, id ); - - if(miicv_get(minc_icv, start, count, im->data) == MI_ERROR) { - fprintf(stderr, "error while reading file\n"); - return 0; - } - (void) miicv_detach( minc_icv ); - (void) miicv_free( minc_icv ); - - ImageIO_closeImage(im); - - /* order data in ZYX */ - (void) ncdiminq(fin, dim[0], name, NULL); - if(!strcmp(name, MIzspace)) { - (void) ncdiminq(fin, dim[1], name, NULL); - /* file is ZYX */ - if(!strcmp(name, MIyspace)) { - miclose(fin); - return 1; - } - } - - (void) ncdiminq(fin, dim[0], name, NULL); - /* file is ZXY */ - if(!strcmp(name, MIzspace)) { - shift[0] = 0; - shift[1] = 2; - shift[2] = 1; - plane2 = im->xdim * im->ydim; - line = im->ydim; - } - /* file is Y first */ - else if(!strcmp(name, MIyspace)) { - shift[0] = 1; - plane2 = im->xdim * im->zdim; - (void) ncdiminq(fin, dim[1], name, NULL); - /* file is YXZ */ - if(!strcmp(name, MIxspace)) { - shift[1] = 2; - shift[2] = 0; - line = im->zdim; - } - /* file is YZX */ - else { - shift[1] = 0; - shift[2] = 2; - line = im->xdim; - } - } - /* file is X first */ - else { - shift[0] = 2; - plane2 = im->ydim * im->zdim; - (void) ncdiminq(fin, dim[1], name, NULL); - /* file is XYZ */ - if(!strcmp(name, MIyspace)) { - shift[1] = 1; - shift[2] = 0; - line = im->zdim; - } - /* file is XZY */ - else { - shift[1] = 0; - shift[2] = 1; - line = im->ydim; - } - } - - data = ImageIO_alloc(im->xdim * im->ydim * im->zdim * im->vdim * im->wdim); - plane1 = im->xdim * im->ydim; - - for(ct[0] = 0; ct[0] < im->zdim; ct[0]++) { - for(ct[1] = 0; ct[1] < im->ydim; ct[1]++) { - for(ct[2] = 0; ct[2] < im->xdim; ct[2]++) { - memcpy((unsigned char *) data + - (ct[0] * plane1 + ct[1] * im->xdim + ct[2]) * im->wdim, - (unsigned char *) im->data + - (ct[shift[0]] * plane2 + ct[shift[1]] * line + ct[shift[2]]) - * im->wdim, im->wdim); - } - } - } - - ImageIO_free(im->data); - im->data = data; - - /* close file */ - miclose(fin); - return 1; -} - - - - -int writeMincFile( const _image* im, const char *filename, - const char *sourcename, - double startx, double starty, double startz, - double stepx, double stepy, double stepz, - const double *Xcosine, const double *Ycosine, - const double *Zcosine, const double *range ) { - int i, cdfid, minc_icv, dim_vars[3], dim_ids[3], img_var_id, min_id, max_id, - src_cdfid, src_img_var, varid, n_excluded, excluded_vars[10], - src_min_id, src_max_id, src_root_id, old_att_length, same_file = 0; - long org[3] = {0, 0, 0}, count[3] = { im->zdim, im->ydim, im->xdim }; - char *dim_names[3] = { MIzspace, MIyspace, MIxspace }, *history, - *excluded_list[3] = { MIxspace, MIyspace, MIzspace }, *newname; - nc_type type, datatype; - double start[3] = { startz, starty, startx }, - vx[3] = { stepz, stepy, stepx }; - - /* if filename is the same file that sourcename, it is needed to create - a temporary file for output */ - newname = (char *) filename; - if(sourcename) { - struct stat out, org; - same_file = (stat(sourcename, &org) == 0) && (stat(filename, &out) == 0) - && (org.st_ino == out.st_ino); - if(same_file) { - newname = (char *) malloc((strlen(filename) + 10) * sizeof(char)); - for(i = strlen(filename) - 1; i > 0 && filename[i] != '/'; i--) ; - if(filename[i] == '/') { - strncpy(newname, filename, i + 1); - newname[i+1] = '\0'; - strcat(newname, "#TMP#"); - strcat(newname, filename + i + 1); - } - else - sprintf(newname, "#TMP#%s", filename); - } - } - - /* open file */ - cdfid = micreate( (char *) newname, NC_CLOBBER ); - if( cdfid == MI_ERROR ) { - fprintf(stderr, "Error: opening MINC file \"%s\".\n", filename ); - return -1; - } - - /* set word type */ - switch(im->wordKind) { - case WK_FIXED: - if(im->wdim == 1) type = NC_BYTE; - else if(im->wdim == 2) type = NC_SHORT; - else type = NC_LONG; - break; - case WK_FLOAT: - default: - if(im->wdim == 4) type = NC_FLOAT; - else type = NC_DOUBLE; - break; - } - - /* create dimensions variables */ - for(i = 0; i < 3; i++) - dim_vars[i] = ncdimdef( cdfid, dim_names[i], count[i] ); - - for(i = 0; i < 3; i++ ) { - dim_ids[i] = micreate_std_variable( cdfid, dim_names[i], NC_DOUBLE, - 0, NULL); - if( dim_ids[i] < 0 ) return -1; - (void) miattputdbl( cdfid, dim_ids[i], MIstep, vx[i]); - (void) miattputdbl( cdfid, dim_ids[i], MIstart, start[i]); - } - if(Zcosine[0] != 0.0 || Zcosine[1] != 0.0 || Zcosine[2] != 0.0) - (void) ncattput(cdfid, dim_ids[0], MIdirection_cosines, NC_DOUBLE, 3, - Zcosine); - if(Ycosine[0] != 0.0 || Ycosine[1] != 0.0 || Ycosine[2] != 0.0) - (void) ncattput(cdfid, dim_ids[1], MIdirection_cosines, NC_DOUBLE, 3, - Ycosine); - if(Xcosine[0] != 0.0 || Xcosine[1] != 0.0 || Xcosine[2] != 0.0) - (void) ncattput(cdfid, dim_ids[2], MIdirection_cosines, NC_DOUBLE, 3, - Xcosine); - - /* create image variable */ - img_var_id = micreate_std_variable( cdfid, MIimage, type, 3, dim_vars ); - (void) miattputstr( cdfid, img_var_id, MIcomplete, MI_FALSE ); - if(im->sign == SGN_SIGNED) - (void) miattputstr( cdfid, img_var_id, MIsigntype, MI_SIGNED ); - else - (void) miattputstr( cdfid, img_var_id, MIsigntype, MI_UNSIGNED ); - (void) ncattput(cdfid, img_var_id, MIvalid_range ,NC_DOUBLE, 2, range); - - min_id = micreate_std_variable(cdfid, MIimagemin, NC_DOUBLE, 0, 0 ); - max_id = micreate_std_variable(cdfid, MIimagemax, NC_DOUBLE, 0, 0 ); - - /* if the image comes from a minc file, duplicate all minc attributes */ - if(sourcename) { - if( (src_cdfid = miopen((char *) sourcename, NC_NOWRITE)) == MI_ERROR ) { - fprintf(stderr, "Error opening %s\n", sourcename ); - return -1; - } - - n_excluded = 0; - for( i = 0; i < 3; i++ ) { - if( (varid = ncvarid(src_cdfid, excluded_list[i] )) != MI_ERROR ) - excluded_vars[n_excluded++] = varid; - } - - if( (src_img_var = ncvarid(src_cdfid, MIimage )) != MI_ERROR ) - excluded_vars[n_excluded++] = src_img_var; - if( (src_max_id = ncvarid(src_cdfid, MIimagemax )) != MI_ERROR ) - excluded_vars[n_excluded++] = src_max_id; - if( (src_min_id = ncvarid(src_cdfid, MIimagemin )) != MI_ERROR ) - excluded_vars[n_excluded++] = src_min_id; - if( (src_root_id = ncvarid(src_cdfid, MIrootvariable )) != MI_ERROR ) - excluded_vars[n_excluded++] = src_root_id; - - (void) micopy_all_var_defs( src_cdfid, cdfid, n_excluded, excluded_vars ); - - if( src_img_var != MI_ERROR ) - (void) micopy_all_atts( src_cdfid, src_img_var, cdfid, img_var_id ); - if( src_root_id != MI_ERROR ) - (void) micopy_all_atts( src_cdfid, src_root_id, cdfid, - ncvarid( cdfid, MIrootvariable) ); - - ncopts = 0; - if(ncattinq(cdfid, NC_GLOBAL, MIhistory, &datatype, &old_att_length) - == MI_ERROR || datatype != NC_CHAR ) - old_att_length = 0; - - history = (char *) malloc(sizeof(char) * old_att_length); - history[0] = (char) 0; - (void) miattgetstr( cdfid, NC_GLOBAL, MIhistory, old_att_length+1, - history ); - ncopts = NC_VERBOSE | NC_FATAL; - (void) miattputstr( cdfid, NC_GLOBAL, MIhistory, history ); - free(history); - - (void) miclose( src_cdfid ); - } - - - /* write data */ - if(ncendef(cdfid) == MI_ERROR) return -1; - minc_icv = miicv_create(); - (void) miicv_setint( minc_icv, MI_ICV_TYPE, type); - (void) miicv_setstr( minc_icv, MI_ICV_SIGN, - (im->sign == SGN_SIGNED) ? MI_SIGNED : MI_UNSIGNED ); - (void) miicv_setint( minc_icv, MI_ICV_DO_NORM, 1 ); - (void) miicv_setint( minc_icv, MI_ICV_USER_NORM, 1 ); - (void) miicv_attach( minc_icv, cdfid, img_var_id ); - - if( miicv_put(minc_icv, org, count, im->data) == MI_ERROR ) - return -1; - - (void) miattputstr( cdfid, img_var_id, MIcomplete, MI_TRUE); - (void) miclose( cdfid ); - (void) miicv_free( minc_icv ); - - if(same_file) { - unlink(filename); - link(newname, filename); - unlink(newname); - free(newname); - } - - return 0; -} - - -#endif +#endif // CGAL_HEADER_ONLY diff --git a/CGAL_ImageIO/src/CGAL_ImageIO/mincio.h b/CGAL_ImageIO/src/CGAL_ImageIO/mincio.h index 53c33a505c0..46cd034abca 100644 --- a/CGAL_ImageIO/src/CGAL_ImageIO/mincio.h +++ b/CGAL_ImageIO/src/CGAL_ImageIO/mincio.h @@ -70,4 +70,8 @@ int writeMincFile( const _image* im, const char *name, const char *sourceName, #endif +#ifdef CGAL_HEADER_ONLY +#include "mincio_impl.h" +#endif // CGAL_HEADER_ONLY + #endif diff --git a/CGAL_ImageIO/src/CGAL_ImageIO/mincio_impl.h b/CGAL_ImageIO/src/CGAL_ImageIO/mincio_impl.h new file mode 100644 index 00000000000..a8ff990ade9 --- /dev/null +++ b/CGAL_ImageIO/src/CGAL_ImageIO/mincio_impl.h @@ -0,0 +1,502 @@ +// Copyright (c) 2005-2008 ASCLEPIOS Project, INRIA Sophia-Antipolis (France) +// All rights reserved. +// +// This file is part of the ImageIO Library, and as been adapted for +// CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the +// GNU Lesser General Public License as published by the Free Software Foundation; +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// These files are provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : ASCLEPIOS Project (INRIA Sophia-Antipolis), Laurent Rineau + +#ifdef CGAL_HEADER_ONLY +#define CGAL_INLINE_FUNCTION inline +#else +#define CGAL_INLINE_FUNCTION +#endif + +#ifdef _MSC_VER +#pragma warning ( disable : 4068 4786 4081 ) +#endif + +#ifdef MINC_FILES + +#include +#include +#include +#include + + + +/** Magic header for MINC (MNI NetCDF) file format */ +#ifdef MINC_FILES +# define MINC_MAGIC "CDF" +#endif + + +CGAL_INLINE_FUNCTION +int readMincHeader(_image *im, const char* filename, + double *startx, double *starty, double *startz, + double *stepx, double *stepy, double *stepz, + double *Xcosines, double *Ycosines, double *Zcosines) { + int fin, varid, minc_icv, id, ndims, i, dim[256], length, + minfound = 0, maxfound = 0; + long start[3] = {0, 0, 0}, count[3]; + unsigned int plane1, plane2, line, shift[3], ct[3]; + nc_type type; + char sign_type[MI_MAX_ATTSTR_LEN], name[MAX_NC_NAME]; + void *data; + double range[2], swap; + + /* open minc file */ + fin = miopen((char *) filename, NC_NOWRITE); + if(fin == MI_ERROR) return 0; + + /* get number of dimensions and type */ + id = ncvarid(fin, MIimage); + (void) ncvarinq(fin, id, NULL, &type, &ndims, dim, NULL); + if(ndims != 3) { + fprintf(stderr, "unsupported %i dimensional minc file\n", ndims); + return 0; + } + + /* get sign */ + if ((miattgetstr(fin, id, MIsigntype, MI_MAX_ATTSTR_LEN, sign_type) + == NULL) || ((strcmp(sign_type, MI_UNSIGNED)!=0) && + (strcmp(sign_type, MI_SIGNED)!=0))) { + if (type == NC_BYTE) + (void) strcpy(sign_type, MI_UNSIGNED); + else + (void) strcpy(sign_type, MI_SIGNED); + } + im->sign = (!strcmp(sign_type, MI_UNSIGNED)) ? SGN_UNSIGNED : SGN_SIGNED; + + switch(type) { + case NC_CHAR: + case NC_BYTE: + im->wdim = 1; + im->wordKind = WK_FIXED; + im->sign = SGN_UNSIGNED; + break; + case NC_SHORT: + im->wdim = 2; + im->wordKind = WK_FIXED; + break; + case NC_LONG: + im->wdim = 4; + im->wordKind = WK_FIXED; + break; + case NC_FLOAT: + im->wdim = 4; + im->wordKind = WK_FLOAT; + im->sign = SGN_UNKNOWN; + break; + case NC_DOUBLE: + im->wdim = 8; + im->wordKind = WK_FLOAT; + im->sign = SGN_UNKNOWN; + break; + default: + fprintf(stderr, "unsupported minc type\n"); + return 0; + break; + } + + + /* get dimensions length and voxel size */ + ncopts = 0; + for(i = 0; i < ndims; i++) { + (void) ncdiminq(fin, dim[i], name, &count[i]); + varid = ncvarid(fin, name); + if(!strcmp(name, MIxspace)) { + im->xdim = count[i]; + if(miattget1(fin, varid, MIstep, NC_DOUBLE, stepx) == MI_ERROR) + *stepx = 1.0; + im->vx = fabs(*stepx); + if(miattget1(fin, varid, MIstart, NC_DOUBLE, startx) == MI_ERROR) + *startx = 0.0; + if(miattget(fin, varid, MIdirection_cosines, NC_DOUBLE, 3, + (void *) Xcosines, (int *) 0) == MI_ERROR ) { + Xcosines[0] = 0.0; + Xcosines[1] = 0.0; + Xcosines[2] = 0.0; + } + } + else if(!strcmp(name, MIyspace)) { + im->ydim = count[i]; + if(miattget1(fin, varid, MIstep, NC_DOUBLE, stepy) == MI_ERROR) + *stepy = 1.0; + im->vy = fabs(*stepy); + if(miattget1(fin, varid, MIstart, NC_DOUBLE, starty) == MI_ERROR) + *starty = 0.0; + if(miattget(fin, varid, MIdirection_cosines, NC_DOUBLE, 3, + (void *) Ycosines, (int *) 0) == MI_ERROR ) { + Ycosines[0] = 0.0; + Ycosines[1] = 0.0; + Ycosines[2] = 0.0; + } + } + else if(!strcmp(name, MIzspace)) { + im->zdim = count[i]; + if(miattget1(fin, varid, MIstep, NC_DOUBLE, stepz) == MI_ERROR) + *stepz = 1.0; + im->vz = fabs(*stepz); + if(miattget1(fin, varid, MIstart, NC_DOUBLE, startz) == MI_ERROR) + *startz = 0.0; + if(miattget(fin, varid, MIdirection_cosines, NC_DOUBLE, 3, + (void *) Zcosines, (int *) 0) == MI_ERROR ) { + Zcosines[0] = 0.0; + Zcosines[1] = 0.0; + Zcosines[2] = 0.0; + } + } else { + fprintf(stderr, "unsupported dimension %s\n", name); + return 0; + } + } + + + if( miattget( fin, id, MIvalid_range, NC_DOUBLE, 2, + (void *) range, &length) == MI_ERROR || length != 2 ) { + if( miattget1( fin, id, MIvalid_min, NC_DOUBLE, + (void *) &range[0] ) != MI_ERROR ) + minfound = 1; + if( miattget1( fin, id, MIvalid_max, NC_DOUBLE, + (void *) &range[1] ) != MI_ERROR ) + maxfound = 1; + } + else { + if( range[0] > range[1] ) { + swap = range[0]; + range[0] = range[1]; + range[1] = swap; + } + minfound = 1; + maxfound = 1; + } + + if(!minfound) { + switch(type) { + case NC_CHAR: + case NC_BYTE: + range[0] = 0; + break; + case NC_SHORT: + if(im->sign == SGN_SIGNED) range[0] = -32768.0; + else range[0] = 0; + break; + case NC_LONG: + if(im->sign == SGN_SIGNED) range[0] = -2147483648.0; + else range[0] = 0; + break; + default: + range[0] = 0.0; + break; + } + } + if(!maxfound) { + switch(type) { + case NC_CHAR: + case NC_BYTE: + range[1] = 255.0; + break; + case NC_SHORT: + if(im->sign == SGN_SIGNED) range[1] = 32767.0; + else range[1] = 65535.0; + break; + case NC_LONG: + if(im->sign == SGN_SIGNED) range[1] = 2147483647.0; + else range[1] = 4294967295.0; + break; + default: + range[1] = 1.0; + break; + } + } + + ncopts = NC_VERBOSE | NC_FATAL; + im->vdim = 1; + + /* read data bytes */ + im->data = ImageIO_alloc(im->xdim * im->ydim * im->zdim * im->vdim * + im->wdim); + + minc_icv = miicv_create(); + (void) miicv_setint( minc_icv, MI_ICV_TYPE, type ); + (void) miicv_setstr( minc_icv, MI_ICV_SIGN, sign_type ); + (void) miicv_setint( minc_icv, MI_ICV_DO_NORM, 1 ); + (void) miicv_setint( minc_icv, MI_ICV_DO_FILLVALUE, 1 ); + (void) miicv_setdbl( minc_icv, MI_ICV_VALID_MIN, range[0] ); + (void) miicv_setdbl( minc_icv, MI_ICV_VALID_MAX, range[1] ); + + (void) miicv_attach( minc_icv, fin, id ); + + if(miicv_get(minc_icv, start, count, im->data) == MI_ERROR) { + fprintf(stderr, "error while reading file\n"); + return 0; + } + (void) miicv_detach( minc_icv ); + (void) miicv_free( minc_icv ); + + ImageIO_closeImage(im); + + /* order data in ZYX */ + (void) ncdiminq(fin, dim[0], name, NULL); + if(!strcmp(name, MIzspace)) { + (void) ncdiminq(fin, dim[1], name, NULL); + /* file is ZYX */ + if(!strcmp(name, MIyspace)) { + miclose(fin); + return 1; + } + } + + (void) ncdiminq(fin, dim[0], name, NULL); + /* file is ZXY */ + if(!strcmp(name, MIzspace)) { + shift[0] = 0; + shift[1] = 2; + shift[2] = 1; + plane2 = im->xdim * im->ydim; + line = im->ydim; + } + /* file is Y first */ + else if(!strcmp(name, MIyspace)) { + shift[0] = 1; + plane2 = im->xdim * im->zdim; + (void) ncdiminq(fin, dim[1], name, NULL); + /* file is YXZ */ + if(!strcmp(name, MIxspace)) { + shift[1] = 2; + shift[2] = 0; + line = im->zdim; + } + /* file is YZX */ + else { + shift[1] = 0; + shift[2] = 2; + line = im->xdim; + } + } + /* file is X first */ + else { + shift[0] = 2; + plane2 = im->ydim * im->zdim; + (void) ncdiminq(fin, dim[1], name, NULL); + /* file is XYZ */ + if(!strcmp(name, MIyspace)) { + shift[1] = 1; + shift[2] = 0; + line = im->zdim; + } + /* file is XZY */ + else { + shift[1] = 0; + shift[2] = 1; + line = im->ydim; + } + } + + data = ImageIO_alloc(im->xdim * im->ydim * im->zdim * im->vdim * im->wdim); + plane1 = im->xdim * im->ydim; + + for(ct[0] = 0; ct[0] < im->zdim; ct[0]++) { + for(ct[1] = 0; ct[1] < im->ydim; ct[1]++) { + for(ct[2] = 0; ct[2] < im->xdim; ct[2]++) { + memcpy((unsigned char *) data + + (ct[0] * plane1 + ct[1] * im->xdim + ct[2]) * im->wdim, + (unsigned char *) im->data + + (ct[shift[0]] * plane2 + ct[shift[1]] * line + ct[shift[2]]) + * im->wdim, im->wdim); + } + } + } + + ImageIO_free(im->data); + im->data = data; + + /* close file */ + miclose(fin); + return 1; +} + + + +CGAL_INLINE_FUNCTION +int writeMincFile( const _image* im, const char *filename, + const char *sourcename, + double startx, double starty, double startz, + double stepx, double stepy, double stepz, + const double *Xcosine, const double *Ycosine, + const double *Zcosine, const double *range ) { + int i, cdfid, minc_icv, dim_vars[3], dim_ids[3], img_var_id, min_id, max_id, + src_cdfid, src_img_var, varid, n_excluded, excluded_vars[10], + src_min_id, src_max_id, src_root_id, old_att_length, same_file = 0; + long org[3] = {0, 0, 0}, count[3] = { im->zdim, im->ydim, im->xdim }; + char *dim_names[3] = { MIzspace, MIyspace, MIxspace }, *history, + *excluded_list[3] = { MIxspace, MIyspace, MIzspace }, *newname; + nc_type type, datatype; + double start[3] = { startz, starty, startx }, + vx[3] = { stepz, stepy, stepx }; + + /* if filename is the same file that sourcename, it is needed to create + a temporary file for output */ + newname = (char *) filename; + if(sourcename) { + struct stat out, org; + same_file = (stat(sourcename, &org) == 0) && (stat(filename, &out) == 0) + && (org.st_ino == out.st_ino); + if(same_file) { + newname = (char *) malloc((strlen(filename) + 10) * sizeof(char)); + for(i = strlen(filename) - 1; i > 0 && filename[i] != '/'; i--) ; + if(filename[i] == '/') { + strncpy(newname, filename, i + 1); + newname[i+1] = '\0'; + strcat(newname, "#TMP#"); + strcat(newname, filename + i + 1); + } + else + sprintf(newname, "#TMP#%s", filename); + } + } + + /* open file */ + cdfid = micreate( (char *) newname, NC_CLOBBER ); + if( cdfid == MI_ERROR ) { + fprintf(stderr, "Error: opening MINC file \"%s\".\n", filename ); + return -1; + } + + /* set word type */ + switch(im->wordKind) { + case WK_FIXED: + if(im->wdim == 1) type = NC_BYTE; + else if(im->wdim == 2) type = NC_SHORT; + else type = NC_LONG; + break; + case WK_FLOAT: + default: + if(im->wdim == 4) type = NC_FLOAT; + else type = NC_DOUBLE; + break; + } + + /* create dimensions variables */ + for(i = 0; i < 3; i++) + dim_vars[i] = ncdimdef( cdfid, dim_names[i], count[i] ); + + for(i = 0; i < 3; i++ ) { + dim_ids[i] = micreate_std_variable( cdfid, dim_names[i], NC_DOUBLE, + 0, NULL); + if( dim_ids[i] < 0 ) return -1; + (void) miattputdbl( cdfid, dim_ids[i], MIstep, vx[i]); + (void) miattputdbl( cdfid, dim_ids[i], MIstart, start[i]); + } + if(Zcosine[0] != 0.0 || Zcosine[1] != 0.0 || Zcosine[2] != 0.0) + (void) ncattput(cdfid, dim_ids[0], MIdirection_cosines, NC_DOUBLE, 3, + Zcosine); + if(Ycosine[0] != 0.0 || Ycosine[1] != 0.0 || Ycosine[2] != 0.0) + (void) ncattput(cdfid, dim_ids[1], MIdirection_cosines, NC_DOUBLE, 3, + Ycosine); + if(Xcosine[0] != 0.0 || Xcosine[1] != 0.0 || Xcosine[2] != 0.0) + (void) ncattput(cdfid, dim_ids[2], MIdirection_cosines, NC_DOUBLE, 3, + Xcosine); + + /* create image variable */ + img_var_id = micreate_std_variable( cdfid, MIimage, type, 3, dim_vars ); + (void) miattputstr( cdfid, img_var_id, MIcomplete, MI_FALSE ); + if(im->sign == SGN_SIGNED) + (void) miattputstr( cdfid, img_var_id, MIsigntype, MI_SIGNED ); + else + (void) miattputstr( cdfid, img_var_id, MIsigntype, MI_UNSIGNED ); + (void) ncattput(cdfid, img_var_id, MIvalid_range ,NC_DOUBLE, 2, range); + + min_id = micreate_std_variable(cdfid, MIimagemin, NC_DOUBLE, 0, 0 ); + max_id = micreate_std_variable(cdfid, MIimagemax, NC_DOUBLE, 0, 0 ); + + /* if the image comes from a minc file, duplicate all minc attributes */ + if(sourcename) { + if( (src_cdfid = miopen((char *) sourcename, NC_NOWRITE)) == MI_ERROR ) { + fprintf(stderr, "Error opening %s\n", sourcename ); + return -1; + } + + n_excluded = 0; + for( i = 0; i < 3; i++ ) { + if( (varid = ncvarid(src_cdfid, excluded_list[i] )) != MI_ERROR ) + excluded_vars[n_excluded++] = varid; + } + + if( (src_img_var = ncvarid(src_cdfid, MIimage )) != MI_ERROR ) + excluded_vars[n_excluded++] = src_img_var; + if( (src_max_id = ncvarid(src_cdfid, MIimagemax )) != MI_ERROR ) + excluded_vars[n_excluded++] = src_max_id; + if( (src_min_id = ncvarid(src_cdfid, MIimagemin )) != MI_ERROR ) + excluded_vars[n_excluded++] = src_min_id; + if( (src_root_id = ncvarid(src_cdfid, MIrootvariable )) != MI_ERROR ) + excluded_vars[n_excluded++] = src_root_id; + + (void) micopy_all_var_defs( src_cdfid, cdfid, n_excluded, excluded_vars ); + + if( src_img_var != MI_ERROR ) + (void) micopy_all_atts( src_cdfid, src_img_var, cdfid, img_var_id ); + if( src_root_id != MI_ERROR ) + (void) micopy_all_atts( src_cdfid, src_root_id, cdfid, + ncvarid( cdfid, MIrootvariable) ); + + ncopts = 0; + if(ncattinq(cdfid, NC_GLOBAL, MIhistory, &datatype, &old_att_length) + == MI_ERROR || datatype != NC_CHAR ) + old_att_length = 0; + + history = (char *) malloc(sizeof(char) * old_att_length); + history[0] = (char) 0; + (void) miattgetstr( cdfid, NC_GLOBAL, MIhistory, old_att_length+1, + history ); + ncopts = NC_VERBOSE | NC_FATAL; + (void) miattputstr( cdfid, NC_GLOBAL, MIhistory, history ); + free(history); + + (void) miclose( src_cdfid ); + } + + + /* write data */ + if(ncendef(cdfid) == MI_ERROR) return -1; + minc_icv = miicv_create(); + (void) miicv_setint( minc_icv, MI_ICV_TYPE, type); + (void) miicv_setstr( minc_icv, MI_ICV_SIGN, + (im->sign == SGN_SIGNED) ? MI_SIGNED : MI_UNSIGNED ); + (void) miicv_setint( minc_icv, MI_ICV_DO_NORM, 1 ); + (void) miicv_setint( minc_icv, MI_ICV_USER_NORM, 1 ); + (void) miicv_attach( minc_icv, cdfid, img_var_id ); + + if( miicv_put(minc_icv, org, count, im->data) == MI_ERROR ) + return -1; + + (void) miattputstr( cdfid, img_var_id, MIcomplete, MI_TRUE); + (void) miclose( cdfid ); + (void) miicv_free( minc_icv ); + + if(same_file) { + unlink(filename); + link(newname, filename); + unlink(newname); + free(newname); + } + + return 0; +} + + +#endif diff --git a/CGAL_ImageIO/src/CGAL_ImageIO/recbuffer.cpp b/CGAL_ImageIO/src/CGAL_ImageIO/recbuffer.cpp index 64f5f5edf9f..e5d6f002a90 100644 --- a/CGAL_ImageIO/src/CGAL_ImageIO/recbuffer.cpp +++ b/CGAL_ImageIO/src/CGAL_ImageIO/recbuffer.cpp @@ -19,1719 +19,9 @@ // // Author(s) : ASCLEPIOS Project (INRIA Sophia-Antipolis), Laurent Rineau -/************************************************************************* - * recbuffer.c - tools for recursive filtering of 3D and 2D image buffers - * - * $Id$ - * - * Copyright©INRIA 1999 - * - * DESCRIPTION: - * - * recursive filtering of a buffer (a [1,2,3]D array) - * according that the filtering is separable - * - * - * - * AUTHOR: - * Gregoire Malandain (greg@sophia.inria.fr) - * - * CREATION DATE: - * June, 9 1998 - * - * Copyright Gregoire Malandain, INRIA - * - * ADDITIONS, CHANGES - * - * * Jul 6 1999 (Gregoire Malandain) - * a bug in RecursiveFilterOnBuffer (*&%^@$^ cut and paste) - * - */ -#include -#include -#include +#ifndef CGAL_HEADER_ONLY - - -#include "convert.h" #include "recbuffer.h" +#include "recbuffer_impl.h" -static int _VERBOSE_ = 0; - -#define EXIT_ON_FAILURE 0 -#define EXIT_ON_SUCCESS 1 - - - - - - - - - - -/* - * - * Gradient modulus - * - * - */ - -int GradientModulus( void *bufferIn, - bufferType typeIn, - void *bufferOut, - bufferType typeOut, - int *bufferDims, - int *borderLengths, - float *filterCoefs, - recursiveFilterType filterType ) -{ - const char *proc = "GradientModulus"; - float *auxBuf = NULL; - float *tmpBuf = NULL, *grdBuf = NULL; - int sizeAuxBuf = 0; - derivativeOrder derivatives[3]; - int i; - - - sizeAuxBuf = bufferDims[0] * bufferDims[1] * bufferDims[2]; - if ( typeOut != CGAL_FLOAT || bufferIn == bufferOut ) - sizeAuxBuf *= 2; - - - /* allocation des buffers de calcul - */ - auxBuf = (float*)malloc( sizeAuxBuf * sizeof(float) ); - if ( auxBuf == NULL ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, "%s: unable to allocate auxiliary buffer\n", proc ); - return( EXIT_ON_FAILURE ); - } - tmpBuf = auxBuf; - if ( typeOut != CGAL_FLOAT || bufferIn == bufferOut ) { - grdBuf = tmpBuf; - grdBuf += bufferDims[0] * bufferDims[1] * bufferDims[2]; - } else { - grdBuf = (float*)bufferOut; - } - - /* cas 2D - */ - if ( bufferDims[2] == 1 ) { - - derivatives[0] = DERIVATIVE_1; - derivatives[1] = DERIVATIVE_0; - derivatives[2] = NODERIVATIVE; - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)grdBuf, CGAL_FLOAT, - bufferDims, borderLengths, derivatives, - filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { - if ( _VERBOSE_ ) - fprintf( stderr, "%s: unable to compute X derivative (2D)\n", proc ); - free( auxBuf ); - return( EXIT_ON_FAILURE ); - } - - derivatives[0] = DERIVATIVE_0; - derivatives[1] = DERIVATIVE_1; - derivatives[2] = NODERIVATIVE; - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)tmpBuf, CGAL_FLOAT, - bufferDims, borderLengths, derivatives, - filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { - if ( _VERBOSE_ ) - fprintf( stderr, "%s: unable to compute Y derivative (2D)\n", proc ); - free( auxBuf ); - return( EXIT_ON_FAILURE ); - } - - sizeAuxBuf = bufferDims[0] * bufferDims[1] * bufferDims[2]; - for ( i = 0; i < sizeAuxBuf; i++ ) - grdBuf[i] = (float)sqrt( grdBuf[i]*grdBuf[i] + tmpBuf[i]*tmpBuf[i] ); - - } else { - - derivatives[0] = NODERIVATIVE; - derivatives[1] = NODERIVATIVE; - derivatives[2] = DERIVATIVE_0; - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)tmpBuf, CGAL_FLOAT, - bufferDims, borderLengths, derivatives, - filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { - if ( _VERBOSE_ ) - fprintf( stderr, "%s: unable to compute Z smoothing (3D)\n", proc ); - free( auxBuf ); - return( EXIT_ON_FAILURE ); - } - - derivatives[0] = DERIVATIVE_1; - derivatives[1] = DERIVATIVE_0; - derivatives[2] = NODERIVATIVE; - if ( RecursiveFilterOnBuffer( (void*)tmpBuf, CGAL_FLOAT, (void*)grdBuf, CGAL_FLOAT, - bufferDims, borderLengths, derivatives, - filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { - if ( _VERBOSE_ ) - fprintf( stderr, "%s: unable to compute X derivative (3D)\n", proc ); - free( auxBuf ); - return( EXIT_ON_FAILURE ); - } - - derivatives[0] = DERIVATIVE_0; - derivatives[1] = DERIVATIVE_1; - derivatives[2] = NODERIVATIVE; - if ( RecursiveFilterOnBuffer( (void*)tmpBuf, CGAL_FLOAT, (void*)tmpBuf, CGAL_FLOAT, - bufferDims, borderLengths, derivatives, - filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { - if ( _VERBOSE_ ) - fprintf( stderr, "%s: unable to compute Y derivative (3D)\n", proc ); - free( auxBuf ); - return( EXIT_ON_FAILURE ); - } - - sizeAuxBuf = bufferDims[0] * bufferDims[1] * bufferDims[2]; - for ( i = 0; i < sizeAuxBuf; i++ ) - grdBuf[i] = grdBuf[i]*grdBuf[i] + tmpBuf[i]*tmpBuf[i]; - - derivatives[0] = DERIVATIVE_0; - derivatives[1] = DERIVATIVE_0; - derivatives[2] = DERIVATIVE_1; - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)tmpBuf, CGAL_FLOAT, - bufferDims, borderLengths, derivatives, - filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { - if ( _VERBOSE_ ) - fprintf( stderr, "%s: unable to compute Z derivative (3D)\n", proc ); - free( auxBuf ); - return( EXIT_ON_FAILURE ); - } - - for ( i = 0; i < sizeAuxBuf; i++ ) - grdBuf[i] = (float)sqrt( grdBuf[i] + tmpBuf[i]*tmpBuf[i] ); - - } - - if ( grdBuf != bufferOut ) - ConvertBuffer( grdBuf, CGAL_FLOAT, bufferOut, typeOut, - bufferDims[0]*bufferDims[1]*bufferDims[2] ); - free( auxBuf ); - return( EXIT_ON_SUCCESS ); -} - - - - - - - - - - - - - - - - - - - - - - - - -/* - * - * Laplacian - * - * - */ -int Laplacian_2D ( void *bufferIn, - bufferType typeIn, - void *bufferOut, - bufferType typeOut, - int *bufferDims, - int *borderLengths, - float *filterCoefs, - recursiveFilterType filterType ) -{ - const char *proc = "Laplacian_2D"; - float *theXX = NULL; - float *theYY = NULL; - - derivativeOrder XXderiv[3] = { DERIVATIVE_2, SMOOTHING, NODERIVATIVE }; - derivativeOrder YYderiv[3] = { SMOOTHING, DERIVATIVE_2, NODERIVATIVE }; - int sliceDims[3]; - int z, i, dimxXdimy; - - void *sliceOut = NULL; - - - - /* - * We check the buffers' dimensions. - */ - if ( (bufferDims[0] <= 0) || (bufferDims[1] <= 0) || (bufferDims[2] <= 0) ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc ); - return( EXIT_ON_FAILURE ); - } - - /* - * test of the coefficients - */ - if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) || - (filterCoefs[2] < 0.0) ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc ); - return( EXIT_ON_FAILURE ); - } - - - /* - * - */ - dimxXdimy = bufferDims[0] * bufferDims[1]; - sliceDims[0] = bufferDims[0]; - sliceDims[1] = bufferDims[1]; - sliceDims[2] = 1; - - - if ( typeOut == CGAL_FLOAT ) { - theXX = (float*)malloc( dimxXdimy * sizeof( float ) ); - } else { - theXX = (float*)malloc( 2 * dimxXdimy * sizeof( float ) ); - } - - if ( theXX == NULL ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to allocate auxiliary buffer.\n" ); - } - return( EXIT_ON_FAILURE ); - } - - if ( typeOut != CGAL_FLOAT ) { - theYY = theXX; - theYY += dimxXdimy; - } - - - - for ( z=0; z 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute X^2 derivative.\n" ); - } - free( theXX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theYY, CGAL_FLOAT, - sliceDims, borderLengths, - YYderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Y^2 derivative.\n" ); - } - free( theXX ); - return( EXIT_ON_FAILURE ); - } - - - for ( i=0; i 0 ) - fprintf( stderr, " Error in %s: such output type not handled.\n", proc ); - free( theXX ); - return( EXIT_ON_FAILURE ); - } - ConvertBuffer( theYY, CGAL_FLOAT, sliceOut, typeOut, dimxXdimy ); - } - } - - return( EXIT_ON_SUCCESS ); -} - - - - - - - - - - - - - - -int Laplacian ( void *bufferIn, - bufferType typeIn, - void *bufferOut, - bufferType typeOut, - int *bufferDims, - int *borderLengths, - float *filterCoefs, - recursiveFilterType filterType ) -{ - const char *proc = "Laplacian"; - float *theSL = NULL; - float *theZZ = NULL; - float *theZ0 = NULL; - - - derivativeOrder XXderiv[3] = { DERIVATIVE_2, SMOOTHING, NODERIVATIVE }; - derivativeOrder YYderiv[3] = { SMOOTHING, DERIVATIVE_2, NODERIVATIVE }; - derivativeOrder Zsmooth[3] = { NODERIVATIVE, NODERIVATIVE, SMOOTHING }; - derivativeOrder ZZderiv[3] = { SMOOTHING, SMOOTHING, DERIVATIVE_2 }; - - int sliceDims[3]; - int z, i, j, dimxXdimy; - - - - - /* - * We check the buffers' dimensions. - */ - if ( bufferDims[2] == 1 ) { - return( Laplacian_2D ( bufferIn, typeIn, bufferOut, typeOut, - bufferDims, borderLengths, filterCoefs, filterType ) ); - } - - if ( (bufferDims[0] <= 0) || (bufferDims[1] <= 0) || (bufferDims[2] <= 0) ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc ); - return( EXIT_ON_FAILURE ); - } - - /* - * test of the coefficients - */ - if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) || - (filterCoefs[2] < 0.0) ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc ); - return( EXIT_ON_FAILURE ); - } - - - /* - * - */ - dimxXdimy = bufferDims[0] * bufferDims[1]; - sliceDims[0] = bufferDims[0]; - sliceDims[1] = bufferDims[1]; - sliceDims[2] = 1; - - - - if ( typeOut == CGAL_FLOAT ) { - theSL = (float*)malloc( (1+bufferDims[2]) * dimxXdimy * sizeof( float ) ); - } else { - theSL = (float*)malloc( (1+2*bufferDims[2]) * dimxXdimy * sizeof( float ) ); - } - - - - if ( theSL == NULL ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to allocate auxiliary buffer.\n" ); - } - return( EXIT_ON_FAILURE ); - } - - - - theZ0 = theSL; - theZ0 += dimxXdimy; - - - - if ( typeOut == CGAL_FLOAT ) { - theZZ = (float*) bufferOut; - } else { - theZZ = theZ0; - theZZ += dimxXdimy * bufferDims[2]; - } - - - - /* - * - * 3D filtering / filtering along Z - * - */ - - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ0, CGAL_FLOAT, - bufferDims, borderLengths, - Zsmooth, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Z^0 derivative.\n" ); - } - free( theSL ); - return( EXIT_ON_FAILURE ); - } - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZZ, CGAL_FLOAT, - bufferDims, borderLengths, - ZZderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Z^2 derivative.\n" ); - } - free( theSL ); - return( EXIT_ON_FAILURE ); - } - - - - - - - for ( z=0; z 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute X^2 derivative.\n" ); - } - free( theSL ); - return( EXIT_ON_FAILURE ); - } - - for ( j=z*dimxXdimy, i=0; i 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Y^2 derivative.\n" ); - } - free( theSL ); - return( EXIT_ON_FAILURE ); - } - - for ( j=z*dimxXdimy, i=0; i 0 ) - fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc ); - return( EXIT_ON_FAILURE ); - } - - /* - * test of the coefficients - */ - if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) || - (filterCoefs[2] < 0.0) ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc ); - return( EXIT_ON_FAILURE ); - } - - - /* - * - */ - dimxXdimy = bufferDims[0] * bufferDims[1]; - sliceDims[0] = bufferDims[0]; - sliceDims[1] = bufferDims[1]; - sliceDims[2] = 1; - - - if ( typeOut == CGAL_FLOAT ) { - theXX = (float*)malloc( 4 * dimxXdimy * sizeof( float ) ); - } else { - theXX = (float*)malloc( 5 * dimxXdimy * sizeof( float ) ); - } - - if ( theXX == NULL ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to allocate auxiliary buffer.\n" ); - } - return( EXIT_ON_FAILURE ); - } - - - - theX = theY = theYY = theXX; - theYY += dimxXdimy; - theX += 2*dimxXdimy; - theY += 3*dimxXdimy; - - - - if ( typeOut != CGAL_FLOAT ) { - theXY = theXX; - theXY += 4*dimxXdimy; - } - - - - for ( z=0; z 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Y^0 derivative.\n" ); - } - free( theXX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( sliceIn, typeIn, theY, CGAL_FLOAT, - sliceDims, borderLengths, - Xsmooth, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute X^0 derivative.\n" ); - } - free( theXX ); - return( EXIT_ON_FAILURE ); - } - - - - - if ( RecursiveFilterOnBuffer( sliceIn, typeIn, theXY, CGAL_FLOAT, - sliceDims, borderLengths, - XYderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute X^1Y^1 derivative.\n" ); - } - free( theXX ); - return( EXIT_ON_FAILURE ); - } - - - - - if ( RecursiveFilterOnBuffer( theX, CGAL_FLOAT, theXX, CGAL_FLOAT, - sliceDims, borderLengths, - XXderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute X^2 derivative.\n" ); - } - free( theXX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( theY, CGAL_FLOAT, theYY, CGAL_FLOAT, - sliceDims, borderLengths, - YYderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Y^2 derivative.\n" ); - } - free( theXX ); - return( EXIT_ON_FAILURE ); - } - - - - - if ( RecursiveFilterOnBuffer( theX, CGAL_FLOAT, theX, CGAL_FLOAT, - sliceDims, borderLengths, - Xderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute X^1 derivative.\n" ); - } - free( theXX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( theY, CGAL_FLOAT, theY, CGAL_FLOAT, - sliceDims, borderLengths, - Yderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Y^1 derivative.\n" ); - } - free( theXX ); - return( EXIT_ON_FAILURE ); - } - - - - - for ( i=0; i 1e-10 ) theXY[i] = (float)(theXY[i] / g); - } - - if ( typeOut != CGAL_FLOAT ) { - switch ( typeOut ) { - case CGAL_UCHAR : - sliceOut = (((u8*)bufferOut) + z * dimxXdimy); - break; - case CGAL_SCHAR : - sliceOut = (((s8*)bufferOut) + z * dimxXdimy); - break; - case CGAL_SSHORT : - sliceOut = (((s16*)bufferOut) + z * dimxXdimy); - break; - case CGAL_DOUBLE : - sliceOut = (((r64*)bufferOut) + z * dimxXdimy); - break; - default : - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Error in %s: such output type not handled.\n", proc ); - free( theXX ); - return( EXIT_ON_FAILURE ); - } - ConvertBuffer( theXY, CGAL_FLOAT, sliceOut, typeOut, dimxXdimy ); - } - } - - return( EXIT_ON_SUCCESS ); -} - - - - - - - - - - - - - - - - - - - - - - -int GradientHessianGradient ( void *bufferIn, - bufferType typeIn, - void *bufferOut, - bufferType typeOut, - int *bufferDims, - int *borderLengths, - float *filterCoefs, - recursiveFilterType filterType ) -{ - const char *proc = "GradientHessianGradient"; - - - - float *theZZ = NULL; - float *theZ = NULL; - float *theZ1 = NULL; - float *theZ0 = NULL; - - float *theXZ = NULL; - float *theYZ = NULL; - - float *theXX = NULL; - float *theYY = NULL; - float *theXY = NULL; - - float *theX = NULL; - float *theY = NULL; - - - derivativeOrder ZZderiv[3] = { SMOOTHING, SMOOTHING, DERIVATIVE_2 }; - derivativeOrder Zderiv[3] = { SMOOTHING, SMOOTHING, DERIVATIVE_1 }; - derivativeOrder Z1deriv[3] = { NODERIVATIVE, NODERIVATIVE, DERIVATIVE_1 }; - derivativeOrder Z0deriv[3] = { NODERIVATIVE, NODERIVATIVE, SMOOTHING }; - - derivativeOrder XZderiv[3] = { DERIVATIVE_1, SMOOTHING, NODERIVATIVE }; - derivativeOrder YZderiv[3] = { SMOOTHING, DERIVATIVE_1, NODERIVATIVE }; - - derivativeOrder XXderiv[3] = { DERIVATIVE_2, SMOOTHING, NODERIVATIVE }; - derivativeOrder YYderiv[3] = { SMOOTHING, DERIVATIVE_2, NODERIVATIVE }; - derivativeOrder XYderiv[3] = { DERIVATIVE_1, DERIVATIVE_1, NODERIVATIVE }; - - derivativeOrder Xderiv[3] = { DERIVATIVE_1, SMOOTHING, NODERIVATIVE }; - derivativeOrder Yderiv[3] = { SMOOTHING, DERIVATIVE_1, NODERIVATIVE }; - - int sliceDims[3]; - int z, i, j, dimxXdimy; - - double gx, gy, gz, g; - - /* - * We check the buffers' dimensions. - */ - if ( bufferDims[2] == 1 ) { - return( GradientHessianGradient_2D ( bufferIn, typeIn, bufferOut, typeOut, - bufferDims, borderLengths, filterCoefs, filterType ) ); - } - - if ( (bufferDims[0] <= 0) || (bufferDims[1] <= 0) || (bufferDims[2] <= 0) ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc ); - return( EXIT_ON_FAILURE ); - } - - /* - * test of the coefficients - */ - if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) || - (filterCoefs[2] < 0.0) ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc ); - return( EXIT_ON_FAILURE ); - } - - - /* - * - */ - dimxXdimy = bufferDims[0] * bufferDims[1]; - sliceDims[0] = bufferDims[0]; - sliceDims[1] = bufferDims[1]; - sliceDims[2] = 1; - - - if ( typeOut == CGAL_FLOAT ) { - theX = (float*)malloc( (7+3*bufferDims[2]) * dimxXdimy * sizeof( float ) ); - } else { - theX = (float*)malloc( (7+4*bufferDims[2]) * dimxXdimy * sizeof( float ) ); - } - - - if ( theX == NULL ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to allocate auxiliary buffer.\n" ); - } - return( EXIT_ON_FAILURE ); - } - - /* - * BUFFERS - * - * slices : theX theY theXY theYY theXX theYZ theXZ - * - * volumes : theZ0 theZ1 theZ theZZ - * - */ - - theY = theXX = theXY = theYY = theYZ = theXZ = theX; - theZ0 = theZ1 = theZ = theX; - - theY += dimxXdimy; - theXY += 2*dimxXdimy; - theYY += 3*dimxXdimy; - theXX += 4*dimxXdimy; - theYZ += 5*dimxXdimy; - theXZ += 6*dimxXdimy; - - theZ0 += 7*dimxXdimy; - theZ1 += 7*dimxXdimy + bufferDims[2]*dimxXdimy; - theZ += 7*dimxXdimy + 2*bufferDims[2]*dimxXdimy; - - if ( typeOut == CGAL_FLOAT ) { - theZZ = (float*)bufferOut; - } else { - theZZ = theX; - theZZ += 7*dimxXdimy + 3*bufferDims[2]*dimxXdimy; - } - - - - /* - * - * 3D filtering / filtering along Z - * - */ - - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ0, CGAL_FLOAT, - bufferDims, borderLengths, - Z0deriv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Z^0 derivative.\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ1, CGAL_FLOAT, - bufferDims, borderLengths, - Z1deriv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Z^1 derivative.\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ, CGAL_FLOAT, - bufferDims, borderLengths, - Zderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Z^1 derivative (edge).\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZZ, CGAL_FLOAT, - bufferDims, borderLengths, - ZZderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Z^2 derivative.\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - - /* - * theZ0 : smoothed along Z - * theZ1 : first derivative along Z - * theZ : first derivative along Z, smoothed along X and Y - * theZZ : second derivative along Z, smoothed along X and Y - */ - - - - for ( z=0; z 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute X^1Z^1 derivative.\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( theZ1+z*dimxXdimy, CGAL_FLOAT, theYZ, CGAL_FLOAT, - sliceDims, borderLengths, - YZderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Y^1Z^1 derivative.\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - - - - if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, CGAL_FLOAT, theXX, CGAL_FLOAT, - sliceDims, borderLengths, - XXderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute X^2 derivative.\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, CGAL_FLOAT, theYY, CGAL_FLOAT, - sliceDims, borderLengths, - YYderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Y^2 derivative.\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, CGAL_FLOAT, theXY, CGAL_FLOAT, - sliceDims, borderLengths, - XYderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute X^1Y^1 derivative.\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - - - if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, CGAL_FLOAT, theX, CGAL_FLOAT, - sliceDims, borderLengths, - Xderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute X^1 derivative (edge).\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, CGAL_FLOAT, theY, CGAL_FLOAT, - sliceDims, borderLengths, - Yderiv, filterCoefs, filterType ) == 0 ) { - if ( _VERBOSE_ > 0 ) { - fprintf( stderr, " Fatal error in %s:", proc ); - fprintf( stderr, " unable to compute Y^1 derivative (edge).\n" ); - } - free( theX ); - return( EXIT_ON_FAILURE ); - } - - - - for ( j=z*dimxXdimy, i=0; i 1e-10 ) theZZ[j] = (float)(theZZ[j] / g); - - } - - } - - if ( typeOut != CGAL_FLOAT ) { - ConvertBuffer( theZZ, CGAL_FLOAT, bufferOut, typeOut, bufferDims[2]*dimxXdimy ); - } - - free( theX ); - - return( EXIT_ON_SUCCESS ); -} - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -/* - * - * - * - * - */ - -int RecursiveFilterOnBuffer( void *bufferIn, - bufferType typeIn, - void *bufferOut, - bufferType typeOut, - int *bufferDims, - int *borderLengths, - derivativeOrder *derivatives, - float *filterCoefs, - recursiveFilterType filterType ) -{ - const char *proc = "RecursiveFilterOnBuffer"; - register int dimx, dimxXdimy; - int dimy, dimz; - register int x, y, z; - /* - *obviously, we need to perform the computation - * with float or double values. For this reason, - * we allocate an auxiliary buffer if the output buffer - * is not of type float or double. - */ - void *bufferToBeProcessed = (void*)NULL; - bufferType typeToBeProcessed = TYPE_UNKNOWN; - void *bufferResult = (void*)NULL; - bufferType typeResult = TYPE_UNKNOWN; - /* - * lines' lengths - */ - int lengthX = 0; - int lengthY = 0; - int lengthZ = 0; - int maxLengthline = 0; - int borderXlength = 0; - int borderYlength = 0; - int borderZlength = 0; - /* - * 1D arrays for computations. - */ - double *theLine = (double*)NULL; - double *resLine = (double*)NULL; - double *tmpLine = (double*)NULL; - /* - * pointers for computations; - */ - register r32 *r32firstPoint = (r32*)NULL; - register r64 *r64firstPoint = (r64*)NULL; - register r32 *r32_pt = (r32*)NULL; - register r64 *r64_pt = (r64*)NULL; - register double *dbl_pt1 = (double*)NULL; - register double *dbl_pt2 = (double*)NULL; - register double dbl_first = 0.0; - register double dbl_last = 0.0; - int offsetLastPoint = 0; - int offsetNextFirstPoint = 0; - register r32 *r32firstPointResult = (r32*)NULL; - register r64 *r64firstPointResult = (r64*)NULL; - double *theLinePlusBorder = (double*)NULL; - double *resLinePlusBorder = (double*)NULL; - - RFcoefficientType *RFC = NULL; - - /* - * We check the buffers' dimensions. - */ - dimx = bufferDims[0]; dimy = bufferDims[1]; dimz = bufferDims[2]; - dimxXdimy = dimx * dimy; - if ( (dimx <= 0) || (dimy <= 0) || (dimz <= 0) ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc ); - return( EXIT_ON_FAILURE ); - } - /* - * We check the pointers. - */ - if ( (bufferIn == (void*)NULL) || (bufferOut == (void*)NULL) ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Fatal error in %s: NULL pointer on buffer.\n", proc ); - return( EXIT_ON_FAILURE ); - } - - /* - * May we use the buffer bufferOut as the bufferResult? - * If its type is CGAL_FLOAT or CGAL_DOUBLE, then yes. - * If not, we have to allocate an auxiliary buffer. - */ - if ( (typeOut == CGAL_FLOAT) || (typeOut == CGAL_DOUBLE) ) { - bufferResult = bufferOut; - typeResult = typeOut; - } else { - bufferResult = (void*)malloc( (dimx*dimy*dimz) * sizeof(r32) ); - if ( bufferResult == (void*)NULL ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Fatal error in %s: unable to allocate auxiliary buffer.\n", proc ); - return( EXIT_ON_FAILURE ); - } - typeResult = CGAL_FLOAT; - } - - /* - * May we consider the buffer bufferIn as the bufferToBeProcessed? - * If its type is CGAL_FLOAT or CGAL_DOUBLE, then yes. - * If not, we convert it into the buffer bufferResult, and this - * last buffer is the bufferToBeProcessed. - */ - if ( (typeIn == CGAL_FLOAT) || (typeIn == CGAL_DOUBLE) ) { - bufferToBeProcessed = bufferIn; - typeToBeProcessed = typeIn; - } else { - ConvertBuffer( bufferIn, typeIn, bufferResult, typeResult, (dimx*dimy*dimz) ); - bufferToBeProcessed = bufferResult; - typeToBeProcessed = typeResult; - } - - /* - * Estimation of the lines' length along each direction. - */ - if ( borderLengths != NULL ) { - borderXlength = borderLengths[0]; - borderYlength = borderLengths[1]; - borderZlength = borderLengths[2]; - if ( borderXlength < 0 ) borderXlength = 0; - if ( borderYlength < 0 ) borderYlength = 0; - if ( borderZlength < 0 ) borderZlength = 0; - } - - /* - * Tue Jul 6 19:15:15 MET DST 1999 (gregoire Malandain) - * changes 3 x dimx -> dimx, dimy, dimz - */ - lengthX = dimx + 2 * borderXlength; - lengthY = dimy + 2 * borderYlength; - lengthZ = dimz + 2 * borderZlength; - maxLengthline = lengthX; - if ( maxLengthline < lengthY ) maxLengthline = lengthY; - if ( maxLengthline < lengthZ ) maxLengthline = lengthZ; - if ( maxLengthline <= 0 ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Error in %s: unable to deal with dimensions = 0.\n", proc ); - if ( (typeOut != CGAL_FLOAT) && (typeOut != CGAL_DOUBLE) ) - free( bufferResult ); - return( EXIT_ON_FAILURE ); - } - /* - * Allocations of work arrays. - * We will use them to process each line. - */ - theLine = (double*)malloc( 3 * maxLengthline * sizeof(double) ); - if ( theLine == (double*)NULL ) { - if ( _VERBOSE_ > 0 ) - fprintf( stderr, " Fatal error in %s: unable to allocate auxiliary work arrays.\n", proc ); - if ( (typeOut != CGAL_FLOAT) && (typeOut != CGAL_DOUBLE) ) - free( bufferResult ); - return( EXIT_ON_FAILURE ); - } - resLine = theLine + maxLengthline; - tmpLine = resLine + maxLengthline; - - /* - * From now, - * typeToBeProcessed is either CGAL_FLOAT or CGAL_DOUBLE - * so is typeResult. - */ - - - /* - * Processing along X. - */ - if ( dimx > 4 ) - if (derivatives[0] != NODERIVATIVE) - if (filterCoefs[0] > 0.0) { - - if ( _VERBOSE_ != 0 ) - fprintf( stderr, " %s: processing along X.\n", proc ); - - RFC = InitRecursiveCoefficients( (double)filterCoefs[0], filterType, derivatives[0] ); - - if ( RFC == NULL ) { - if ( _VERBOSE_ != 0 ) - fprintf( stderr, " %s: unable to allocate coefficients\n", proc ); - if ( (typeOut != CGAL_FLOAT) && (typeOut != CGAL_DOUBLE) ) - free( bufferResult ); - return( EXIT_ON_FAILURE ); - } - - r64firstPoint = (r64*)bufferToBeProcessed; - r32firstPoint = (r32*)bufferToBeProcessed; - - r64firstPointResult = (r64*)bufferResult; - r32firstPointResult = (r32*)bufferResult; - - offsetLastPoint = borderXlength + dimx - 1; - - theLinePlusBorder = theLine + borderXlength; - resLinePlusBorder = resLine + borderXlength; - - /* - * There are dimz*dimy X lines to be processed. - */ - for ( z=0; z 0 ) { - dbl_pt1 = theLine + borderXlength; dbl_first = *dbl_pt1; - dbl_pt2 = theLine + offsetLastPoint; dbl_last = *dbl_pt2; - for ( x=0; x 4 ) - if (derivatives[1] != NODERIVATIVE) - if (filterCoefs[1] > 0.0) { - - if ( _VERBOSE_ != 0 ) - fprintf( stderr, " %s: processing along Y.\n", proc ); - - RFC = InitRecursiveCoefficients( (double)filterCoefs[1], filterType, derivatives[1] ); - - if ( RFC == NULL ) { - if ( _VERBOSE_ != 0 ) - fprintf( stderr, " %s: unable to allocate coefficients\n", proc ); - if ( (typeOut != CGAL_FLOAT) && (typeOut != CGAL_DOUBLE) ) - free( bufferResult ); - return( EXIT_ON_FAILURE ); - } - - r64firstPoint = (r64*)bufferToBeProcessed; - r32firstPoint = (r32*)bufferToBeProcessed; - - r64firstPointResult = (r64*)bufferResult; - r32firstPointResult = (r32*)bufferResult; - - offsetLastPoint = borderYlength + dimy - 1; - offsetNextFirstPoint = dimx * dimy - dimx; - - theLinePlusBorder = theLine + borderYlength; - resLinePlusBorder = resLine + borderYlength; - - /* - * There are dimz*dimx Y lines to be processed. - */ - for ( z=0; z 0 ) { - dbl_pt1 = theLine + borderYlength; dbl_first = *dbl_pt1; - dbl_pt2 = theLine + offsetLastPoint; dbl_last = *dbl_pt2; - for ( y=0; y 4 ) - if (derivatives[2] != NODERIVATIVE) - if (filterCoefs[2] > 0.0) { - - if ( _VERBOSE_ != 0 ) - fprintf( stderr, " %s: processing along Z.\n", proc ); - - RFC = InitRecursiveCoefficients( (double)filterCoefs[2], filterType, derivatives[2] ); - - if ( RFC == NULL ) { - if ( _VERBOSE_ != 0 ) - fprintf( stderr, " %s: unable to allocate coefficients\n", proc ); - if ( (typeOut != CGAL_FLOAT) && (typeOut != CGAL_DOUBLE) ) - free( bufferResult ); - return( EXIT_ON_FAILURE ); - } - - r64firstPoint = (r64*)bufferToBeProcessed; - r32firstPoint = (r32*)bufferToBeProcessed; - - offsetLastPoint = borderZlength + dimz - 1; - - r64firstPointResult = (r64*)bufferResult; - r32firstPointResult = (r32*)bufferResult; - - offsetLastPoint = borderZlength + dimz - 1; - - theLinePlusBorder = theLine + borderYlength; - resLinePlusBorder = resLine + borderYlength; - - /* - * There are dimy*dimx Z lines to be processed. - */ - for ( y=0; y 0 ) { - dbl_pt1 = theLine + borderZlength; dbl_first = *dbl_pt1; - dbl_pt2 = theLine + offsetLastPoint; dbl_last = *dbl_pt2; - for ( z=0; z +#include +#include + +#include "convert.h" +#include "recbuffer.h" + +#ifdef CGAL_HEADER_ONLY + +inline int& get_static_verbose_recbuffer() +{ + static int _VERBOSE_ = 0; + return _VERBOSE_; +} + +#else // CGAL_HEADER_ONLY + +static int _VERBOSE_ = 0; + +inline int& get_static_verbose_recbuffer() +{ return _VERBOSE_; } + +#endif // CGAL_HEADER_ONLY + + +#define EXIT_ON_FAILURE 0 +#define EXIT_ON_SUCCESS 1 + +/* + * + * Gradient modulus + * + * + */ +CGAL_INLINE_FUNCTION +int GradientModulus( void *bufferIn, + bufferType typeIn, + void *bufferOut, + bufferType typeOut, + int *bufferDims, + int *borderLengths, + float *filterCoefs, + recursiveFilterType filterType ) +{ + const char *proc = "GradientModulus"; + float *auxBuf = NULL; + float *tmpBuf = NULL, *grdBuf = NULL; + int sizeAuxBuf = 0; + derivativeOrder derivatives[3]; + int i; + + + sizeAuxBuf = bufferDims[0] * bufferDims[1] * bufferDims[2]; + if ( typeOut != CGAL_FLOAT || bufferIn == bufferOut ) + sizeAuxBuf *= 2; + + + /* allocation des buffers de calcul + */ + auxBuf = (float*)malloc( sizeAuxBuf * sizeof(float) ); + if ( auxBuf == NULL ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, "%s: unable to allocate auxiliary buffer\n", proc ); + return( EXIT_ON_FAILURE ); + } + tmpBuf = auxBuf; + if ( typeOut != CGAL_FLOAT || bufferIn == bufferOut ) { + grdBuf = tmpBuf; + grdBuf += bufferDims[0] * bufferDims[1] * bufferDims[2]; + } else { + grdBuf = (float*)bufferOut; + } + + /* cas 2D + */ + if ( bufferDims[2] == 1 ) { + + derivatives[0] = DERIVATIVE_1; + derivatives[1] = DERIVATIVE_0; + derivatives[2] = NODERIVATIVE; + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)grdBuf, CGAL_FLOAT, + bufferDims, borderLengths, derivatives, + filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { + if ( _VERBOSE_ ) + fprintf( stderr, "%s: unable to compute X derivative (2D)\n", proc ); + free( auxBuf ); + return( EXIT_ON_FAILURE ); + } + + derivatives[0] = DERIVATIVE_0; + derivatives[1] = DERIVATIVE_1; + derivatives[2] = NODERIVATIVE; + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)tmpBuf, CGAL_FLOAT, + bufferDims, borderLengths, derivatives, + filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { + if ( _VERBOSE_ ) + fprintf( stderr, "%s: unable to compute Y derivative (2D)\n", proc ); + free( auxBuf ); + return( EXIT_ON_FAILURE ); + } + + sizeAuxBuf = bufferDims[0] * bufferDims[1] * bufferDims[2]; + for ( i = 0; i < sizeAuxBuf; i++ ) + grdBuf[i] = (float)sqrt( grdBuf[i]*grdBuf[i] + tmpBuf[i]*tmpBuf[i] ); + + } else { + + derivatives[0] = NODERIVATIVE; + derivatives[1] = NODERIVATIVE; + derivatives[2] = DERIVATIVE_0; + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)tmpBuf, CGAL_FLOAT, + bufferDims, borderLengths, derivatives, + filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { + if ( _VERBOSE_ ) + fprintf( stderr, "%s: unable to compute Z smoothing (3D)\n", proc ); + free( auxBuf ); + return( EXIT_ON_FAILURE ); + } + + derivatives[0] = DERIVATIVE_1; + derivatives[1] = DERIVATIVE_0; + derivatives[2] = NODERIVATIVE; + if ( RecursiveFilterOnBuffer( (void*)tmpBuf, CGAL_FLOAT, (void*)grdBuf, CGAL_FLOAT, + bufferDims, borderLengths, derivatives, + filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { + if ( _VERBOSE_ ) + fprintf( stderr, "%s: unable to compute X derivative (3D)\n", proc ); + free( auxBuf ); + return( EXIT_ON_FAILURE ); + } + + derivatives[0] = DERIVATIVE_0; + derivatives[1] = DERIVATIVE_1; + derivatives[2] = NODERIVATIVE; + if ( RecursiveFilterOnBuffer( (void*)tmpBuf, CGAL_FLOAT, (void*)tmpBuf, CGAL_FLOAT, + bufferDims, borderLengths, derivatives, + filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { + if ( _VERBOSE_ ) + fprintf( stderr, "%s: unable to compute Y derivative (3D)\n", proc ); + free( auxBuf ); + return( EXIT_ON_FAILURE ); + } + + sizeAuxBuf = bufferDims[0] * bufferDims[1] * bufferDims[2]; + for ( i = 0; i < sizeAuxBuf; i++ ) + grdBuf[i] = grdBuf[i]*grdBuf[i] + tmpBuf[i]*tmpBuf[i]; + + derivatives[0] = DERIVATIVE_0; + derivatives[1] = DERIVATIVE_0; + derivatives[2] = DERIVATIVE_1; + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)tmpBuf, CGAL_FLOAT, + bufferDims, borderLengths, derivatives, + filterCoefs, filterType ) != EXIT_ON_SUCCESS ) { + if ( _VERBOSE_ ) + fprintf( stderr, "%s: unable to compute Z derivative (3D)\n", proc ); + free( auxBuf ); + return( EXIT_ON_FAILURE ); + } + + for ( i = 0; i < sizeAuxBuf; i++ ) + grdBuf[i] = (float)sqrt( grdBuf[i] + tmpBuf[i]*tmpBuf[i] ); + + } + + if ( grdBuf != bufferOut ) + ConvertBuffer( grdBuf, CGAL_FLOAT, bufferOut, typeOut, + bufferDims[0]*bufferDims[1]*bufferDims[2] ); + free( auxBuf ); + return( EXIT_ON_SUCCESS ); +} + + + + + + + + + + + + + + + + + + + + + + + + +/* + * + * Laplacian + * + * + */ +CGAL_INLINE_FUNCTION +int Laplacian_2D ( void *bufferIn, + bufferType typeIn, + void *bufferOut, + bufferType typeOut, + int *bufferDims, + int *borderLengths, + float *filterCoefs, + recursiveFilterType filterType ) +{ + const char *proc = "Laplacian_2D"; + float *theXX = NULL; + float *theYY = NULL; + + derivativeOrder XXderiv[3] = { DERIVATIVE_2, SMOOTHING, NODERIVATIVE }; + derivativeOrder YYderiv[3] = { SMOOTHING, DERIVATIVE_2, NODERIVATIVE }; + int sliceDims[3]; + int z, i, dimxXdimy; + + void *sliceOut = NULL; + + + + /* + * We check the buffers' dimensions. + */ + if ( (bufferDims[0] <= 0) || (bufferDims[1] <= 0) || (bufferDims[2] <= 0) ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc ); + return( EXIT_ON_FAILURE ); + } + + /* + * test of the coefficients + */ + if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) || + (filterCoefs[2] < 0.0) ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc ); + return( EXIT_ON_FAILURE ); + } + + + /* + * + */ + dimxXdimy = bufferDims[0] * bufferDims[1]; + sliceDims[0] = bufferDims[0]; + sliceDims[1] = bufferDims[1]; + sliceDims[2] = 1; + + + if ( typeOut == CGAL_FLOAT ) { + theXX = (float*)malloc( dimxXdimy * sizeof( float ) ); + } else { + theXX = (float*)malloc( 2 * dimxXdimy * sizeof( float ) ); + } + + if ( theXX == NULL ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to allocate auxiliary buffer.\n" ); + } + return( EXIT_ON_FAILURE ); + } + + if ( typeOut != CGAL_FLOAT ) { + theYY = theXX; + theYY += dimxXdimy; + } + + + + for ( z=0; z 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute X^2 derivative.\n" ); + } + free( theXX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theYY, CGAL_FLOAT, + sliceDims, borderLengths, + YYderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Y^2 derivative.\n" ); + } + free( theXX ); + return( EXIT_ON_FAILURE ); + } + + + for ( i=0; i 0 ) + fprintf( stderr, " Error in %s: such output type not handled.\n", proc ); + free( theXX ); + return( EXIT_ON_FAILURE ); + } + ConvertBuffer( theYY, CGAL_FLOAT, sliceOut, typeOut, dimxXdimy ); + } + } + + return( EXIT_ON_SUCCESS ); +} + + + + + + + + + + + + + +CGAL_INLINE_FUNCTION +int Laplacian ( void *bufferIn, + bufferType typeIn, + void *bufferOut, + bufferType typeOut, + int *bufferDims, + int *borderLengths, + float *filterCoefs, + recursiveFilterType filterType ) +{ + const char *proc = "Laplacian"; + float *theSL = NULL; + float *theZZ = NULL; + float *theZ0 = NULL; + + + derivativeOrder XXderiv[3] = { DERIVATIVE_2, SMOOTHING, NODERIVATIVE }; + derivativeOrder YYderiv[3] = { SMOOTHING, DERIVATIVE_2, NODERIVATIVE }; + derivativeOrder Zsmooth[3] = { NODERIVATIVE, NODERIVATIVE, SMOOTHING }; + derivativeOrder ZZderiv[3] = { SMOOTHING, SMOOTHING, DERIVATIVE_2 }; + + int sliceDims[3]; + int z, i, j, dimxXdimy; + + + + + /* + * We check the buffers' dimensions. + */ + if ( bufferDims[2] == 1 ) { + return( Laplacian_2D ( bufferIn, typeIn, bufferOut, typeOut, + bufferDims, borderLengths, filterCoefs, filterType ) ); + } + + if ( (bufferDims[0] <= 0) || (bufferDims[1] <= 0) || (bufferDims[2] <= 0) ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc ); + return( EXIT_ON_FAILURE ); + } + + /* + * test of the coefficients + */ + if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) || + (filterCoefs[2] < 0.0) ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc ); + return( EXIT_ON_FAILURE ); + } + + + /* + * + */ + dimxXdimy = bufferDims[0] * bufferDims[1]; + sliceDims[0] = bufferDims[0]; + sliceDims[1] = bufferDims[1]; + sliceDims[2] = 1; + + + + if ( typeOut == CGAL_FLOAT ) { + theSL = (float*)malloc( (1+bufferDims[2]) * dimxXdimy * sizeof( float ) ); + } else { + theSL = (float*)malloc( (1+2*bufferDims[2]) * dimxXdimy * sizeof( float ) ); + } + + + + if ( theSL == NULL ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to allocate auxiliary buffer.\n" ); + } + return( EXIT_ON_FAILURE ); + } + + + + theZ0 = theSL; + theZ0 += dimxXdimy; + + + + if ( typeOut == CGAL_FLOAT ) { + theZZ = (float*) bufferOut; + } else { + theZZ = theZ0; + theZZ += dimxXdimy * bufferDims[2]; + } + + + + /* + * + * 3D filtering / filtering along Z + * + */ + + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ0, CGAL_FLOAT, + bufferDims, borderLengths, + Zsmooth, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Z^0 derivative.\n" ); + } + free( theSL ); + return( EXIT_ON_FAILURE ); + } + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZZ, CGAL_FLOAT, + bufferDims, borderLengths, + ZZderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Z^2 derivative.\n" ); + } + free( theSL ); + return( EXIT_ON_FAILURE ); + } + + + + + + + for ( z=0; z 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute X^2 derivative.\n" ); + } + free( theSL ); + return( EXIT_ON_FAILURE ); + } + + for ( j=z*dimxXdimy, i=0; i 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Y^2 derivative.\n" ); + } + free( theSL ); + return( EXIT_ON_FAILURE ); + } + + for ( j=z*dimxXdimy, i=0; i 0 ) + fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc ); + return( EXIT_ON_FAILURE ); + } + + /* + * test of the coefficients + */ + if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) || + (filterCoefs[2] < 0.0) ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc ); + return( EXIT_ON_FAILURE ); + } + + + /* + * + */ + dimxXdimy = bufferDims[0] * bufferDims[1]; + sliceDims[0] = bufferDims[0]; + sliceDims[1] = bufferDims[1]; + sliceDims[2] = 1; + + + if ( typeOut == CGAL_FLOAT ) { + theXX = (float*)malloc( 4 * dimxXdimy * sizeof( float ) ); + } else { + theXX = (float*)malloc( 5 * dimxXdimy * sizeof( float ) ); + } + + if ( theXX == NULL ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to allocate auxiliary buffer.\n" ); + } + return( EXIT_ON_FAILURE ); + } + + + + theX = theY = theYY = theXX; + theYY += dimxXdimy; + theX += 2*dimxXdimy; + theY += 3*dimxXdimy; + + + + if ( typeOut != CGAL_FLOAT ) { + theXY = theXX; + theXY += 4*dimxXdimy; + } + + + + for ( z=0; z 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Y^0 derivative.\n" ); + } + free( theXX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( sliceIn, typeIn, theY, CGAL_FLOAT, + sliceDims, borderLengths, + Xsmooth, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute X^0 derivative.\n" ); + } + free( theXX ); + return( EXIT_ON_FAILURE ); + } + + + + + if ( RecursiveFilterOnBuffer( sliceIn, typeIn, theXY, CGAL_FLOAT, + sliceDims, borderLengths, + XYderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute X^1Y^1 derivative.\n" ); + } + free( theXX ); + return( EXIT_ON_FAILURE ); + } + + + + + if ( RecursiveFilterOnBuffer( theX, CGAL_FLOAT, theXX, CGAL_FLOAT, + sliceDims, borderLengths, + XXderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute X^2 derivative.\n" ); + } + free( theXX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( theY, CGAL_FLOAT, theYY, CGAL_FLOAT, + sliceDims, borderLengths, + YYderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Y^2 derivative.\n" ); + } + free( theXX ); + return( EXIT_ON_FAILURE ); + } + + + + + if ( RecursiveFilterOnBuffer( theX, CGAL_FLOAT, theX, CGAL_FLOAT, + sliceDims, borderLengths, + Xderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute X^1 derivative.\n" ); + } + free( theXX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( theY, CGAL_FLOAT, theY, CGAL_FLOAT, + sliceDims, borderLengths, + Yderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Y^1 derivative.\n" ); + } + free( theXX ); + return( EXIT_ON_FAILURE ); + } + + + + + for ( i=0; i 1e-10 ) theXY[i] = (float)(theXY[i] / g); + } + + if ( typeOut != CGAL_FLOAT ) { + switch ( typeOut ) { + case CGAL_UCHAR : + sliceOut = (((u8*)bufferOut) + z * dimxXdimy); + break; + case CGAL_SCHAR : + sliceOut = (((s8*)bufferOut) + z * dimxXdimy); + break; + case CGAL_SSHORT : + sliceOut = (((s16*)bufferOut) + z * dimxXdimy); + break; + case CGAL_DOUBLE : + sliceOut = (((r64*)bufferOut) + z * dimxXdimy); + break; + default : + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Error in %s: such output type not handled.\n", proc ); + free( theXX ); + return( EXIT_ON_FAILURE ); + } + ConvertBuffer( theXY, CGAL_FLOAT, sliceOut, typeOut, dimxXdimy ); + } + } + + return( EXIT_ON_SUCCESS ); +} + + + + + + + + + + + + + + + + + + + + + +CGAL_INLINE_FUNCTION +int GradientHessianGradient ( void *bufferIn, + bufferType typeIn, + void *bufferOut, + bufferType typeOut, + int *bufferDims, + int *borderLengths, + float *filterCoefs, + recursiveFilterType filterType ) +{ + const char *proc = "GradientHessianGradient"; + + + + float *theZZ = NULL; + float *theZ = NULL; + float *theZ1 = NULL; + float *theZ0 = NULL; + + float *theXZ = NULL; + float *theYZ = NULL; + + float *theXX = NULL; + float *theYY = NULL; + float *theXY = NULL; + + float *theX = NULL; + float *theY = NULL; + + + derivativeOrder ZZderiv[3] = { SMOOTHING, SMOOTHING, DERIVATIVE_2 }; + derivativeOrder Zderiv[3] = { SMOOTHING, SMOOTHING, DERIVATIVE_1 }; + derivativeOrder Z1deriv[3] = { NODERIVATIVE, NODERIVATIVE, DERIVATIVE_1 }; + derivativeOrder Z0deriv[3] = { NODERIVATIVE, NODERIVATIVE, SMOOTHING }; + + derivativeOrder XZderiv[3] = { DERIVATIVE_1, SMOOTHING, NODERIVATIVE }; + derivativeOrder YZderiv[3] = { SMOOTHING, DERIVATIVE_1, NODERIVATIVE }; + + derivativeOrder XXderiv[3] = { DERIVATIVE_2, SMOOTHING, NODERIVATIVE }; + derivativeOrder YYderiv[3] = { SMOOTHING, DERIVATIVE_2, NODERIVATIVE }; + derivativeOrder XYderiv[3] = { DERIVATIVE_1, DERIVATIVE_1, NODERIVATIVE }; + + derivativeOrder Xderiv[3] = { DERIVATIVE_1, SMOOTHING, NODERIVATIVE }; + derivativeOrder Yderiv[3] = { SMOOTHING, DERIVATIVE_1, NODERIVATIVE }; + + int sliceDims[3]; + int z, i, j, dimxXdimy; + + double gx, gy, gz, g; + + /* + * We check the buffers' dimensions. + */ + if ( bufferDims[2] == 1 ) { + return( GradientHessianGradient_2D ( bufferIn, typeIn, bufferOut, typeOut, + bufferDims, borderLengths, filterCoefs, filterType ) ); + } + + if ( (bufferDims[0] <= 0) || (bufferDims[1] <= 0) || (bufferDims[2] <= 0) ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc ); + return( EXIT_ON_FAILURE ); + } + + /* + * test of the coefficients + */ + if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) || + (filterCoefs[2] < 0.0) ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc ); + return( EXIT_ON_FAILURE ); + } + + + /* + * + */ + dimxXdimy = bufferDims[0] * bufferDims[1]; + sliceDims[0] = bufferDims[0]; + sliceDims[1] = bufferDims[1]; + sliceDims[2] = 1; + + + if ( typeOut == CGAL_FLOAT ) { + theX = (float*)malloc( (7+3*bufferDims[2]) * dimxXdimy * sizeof( float ) ); + } else { + theX = (float*)malloc( (7+4*bufferDims[2]) * dimxXdimy * sizeof( float ) ); + } + + + if ( theX == NULL ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to allocate auxiliary buffer.\n" ); + } + return( EXIT_ON_FAILURE ); + } + + /* + * BUFFERS + * + * slices : theX theY theXY theYY theXX theYZ theXZ + * + * volumes : theZ0 theZ1 theZ theZZ + * + */ + + theY = theXX = theXY = theYY = theYZ = theXZ = theX; + theZ0 = theZ1 = theZ = theX; + + theY += dimxXdimy; + theXY += 2*dimxXdimy; + theYY += 3*dimxXdimy; + theXX += 4*dimxXdimy; + theYZ += 5*dimxXdimy; + theXZ += 6*dimxXdimy; + + theZ0 += 7*dimxXdimy; + theZ1 += 7*dimxXdimy + bufferDims[2]*dimxXdimy; + theZ += 7*dimxXdimy + 2*bufferDims[2]*dimxXdimy; + + if ( typeOut == CGAL_FLOAT ) { + theZZ = (float*)bufferOut; + } else { + theZZ = theX; + theZZ += 7*dimxXdimy + 3*bufferDims[2]*dimxXdimy; + } + + + + /* + * + * 3D filtering / filtering along Z + * + */ + + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ0, CGAL_FLOAT, + bufferDims, borderLengths, + Z0deriv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Z^0 derivative.\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ1, CGAL_FLOAT, + bufferDims, borderLengths, + Z1deriv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Z^1 derivative.\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ, CGAL_FLOAT, + bufferDims, borderLengths, + Zderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Z^1 derivative (edge).\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZZ, CGAL_FLOAT, + bufferDims, borderLengths, + ZZderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Z^2 derivative.\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + + /* + * theZ0 : smoothed along Z + * theZ1 : first derivative along Z + * theZ : first derivative along Z, smoothed along X and Y + * theZZ : second derivative along Z, smoothed along X and Y + */ + + + + for ( z=0; z 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute X^1Z^1 derivative.\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( theZ1+z*dimxXdimy, CGAL_FLOAT, theYZ, CGAL_FLOAT, + sliceDims, borderLengths, + YZderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Y^1Z^1 derivative.\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + + + + if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, CGAL_FLOAT, theXX, CGAL_FLOAT, + sliceDims, borderLengths, + XXderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute X^2 derivative.\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, CGAL_FLOAT, theYY, CGAL_FLOAT, + sliceDims, borderLengths, + YYderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Y^2 derivative.\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, CGAL_FLOAT, theXY, CGAL_FLOAT, + sliceDims, borderLengths, + XYderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute X^1Y^1 derivative.\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + + + if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, CGAL_FLOAT, theX, CGAL_FLOAT, + sliceDims, borderLengths, + Xderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute X^1 derivative (edge).\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, CGAL_FLOAT, theY, CGAL_FLOAT, + sliceDims, borderLengths, + Yderiv, filterCoefs, filterType ) == 0 ) { + if ( _VERBOSE_ > 0 ) { + fprintf( stderr, " Fatal error in %s:", proc ); + fprintf( stderr, " unable to compute Y^1 derivative (edge).\n" ); + } + free( theX ); + return( EXIT_ON_FAILURE ); + } + + + + for ( j=z*dimxXdimy, i=0; i 1e-10 ) theZZ[j] = (float)(theZZ[j] / g); + + } + + } + + if ( typeOut != CGAL_FLOAT ) { + ConvertBuffer( theZZ, CGAL_FLOAT, bufferOut, typeOut, bufferDims[2]*dimxXdimy ); + } + + free( theX ); + + return( EXIT_ON_SUCCESS ); +} + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +/* + * + * + * + * + */ +CGAL_INLINE_FUNCTION +int RecursiveFilterOnBuffer( void *bufferIn, + bufferType typeIn, + void *bufferOut, + bufferType typeOut, + int *bufferDims, + int *borderLengths, + derivativeOrder *derivatives, + float *filterCoefs, + recursiveFilterType filterType ) +{ + const char *proc = "RecursiveFilterOnBuffer"; + register int dimx, dimxXdimy; + int dimy, dimz; + register int x, y, z; + /* + *obviously, we need to perform the computation + * with float or double values. For this reason, + * we allocate an auxiliary buffer if the output buffer + * is not of type float or double. + */ + void *bufferToBeProcessed = (void*)NULL; + bufferType typeToBeProcessed = TYPE_UNKNOWN; + void *bufferResult = (void*)NULL; + bufferType typeResult = TYPE_UNKNOWN; + /* + * lines' lengths + */ + int lengthX = 0; + int lengthY = 0; + int lengthZ = 0; + int maxLengthline = 0; + int borderXlength = 0; + int borderYlength = 0; + int borderZlength = 0; + /* + * 1D arrays for computations. + */ + double *theLine = (double*)NULL; + double *resLine = (double*)NULL; + double *tmpLine = (double*)NULL; + /* + * pointers for computations; + */ + register r32 *r32firstPoint = (r32*)NULL; + register r64 *r64firstPoint = (r64*)NULL; + register r32 *r32_pt = (r32*)NULL; + register r64 *r64_pt = (r64*)NULL; + register double *dbl_pt1 = (double*)NULL; + register double *dbl_pt2 = (double*)NULL; + register double dbl_first = 0.0; + register double dbl_last = 0.0; + int offsetLastPoint = 0; + int offsetNextFirstPoint = 0; + register r32 *r32firstPointResult = (r32*)NULL; + register r64 *r64firstPointResult = (r64*)NULL; + double *theLinePlusBorder = (double*)NULL; + double *resLinePlusBorder = (double*)NULL; + + RFcoefficientType *RFC = NULL; + + /* + * We check the buffers' dimensions. + */ + dimx = bufferDims[0]; dimy = bufferDims[1]; dimz = bufferDims[2]; + dimxXdimy = dimx * dimy; + if ( (dimx <= 0) || (dimy <= 0) || (dimz <= 0) ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc ); + return( EXIT_ON_FAILURE ); + } + /* + * We check the pointers. + */ + if ( (bufferIn == (void*)NULL) || (bufferOut == (void*)NULL) ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Fatal error in %s: NULL pointer on buffer.\n", proc ); + return( EXIT_ON_FAILURE ); + } + + /* + * May we use the buffer bufferOut as the bufferResult? + * If its type is CGAL_FLOAT or CGAL_DOUBLE, then yes. + * If not, we have to allocate an auxiliary buffer. + */ + if ( (typeOut == CGAL_FLOAT) || (typeOut == CGAL_DOUBLE) ) { + bufferResult = bufferOut; + typeResult = typeOut; + } else { + bufferResult = (void*)malloc( (dimx*dimy*dimz) * sizeof(r32) ); + if ( bufferResult == (void*)NULL ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Fatal error in %s: unable to allocate auxiliary buffer.\n", proc ); + return( EXIT_ON_FAILURE ); + } + typeResult = CGAL_FLOAT; + } + + /* + * May we consider the buffer bufferIn as the bufferToBeProcessed? + * If its type is CGAL_FLOAT or CGAL_DOUBLE, then yes. + * If not, we convert it into the buffer bufferResult, and this + * last buffer is the bufferToBeProcessed. + */ + if ( (typeIn == CGAL_FLOAT) || (typeIn == CGAL_DOUBLE) ) { + bufferToBeProcessed = bufferIn; + typeToBeProcessed = typeIn; + } else { + ConvertBuffer( bufferIn, typeIn, bufferResult, typeResult, (dimx*dimy*dimz) ); + bufferToBeProcessed = bufferResult; + typeToBeProcessed = typeResult; + } + + /* + * Estimation of the lines' length along each direction. + */ + if ( borderLengths != NULL ) { + borderXlength = borderLengths[0]; + borderYlength = borderLengths[1]; + borderZlength = borderLengths[2]; + if ( borderXlength < 0 ) borderXlength = 0; + if ( borderYlength < 0 ) borderYlength = 0; + if ( borderZlength < 0 ) borderZlength = 0; + } + + /* + * Tue Jul 6 19:15:15 MET DST 1999 (gregoire Malandain) + * changes 3 x dimx -> dimx, dimy, dimz + */ + lengthX = dimx + 2 * borderXlength; + lengthY = dimy + 2 * borderYlength; + lengthZ = dimz + 2 * borderZlength; + maxLengthline = lengthX; + if ( maxLengthline < lengthY ) maxLengthline = lengthY; + if ( maxLengthline < lengthZ ) maxLengthline = lengthZ; + if ( maxLengthline <= 0 ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Error in %s: unable to deal with dimensions = 0.\n", proc ); + if ( (typeOut != CGAL_FLOAT) && (typeOut != CGAL_DOUBLE) ) + free( bufferResult ); + return( EXIT_ON_FAILURE ); + } + /* + * Allocations of work arrays. + * We will use them to process each line. + */ + theLine = (double*)malloc( 3 * maxLengthline * sizeof(double) ); + if ( theLine == (double*)NULL ) { + if ( _VERBOSE_ > 0 ) + fprintf( stderr, " Fatal error in %s: unable to allocate auxiliary work arrays.\n", proc ); + if ( (typeOut != CGAL_FLOAT) && (typeOut != CGAL_DOUBLE) ) + free( bufferResult ); + return( EXIT_ON_FAILURE ); + } + resLine = theLine + maxLengthline; + tmpLine = resLine + maxLengthline; + + /* + * From now, + * typeToBeProcessed is either CGAL_FLOAT or CGAL_DOUBLE + * so is typeResult. + */ + + + /* + * Processing along X. + */ + if ( dimx > 4 ) + if (derivatives[0] != NODERIVATIVE) + if (filterCoefs[0] > 0.0) { + + if ( _VERBOSE_ != 0 ) + fprintf( stderr, " %s: processing along X.\n", proc ); + + RFC = InitRecursiveCoefficients( (double)filterCoefs[0], filterType, derivatives[0] ); + + if ( RFC == NULL ) { + if ( _VERBOSE_ != 0 ) + fprintf( stderr, " %s: unable to allocate coefficients\n", proc ); + if ( (typeOut != CGAL_FLOAT) && (typeOut != CGAL_DOUBLE) ) + free( bufferResult ); + return( EXIT_ON_FAILURE ); + } + + r64firstPoint = (r64*)bufferToBeProcessed; + r32firstPoint = (r32*)bufferToBeProcessed; + + r64firstPointResult = (r64*)bufferResult; + r32firstPointResult = (r32*)bufferResult; + + offsetLastPoint = borderXlength + dimx - 1; + + theLinePlusBorder = theLine + borderXlength; + resLinePlusBorder = resLine + borderXlength; + + /* + * There are dimz*dimy X lines to be processed. + */ + for ( z=0; z 0 ) { + dbl_pt1 = theLine + borderXlength; dbl_first = *dbl_pt1; + dbl_pt2 = theLine + offsetLastPoint; dbl_last = *dbl_pt2; + for ( x=0; x 4 ) + if (derivatives[1] != NODERIVATIVE) + if (filterCoefs[1] > 0.0) { + + if ( _VERBOSE_ != 0 ) + fprintf( stderr, " %s: processing along Y.\n", proc ); + + RFC = InitRecursiveCoefficients( (double)filterCoefs[1], filterType, derivatives[1] ); + + if ( RFC == NULL ) { + if ( _VERBOSE_ != 0 ) + fprintf( stderr, " %s: unable to allocate coefficients\n", proc ); + if ( (typeOut != CGAL_FLOAT) && (typeOut != CGAL_DOUBLE) ) + free( bufferResult ); + return( EXIT_ON_FAILURE ); + } + + r64firstPoint = (r64*)bufferToBeProcessed; + r32firstPoint = (r32*)bufferToBeProcessed; + + r64firstPointResult = (r64*)bufferResult; + r32firstPointResult = (r32*)bufferResult; + + offsetLastPoint = borderYlength + dimy - 1; + offsetNextFirstPoint = dimx * dimy - dimx; + + theLinePlusBorder = theLine + borderYlength; + resLinePlusBorder = resLine + borderYlength; + + /* + * There are dimz*dimx Y lines to be processed. + */ + for ( z=0; z 0 ) { + dbl_pt1 = theLine + borderYlength; dbl_first = *dbl_pt1; + dbl_pt2 = theLine + offsetLastPoint; dbl_last = *dbl_pt2; + for ( y=0; y 4 ) + if (derivatives[2] != NODERIVATIVE) + if (filterCoefs[2] > 0.0) { + + if ( _VERBOSE_ != 0 ) + fprintf( stderr, " %s: processing along Z.\n", proc ); + + RFC = InitRecursiveCoefficients( (double)filterCoefs[2], filterType, derivatives[2] ); + + if ( RFC == NULL ) { + if ( _VERBOSE_ != 0 ) + fprintf( stderr, " %s: unable to allocate coefficients\n", proc ); + if ( (typeOut != CGAL_FLOAT) && (typeOut != CGAL_DOUBLE) ) + free( bufferResult ); + return( EXIT_ON_FAILURE ); + } + + r64firstPoint = (r64*)bufferToBeProcessed; + r32firstPoint = (r32*)bufferToBeProcessed; + + offsetLastPoint = borderZlength + dimz - 1; + + r64firstPointResult = (r64*)bufferResult; + r32firstPointResult = (r32*)bufferResult; + + offsetLastPoint = borderZlength + dimz - 1; + + theLinePlusBorder = theLine + borderYlength; + resLinePlusBorder = resLine + borderYlength; + + /* + * There are dimy*dimx Z lines to be processed. + */ + for ( y=0; y 0 ) { + dbl_pt1 = theLine + borderZlength; dbl_first = *dbl_pt1; + dbl_pt2 = theLine + offsetLastPoint; dbl_last = *dbl_pt2; + for ( z=0; z