/*
* stickkit.c
*
* copyright 2007,8 Mark J. Stock mstock@umich.edu
*
* a program to read, manipulate, and write string, segmented, or
* graph-like data in fewer than 16 space dimensions
*
* compile with
*
* cc -O2 -std=c99 -o stickkit stickkit.c -lm
* /usr/local/bin/i386-mingw32-gcc -O2 -o stickkit.exe stickkit.c -lm
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
*/
// user-changeable parameters
#define BUCKET 100
#define DOTPER 10000
#define MAXDIM 16
#define MAX_ACTIONS 100
// necessary defines
#include
#include
#include
#include
#include
//#include
#define TRUE 1
#define FALSE 0
#define EPSILON 1.0e-6
#define EPSILONSQRD 1.0e-12
#define M_PI 3.14159265358979323846
#define MAXSTR 160
//#define VERSION ("stickkit, version 0.1 2007-06-23 MJS")
//#define VERSION ("stickkit, version 0.2 2007-07-22 MJS")
//#define VERSION ("stickkit, version 0.3 2007-07-28 MJS")
// Version 0.4: added -info and -dots options
//#define VERSION ("stickkit, version 0.4 2008-04-29 MJS")
// Version 0.5: added -split option
#define VERSION ("stickkit, version 0.5 2009-04-05 MJS")
// begin with the data structures
// predefine some structure pointers
typedef struct node *node_ptr;
typedef struct node_group *node_group_ptr;
typedef struct segment *seg_ptr;
typedef struct seg_group *seg_group_ptr;
typedef struct tangent *tan_ptr;
typedef struct tan_group *tan_group_ptr;
typedef struct radius *rad_ptr;
typedef struct rad_group *rad_group_ptr;
// a node
typedef struct node {
unsigned long int index;
char flag;
double *x;
unsigned int numconn0;
seg_ptr *conn0;
unsigned int numconn1;
seg_ptr *conn1;
node_group_ptr parent;
node_ptr prev;
node_ptr next;
} NODE;
// a collection of nodes, like a cell in a tree
typedef struct node_group {
unsigned int num;
node_ptr first;
node_group_ptr parent;
node_group_ptr child[2];
unsigned char axis;
double *min;
double *max;
} NODE_GROUP;
// a segment
typedef struct segment {
unsigned long int index; // id of individual segment
unsigned int block; // id of disconnected block, start at 1
char flag; // temporary flag
node_ptr n[2];
rad_ptr r[2];
tan_ptr t[2];
seg_group_ptr parent;
seg_ptr prev;
seg_ptr next;
} SEGMENT;
// a collection of segments
typedef struct seg_group {
unsigned char dim; // dimensionality of system
double radius; // global radius, -1 if none defined
unsigned long int num; // number of segments in this group
unsigned int numblock; // number of segment blocks in this group
seg_ptr first; // pointer to first segment in linked list
node_group_ptr nodes; // pointer to group of nodes
rad_group_ptr radii; // pointer to group of radius objects
tan_group_ptr tangents; // pointer to group of tangent vectors
} SEGMENT_GROUP;
// a tangent vector
typedef struct tangent {
unsigned long int index;
char flag;
double *x;
tan_group_ptr parent;
tan_ptr prev;
tan_ptr next;
} TANGENT;
// a collection of tangent vectors
typedef struct tan_group {
unsigned int num;
tan_ptr first;
tan_group_ptr parent;
tan_group_ptr child[2];
unsigned char axis;
double *min;
double *max;
} TANGENT_GROUP;
// a radius
typedef struct radius {
unsigned long int index;
char flag;
double r;
rad_group_ptr parent;
rad_ptr prev;
rad_ptr next;
} RADIUS;
// a collection of radii
typedef struct rad_group {
unsigned int num;
rad_ptr first;
rad_group_ptr parent;
rad_group_ptr child[2];
double min;
double max;
} RADIUS_GROUP;
// define each action possible on the data
typedef enum action_type {
none, // no value, no action (default)
read, // read a file
write, // write a file
info, // dump min/max bounds, number of nodes, segments
coarsen, // remove nodes and merge segments
refine, // split segments and add nodes, assume straight-line segments
splrefine, // split segments and add nodes, use spline interpolation
roughen, // nudge nodes to roughen the curve
treeradius, // set radii according to principle of constant stress
translate, // translate all nodes by x,y,z
scale, // scale all nodes by x,y,z
split // split all segments along the longest axis, write both to .rad files
//resample, // resample a strand to a fixed segment length
//smooth, // nudge nodes to smooth the curve
} ACTION_NAME;
typedef struct an_action {
ACTION_NAME type;
unsigned char nc;
char **carg;
unsigned char ni;
int *iarg;
unsigned char nd;
double *darg;
} ACTION;
// define the possible output file types
typedef enum output_format_type {
noout, // default is no output
seg, // Wavefront .obj-like .seg file
rad, // radiance cylinders
vtk, // Mayavi-readable
dots, // ASCII space-delimited, equidistant dots
png // image of segments, for debugging only
} OUT_FORMAT;
// prototypes
// processing
void allocate_action (ACTION*, ACTION_NAME);
seg_group_ptr make_seg_structure (int);
int update_all_stats (seg_group_ptr);
int update_node_group_stats (node_group_ptr, int);
int update_radius_group_stats (rad_group_ptr);
int update_tangent_group_stats (tan_group_ptr, int);
int dump_stats (seg_group_ptr);
int merge_close_nodes (node_group_ptr, int, double);
void delete_segment (node_ptr, node_ptr, seg_ptr, int);
int split_long_segs (seg_group_ptr, double, int, int);
void split_segment (seg_ptr, int, int);
void find_seg_midpt_using_midpoint (double*, double*, seg_ptr, int);
void find_seg_midpt_using_spline (double*, double*, seg_ptr, int);
void find_seg_tangents (seg_ptr, int);
int roughen_nodes (node_group_ptr, int, double);
int identify_separate_strands (seg_group_ptr);
int identify_connected_segments (seg_ptr, unsigned int);
int find_root_of_each_strand (seg_group_ptr, ACTION*);
//int set_treelike_radii_2d (seg_group_ptr, ACTION*);
int set_treelike_radii_3d (seg_group_ptr, ACTION*);
int translate_nodes (node_group_ptr, int, double*);
int scale_nodes (node_group_ptr, int, double*);
int split_and_write_rad (seg_group_ptr, char*);
int split_into_two_seg_groups (node_group_ptr, int, double, seg_group_ptr, seg_group_ptr);
// input/output
int read_seg (char*, seg_group_ptr);
int read_rad (char*, seg_group_ptr);
int read_radiance_floats (FILE*, double*);
seg_ptr add_segment (seg_group_ptr, node_ptr, node_ptr);
node_ptr add_node (node_group_ptr, unsigned char, double*, int);
rad_ptr add_radius (rad_group_ptr, double, int);
tan_ptr add_tangent (tan_group_ptr, unsigned char, double*, int);
node_ptr find_closest_node (node_group_ptr, unsigned char, double *,
node_ptr, double *);
int write_seg (FILE*, seg_group_ptr, int, char**);
int write_vtk (FILE*, seg_group_ptr);
int write_vtk_nodes_and_set_indexes (FILE*, node_group_ptr, int, int);
int write_vtk_nodes_radii (FILE*, node_group_ptr, double);
int write_rad (FILE*, seg_group_ptr);
int write_radiance_nodes (FILE*, node_group_ptr, int, double, int*, int*);
int write_dots (FILE*, seg_group_ptr, double);
int write_dots_1 (node_group_ptr, node_group_ptr, int, double *);
int write_dots_2 (seg_group_ptr, node_group_ptr, int, double *);
int write_dots_3 (FILE *, node_group_ptr, int);
// utility
int set_seg_flags (seg_group_ptr, int);
int set_node_flags (node_group_ptr, int);
int set_radius_flags (rad_group_ptr, int);
void perturb_gaussian (double*, double, int);
double seg_length (seg_ptr, int);
double seg_length_sqrd (seg_ptr, int);
double vec_dist_sqrd (double*, double*, int);
double vec_inv_sqrt (double*, int);
void vec_normalize (double*, double*, int);
void vec_cross (double*, double*, double*);
int free_nodes (node_group_ptr);
int free_radii (rad_group_ptr);
int free_tangents (tan_group_ptr);
int free_segs (seg_group_ptr);
int Usage(char[MAXSTR],int);
// dj don't bore us, get back to the chorus!
int main (int argc, char **argv) {
// main non-global definitions
int i,j,act;
int have_input_file = FALSE;
double dot_scale = 0.3;
seg_group_ptr segs;
//seg_group_ptr segs = (SEGMENT_GROUP*) malloc (sizeof(SEGMENT_GROUP));
char progname[MAXSTR],infile[MAXSTR],root[MAXSTR],extension[4];
OUT_FORMAT out_format = noout;
int nactions = 0;
ACTION action[MAX_ACTIONS];
FILE *outptr = stdout;
// parse command-line arguments
// input file can appear anywhere
// processing options are performed in the order they appear
(void) strcpy(progname,argv[0]);
if (argc < 2) (void) Usage(progname,0);
for (i=1; i i+1) {
if (!isalpha(argv[i+1][1])) {
dot_scale = atof(argv[++i]);
}
}
} else if (strncmp(argv[i], "-png", 4) == 0) {
out_format = png;
} else if (strncmp(argv[i], "-noop", 5) == 0) {
(void) allocate_action (&action[nactions], none);
nactions++;
} else if (strncmp(argv[i], "-coarsen", 5) == 0) {
(void) allocate_action (&action[nactions], coarsen);
if (isalpha(argv[i+1][1])) {
action[nactions].darg[0] = 1.;
} else {
action[nactions].darg[0] = atof(argv[++i]);
}
nactions++;
} else if (strncmp(argv[i], "-srefine", 5) == 0) {
(void) allocate_action (&action[nactions], splrefine);
if (isalpha(argv[i+1][1])) {
action[nactions].darg[0] = 1.;
} else {
action[nactions].darg[0] = atof(argv[++i]);
}
nactions++;
} else if (strncmp(argv[i], "-refine", 5) == 0) {
(void) allocate_action (&action[nactions], refine);
if (isalpha(argv[i+1][1])) {
action[nactions].darg[0] = 1.;
} else {
action[nactions].darg[0] = atof(argv[++i]);
}
nactions++;
} else if (strncmp(argv[i], "-roughen", 5) == 0) {
(void) allocate_action (&action[nactions], roughen);
if (isalpha(argv[i+1][1])) {
action[nactions].darg[0] = 1.;
} else {
action[nactions].darg[0] = atof(argv[++i]);
}
nactions++;
} else if (strncmp(argv[i], "-treeradius", 4) == 0) {
(void) allocate_action (&action[nactions], treeradius);
if (isalpha(argv[i+1][1])) {
action[nactions].darg[0] = 1.e+3;
} else {
action[nactions].darg[0] = atof(argv[++i]);
if (isalpha(argv[i+1][1])) {
action[nactions].iarg[0] = 0;
} else {
action[nactions].iarg[0] = atoi(argv[++i]);
if (isalpha(argv[i+1][1])) {
action[nactions].darg[1] = 0.;
} else {
action[nactions].darg[1] = atof(argv[++i]);
}
}
}
nactions++;
} else if (strncmp(argv[i], "-translate", 4) == 0) {
(void) allocate_action (&action[nactions], translate);
if (isalpha(argv[i+1][1])) {
action[nactions].darg[0] = 0.;
} else {
action[nactions].darg[0] = atof(argv[++i]);
if (isalpha(argv[i+1][1])) {
action[nactions].darg[1] = 0.;
} else {
action[nactions].darg[1] = atof(argv[++i]);
if (isalpha(argv[i+1][1])) {
action[nactions].darg[2] = 0.;
} else {
action[nactions].darg[2] = atof(argv[++i]);
}
}
}
nactions++;
} else if (strncmp(argv[i], "-scale", 3) == 0) {
(void) allocate_action (&action[nactions], scale);
for (j=0;j i+1) {
if (!isalpha(argv[i+1][1])) {
i++;
for (j=0;j i+1) {
if (!isalpha(argv[i+1][1])) {
i++;
for (j=1;j i+1) {
if (!isalpha(argv[i+1][1])) {
i++;
for (j=2;jnum,segs->nodes->num);
// do the merge, interate until all are done
i = 1;
while (i > 0) {
(void) set_node_flags (segs->nodes, FALSE);
i = merge_close_nodes (segs->nodes, segs->dim,
pow(action[act].darg[0],2));
fprintf (stderr," performed %d merges\n",i);
}
// reset the counters and bounds
(void) update_all_stats (segs);
fprintf (stderr," end with %d segs, %d nodes\n",
segs->num,segs->nodes->num);
// remove nodes closer than a threshold distance
} else if (action[act].type == splrefine) {
fprintf (stderr,"spline refine\n");
fprintf (stderr," begin with %d segs, %d nodes\n",
segs->num,segs->nodes->num);
// should we do the normalized version? if thresh is negative
if (action[act].darg[0] > 0.) j = FALSE;
else j = TRUE;
//fprintf (stderr," normalization is %d\n",j);
// do the splits, interate until all are done
i = 1;
while (i > 0) {
(void) set_seg_flags (segs, FALSE);
i = split_long_segs (segs, pow(action[act].darg[0],2), TRUE, j);
fprintf (stderr," performed %d splits\n",i);
}
// reset the counters and bounds
(void) update_all_stats (segs);
fprintf (stderr," end with %d segs, %d nodes\n",
segs->num,segs->nodes->num);
// remove nodes closer than a threshold distance
} else if (action[act].type == refine) {
fprintf (stderr,"refine\n");
fprintf (stderr," begin with %d segs, %d nodes\n",
segs->num,segs->nodes->num);
// should we do the normalized version? if thresh is negative
if (action[act].darg[0] > 0.) j = FALSE;
else j = TRUE;
//fprintf (stderr," normalization is %d\n",j);
// do the splits, interate until all are done
i = 1;
while (i > 0) {
(void) set_seg_flags (segs, FALSE);
i = split_long_segs (segs, pow(action[act].darg[0],2), FALSE, j);
fprintf (stderr," performed %d splits\n",i);
}
// reset the counters and bounds
(void) update_all_stats (segs);
fprintf (stderr," end with %d segs, %d nodes\n",
segs->num,segs->nodes->num);
// jiggle nodes around
} else if (action[act].type == roughen) {
fprintf (stderr,"roughen\n");
(void) roughen_nodes (segs->nodes, segs->dim, action[act].darg[0]);
// set radii according to principle of constant stress
} else if (action[act].type == treeradius) {
fprintf (stderr,"treeradius\n");
// first, identify separate disconnected strands
i = identify_separate_strands (segs);
fprintf (stderr," found %d separate strands\n",i);
// then, find the root node of each strand
i = find_root_of_each_strand (segs, &action[act]);
fprintf (stderr," found %d strand roots\n",i);
// finally, set the radii
(void) set_seg_flags (segs, FALSE);
i = 1;
while (i > 0) {
if (segs->dim == 2) {
//i = set_treelike_radii_2d (segs, &action[act]);
} else if (segs->dim == 3) {
i = set_treelike_radii_3d (segs, &action[act]);
} else {
fprintf (stderr," Error: treeradius cannot work on systems with");
fprintf (stderr," %d dimensions\n",segs->dim);
fprintf (stderr," works on dim 2 and 3\n");
fprintf (stderr," skipping this step\n");
i = 0;
}
fprintf (stderr,".");
}
fprintf (stderr,"\n");
// remove nodes closer than a threshold distance
} else if (action[act].type == translate) {
fprintf (stderr,"translate\n");
(void) translate_nodes (segs->nodes, segs->dim, action[act].darg);
// remove nodes closer than a threshold distance
} else if (action[act].type == scale) {
fprintf (stderr,"scale\n");
(void) scale_nodes (segs->nodes, segs->dim, action[act].darg);
// split and write .rad files
} else if (action[act].type == split) {
// find input filename root
if (strrchr(infile,'.')) {
strncpy(root,infile,strrchr(infile,'.')-infile);
root[strrchr(infile,'.')-infile] = '_';
root[strrchr(infile,'.')-infile+1] = '\0';
} else {
// there is no "." in the filename?
sprintf(root,"%s_",infile);
}
//fprintf(stderr,"root is (%s)\n",root);
fprintf (stderr,"split\n");
(void) split_and_write_rad (segs, root);
// dump statistics and quit
} else if (action[act].type == info) {
fprintf (stderr,"info\n");
(void) dump_stats (segs);
// no action
} else {
fprintf(stderr,"none\n");
}
}
exit(0);
}
/*
* malloc the necessary space for this type of action
*/
void allocate_action (ACTION *action, ACTION_NAME type) {
int i;
action->type = type;
if (type == read) {
// a file read, allocate one character variable, none others
action->nc = 1;
action->ni = 0;
action->nd = 0;
} else if (type == write) {
// a file write, allocate one character variable, none others
action->nc = 1;
action->ni = 0;
action->nd = 0;
} else if (type == info) {
// dump current segments information
action->nc = 0;
action->ni = 0;
action->nd = 0;
} else if (type == coarsen) {
// need a threshold
action->nc = 0;
action->ni = 0;
action->nd = 1;
} else if (type == splrefine) {
// need a threshold
action->nc = 0;
action->ni = 0;
action->nd = 1;
} else if (type == refine) {
// need a threshold
action->nc = 0;
action->ni = 0;
action->nd = 1;
} else if (type == roughen) {
// need a scale
action->nc = 0;
action->ni = 0;
action->nd = 1;
} else if (type == treeradius) {
// need a vector
action->nc = 0;
action->ni = 1;
action->nd = 2;
} else if (type == translate) {
// need a vector
action->nc = 0;
action->ni = 0;
action->nd = MAXDIM;
} else if (type == scale) {
// need a vector
action->nc = 0;
action->ni = 0;
action->nd = MAXDIM;
} else if (type == split) {
// no need for arguments
action->nc = 0;
action->ni = 0;
action->nd = 0;
} else if (type == none) {
// no operation!
action->nc = 0;
action->ni = 0;
action->nd = 0;
} else {
fprintf (stderr,"ERROR (allocate_action): not supposed to get here.\n");
fprintf (stderr,"Quitting.\n");
exit(0);
}
// now, malloc the space
if (action->nc > 0) {
action->carg = (char**) malloc (action->nc * sizeof(char*));
for (i=0; inc; i++)
action->carg[i] = (char*) malloc (MAXSTR * sizeof(char));
} else {
action->carg = NULL;
}
if (action->ni > 0) {
action->iarg = (int*) malloc (action->ni * sizeof(int));
} else {
action->iarg = NULL;
}
if (action->nd > 0) {
action->darg = (double*) malloc (action->nd * sizeof(double));
} else {
action->darg = NULL;
}
return;
}
/*
* Create an independent set of related nodes and segments
*
* Should we support having multiple segment groups active?
*/
seg_group_ptr make_seg_structure (int dim) {
seg_group_ptr segs = (SEGMENT_GROUP*) malloc (sizeof(SEGMENT_GROUP));
segs->dim = dim;
segs->num = 0;
segs->first = NULL;
segs->nodes = (NODE_GROUP*) malloc (sizeof(NODE_GROUP));
segs->nodes->num = 0;
segs->nodes->first = NULL;
segs->nodes->parent = NULL;
segs->nodes->child[0] = NULL;
segs->nodes->child[1] = NULL;
segs->nodes->axis = -1;
segs->nodes->min = NULL;
segs->nodes->max = NULL;
segs->radii = (RADIUS_GROUP*) malloc (sizeof(RADIUS_GROUP));
segs->radii->num = 0;
segs->radii->first = NULL;
segs->radii->parent = NULL;
segs->radii->child[0] = NULL;
segs->radii->child[1] = NULL;
segs->radii->min = 9.9e+9;
segs->radii->max = -9.9e+9;
segs->tangents = (TANGENT_GROUP*) malloc (sizeof(TANGENT_GROUP));
segs->tangents->num = 0;
segs->tangents->first = NULL;
segs->tangents->parent = NULL;
segs->tangents->child[0] = NULL;
segs->tangents->child[1] = NULL;
segs->tangents->axis = -1;
segs->tangents->min = NULL;
segs->tangents->max = NULL;
return (segs);
}
// ==========================================================================
// process the data
//
/*
* Update numbers, bounds of all data
*/
int update_all_stats (seg_group_ptr this) {
int nrad, nnode, nseg; // ntan
seg_ptr curr;
// do segments first (they're easy)
nseg = 0;
curr = this->first;
while (curr) {
curr->parent = this;
nseg++;
curr = curr->next;
}
this->num = nseg;
// then do nodes
nnode = update_node_group_stats (this->nodes, this->dim);
// do radii next
nrad = update_radius_group_stats (this->radii);
// do tangents last
//ntan = update_tangent_group_stats (this->tangents, this->dim);
// if things bomb, send back a nonzero
return (0);
}
/*
* Print the system statistics
*/
int dump_stats (seg_group_ptr this) {
int i;
fprintf(stdout,"# number of dimensions %d\n",this->dim);
fprintf(stdout,"# number of segment blocks %d\n",this->numblock);
fprintf(stdout,"# number of segments %ld\n",this->num);
fprintf(stdout,"# number of nodes %ld\n",this->nodes->num);
fprintf(stdout,"# node minima");
for (i=0; idim; i++) fprintf(stdout," %g",this->nodes->min[i]);
fprintf(stdout,"\n");
fprintf(stdout,"# node maxima");
for (i=0; idim; i++) fprintf(stdout," %g",this->nodes->max[i]);
fprintf(stdout,"\n");
fprintf(stdout,"# number of radii %ld\n",this->radii->num);
fprintf(stdout,"# radius minimum %g\n",this->radii->min);
fprintf(stdout,"# radius maximum %g\n",this->radii->max);
fprintf(stdout,"# number of tangents %ld\n",this->tangents->num);
return(0);
}
/*
* Update the num, min, max of a group of nodes
*/
int update_node_group_stats (node_group_ptr this, int dim) {
int i,j;
int cnt = 0;
node_ptr curr;
// reset the min/max
for (i=0; imin[i] = 9.9e+9;
this->max[i] = -9.9e+9;
}
// check all children
if (this->child[0]) {
// update counts and bounds of children
cnt += update_node_group_stats (this->child[0], dim);
cnt += update_node_group_stats (this->child[1], dim);
// use those bounds to update this bounds
for (j=0; j<2; j++) {
for (i=0; ichild[j]->min[i] < this->min[i])
this->min[i] = this->child[j]->min[i];
if (this->child[j]->max[i] > this->max[i])
this->max[i] = this->child[j]->max[i];
}
}
} else {
// check all local nodes
curr = this->first;
while (curr) {
// update parent
curr->parent = this;
// update count
cnt++;
// update bounds
for (i=0; ix[i] < this->min[i])
this->min[i] = curr->x[i];
if (curr->x[i] > this->max[i])
this->max[i] = curr->x[i];
}
curr = curr->next;
}
}
this->num = cnt;
return (cnt);
}
/*
* Update the num, min, max of a group of radii
*/
int update_radius_group_stats (rad_group_ptr this) {
int cnt = 0;
rad_ptr curr;
// if this group is NULL, skip and return
if (this) {
// reset the min/max
this->min = 9.9e+9;
this->max = -9.9e+9;
// first, count the local entries and compute min/max
curr = this->first;
while (curr) {
// update parent
curr->parent = this;
// update count
cnt++;
// update bounds
if (curr->r < this->min) this->min = curr->r;
if (curr->r > this->max) this->max = curr->r;
curr = curr->next;
}
// then, add on all children
if (this->child[0]) {
cnt += update_radius_group_stats (this->child[0]);
// and check the children's bounds
if (this->child[0]->min < this->min) this->min = this->child[0]->min;
if (this->child[0]->max > this->max) this->max = this->child[0]->max;
}
if (this->child[1]) {
cnt += update_radius_group_stats (this->child[1]);
if (this->child[1]->min < this->min) this->min = this->child[1]->min;
if (this->child[1]->max > this->max) this->max = this->child[1]->max;
}
// set the new count
this->num = cnt;
}
return (cnt);
}
/*
* Merge nodes that share a common segment if their distance is under
* a threshold (the threshold is squared already!)
*
* cnt is the number of merges that took place
*/
int merge_close_nodes (node_group_ptr this, int dim, double thresh) {
int i;
int cnt = 0;
double dist,mindist;
node_ptr curr,delnode;
seg_ptr delseg;
// check all children
if (this->child[0]) {
cnt += merge_close_nodes (this->child[0], dim, thresh);
cnt += merge_close_nodes (this->child[1], dim, thresh);
} else {
// check all connected segments
//fprintf (stderr,"group with %d nodes, thresh %g\n",this->num,thresh);
curr = this->first;
while (curr) {
// loop through all connected segments
//fprintf (stderr, "node %ld at %g %g %g has neighbors at\n",
// curr->index, curr->x[0], curr->x[1], curr->x[2]);
// loop until there are no more close nodes
//mindist = 0.;
//while (mindist < thresh) {
// no, do looping higher up!
// find the closest adjacent connected node
mindist = 9.9e+9;
delnode = NULL;
delseg = NULL;
for (i=0; inumconn0; i++) {
// how close is the connected node?
dist = vec_dist_sqrd (curr->x, curr->conn0[i]->n[1]->x, dim);
//fprintf (stderr, " %ld at %g %g %g dist %g\n",
// curr->conn0[i]->n[1]->index,curr->conn0[i]->n[1]->x[0],
// curr->conn0[i]->n[1]->x[1],curr->conn0[i]->n[1]->x[2],
// sqrt(dist));
if (dist < mindist) {
mindist = dist;
delnode = curr->conn0[i]->n[1];
delseg = curr->conn0[i];
}
}
for (i=0; inumconn1; i++) {
dist = vec_dist_sqrd (curr->x, curr->conn1[i]->n[0]->x, dim);
//fprintf (stderr, " %ld at %g %g %g dist %g\n",
// curr->conn1[i]->n[0]->index,curr->conn1[i]->n[0]->x[0],
// curr->conn1[i]->n[0]->x[1],curr->conn1[i]->n[0]->x[2],
// sqrt(dist));
if (dist < mindist) {
mindist = dist;
delnode = curr->conn1[i]->n[0];
delseg = curr->conn1[i];
}
}
// delete them both
if (mindist < thresh && !curr->flag && !delnode->flag) {
delete_segment (curr, delnode, delseg, dim);
cnt++;
}
//}
curr = curr->next;
}
}
return (cnt);
}
/*
* General routine for deleting a segment (and joining two nodes)
*/
void delete_segment (node_ptr keepnode, node_ptr delnode, seg_ptr delseg,
int dim) {
int i,num_link_keep,num_link_del,copyto;
int keepis0 = -1;
double weight,newloc[dim];
seg_ptr *newconnlist;
rad_ptr newrad,keeprad,delrad;
//fprintf (stderr, " should we delete?\n");
// first, make sure that the two given nodes are actually joined
// by the segment in question!
if (delseg->n[0] == keepnode) {
if (delseg->n[1] == delnode) {
keepis0 = TRUE;
} else {
keepis0 = -1; // bad!
}
} else if (delseg->n[1] == keepnode) {
if (delseg->n[0] == delnode) {
keepis0 = FALSE;
} else {
keepis0 = -1; // bad!
}
} else {
keepis0 = -1; // bad!
}
// if they're not right, dump and quit
if (keepis0 == -1) {
fprintf (stderr,"ERROR (delete_segment): segment %ld was given nodes\n",delseg->index);
fprintf (stderr," (%ld) and (%ld), but that segment's end nodes are\n",keepnode->index,delnode->index);
fprintf (stderr," actually (%ld) and (%ld). Quitting.\n",delseg->n[0]->index,delseg->n[1]->index);
exit(1);
}
//fprintf (stderr, " deleting node %ld and seg %ld, keeping node %ld\n",delnode->index,delseg->index,keepnode->index);
// throw away the index and the flag of the old node
// and set the new flag to indicate that it's participated in a merge
keepnode->flag = TRUE;
// relocate keepnode, weight by number of connections
num_link_keep = keepnode->numconn0 + keepnode->numconn1;
num_link_del = delnode->numconn0 + delnode->numconn1;
//fprintf (stderr, " num_link_keep %d, num_link_del %d\n",num_link_keep,num_link_del);
if (num_link_keep == 1 && num_link_del == 2) {
weight = 1.0;
} else if (num_link_keep == 2 && num_link_del == 1) {
weight = 0.0;
} else {
weight = (double)num_link_keep/(double)(num_link_keep+num_link_del);
}
for (i=0; ix[i] = weight*keepnode->x[i] + (1.-weight)*delnode->x[i];
// now, find the new radius
if (keepis0) {
keeprad = delseg->r[0];
delrad = delseg->r[1];
newrad = add_radius (delseg->parent->radii, weight*delseg->r[0]->r +
(1.-weight)*delseg->r[1]->r, 0);
} else {
keeprad = delseg->r[1];
delrad = delseg->r[0];
newrad = add_radius (delseg->parent->radii, weight*delseg->r[1]->r +
(1.-weight)*delseg->r[0]->r, 0);
}
// do something with this!
// assign it to all still-connected segments
for (i=0; inumconn0; i++)
if (keepnode->conn0[i]->r[0] == keeprad)
keepnode->conn0[i]->r[0] = newrad;
for (i=0; inumconn1; i++)
if (keepnode->conn1[i]->r[1] == keeprad)
keepnode->conn1[i]->r[1] = newrad;
for (i=0; inumconn0; i++)
if (delnode->conn0[i]->r[0] == delrad)
delnode->conn0[i]->r[0] = newrad;
for (i=0; inumconn1; i++)
if (delnode->conn1[i]->r[1] == delrad)
delnode->conn1[i]->r[1] = newrad;
// copy all conn0 linkages from delnode to keepnode
//fprintf (stderr, " numconn0: %d + %d",keepnode->numconn0,delnode->numconn0);
newconnlist = (seg_ptr*) malloc (
(keepnode->numconn0 + delnode->numconn0) * sizeof(seg_ptr));
for (i=0; inumconn0; i++) newconnlist[i] = keepnode->conn0[i];
for (i=0; inumconn0; i++) {
delnode->conn0[i]->n[0] = keepnode;
// add these to keepnode->conn0
newconnlist[i + keepnode->numconn0] = delnode->conn0[i];
}
free (keepnode->conn0);
keepnode->conn0 = newconnlist;
keepnode->numconn0 += delnode->numconn0;
//fprintf (stderr, " = %d\n",keepnode->numconn0);
// copy all conn1 linkages from delnode to keepnode
//fprintf (stderr, " numconn1: %d + %d",keepnode->numconn1,delnode->numconn1);
newconnlist = (seg_ptr*) malloc (
(keepnode->numconn1 + delnode->numconn1) * sizeof(seg_ptr));
for (i=0; inumconn1; i++) newconnlist[i] = keepnode->conn1[i];
for (i=0; inumconn1; i++) {
delnode->conn1[i]->n[1] = keepnode;
// add these to keepnode->conn1
newconnlist[i + keepnode->numconn1] = delnode->conn1[i];
}
free (keepnode->conn1);
keepnode->conn1 = newconnlist;
keepnode->numconn1 += delnode->numconn1;
//fprintf (stderr, " = %d\n",keepnode->numconn1);
// if it's the first node, reset the first pointer
if (delnode->parent->first == delnode)
delnode->parent->first = delnode->next;
delnode->parent->num--;
// reset pointers
if (delnode->prev) delnode->prev->next = delnode->next;
if (delnode->next) delnode->next->prev = delnode->prev;
// remove delnode
free (delnode->conn0);
free (delnode->conn1);
free (delnode);
// now, remove delseg
// if it's the first segment, reset the first pointer
if (delseg->parent->first == delseg)
delseg->parent->first = delseg->next;
delseg->parent->num--;
// lose everything except the pointers
if (delseg->prev) delseg->prev->next = delseg->next;
if (delseg->next) delseg->next->prev = delseg->prev;
// remove *both* entries to delseg in keepnode's conn array
// locate delseg in the conn0 array and compress the array
copyto = 0;
for (i=0; inumconn0; i++) {
keepnode->conn0[copyto] = keepnode->conn0[i];
if (keepnode->conn0[i] != delseg) copyto++;
}
keepnode->numconn0 = copyto;
copyto = 0;
// compress the conn1 array instead
for (i=0; inumconn1; i++) {
keepnode->conn1[copyto] = keepnode->conn1[i];
if (keepnode->conn1[i] != delseg) copyto++;
}
keepnode->numconn1 = copyto;
// don't worry about reallocating memory here
// debug print the keepnode's connectivity list
//fprintf (stderr, " numconn0 are\n");
//for (i=0; inumconn0; i++) {
// fprintf (stderr, " seg %ld, nodes %ld %ld\n",keepnode->conn0[i]->index,keepnode->conn0[i]->n[0]->index,keepnode->conn0[i]->n[1]->index);
//}
//fprintf (stderr, " numconn1 are\n");
//for (i=0; inumconn1; i++) {
// fprintf (stderr, " seg %ld, nodes %ld %ld\n",keepnode->conn1[i]->index,keepnode->conn1[i]->n[0]->index,keepnode->conn1[i]->n[1]->index);
//}
// finally, axe the segment
free (delseg);
//fprintf (stderr, " done deleting\n");
return;
}
/*
* Split nodes that share a common segment if their distance is over
* a threshold (the threshold is squared already!)
*
* cnt is the number of splits that took place
* thresh is the *square* of the input threshold length
*/
int split_long_segs (seg_group_ptr this, double thresh, int use_splines, int use_normalized_lengths) {
int cnt = 0;
double dist,rad;
seg_ptr curr;
curr = this->first;
while (curr) {
// how long is this segment
dist = vec_dist_sqrd (curr->n[0]->x, curr->n[1]->x, this->dim);
// scale the length by the radius
if (use_normalized_lengths) {
rad = 0.5*(curr->r[0]->r + curr->r[1]->r);
dist /= pow(rad,2);
//fprintf(stderr,"dist %g rad %g thresh %g doit %d\n",sqrt(dist),rad,sqrt(thresh),dist>thresh);
}
// split this segment
if (dist > thresh && !curr->flag) {
split_segment (curr, this->dim, use_splines);
cnt++;
}
curr = curr->next;
}
return (cnt);
}
/*
* General routine for splitting a segment
*/
void split_segment (seg_ptr this, int dim, int use_spline) {
int i,copyto;
//int use_spline = TRUE;
double radval,newloc[dim];
seg_ptr newseg = NULL;
node_ptr newnode = NULL;
rad_ptr newrad = NULL;
//fprintf (stderr, " splitting seg %ld, nodes %ld %ld\n",this->index,this->n[0]->index,this->n[1]->index);
// locate the new node
if (use_spline) {
// spline-fit
find_seg_midpt_using_spline (newloc, &radval, this, dim);
} else {
// naive method (midpoint)
find_seg_midpt_using_midpoint (newloc, &radval, this, dim);
}
// make the new node
newnode = add_node (this->parent->nodes, dim, newloc, 0);
// make the segment
newseg = add_segment (this->parent, newnode, this->n[1]);
// now, set the new radius
newrad = add_radius (this->parent->radii, radval, 0);
// assign a bunch of pointers!
// first, the new node
newnode->numconn0 = 1;
newnode->conn0 = (seg_ptr*) malloc (newnode->numconn0 * sizeof(seg_ptr));
newnode->conn0[0] = newseg;
newnode->numconn1 = 1;
newnode->conn1 = (seg_ptr*) malloc (newnode->numconn1 * sizeof(seg_ptr));
newnode->conn1[0] = this;
// but now this segment is no longer connected to its n[1]
copyto = 0;
for (i=0; in[1]->numconn1; i++) {
this->n[1]->conn1[copyto] = this->n[1]->conn1[i];
if (this->n[1]->conn1[i] != this) copyto++;
}
this->n[1]->numconn1 = copyto;
// then the new segment
newseg->r[0] = newrad;
newseg->r[1] = this->r[1];
newseg->t[1] = this->t[1];
// and the old segment
this->n[1] = newnode;
this->r[1] = newrad;
this->t[1] = NULL;
// set the two seg flags to indicate that they've participated in a split
this->flag = TRUE;
newseg->flag = TRUE;
return;
}
/*
* use a naive method to determine the new node location - the midpoint
*/
void find_seg_midpt_using_midpoint (double *loc, double *rad,
seg_ptr this, int dim) {
int i;
for (i=0;in[0]->x[i] + this->n[1]->x[i]);
*rad = 0.5*this->r[0]->r + 0.5*this->r[1]->r;
return;
}
/*
* Make the new location at the midpoint of a cubic spline threaded
* through the two endpoints
*/
void find_seg_midpt_using_spline (double *loc, double *rad,
seg_ptr this, int dim) {
int i,j;
//double dl[dim],len,fp[2][3][3],p1[dim],p2[dim];
double a[4];
double dl,tan0[dim],tan1[dim];
double r0,r1,r2;
//double difflen = 0.;
//static double maxdiff = 0.;
//node_ptr n1 = this->n[0];
//node_ptr n2 = this->n[1];
// test print
//fprintf(stderr,"\nnode0 at %g %g %g\n",this->n[0]->x[0],this->n[0]->x[1],this->n[0]->x[2]);
//fprintf(stderr,"node1 at %g %g %g\n",this->n[1]->x[0],this->n[1]->x[1],this->n[1]->x[2]);
// is a tangent defined here?
if (!this->t[0] || !this->t[1]) find_seg_tangents (this, dim);
//fprintf(stderr," tangent0 %g %g %g\n",this->t[0]->x[0],this->t[0]->x[1],this->t[0]->x[2]);
//fprintf(stderr," tangent1 %g %g %g\n",this->t[1]->x[0],this->t[1]->x[1],this->t[1]->x[2]);
// compute the length of the edge
/*
len = 0.;
for (i=0; ix[i] - n1->x[i];
len += dl[i]*dl[i];
}
len = sqrt(len);
// compute the tangential operator for each node (P = I - nn^T)
for (i=0; inorm[i]*n1->norm[j];
fp[1][i][j] -= n2->norm[i]*n2->norm[j];
}
}
// find the vector product of each of these with dl
for (i=0; in[0]->x, this->n[1]->x, dim));
for (i=0; it[0]->x[i];
for (i=0; it[1]->x[i];
// compute a spline for each coordinate axis
for (i=0; in[0]->x[i];
// a[1] = f'_0
a[1] = tan0[i];
// a[2] = 3 (f_1 - f_0) / h^2 - (f'_1 + 2 f'_0) / h
// but h=1 for us
a[2] = 3. * (this->n[1]->x[i] - this->n[0]->x[i])
- (tan1[i] + 2. * tan0[i]);
// a[3] = 2 (f_0 - f_1) / h^3 + (f'_1 + f'_0) / h^2
a[3] = 2. * (this->n[0]->x[i] - this->n[1]->x[i])
+ (tan0[i] + tan1[i]);
// and the node that we want is 1/2 along the curve
loc[i] = a[0] + 0.5*a[1] + 0.25*a[2] + 0.125*a[3];
//fprintf(stdout," a %g %g %g %g\n",a[0],a[1],a[2],a[3]);
}
//fprintf(stderr," real node loc %g %g %g\n",loc[0],loc[1],loc[2]);
//exit(0);
// Now, do the radius (must find slope of radius?) -------------------
// find the end radii of this segment
if (this->r[0]) r0 = this->r[0]->r;
else r0 = this->parent->radius;
if (this->r[1]) r1 = this->r[1]->r;
else r1 = this->parent->radius;
// find the length of this segment
dl = sqrt(vec_dist_sqrd (this->n[0]->x, this->n[1]->x, dim));
//fprintf(stderr,"\ndoing radius between %g and %g, length %g\n",r0,r1,dl);
if (r0 < -EPSILON || r1 < -EPSILON) exit(0);
// find the radius of the next node past n[0]
j = this->n[0]->numconn0 + this->n[0]->numconn1;
if (j == 1) {
// n0 is an end, what to do?
tan0[0] = r1-r0;
} else if (j==2) {
// there's one node past n[0], find its radius
for (i=0; in[0]->numconn0; i++) {
if (this->n[0]->conn0[i] != this) {
// find the length of that segment
tan0[1] = sqrt (vec_dist_sqrd (this->n[0]->conn0[i]->n[0]->x,
this->n[0]->conn0[i]->n[1]->x, dim));
// find the radius of that segment
if (this->n[0]->conn0[i]->r[1]) r2 = this->n[0]->conn0[i]->r[1]->r;
else r2 = this->parent->radius;
// find the slope of the tangent
//tan0[0] = (r1-r2) / (dl+tan0[1]);
tan0[0] = 0.5*(r0-r2)*dl/tan0[1] + 0.5*(r1-r0);
}
}
for (i=0; in[0]->numconn1; i++) {
if (this->n[0]->conn1[i] != this) {
// find the length of that segment
tan0[1] = sqrt (vec_dist_sqrd (this->n[0]->conn1[i]->n[0]->x,
this->n[0]->conn1[i]->n[1]->x, dim));
// find the radius of that segment
if (this->n[0]->conn1[i]->r[0]) r2 = this->n[0]->conn1[i]->r[0]->r;
else r2 = this->parent->radius;
// find the slope of the tangent
//tan0[0] = (r1-r2) / (dl+tan0[1]);
tan0[0] = 0.5*(r0-r2)*dl/tan0[1] + 0.5*(r1-r0);
}
}
//fprintf(stderr," other seg ends with rad %g and len %g\n",r2,tan0[1]);
} else {
// if there are 3 segments connected to n[0], then just assume linear
tan0[0] = r1-r0;
}
// find the radius of the next node past n[1]
j = this->n[1]->numconn0 + this->n[1]->numconn1;
if (j == 1) {
// n0 is an end, what to do?
tan1[0] = r1-r0;
} else if (j==2) {
// there's one node past n[1], find its radius
for (i=0; in[1]->numconn0; i++) {
if (this->n[1]->conn0[i] != this) {
// find the length of that segment
tan1[1] = sqrt (vec_dist_sqrd (this->n[1]->conn0[i]->n[0]->x,
this->n[1]->conn0[i]->n[1]->x, dim));
// find the radius of that segment
if (this->n[1]->conn0[i]->r[1]) r2 = this->n[1]->conn0[i]->r[1]->r;
else r2 = this->parent->radius;
// find the slope of the tangent
//tan1[0] = (r2-r0) / (dl+tan1[1]);
tan1[0] = 0.5*(r2-r1)*dl/tan1[1] + 0.5*(r1-r0);
}
}
for (i=0; in[1]->numconn1; i++) {
if (this->n[1]->conn1[i] != this) {
// find the length of that segment
tan1[1] = sqrt (vec_dist_sqrd (this->n[1]->conn1[i]->n[0]->x,
this->n[1]->conn1[i]->n[1]->x, dim));
// find the radius of that segment
if (this->n[1]->conn1[i]->r[0]) r2 = this->n[1]->conn1[i]->r[0]->r;
else r2 = this->parent->radius;
// find the slope of the tangent
//tan1[0] = (r2-r0) / (dl+tan1[1]);
tan1[0] = 0.5*(r2-r1)*dl/tan1[1] + 0.5*(r1-r0);
}
}
//fprintf(stderr," other seg ends with rad %g and len %g\n",r2,tan1[1]);
} else {
// if there are 3 segments connected to n[1], then just assume linear
tan1[0] = r1-r0;
}
//fprintf(stderr," tangents %g and %g\n",tan0[0],tan1[0]);
// do the spline interpolation
a[0] = r0;
a[1] = tan0[0];
a[2] = 3. * (r1 - r0) - (tan1[0] + 2. * tan0[0]);
a[3] = 2. * (r0 - r1) + (tan0[0] + tan1[0]);
*rad = a[0] + 0.5*a[1] + 0.25*a[2] + 0.125*a[3];
//fprintf(stderr," new radius %g\n",*rad);
if (*rad < EPSILON) {
//fprintf(stderr," new radius %g\n",*rad);
*rad = 0.;
//exit(0);
}
return;
}
/*
* Find the tangent vectors for a segment
*/
void find_seg_tangents (seg_ptr this, int dim) {
int i,j;
int do_later0 = FALSE;
int do_later1 = FALSE;
double dotp,thisrad,otherrad,radsum;
double x[dim],dl[dim],dlnorm[dim],other[dim];
node_ptr n0 = this->n[0];
node_ptr n1 = this->n[1];
// we may need dl
for (i=0; ix[i] - n0->x[i];
(void) vec_normalize (dl, dlnorm, dim);
//fprintf(stderr," dlnorm %g %g %g\n",dlnorm[0],dlnorm[1],dlnorm[2]);
// do we have a tangent vector already?
if (!this->t[0]) {
// create the tangent
this->t[0] = add_tangent (this->parent->tangents, dim, dl, 0);
// if this is an end, then mirror the other tangent
if (n0->numconn0 + n0->numconn1 == 1) {
do_later0 = TRUE;
// otherwise, use summations
} else {
// now, loop over connected segments, counting this one
for (i=0; inumconn0,n0->numconn1);
for (j=0; jnumconn0; j++) {
// sum the radius-weighted, normalized tangent vectors
for (i=0; iconn0[j]->n[1]->x[i] -
n0->conn0[j]->n[0]->x[i];
(void) vec_normalize (other, other, dim);
if (n0->conn0[j]->r[0])
otherrad = pow (n0->conn0[j]->r[0]->r, 2);
else
otherrad = pow (this->parent->radius, 2);
for (i=0; inumconn1; j++) {
for (i=0; iconn1[j]->n[1]->x[i] -
n0->conn1[j]->n[0]->x[i];
(void) vec_normalize (other, other, dim);
if (n0->conn1[j]->r[1])
otherrad = pow (n0->conn1[j]->r[1]->r, 2);
else
otherrad = pow (this->parent->radius, 2);
for (i=0; ir[0]) thisrad = pow (this->r[0]->r, 2);
else thisrad = pow (this->parent->radius, 2);
// how important is this segment compared to all other segs?
thisrad = 2.*thisrad/radsum;
// finally, if this node just has 2 segments, use the same
// tangent on the other end
if (n0->numconn0 + n0->numconn1 == 2) {
// just use the weighted norm
for (i=0; it[0]->x[i] = x[i];
(void) vec_normalize (this->t[0]->x, this->t[0]->x, dim);
// and copy it to the other segment
for (j=0; jnumconn0; j++) {
if (n0->conn0[j] != this) {
if (!n0->conn0[j]->t[0]) {
n0->conn0[j]->t[0] = add_tangent (this->parent->tangents,
dim, this->t[0]->x, 0);
} } }
for (j=0; jnumconn1; j++) {
if (n0->conn1[j] != this) {
if (!n0->conn1[j]->t[1]) {
n0->conn1[j]->t[1] = add_tangent (this->parent->tangents,
dim, this->t[0]->x, 0);
} } }
} else {
// if this segment is tiny, then just plug it in straight
// if its half, then just take the mean
// if its a lot more, then, um,
for (i=0; it[0]->x[i] = thisrad*x[i] + (1.-thisrad)*dlnorm[i];
(void) vec_normalize (this->t[0]->x, this->t[0]->x, dim);
}
}
}
// is there one on this side?
if (!this->t[1]) {
this->t[1] = add_tangent (this->parent->tangents, dim, dl, 0);
// if this is an end, then mirror the other tangent
if (n1->numconn0 + n1->numconn1 == 1) {
do_later1 = TRUE;
// otherwise, use summations
} else {
// now, loop over connected segments, counting this one
for (i=0; inumconn0; j++) {
// sum the radius-weighted, normalized tangent vectors
for (i=0; iconn0[j]->n[1]->x[i] -
n1->conn0[j]->n[0]->x[i];
(void) vec_normalize (other, other, dim);
if (n1->conn0[j]->r[0])
otherrad = pow (n1->conn0[j]->r[0]->r, 2);
else
otherrad = pow (this->parent->radius, 2);
for (i=0; inumconn1; j++) {
for (i=0; iconn1[j]->n[1]->x[i] -
n1->conn1[j]->n[0]->x[i];
(void) vec_normalize (other, other, dim);
if (n1->conn1[j]->r[1])
otherrad = pow (n1->conn1[j]->r[1]->r, 2);
else
otherrad = pow (this->parent->radius, 2);
for (i=0; ir[1]) thisrad = pow (this->r[1]->r, 2);
else thisrad = pow (this->parent->radius, 2);
// how important is this segment compared to all other segs?
thisrad = 2.*thisrad/radsum;
// finally, if this node just has 2 segments, use the same
// tangent on the other end
if (n1->numconn0 + n1->numconn1 == 2) {
// just use the weighted norm
for (i=0; it[1]->x[i] = x[i];
(void) vec_normalize (this->t[1]->x, this->t[1]->x, dim);
// and copy it to the other segment
for (j=0; jnumconn0; j++) {
if (n1->conn0[j] != this) {
if (!n1->conn0[j]->t[0]) {
n1->conn0[j]->t[0] = add_tangent (this->parent->tangents,
dim, this->t[1]->x, 0);
} } }
for (j=0; jnumconn1; j++) {
if (n1->conn1[j] != this) {
if (!n1->conn1[j]->t[1]) {
n1->conn1[j]->t[1] = add_tangent (this->parent->tangents,
dim, this->t[1]->x, 0);
} } }
} else {
// if this segment is tiny, then just plug it in straight
// if its half, then just take the mean
// if its a lot more, then, um,
for (i=0; it[1]->x[i] = thisrad*x[i] + (1.-thisrad)*dlnorm[i];
(void) vec_normalize (this->t[1]->x, this->t[1]->x, dim);
}
}
}
// if both nodes are ends, then the segment is straight
if (do_later0 && do_later1) {
for (i=0; it[0]->x[i] = dl[i];
for (i=0; it[1]->x[i] = this->t[0]->x[i];
//fprintf(stderr," special case 1\n");
// if just one is an end, copy the other end's tangent, but mirror it
} else if (do_later0) {
// use the mirror image of tangent 1
dotp = 0.;
for (i=0; it[1]->x[i];
for (i=0; it[0]->x[i] = 2.*dlnorm[i]*dotp - this->t[1]->x[i];
//fprintf(stderr," special case 2\n");
} else if (do_later1) {
// use the mirror image of tangent 0
dotp = 0.;
for (i=0; it[0]->x[i];
for (i=0; it[1]->x[i] = 2.*dlnorm[i]*dotp - this->t[0]->x[i];
//fprintf(stderr," special case 3\n");
}
return;
}
/*
* Perturb all nodes according to Gaussian, scale by mean segment length
*/
int roughen_nodes (node_group_ptr this, int dim, double scale) {
int i;
int cnt = 0;
double thislen;
double len = 9.9e+9;
double pert[dim];
node_ptr curr;
// check all children
if (this->child[0]) {
cnt += roughen_nodes (this->child[0], dim, scale);
cnt += roughen_nodes (this->child[1], dim, scale);
} else {
// check all local nodes
curr = this->first;
while (curr) {
// check all connected segments, find shortest length
for (i=0; inumconn0; i++) {
thislen = seg_length (curr->conn0[i],dim);
if (thislen < len) len = thislen;
}
for (i=0; inumconn1; i++) {
thislen = seg_length (curr->conn1[i],dim);
if (thislen < len) len = thislen;
}
(void) perturb_gaussian (pert, scale*len, dim);
for (i=0; ix[i] += pert[i];
curr = curr->next;
}
}
return (cnt);
}
/*
* Set the block id on each segment, a block is a set of connected segs
*/
int identify_separate_strands (seg_group_ptr this) {
int i;
int cnt = 0;
int blockcnt = 0;
seg_ptr curr;
// reset the counters
this->numblock = 0;
curr = this->first;
while (curr) {
curr->flag = FALSE;
curr->block = 0;
curr = curr->next;
}
// check all elements
curr = this->first;
while (curr) {
// has this segment been done yet?
if (!curr->flag) {
// reset the counter
cnt = 0;
// set the ID number for this block
blockcnt++;
// recursively look through the segments, setting the block and flag
cnt = identify_connected_segments (curr, blockcnt);
// what happened?
//fprintf (stderr,"block %d has %d segments\n",blockcnt,cnt);
}
curr = curr->next;
}
this->numblock = blockcnt;
return (blockcnt);
}
int identify_connected_segments (seg_ptr, unsigned int);
/*
* Recursively look at all connected segments! Set block
*/
int identify_connected_segments (seg_ptr this, unsigned int blockid) {
int i;
int cnt = 1;
// if this segment has been done, return immediately
if (this->flag) return (0);
// set ID
this->block = blockid;
// set flag
this->flag = TRUE;
// recursively call all connected elems
for (i=0; in[0]->numconn0; i++)
cnt += identify_connected_segments (this->n[0]->conn0[i], blockid);
for (i=0; in[0]->numconn1; i++)
cnt += identify_connected_segments (this->n[0]->conn1[i], blockid);
for (i=0; in[1]->numconn0; i++)
cnt += identify_connected_segments (this->n[1]->conn0[i], blockid);
for (i=0; in[1]->numconn1; i++)
cnt += identify_connected_segments (this->n[1]->conn1[i], blockid);
return (cnt);
}
/*
* Set the flag for each node that is a root of its strand!
*
* Make sure each strand has *at least* one root. It may have more.
*
* A strand will have one root for each segment that *crosses*
* the root plane. If no segments cross the root plane, it will pick
* the node that is closest to the root plane.
*
* NO! To make the solution possible, we must only have one root per strand!
*/
int find_root_of_each_strand (seg_group_ptr this, ACTION *args) {
int i;
int cnt = 0;
int axis = args->iarg[0];
double intercept = args->darg[1];
double dist;
unsigned char *found_root;
double *close_dist;
node_ptr *close_node;
seg_ptr curr;
// set appropriate flags
(void) set_node_flags (this->nodes, FALSE);
// create array for each block
found_root = (unsigned char*) malloc ((int)(this->numblock+1)*sizeof(char));
close_dist = (double*) malloc ((int)(this->numblock+1) * sizeof(double));
close_node = (node_ptr*) malloc ((int)(this->numblock+1) * sizeof(node_ptr));
// initialize those arrays
for (i=0; inumblock+1; i++) {
found_root[i] = FALSE;
close_dist[i] = 9.9e+9;
close_node[i] = NULL;
}
// check all elements
curr = this->first;
while (curr) {
// have we found a root for this strand/block?
if (!found_root[curr->block]) {
// does this segment intersect the root plane?
if ((curr->n[0]->x[axis] - intercept)*
(curr->n[1]->x[axis] - intercept) < 0.) {
found_root[curr->block] = TRUE;
// one of the two nodes will be the root
if (fabs(curr->n[0]->x[axis] - intercept) >
fabs(curr->n[1]->x[axis] - intercept)) {
curr->n[0]->flag = TRUE;
} else {
curr->n[1]->flag = TRUE;
}
cnt++;
// if not, how close to the root plane does it get?
} else {
dist = fabs(curr->n[0]->x[axis] - intercept);
if (dist < close_dist[curr->block]) {
close_dist[curr->block] = dist;
close_node[curr->block] = curr->n[0];
}
dist = fabs(curr->n[1]->x[axis] - intercept);
if (dist < close_dist[curr->block]) {
close_dist[curr->block] = dist;
close_node[curr->block] = curr->n[1];
}
}
}
curr = curr->next;
}
// check for correctness of results
for (i=1; inumblock+1; i++) {
// if there's no intersection, there's no root yet
if (!found_root[i]) {
if (close_node[i]) {
close_node[i]->flag = TRUE;
//fprintf (stderr,"seg %d closest node at %g\n",i,close_dist[i]);
cnt++;
} else {
fprintf (stderr,"ERROR (f_r_o_e_s): no closest node found\n");
exit(0);
}
}
}
return (cnt);
}
/*
* Set as many segment radii as possible
*
* Only works for dimensions 2 and 3 (or first 3 dimensions of d>3)
*
* A node flag of TRUE means that that node is a root!
*/
int set_treelike_radii_3d (seg_group_ptr this, ACTION *args) {
int i,j,n0conn,n1conn,rindex,tindex;
int cnt = 0;
double lensq,armsq,rad,thisforce,dl[3];
double **moment;
double **force;
double stress = args->darg[0];
int axis = args->iarg[0];
seg_ptr curr;
node_ptr rootnode,tipnode;
// allocate space for a moment on each node
moment = (double**) malloc ((int)this->nodes->num * sizeof(double*));
for (i=0; inodes->num; i++)
moment[i] = (double*) malloc ((int)this->dim * sizeof(double));
force = (double**) malloc ((int)this->nodes->num * sizeof(double*));
for (i=0; inodes->num; i++)
force[i] = (double*) malloc ((int)this->dim * sizeof(double));
// zero those moments
for (i=0; inodes->num; i++) {
moment[i][0] = 0.;
moment[i][1] = 0.;
moment[i][2] = 0.;
force[i][0] = 0.;
force[i][1] = 0.;
force[i][2] = 0.;
}
// check all elements
curr = this->first;
while (curr) {
// has this segment been done yet?
if (!curr->flag) {
// how many connected segments does it have
n0conn = curr->n[0]->numconn0 + curr->n[0]->numconn1;
n1conn = curr->n[1]->numconn0 + curr->n[1]->numconn1;
//fprintf (stderr,"seg %d numconn %d\n",curr->index,numconn);
// it's independent
if (n0conn == 1 && n1conn == 1) {
// it doesn't matter which end is the root
// segment length
lensq = seg_length_sqrd (curr, this->dim);
// moment arm
armsq = lensq - pow(curr->n[1]->x[axis] - curr->n[0]->x[axis],2);
// here's the calculation for a constant-radius segment
// with no tip force or moment
rad = 2.*sqrt(lensq*armsq)/stress;
curr->r[0] = add_radius (this->radii, rad, 0);
curr->r[1] = add_radius (this->radii, rad, 0);
fprintf (stderr," independent, has radius %g\n",rad);
// it's an end, do it
} else if (n0conn == 1 || n1conn == 1) {
// define our directions
if (n1conn == 1) {
rootnode = curr->n[0];
rindex = curr->n[0]->index;
tipnode = curr->n[1];
tindex = curr->n[1]->index;
} else {
rootnode = curr->n[1];
rindex = curr->n[1]->index;
tipnode = curr->n[0];
tindex = curr->n[0]->index;
}
// segment length
for (j=0; j<3; j++) dl[j] = tipnode->x[j] - rootnode->x[j];
lensq = dl[0]*dl[0] + dl[1]*dl[1] + dl[2]*dl[2];
// moment arm
armsq = lensq - pow(curr->n[1]->x[axis] - curr->n[0]->x[axis],2);
// here's the calculation for a constant-radius segment
// with no tip force or moment
rad = 2.*sqrt(lensq*armsq)/stress;
curr->r[0] = add_radius (this->radii, rad, 0);
curr->r[1] = add_radius (this->radii, rad, 0);
fprintf (stderr,"len %g, arm\n",sqrt(lensq),sqrt(armsq));
fprintf (stderr," tip has radius %g\n",rad);
// set the force and moment for this tip segment
force[rindex][axis] = M_PI*sqrt(lensq)*rad*rad;
vec_cross (dl, force[rindex], moment[rindex]);
fprintf (stderr," force %g\n",force[rindex][2]);
fprintf (stderr," moment %g %g\n",moment[rindex][0],moment[rindex][1]);
// it's not an end, check all connected segs for doneness
} else {
// only continue if ONE connected adjacent seg is not done
// define our directions
if (n1conn == 1) {
rootnode = curr->n[0];
rindex = curr->n[0]->index;
tipnode = curr->n[1];
tindex = curr->n[1]->index;
} else {
rootnode = curr->n[1];
rindex = curr->n[1]->index;
tipnode = curr->n[0];
tindex = curr->n[0]->index;
}
// carry the force from the tipward segments
//for (j=0; j<3; j++)
//force[rootnode->index][j] += force[tipnode->index][j];
}
// for testing, force all segments to be done
curr->flag = TRUE;
}
curr = curr->next;
}
return (cnt);
}
/*
* Translate all nodes by a fixed vector
*/
int translate_nodes (node_group_ptr this, int dim, double *vec) {
int i;
int cnt = 0;
node_ptr curr;
// check all children
if (this->child[0]) {
cnt += translate_nodes (this->child[0], dim, vec);
cnt += translate_nodes (this->child[1], dim, vec);
} else {
// check all local nodes
curr = this->first;
while (curr) {
for (i=0; ix[i] += vec[i];
curr = curr->next;
}
}
return (cnt);
}
/*
* Scale all nodes by a fixed vector
*/
int scale_nodes (node_group_ptr this, int dim, double *vec) {
int i;
int cnt = 0;
node_ptr curr;
// check all children
if (this->child[0]) {
cnt += scale_nodes (this->child[0], dim, vec);
cnt += scale_nodes (this->child[1], dim, vec);
} else {
// check all local nodes
curr = this->first;
while (curr) {
for (i=0; ix[i] *= vec[i];
curr = curr->next;
}
}
// finally, scale the min/max bounds
for (i=0; imin[i] *= vec[i];
for (i=0; imax[i] *= vec[i];
return (cnt);
}
/*
* Scale all nodes by a fixed vector
*/
int split_and_write_rad (seg_group_ptr this, char root[MAXSTR]) {
int i;
int long_axis = -1;
double long_length = -1.;
double this_length = -1.;
char leftroot[MAXSTR],rightroot[MAXSTR];
FILE *outptr = stdout;
// two child segment groups
seg_group_ptr left,right;
//fprintf(stderr,"splitting (%d)\n",this->num);
fprintf(stderr,"splitting (%s)\n",root);
// is this group small enough to just print?
// construct the filename
//if (TRUE) {
if (this->num < 100000) {
sprintf(leftroot,"%s.rad",root);
outptr = fopen(leftroot,"w");
write_rad(outptr,this);
fclose(outptr);
// don't delete the structure here because it may be the original segs
return(-1);
}
// if we're still here, then we need to split this seg group into two child seg groups
// prepare both groups
left = make_seg_structure (this->dim);
right = make_seg_structure (this->dim);
//dump_stats(this);
// find the longest edge
long_length = -1.;
for (i=0; idim; i++) {
this_length = pow(this->nodes->max[i]-this->nodes->min[i],2);
if (this_length > long_length) {
long_length = this_length;
long_axis = i;
}
}
long_length = this->nodes->min[long_axis] + 0.5*sqrt(long_length);
fprintf(stderr," axis %d, split position %g\n",long_axis,long_length);
// recurse over all segments in "this" and add them to left, right, or both
split_into_two_seg_groups (this->nodes, long_axis, long_length, left, right);
fprintf(stderr," left segments (%d)\n",left->num);
fprintf(stderr," right segments (%d)\n",right->num);
//exit(0);
// now, call this routine on both
// and remove the structure when we're done with it!
sprintf(leftroot,"%sl",root);
(void) split_and_write_rad(left, leftroot);
(void) free_segs (left);
sprintf(rightroot,"%sr",root);
(void) split_and_write_rad(right, rightroot);
(void) free_segs (right);
return(0);
}
/*
* Put all segments, nodes, and tangents into either left, right, or both;
* cnt is number of segments put on the left.
*/
int split_into_two_seg_groups (node_group_ptr this, int split_dim, double split_val,
seg_group_ptr left, seg_group_ptr right) {
int i;
int cnt = 0;
node_ptr curr,n0,n1;
seg_ptr this_seg,newseg;
rad_ptr r0,r1;
//tan_ptr t0,t1;
enum { leftonly,
rightonly,
split } which_side;
// check all children
if (this->child[0]) {
cnt += split_into_two_seg_groups (this->child[0], split_dim, split_val, left, right);
cnt += split_into_two_seg_groups (this->child[1], split_dim, split_val, left, right);
} else {
// check all local nodes
curr = this->first;
while (curr) {
// loop over all connected segments for which this is node 0
for (i=0; inumconn0; i++) {
this_seg = curr->conn0[i];
// figure out which side to put it on
// for now, make it totally simple
if (curr->x[split_dim] < split_val) {
which_side = leftonly;
} else {
which_side = rightonly;
}
// now, move it there
if (which_side == leftonly) {
n0 = add_node (left->nodes, left->dim, curr->x, 0);
n1 = add_node (left->nodes, left->dim, this_seg->n[1]->x, 0);
r0 = add_radius (left->radii, this_seg->r[0]->r, 0);
r1 = add_radius (left->radii, this_seg->r[1]->r, 0);
//t0 = add_tangent (left->tangents, this_seg->t[0]->x, 0);
//t1 = add_tangent (left->tangents, this_seg->t[1]->x, 0);
newseg = add_segment (left, n0, n1);
if (r0) newseg->r[0] = r0;
if (r1) newseg->r[1] = r1;
//if (t0) newseg->t[0] = t0;
//if (t1) newseg->t[1] = t1;
} else if (which_side == rightonly) {
n0 = add_node (right->nodes, right->dim, curr->x, 0);
n1 = add_node (right->nodes, right->dim, this_seg->n[1]->x, 0);
r0 = add_radius (right->radii, this_seg->r[0]->r, 0);
r1 = add_radius (right->radii, this_seg->r[1]->r, 0);
newseg = add_segment (right, n0, n1);
if (r0) newseg->r[0] = r0;
if (r1) newseg->r[1] = r1;
} else {
fprintf(stderr,"ERROR (split_into_two_seg_groups): cannot split\n");
exit(0);
}
}
curr = curr->next;
}
}
return (cnt);
}
// ==========================================================================
// input/output
//
/*
* Read a .seg file (syntax very much like .obj)
*/
int read_seg (char *infile, seg_group_ptr this) {
int i,j,k,nnode,nrad,ntan,nlines;
int n0index,n1index,n0rad,n1rad,n0tan,n1tan;
char onechar,anotherchar;
char twochar[2],sbuf[128],newval[32],sub[10];
double tempd[MAXDIM];
node_ptr newnode,newnode0,newnode1;
node_ptr *thenodes;
rad_ptr newrad;
rad_ptr *therads;
tan_ptr newtan;
tan_ptr *thetans;
seg_ptr newseg;
FILE *infp;
long int fpos;
// before proceeding, prepare the seg_group struct
this->dim = 3; // provide a default
this->radius = -1.; // negative == "unset"
this->num = 0; // gotta start somewhere
// first, scan the file and record the dimensionality and bounds =======
// open the file for reading
infp = fopen (infile,"r");
if (infp==NULL) {
fprintf (stderr,"Could not open input file %s\n",infile);
fprintf (stderr,"Exiting\n");
exit (0);
}
fprintf (stderr,"Opening file %s\n",infile);
fprintf (stderr,"Prescanning");
fflush (stderr);
// read the first character after the newline
nnode = 0;
nrad = 0;
nlines = 0;
while (fread (&onechar,sizeof(char),1,infp) == 1) {
if (onechar == 'v') {
// this is some form of vertex information
fread (&anotherchar,sizeof(char),1,infp);
if (isspace (anotherchar)) {
// this is an actual vertex
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
// count it
nnode++;
} else if (anotherchar == 'r') {
// this is a radius
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
// count it
nrad++;
} else if (anotherchar == 'n') {
// this is a tangent
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
// count it
ntan++;
} else {
// skip it
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
}
} else {
// if its not identifiable, skip it
fscanf (infp,"%[^\n]",sbuf); // read comment beyond '#'
fscanf (infp,"%[\n]",twochar); // read newline
}
if (++nlines%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
}
// close it
fclose(infp);
fprintf(stderr,"%d nodes, %d radii\n",nnode,nrad);
fflush(stderr);
// malloc space for an array of pointers
thenodes = (node_ptr*) malloc ((nnode+1) * sizeof(node_ptr));
therads = (rad_ptr*) malloc ((nrad+1) * sizeof(rad_ptr));
thetans = (tan_ptr*) malloc ((ntan+1) * sizeof(tan_ptr));
// then, reopen the file and read in the vertexes and segments =========
// open the file for reading
infp = fopen (infile,"r");
if (infp==NULL) {
fprintf (stderr,"Could not open input file %s\n",infile);
fprintf (stderr,"Exiting\n");
exit (0);
}
//fprintf (stderr,"Opening file %s\n",infile);
fprintf (stderr,"Reading");
fflush (stderr);
// read the first character after the newline
nnode = 0;
nrad = 0;
nlines = 0;
while (fread (&onechar,sizeof(char),1,infp) == 1) {
// split on first character
if (onechar == '#') {
// read a comment line
fscanf (infp,"%[^\n]",sbuf); // read comment beyond '#'
fscanf (infp,"%[\n]",twochar); // read newline
} else if (onechar == 'g') {
// this is some form of global variable
fread (&anotherchar,sizeof(char),1,infp);
if (anotherchar == 'r') {
// the global radius
fscanf (infp,"%s",newval);
this->radius = atof(newval);
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
//fprintf (stderr,"set radius to %g\n",this->radius);
} else {
// nothing important
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
}
} else if (onechar == 'd') {
// the dimensionality
fscanf (infp,"%s",newval);
this->dim = (char)atoi(newval);
if (this->dim > MAXDIM) {
fprintf (stderr,"ERROR (read_seg): input file's dim is %d, but\n",
this->dim);
fprintf (stderr," the maximum spatial dimensions allowed is %d.\n",
(MAXDIM));
fprintf (stderr," Quitting.\n");
exit(1);
}
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
//fprintf (stderr,"set dimensionality to %d\n",(int)this->dim);
} else if (onechar == 'v') {
// this is some form of vertex information
fread (&anotherchar,sizeof(char),1,infp);
if (isspace(anotherchar)) {
// this is an actual vertex
//fprintf (stderr,"adding node\n"); fflush (stderr);
for (i=0; idim; i++) {
fpos = ftell (infp);
fscanf (infp,"%s",newval);
// check this value for \n or non-math character!
//fprintf (stderr," (%s)\n",newval); fflush (stderr);
if (newval[0] < 43 || newval[0] == 47 | newval[0] > 57) {
// if it is, fill the rest of the data with zeros
for (j=i; jdim; j++) tempd[j] = 0.;
// rewind the file
fseek (infp, fpos, SEEK_SET);
// and break out of this loop
break;
} else {
tempd[i] = atof(newval);
}
}
newnode = add_node (this->nodes, this->dim, tempd, 0);
//for (i=0; idim; i++) newnode->x[i] = tempd[i];
//fprintf (stderr," node index %ld\n",newnode->index); fflush (stderr);
//fprintf (stderr," node loc %g %g\n",newnode->x[0],newnode->x[1]); fflush (stderr);
// put it in the list
//thenodes[this->nodes->num] = newnode;
thenodes[++nnode] = newnode;
// and finish reading
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
} else if (anotherchar == 'r') {
// this is a radius
//fprintf (stderr,"adding radius\n"); fflush (stderr);
fscanf (infp,"%s",newval);
newrad = add_radius (this->radii, atof(newval), 0);
therads[++nrad] = newrad;
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
} else if (anotherchar == 'n') {
// this is a tangent
for (i=0; idim; i++) {
fscanf (infp,"%s",newval);
tempd[i] = atof(newval);
}
//newtan = add_tangent (this->tangents, atof(newval), 0);
//thetans[++ntan] = newtan;
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
} else {
// skip it
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
}
} else if (onechar == 's') {
// this is a segment
//fprintf (stderr,"adding segment\n"); fflush (stderr);
n0index = -1;
n1index = -1;
n0rad = -1;
n1rad = -1;
n0tan = -1;
n1tan = -1;
// read the first set of indices
fscanf (infp,"%s",newval);
// look for a slash
for (i=0; inode indicies incorrect\n");
fprintf (stderr," n0index=%d, n1index=%d\n",n0index,n1index);
fprintf (stderr," Quitting.\n");
exit(1);
}
if (n0index > nnode || n1index > nnode) {
fprintf (stderr,"ERROR (read_seg): segment->node indicies incorrect\n");
fprintf (stderr," n0index=%d, n1index=%d\n",n0index,n1index);
fprintf (stderr," only %d nodes have been defined so far\n",nnode);
fprintf (stderr," Quitting.\n");
exit(1);
}
if (n0rad == 0 || n1rad == 0) {
fprintf (stderr,"ERROR (read_seg): segment->rad indicies incorrect\n");
fprintf (stderr," n0rad=%d, n1rad=%d\n",n0rad,n1rad);
fprintf (stderr," Quitting.\n");
exit(1);
}
if (n0rad > nrad || n1rad > nrad) {
fprintf (stderr,"ERROR (read_seg): segment->rad indicies incorrect\n");
fprintf (stderr," n0rad=%d, n1rad=%d\n",n0rad,n1rad);
fprintf (stderr," only %d rads have been defined so far\n",nrad);
fprintf (stderr," Quitting.\n");
exit(1);
}
// create the segment
newseg = add_segment (this, thenodes[n0index], thenodes[n1index]);
// assign the radius and texture/tangent
if (n0rad > -1) newseg->r[0] = therads[n0rad];
if (n1rad > -1) newseg->r[1] = therads[n1rad];
if (n0tan > -1) newseg->t[0] = thetans[n0tan];
if (n1tan > -1) newseg->t[1] = thetans[n1tan];
// finish reading the line
fscanf (infp,"%[^\n]",sbuf); // read comment beyond '#'
fscanf (infp,"%[\n]",twochar); // read newline
} else {
// if its not identifiable, skip it
fscanf (infp,"%[^\n]",sbuf); // read comment beyond '#'
fscanf (infp,"%[\n]",twochar); // read newline
}
if (++nlines%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
}
// close it
fclose (infp);
fprintf (stderr,"%d segments\n",this->num);
fflush (stderr);
// remove the temporary array
free (thenodes);
// if things bomb, send back a nonzero
return (0);
}
/*
* Read a .rad file (detects Radiance cones and cylinders only)
*
* For now, naive read: will maintain copies of nodes
*/
int read_rad (char *infile, seg_group_ptr this) {
int i,nlines,nfloats;
char onechar,anotherchar;
char twochar[2],sbuf[512],token[14][32];
double tempd[8];
node_ptr newnode0,newnode1;
rad_ptr newrad;
seg_ptr newseg;
FILE *infp;
// before proceeding, prepare the seg_group struct
this->dim = 3; // Radiance files are always 3D
this->radius = -1.; // negative == "unset"
this->num = 0; // gotta start somewhere
// first, scan the file and record the dimensionality and bounds =======
// no need!
// then, reopen the file and read in the vertexes and segments =========
// open the file for reading
infp = fopen (infile,"r");
if (infp==NULL) {
fprintf (stderr,"Could not open input file %s\n",infile);
fprintf (stderr,"Exiting\n");
exit (0);
}
//fprintf (stderr,"Opening file %s\n",infile);
fprintf (stderr,"Reading");
fflush (stderr);
// read a whole line from the input file
nlines = 0;
//while (fscanf(infp,"%[^\n]",sbuf) != EOF) {
// naw, read only the first three tokens
while (fscanf(infp,"%s %s %s",token[0],token[1],token[2]) != EOF) {
// clear out the token buffers!
//for (i=0; i<3; i++) token[i][0] = '\0';
// grab the first three words (mat key name)
//sscanf(sbuf,"%s %s %s",token[0],token[1],token[2]);
//fscanf(infp,"%s %s %s",token[0],token[1],token[2]);
// grab the first six words (mat key name nc ni nf)
//sscanf(sbuf,"%s %s %s %s %s %s",token[0],token[1],token[2],
// token[3],token[4],token[5]);
//sscanf(sbuf,"%s %s %s %s %s %s %s %s %s %s %s %s %s %s",token[0],
// token[1],token[2],token[3],token[4],token[5],token[6],token[7],
// token[8],token[9],token[10],token[11],token[12],token[13]);
//fprintf (stdout," first token is %s\n",token[0]); fflush(stdout);
//for (i=0; i<3; i++) fprintf (stdout,"%d %s\n",i,token[i]); fflush(stdout);
//exit(0);
// split on those three tokens only
if (token[0][0] == '#' || token[0][0] == '!') {
// read a comment line, or some other descriptor
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read up to newline
// fprintf (stdout,"%s\n",sbuf); // write comment
} else if (strncmp(token[1],"cylinder",8) == 0) {
// read in and process a cylinder
// read in the next 7 tokens
//sscanf(sbuf,"%s %s %s %s %s %s %s",token[0],token[1],token[2],
// token[3],token[4],token[5],token[6]);
//sscanf(sbuf,"%s %s %s %s %s %s %s %s %s %s %s %s %s",token[0],
// token[1],token[2],token[3],token[4],token[5],token[6],token[7],
// token[8],token[9],token[10],token[11],token[12]);
// call subroutine to find 7 float values
nfloats = read_radiance_floats (infp, tempd);
if (nfloats != 7) {
fprintf (stderr,"ERROR (read_rad): too many arguments for\n");
fprintf (stderr," cylinder (%d).\n",nfloats);
fprintf (stderr,"Quitting.\n");
exit(1);
}
//for (i=0; i<7; i++) fprintf (stdout,"%d %g\n",i,tempd[i]); fflush(stdout);
// create nodes at the two endpoints
newnode0 = add_node (this->nodes, this->dim, &tempd[0], 0);
//for (i=0; idim; i++) newnode0->x[i] = tempd[i];
newnode1 = add_node (this->nodes, this->dim, &tempd[3], 0);
//for (i=0; idim; i++) newnode1->x[i] = tempd[3+i];
// create the segment
newseg = add_segment (this, newnode0, newnode1);
// and create one radius
newrad = add_radius (this->radii, tempd[6], 0);
newseg->r[0] = newrad;
newseg->r[1] = newrad;
// and finish reading
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
} else if (strncmp(token[1],"cone",4) == 0) {
// read in and process a cone
// read in the next 8 tokens
//sscanf(sbuf,"%s %s %s %s %s %s %s %s",token[0],token[1],token[2],
// token[3],token[4],token[5],token[6],token[7]);
//sscanf(sbuf,"%s %s %s %s %s %s %s %s %s %s %s %s %s %s",token[0],
// token[1],token[2],token[3],token[4],token[5],token[6],token[7],
// token[8],token[9],token[10],token[11],token[12],token[13]);
// call subroutine to find 8 float values
nfloats = read_radiance_floats (infp, tempd);
if (nfloats != 8) {
fprintf (stderr,"ERROR (read_rad): too many arguments for\n");
fprintf (stderr," cone (%d).\n",nfloats);
fprintf (stderr,"Quitting.\n");
exit(1);
}
//for (i=0; i<8; i++) fprintf (stdout,"%d %g\n",i,tempd[i]); fflush(stdout);
// create nodes at the two endpoints
newnode0 = add_node (this->nodes, this->dim, &tempd[0], 0);
//for (i=0; idim; i++) newnode0->x[i] = tempd[i];
newnode1 = add_node (this->nodes, this->dim, &tempd[3], 0);
//for (i=0; idim; i++) newnode1->x[i] = tempd[3+i];
// create the segment
newseg = add_segment (this, newnode0, newnode1);
// and create two radii
newrad = add_radius (this->radii, tempd[6], 0);
newseg->r[0] = newrad;
newrad = add_radius (this->radii, tempd[7], 0);
newseg->r[1] = newrad;
// and finish reading
fscanf (infp,"%[^\n]",sbuf); // read line up to newline
fscanf (infp,"%[\n]",twochar); // read newline
} else {
// ignore anything else
fscanf (infp,"%[^\n]",sbuf); // read comment beyond '#'
fscanf (infp,"%[\n]",twochar); // read newline
}
// write the counter
if (++nlines%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
// clear out the token buffers!
for (i=0; i<3; i++) token[i][0] = '\0';
}
// close it
fclose (infp);
fprintf (stderr,"%d segments\n",this->num);
fflush (stderr);
// if things bomb, send back a nonzero
return (0);
}
/*
* Parse a Radiance file, passing over whitespace, to get to the
* floating-point parameters
*/
int read_radiance_floats (FILE *infp, double *x) {
char token[32];
int num = 0;
int i;
// read the character args
(void) fscanf(infp,"%s",token);
num = atoi(token);
for (i=0; inext = this->first;
newseg->prev = NULL;
if (newseg->next) newseg->next->prev = newseg;
this->first = newseg;
newseg->parent = this;
//if (n0) fprintf (stderr, " node 0 exists, index %d\n",n0->index);
//else fprintf (stderr, " node 0 does not exist, why?\n");
//if (n1) fprintf (stderr, " node 1 exists, index %d\n",n1->index);
//else fprintf (stderr, " node 1 does not exist, why?\n");
// null some pointers
newseg->n[0] = n0;
newseg->n[1] = n1;
newseg->t[0] = NULL;
newseg->t[1] = NULL;
newseg->r[0] = NULL;
newseg->r[1] = NULL;
// first node is index 1
newseg->index = ++this->num;
//fprintf (stderr,"\ncreated segment %ld\n",newseg->index);
// set the nodes' segment pointers
//fprintf (stderr," n0 had %d conn0\n",n0->numconn0);
n0->numconn0++;
// alloc space for new conn list
newconnlist = (seg_ptr*) malloc (n0->numconn0 * sizeof(seg_ptr));
// copy data from old conn list
for (i=0; inumconn0-1; i++) {
newconnlist[i] = n0->conn0[i];
//fprintf (stderr," one is %ld\n",n0->conn0[i]->index);
}
newconnlist[n0->numconn0-1] = newseg;
// delete old list
free (n0->conn0);
// copy new list
n0->conn0 = newconnlist;
//fprintf (stderr," now has %d conn0\n",n0->numconn0);
//fprintf (stderr," new one is %ld\n",n0->conn0[n0->numconn0-1]->index);
// now do it for node 1
//fprintf (stderr," n1 had %d conn1\n",n1->numconn1);
n1->numconn1++;
// alloc space for new conn list
newconnlist = (seg_ptr*) malloc (n1->numconn1 * sizeof(seg_ptr));
// copy data from old conn list
for (i=0; inumconn1-1; i++) {
newconnlist[i] = n1->conn1[i];
//fprintf (stderr," one is %ld\n",n1->conn1[i]->index);
}
newconnlist[n1->numconn1-1] = newseg;
// delete old list
free (n1->conn1);
// copy new list
n1->conn1 = newconnlist;
// return the pointer
return (newseg);
}
/*
* Generic procedure for creating a new node in a node group
*
* Does not do any tree stuff yet
*/
/*
node_ptr add_node (node_group_ptr this, unsigned char dim) {
int i = 0;
node_ptr newnode = NULL;
// if there no nodes, set the bounds
if (!this->first) {
this->min = (double*) malloc ((int)dim * sizeof(double));
this->max = (double*) malloc ((int)dim * sizeof(double));
for (i=0; imin[i] = 9.9e+9;
for (i=0; imax[i] = -9.9e+9;
}
// get memory for the node
newnode = (NODE*) malloc (sizeof(NODE));
// add it to the head of the list
newnode->next = this->first;
newnode->prev = NULL;
if (newnode->next) newnode->next->prev = newnode;
this->first = newnode;
// assign some shit
newnode->numconn0 = 0;
newnode->conn0 = NULL;
newnode->numconn1 = 0;
newnode->conn1 = NULL;
newnode->x = (double*) malloc ((int)dim * sizeof(double));
// first node is index 1
newnode->index = ++this->num;
// return the pointer
return (newnode);
}
*/
/*
* Generic procedure for creating a new node in a node group
*/
node_ptr add_node (node_group_ptr this, unsigned char dim, double *x,
int totnode) {
int i = 0;
int axis;
double dist;
node_ptr newnode = NULL;
node_ptr curr = NULL;
node_ptr currnext = NULL;
// set the index of the new entry, if we use it
if (totnode == 0) totnode = this->num;
// first, see if we have this entry already
if (this->child[0]) {
// there are children, check in the appropriate one
if (2.*x[this->axis] <
this->child[0]->max[this->axis] + this->child[1]->min[this->axis]) {
newnode = add_node (this->child[0], dim, x, totnode);
// reset the bounds
for (i=0; ichild[0]->min[i] < this->min[i])
this->min[i] = this->child[0]->min[i];
if (this->child[0]->max[i] > this->max[i])
this->max[i] = this->child[0]->max[i];
}
} else {
newnode = add_node (this->child[1], dim, x, totnode);
// reset the bounds
for (i=0; ichild[1]->min[i] < this->min[i])
this->min[i] = this->child[1]->min[i];
if (this->child[1]->max[i] > this->max[i])
this->max[i] = this->child[1]->max[i];
}
}
// reset the count
this->num = this->child[0]->num + this->child[1]->num;
} else {
// there are no children, check all local entries
curr = this->first;
while (curr) {
//fprintf (stderr," checking vs %g\n",curr->r);
dist = 0;
for (i=0; ix[i],2);
if (fabs(dist) < EPSILONSQRD) {
// these are the same node!
newnode = curr;
//fprintf (stderr," found match %g\n",curr->r);
}
curr = curr->next;
}
}
// if we couldn't find a match, make a new entry
if (!newnode) {
//fprintf (stderr," no match, creating %g\n",rad);
// get memory for the node
newnode = (NODE*) malloc (sizeof(NODE));
// assign some shit
newnode->numconn0 = 0;
newnode->conn0 = NULL;
newnode->numconn1 = 0;
newnode->conn1 = NULL;
newnode->x = (double*) malloc ((int)dim * sizeof(double));
// set the value
for (i=0; ix[i] = x[i];
// if this is the first node in this group, malloc the space
if (!this->first) {
this->axis = -1;
this->min = (double*) malloc ((int)dim * sizeof(double));
this->max = (double*) malloc ((int)dim * sizeof(double));
for (i=0; imin[i] = x[i];
for (i=0; imax[i] = x[i];
}
// add it to the head of the list
newnode->next = this->first;
this->first = newnode;
newnode->prev = NULL;
if (newnode->next) newnode->next->prev = newnode;
newnode->parent = this;
// set the index
newnode->index = totnode + 1;
// and increment the counter for this group
this->num++;
// reset the bounds
for (i=0; imin[i]) this->min[i] = x[i];
if (x[i] > this->max[i]) this->max[i] = x[i];
}
// finally, if we surpassed the bucket size, split this group!
if (this->num > BUCKET) {
// make two children
this->child[0] = (NODE_GROUP*) malloc (sizeof(NODE_GROUP));
this->child[1] = (NODE_GROUP*) malloc (sizeof(NODE_GROUP));
// choose which axis to split (use dist as the max size here)
dist = 0.;
for (i=0; imax[i] - this->min[i] > dist) {
dist = this->max[i] - this->min[i];
//fprintf (stderr," %d %g\n",i,dist);
axis = i;
}
this->axis = axis;
// dump a sentence
//fprintf (stderr,"\noriginal: num=%d, bounds %g %g, split axis %d\n",
// this->num,this->min[axis],this->max[axis],axis);
// children point to parent
this->child[0]->num = 0;
this->child[0]->first = NULL;
this->child[0]->parent = this;
this->child[0]->child[0] = NULL;
this->child[0]->child[1] = NULL;
this->child[0]->axis = -1;
this->child[0]->min = (double*) malloc ((int)dim * sizeof(double));
this->child[0]->max = (double*) malloc ((int)dim * sizeof(double));
for (i=0; ichild[0]->min[i] = 9.9e+9;
for (i=0; ichild[0]->max[i] = -9.9e+9;
this->child[1]->num = 0;
this->child[1]->first = NULL;
this->child[1]->parent = this;
this->child[1]->child[0] = NULL;
this->child[1]->child[1] = NULL;
this->child[1]->axis = -1;
this->child[1]->min = (double*) malloc ((int)dim * sizeof(double));
this->child[1]->max = (double*) malloc ((int)dim * sizeof(double));
for (i=0; ichild[1]->min[i] = 9.9e+9;
for (i=0; ichild[1]->max[i] = -9.9e+9;
// split the cell along this edge
dist = 0.5 * (this->max[axis] + this->min[axis]);
// move all entries to one of the two children
curr = this->first;
while (curr) {
currnext = curr->next;
if (curr->x[axis] < dist) {
curr->next = this->child[0]->first;
this->child[0]->first = curr;
if (curr->next) curr->next->prev = curr;
curr->prev = NULL;
this->child[0]->num++;
curr->parent = this->child[0];
for (i=0; ix[i] < this->child[0]->min[i])
this->child[0]->min[i] = curr->x[i];
for (i=0; ix[i] > this->child[0]->max[i])
this->child[0]->max[i] = curr->x[i];
} else {
curr->next = this->child[1]->first;
this->child[1]->first = curr;
if (curr->next) curr->next->prev = curr;
curr->prev = NULL;
this->child[1]->num++;
curr->parent = this->child[1];
for (i=0; ix[i] < this->child[1]->min[i])
this->child[1]->min[i] = curr->x[i];
for (i=0; ix[i] > this->child[1]->max[i])
this->child[1]->max[i] = curr->x[i];
}
curr = currnext;
}
// blank out the entry list
this->first = NULL;
// update counter and min/max
//fprintf (stderr," makes children with num=%d, bounds %g %g\n",
// this->child[0]->num,this->child[0]->min[axis],
// this->child[0]->max[axis]);
//fprintf (stderr," and num=%d, bounds %g %g\n",
// this->child[1]->num,this->child[1]->min[axis],
// this->child[1]->max[axis]);
}
}
//if (newnode) {
//fprintf (stderr,"Added a node %d/%d at %g %g %g\n",totnode,newnode->index,x[0],x[1],x[2]);
//fprintf (stderr,"Added node %d\n",newnode->index);
//} else {
// fprintf (stderr,"Did not add a node? %d\n",totnode);
//}
// return the pointer
return (newnode);
}
/*
* Generic procedure for creating a new radius in a radius group
*/
rad_ptr add_radius (rad_group_ptr this, double rad, int totrad) {
double midpt;
rad_ptr newrad = NULL;
rad_ptr curr = NULL;
rad_ptr currnext = NULL;
// set the index of the new entry, if we use it
if (totrad == 0) totrad = this->num;
// first, see if we have this entry already
if (this->child[0]) {
// there are children, check in the appropriate one
if (2.*rad < this->child[0]->max + this->child[1]->min) {
newrad = add_radius (this->child[0], rad, totrad);
// reset the bounds
if (this->child[0]->min < this->min) this->min = this->child[0]->min;
if (this->child[0]->max > this->max) this->max = this->child[0]->max;
} else {
newrad = add_radius (this->child[1], rad, totrad);
// reset the bounds
if (this->child[1]->min < this->min) this->min = this->child[1]->min;
if (this->child[1]->max > this->max) this->max = this->child[1]->max;
}
// reset the count
this->num = this->child[0]->num + this->child[1]->num;
} else {
// there are no children, check all local entries
curr = this->first;
while (curr) {
//fprintf (stderr," checking vs %g\n",curr->r);
if (fabs(rad-curr->r) < EPSILON) {
// these are the same node!
newrad = curr;
//fprintf (stderr," found match %g\n",curr->r);
}
curr = curr->next;
}
}
// if we couldn't find a match, make a new entry
if (!newrad) {
//fprintf (stderr," no match, creating %g\n",rad);
// get memory for the node
newrad = (RADIUS*) malloc (sizeof(RADIUS));
// set the value
newrad->r = rad;
// add it to the head of the list
newrad->next = this->first;
this->first = newrad;
if (newrad->next) newrad->next->prev = newrad;
newrad->prev = NULL;
newrad->parent = this;
// set the index
newrad->index = totrad + 1;
// and increment the counter for this group
this->num++;
// reset the bounds
if (rad < this->min) this->min = rad;
if (rad > this->max) this->max = rad;
// finally, if we surpassed the bucket size, split this group!
if (this->num > BUCKET) {
//fprintf (stderr,"\noriginal box with num=%d, bounds %g %g\n",
// this->num,this->min,this->max);
// make two children
this->child[0] = (RADIUS_GROUP*) malloc (sizeof(RADIUS_GROUP));
this->child[1] = (RADIUS_GROUP*) malloc (sizeof(RADIUS_GROUP));
// children point to parent
this->child[0]->num = 0;
this->child[0]->first = NULL;
this->child[0]->parent = this;
this->child[0]->child[0] = NULL;
this->child[0]->child[1] = NULL;
this->child[0]->min = this->min;
this->child[0]->max = this->min;
this->child[1]->num = 0;
this->child[1]->first = NULL;
this->child[1]->parent = this;
this->child[1]->child[0] = NULL;
this->child[1]->child[1] = NULL;
this->child[1]->min = this->max;
this->child[1]->max = this->max;
// split the cell along this edge
midpt = 0.5 * (this->max + this->min);
// move all entries to one of the two children
curr = this->first;
while (curr) {
currnext = curr->next;
if (curr->r < midpt) {
curr->next = this->child[0]->first;
if (curr->next) curr->next->prev = curr;
curr->prev = NULL;
this->child[0]->first = curr;
this->child[0]->num++;
curr->parent = this->child[0];
if (curr->r > this->child[0]->max) this->child[0]->max = curr->r;
} else {
curr->next = this->child[1]->first;
if (curr->next) curr->next->prev = curr;
curr->prev = NULL;
this->child[1]->first = curr;
this->child[1]->num++;
curr->parent = this->child[1];
if (curr->r < this->child[1]->min) this->child[1]->min = curr->r;
}
curr = currnext;
}
// blank out the entry list
this->first = NULL;
// update counter and min/max
//fprintf (stderr," makes children with num=%d, bounds %g %g\n",
// this->child[0]->num,this->child[0]->min,this->child[0]->max);
//fprintf (stderr," and num=%d, bounds %g %g\n",
// this->child[1]->num,this->child[1]->min,this->child[1]->max);
}
}
// return the pointer
return (newrad);
}
/*
* Generic procedure for creating a new tangent in a tangent group
*/
tan_ptr add_tangent (tan_group_ptr this, unsigned char dim, double *x,
int tottan) {
int i = 0;
int axis;
double dist;
tan_ptr newtan = NULL;
tan_ptr curr = NULL;
tan_ptr currnext = NULL;
// set the index of the new entry, if we use it
if (tottan == 0) tottan = this->num;
// first, see if we have this entry already
if (this->child[0]) {
// there are children, check in the appropriate one
if (2.*x[this->axis] <
this->child[0]->max[this->axis] + this->child[1]->min[this->axis]) {
newtan = add_tangent (this->child[0], dim, x, tottan);
// reset the bounds
for (i=0; ichild[0]->min[i] < this->min[i])
this->min[i] = this->child[0]->min[i];
if (this->child[0]->max[i] > this->max[i])
this->max[i] = this->child[0]->max[i];
}
} else {
newtan = add_tangent (this->child[1], dim, x, tottan);
// reset the bounds
for (i=0; ichild[1]->min[i] < this->min[i])
this->min[i] = this->child[1]->min[i];
if (this->child[1]->max[i] > this->max[i])
this->max[i] = this->child[1]->max[i];
}
}
// reset the count
this->num = this->child[0]->num + this->child[1]->num;
} else {
// there are no children, check all local entries
curr = this->first;
while (curr) {
//fprintf (stderr," checking vs %g\n",curr->r);
dist = 0;
for (i=0; ix[i],2);
if (fabs(dist) < EPSILONSQRD) {
// these are the same tan!
newtan = curr;
//fprintf (stderr," found match %g\n",curr->r);
}
curr = curr->next;
}
}
// if we couldn't find a match, make a new entry
if (!newtan) {
//fprintf (stderr," no match, creating %g\n",rad);
// get memory for the tan
newtan = (TANGENT*) malloc (sizeof(TANGENT));
// assign some shit
newtan->x = (double*) malloc ((int)dim * sizeof(double));
// set the value
for (i=0; ix[i] = x[i];
// if this is the first tan in this group, malloc the space
if (!this->first) {
this->axis = -1;
this->min = (double*) malloc ((int)dim * sizeof(double));
this->max = (double*) malloc ((int)dim * sizeof(double));
for (i=0; imin[i] = x[i];
for (i=0; imax[i] = x[i];
}
// add it to the head of the list
newtan->next = this->first;
this->first = newtan;
newtan->prev = NULL;
if (newtan->next) newtan->next->prev = newtan;
newtan->parent = this;
// set the index
newtan->index = tottan + 1;
// and increment the counter for this group
this->num++;
// reset the bounds
for (i=0; imin[i]) this->min[i] = x[i];
if (x[i] > this->max[i]) this->max[i] = x[i];
}
// finally, if we surpassed the bucket size, split this group!
if (this->num > BUCKET) {
// make two children
this->child[0] = (TANGENT_GROUP*) malloc (sizeof(TANGENT_GROUP));
this->child[1] = (TANGENT_GROUP*) malloc (sizeof(TANGENT_GROUP));
// choose which axis to split (use dist as the max size here)
dist = 0.;
for (i=0; imax[i] - this->min[i] > dist) {
dist = this->max[i] - this->min[i];
//fprintf (stderr," %d %g\n",i,dist);
axis = i;
}
this->axis = axis;
// dump a sentence
//fprintf (stderr,"\noriginal: num=%d, bounds %g %g, split axis %d\n",
// this->num,this->min[axis],this->max[axis],axis);
// children point to parent
this->child[0]->num = 0;
this->child[0]->first = NULL;
this->child[0]->parent = this;
this->child[0]->child[0] = NULL;
this->child[0]->child[1] = NULL;
this->child[0]->axis = -1;
this->child[0]->min = (double*) malloc ((int)dim * sizeof(double));
this->child[0]->max = (double*) malloc ((int)dim * sizeof(double));
for (i=0; ichild[0]->min[i] = 9.9e+9;
for (i=0; ichild[0]->max[i] = -9.9e+9;
this->child[1]->num = 0;
this->child[1]->first = NULL;
this->child[1]->parent = this;
this->child[1]->child[0] = NULL;
this->child[1]->child[1] = NULL;
this->child[1]->axis = -1;
this->child[1]->min = (double*) malloc ((int)dim * sizeof(double));
this->child[1]->max = (double*) malloc ((int)dim * sizeof(double));
for (i=0; ichild[1]->min[i] = 9.9e+9;
for (i=0; ichild[1]->max[i] = -9.9e+9;
// split the cell along this edge
dist = 0.5 * (this->max[axis] + this->min[axis]);
// move all entries to one of the two children
curr = this->first;
while (curr) {
currnext = curr->next;
if (curr->x[axis] < dist) {
curr->next = this->child[0]->first;
this->child[0]->first = curr;
if (curr->next) curr->next->prev = curr;
curr->prev = NULL;
this->child[0]->num++;
curr->parent = this->child[0];
for (i=0; ix[i] < this->child[0]->min[i])
this->child[0]->min[i] = curr->x[i];
for (i=0; ix[i] > this->child[0]->max[i])
this->child[0]->max[i] = curr->x[i];
} else {
curr->next = this->child[1]->first;
this->child[1]->first = curr;
if (curr->next) curr->next->prev = curr;
curr->prev = NULL;
this->child[1]->num++;
curr->parent = this->child[1];
for (i=0; ix[i] < this->child[1]->min[i])
this->child[1]->min[i] = curr->x[i];
for (i=0; ix[i] > this->child[1]->max[i])
this->child[1]->max[i] = curr->x[i];
}
curr = currnext;
}
// blank out the entry list
this->first = NULL;
// update counter and min/max
//fprintf (stderr," makes children with num=%d, bounds %g %g\n",
// this->child[0]->num,this->child[0]->min[axis],
// this->child[0]->max[axis]);
//fprintf (stderr," and num=%d, bounds %g %g\n",
// this->child[1]->num,this->child[1]->min[axis],
// this->child[1]->max[axis]);
}
}
//if (newtan) {
//fprintf (stderr,"Added a tan %d/%d at %g %g %g\n",tottan,newtan->index,x[0],x[1],x[2]);
//fprintf (stderr,"Added tan %d\n",newtan->index);
//} else {
// fprintf (stderr,"Did not add a tan? %d\n",tottan);
//}
// return the pointer
return (newtan);
}
/*
* Search for the closest node to a point, or return NULL if there are
* no nodes within the value of dist given in the initial call.
*/
node_ptr find_closest_node (node_group_ptr this, unsigned char dim, double *x,
node_ptr closenode, double *dist) {
int i;
double thisdistsq;
double distsq = pow(*dist,2);
//node_ptr closenode = NULL;
node_ptr curr = NULL;
// first, can any nodes in this group be close enough?
if (this->first != NULL || this->child[0] != NULL) {
thisdistsq = 0;
for (i=0; imax[i], this->min[i]-x[i]) ) ,2);
// if no nodes can possibly be closer, return
if (thisdistsq > distsq) return (closenode);
}
// if this group could contain a closer node, check the contents
if (this->child[0]) {
// there are children, check both
closenode = find_closest_node (this->child[0], dim, x, closenode, dist);
closenode = find_closest_node (this->child[1], dim, x, closenode, dist);
} else {
// there are no children, check all local entries
curr = this->first;
while (curr) {
thisdistsq = 0;
for (i=0; ix[i],2);
if (thisdistsq < distsq) {
// this is the new closest node!
closenode = curr;
distsq = thisdistsq;
}
curr = curr->next;
}
// reset the distance
*dist = sqrt(distsq);
}
// return the pointer
return (closenode);
}
/*
* Write a .seg file (syntax very much like .obj)
*/
int write_seg (FILE *out, seg_group_ptr this, int argc, char **argv) {
int i,j;
int nnode = 0;
int nrad = 0;
int ntan = 0;
int nlines = 0;
node_ptr currn;
rad_ptr currr;
tan_ptr currt;
seg_ptr curr;
fprintf (stderr,"Writing .seg file");
fflush (stderr);
// set all flags to TRUE (true==not been printed yet)
(void) set_radius_flags (this->radii, TRUE);
(void) set_node_flags (this->nodes, TRUE);
// for this first take, blindly dump all segments and nodes
fprintf (out,"# %s\n",VERSION);
fprintf (out,"#");
for (i=0; idim);
fprintf (out,"gr %g\n",this->radius);
// first, the nodes
/*
currn = this->nodes->first;
while (currn) {
fprintf (out,"v");
for (i=0; idim; i++) fprintf (out," %g",currn->x[i]);
fprintf (out,"\n");
currn->index = ++nnode;
currn = currn->next;
if (++nlines%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
}
// then, all the radii
currr = this->radii->first;
while (currr) {
fprintf (out,"vr %g\n",currr->r);
currr->index = ++nrad;
currr = currr->next;
if (++nlines%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
}
*/
// then, all the tangents
currt = this->tangents->first;
while (currt) {
fprintf (out,"vn");
for (i=0; idim; i++) fprintf (out," %g",currt->x[i]);
fprintf (out,"\n");
currt->index = ++ntan;
currt = currt->next;
if (++nlines%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
}
// then all of the segments
curr = this->first;
while (curr) {
// print the nodes, if they've not yet been printed
for (i=0; i<2; i++) {
if (curr->n[i]->flag) {
fprintf (out,"v");
for (j=0; jdim; j++) fprintf (out," %g",curr->n[i]->x[j]);
fprintf (out,"\n");
curr->n[i]->index = ++nnode;
curr->n[i]->flag = FALSE;
}
}
// print the radius, if it hasn't been printed yet
for (i=0; i<2; i++) {
if (curr->r[i]) {
if (curr->r[i]->flag) {
fprintf (out,"vr %g\n",curr->r[i]->r);
curr->r[i]->index = ++nrad;
curr->r[i]->flag = FALSE;
}
}
}
// then, print the segment
fprintf (out,"s");
for (i=0; i<2; i++) {
if (curr->r[i]) {
if (curr->t[i]) {
fprintf (out," %ld/%ld/%ld",curr->n[i]->index,curr->r[i]->index,
curr->t[i]->index);
} else {
fprintf (out," %ld/%ld",curr->n[i]->index,curr->r[i]->index);
}
} else {
if (curr->t[i]) {
fprintf (out," %ld//%ld",curr->n[i]->index,curr->t[i]->index);
} else {
fprintf (out," %ld",curr->n[i]->index);
}
}
}
fprintf (out,"\n");
curr = curr->next;
if (++nlines%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
}
fprintf(stderr,"%d elements\n",nlines);
fflush(stderr);
// if things bomb, send back a nonzero
return (0);
}
/*
* Write a .vtk file (steal syntax from fbbconvert.c)
*/
int write_vtk (FILE *out, seg_group_ptr this) {
int i,j;
int thisdim = 3;
int nnode = 0;
int nrad = 0;
int ntan = 0;
int nlines = 0;
seg_ptr curr;
fprintf (stderr,"Writing .vtk file");
fflush (stderr);
if (this->dim < thisdim) thisdim = this->dim;
// for this first take, blindly dump all segments and nodes
fprintf (out,"# vtk DataFile Version 2.0\n");
fprintf (out,"%s\n",VERSION);
fprintf (out,"ASCII\n");
fprintf (out,"DATASET UNSTRUCTURED_GRID\n");
// note: we assume that the node count is correct
fprintf (out,"POINTS %d float\n",this->nodes->num);
fprintf(stderr,"."); fflush(stderr);
// first, print the nodes
nnode = write_vtk_nodes_and_set_indexes (out, this->nodes, thisdim, 0);
fprintf(stderr,"."); fflush(stderr);
/*
// then, all the radii
currr = this->radii->first;
while (currr) {
fprintf (out,"vr %g\n",currr->r);
currr->index = ++nrad;
currr = currr->next;
if (++nlines%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
}
// then, all the tangents
currt = this->tangents->first;
while (currt) {
fprintf (out,"vn");
for (i=0; idim; i++) fprintf (out," %g",currt->t[i]);
fprintf (out,"\n");
currt->index = ++ntan;
currt = currt->next;
if (++nlines%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
}
*/
// then all of the segments
fprintf (out,"CELLS %d %d\n",this->num,this->num*3);
curr = this->first;
while (curr) {
// then, print the segment
fprintf (out,"2 %ld %ld\n",curr->n[0]->index,curr->n[1]->index);
curr = curr->next;
}
fprintf(stderr,"."); fflush(stderr);
fprintf (out,"CELL_TYPES %d\n",this->num);
for (i=0; inum; i++) fprintf (out,"3\n");
fprintf(stderr,"."); fflush(stderr);
// then, scalar values, namely, the radii
fprintf (out,"POINT_DATA %d\n",this->nodes->num);
fprintf (out,"SCALARS radius float 1\n");
fprintf (out,"LOOKUP_TABLE default\n");
(void) write_vtk_nodes_radii (out, this->nodes, this->radius);
fprintf(stderr,".\n"); fflush(stderr);
//fprintf(stderr,"%d elements\n",nlines);
//fflush(stderr);
// if things bomb, send back a nonzero
return (0);
}
/*
* Recursively index and print nodes to .vtk file
*
* Note: indexes start at 0 !
*/
int write_vtk_nodes_and_set_indexes (FILE *outfp, node_group_ptr this,
int dim, int sum) {
int i;
node_ptr curr;
// first, see if we have this entry already
if (this->child[0]) {
// there are children, call each of them
sum = write_vtk_nodes_and_set_indexes (outfp, this->child[0], dim, sum);
sum = write_vtk_nodes_and_set_indexes (outfp, this->child[1], dim, sum);
} else {
// there are no children, check all local entries
curr = this->first;
while (curr) {
fprintf (outfp,"%g",curr->x[0]);
for (i=1; ix[i]);
fprintf (outfp,"\n");
curr->index = sum++;
curr = curr->next;
}
}
return (sum);
}
/*
* Recursively index and print nodes to .vtk file
*
* Note: indexes start at 0 !
*/
int write_vtk_nodes_radii (FILE *outfp, node_group_ptr this, double rad) {
int i;
double bigrad;
node_ptr curr;
// first, see if we have this entry already
if (this->child[0]) {
// there are children, call each of them
(void) write_vtk_nodes_radii (outfp, this->child[0], rad);
(void) write_vtk_nodes_radii (outfp, this->child[1], rad);
} else {
// there are no children, check all local entries
curr = this->first;
while (curr) {
//fprintf (stderr, "node at %g %g %g\n",curr->x[0],curr->x[1],curr->x[2]);
//for (i=0; inumconn0; i++)
//fprintf (stderr, " 0 conn %d has nodes %ld %ld\n",i,
// curr->conn0[i]->n[0]->index,curr->conn0[i]->n[1]->index);
//for (i=0; inumconn1; i++)
//fprintf (stderr, " 1 conn %d has nodes %ld %ld\n",i,
// curr->conn1[i]->n[0]->index,curr->conn1[i]->n[1]->index);
// choose the largest radius at this node
bigrad = -1.;
for (i=0; inumconn0; i++)
if (curr->conn0[i]->r[0])
if (curr->conn0[i]->r[0]->r > bigrad)
bigrad = curr->conn0[i]->r[0]->r;
//fprintf (stderr, " rad %g\n",bigrad);
for (i=0; inumconn1; i++)
if (curr->conn1[i]->r[1])
if (curr->conn1[i]->r[1]->r > bigrad)
bigrad = curr->conn1[i]->r[1]->r;
//fprintf (stderr, " rad %g\n",bigrad);
// print either the largest radius, or the default radius
if (bigrad > -EPSILON) fprintf (outfp,"%g\n",bigrad);
else fprintf (outfp,"%g\n",rad);
curr = curr->next;
}
}
return (0);
}
/*
* Recursively set all flags
*/
int set_seg_flags (seg_group_ptr this, int val) {
seg_ptr curr;
// there are no children, check all local entries
curr = this->first;
while (curr) {
curr->flag = val;
curr = curr->next;
}
return (0);
}
/*
* Recursively set all flags
*/
int set_node_flags (node_group_ptr this, int val) {
node_ptr curr;
// first, see if we have this entry already
if (this->child[0]) {
// there are children, call each of them
(void) set_node_flags (this->child[0], val);
(void) set_node_flags (this->child[1], val);
} else {
// there are no children, check all local entries
curr = this->first;
while (curr) {
curr->flag = val;
curr = curr->next;
}
}
return (0);
}
/*
* Recursively set all flags
*/
int set_radius_flags (rad_group_ptr this, int val) {
rad_ptr curr;
// first, see if we have this entry already
if (this->child[0]) {
// there are children, call each of them
(void) set_radius_flags (this->child[0], val);
(void) set_radius_flags (this->child[1], val);
} else {
// there are no children, check all local entries
curr = this->first;
while (curr) {
curr->flag = val;
curr = curr->next;
}
}
return (0);
}
/*
* Write a Radiance file (fat)
*/
int write_rad (FILE *out, seg_group_ptr this) {
int write_no_cones = FALSE;
int write_many_colors = FALSE;
int i;
int count = 0;
int nlines = 0;
node_ptr currn;
seg_ptr curr;
int thisdim = this->dim;
double thisrad = -1.;
char color[5];
fprintf (stderr,"Writing .rad file");
fflush (stderr);
strcpy(color,"def");
// make sure we dump 3-dimensional data
thisdim = this->dim;
if (thisdim > 3) thisdim = 3;
// for this first take, blindly dump all segments and nodes
fprintf (out,"# %s\n",VERSION);
// first, the nodes
(void) write_radiance_nodes (out, this->nodes, thisdim, this->radius,
&count, &nlines);
// then, all of the segments
count = 0;
curr = this->first;
while (curr) {
if (write_many_colors) {
//sprintf(color,"def%1d",(int)(10.*rand()/(double)(RAND_MAX)));
sprintf(color,"def%1d",(int)(-0.5*log(rand()/(double)(RAND_MAX))));
}
// do we write a cylinder or a cone?
if (curr->r[0] && curr->r[1]) {
// if we have radii for both nodes...
if (write_no_cones) {
// find smallest radius
if (curr->r[0]->r > curr->r[1]->r) thisrad = curr->r[0]->r;
else thisrad = curr->r[1]->r;
fprintf (out,"%s cylinder c%ld 0 0 7",color,++count);
for (i=0; in[0]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
for (i=0; in[1]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e\n",thisrad);
} else if (fabs(curr->r[0]->r - curr->r[1]->r) < EPSILON) {
// and they're the same radius, write a cylinder
fprintf (out,"%s cylinder c%ld 0 0 7",color,++count);
for (i=0; in[0]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
for (i=0; in[1]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e\n",curr->r[0]->r);
} else {
// different radius, write a cone
fprintf (out,"%s cone c%ld 0 0 8",color,++count);
for (i=0; in[0]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
for (i=0; in[1]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e %12.8e\n",curr->r[0]->r,curr->r[1]->r);
}
} else if (curr->r[0]) {
// if we have radii for only one node...
if (write_no_cones) {
fprintf (out,"%s cylinder c%ld 0 0 7",color,++count);
for (i=0; in[0]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
for (i=0; in[1]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e\n",curr->r[0]->r);
} else if (fabs(curr->r[0]->r - this->radius) < EPSILON) {
// same radius, write a cylinder
fprintf (out,"%s cylinder c%ld 0 0 7",color,++count);
for (i=0; in[0]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
for (i=0; in[1]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e\n",this->radius);
} else {
// different radius, write a cone
fprintf (out,"%s cone c%ld 0 0 8",color,++count);
for (i=0; in[0]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
for (i=0; in[1]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e %12.8e\n",curr->r[0]->r,this->radius);
}
} else if (curr->r[1]) {
// if we have radii for only one node...
if (write_no_cones) {
fprintf (out,"%s cylinder c%ld 0 0 7",color,++count);
for (i=0; in[0]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
for (i=0; in[1]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e\n",curr->r[1]->r);
} else if (fabs(curr->r[1]->r - this->radius) < EPSILON) {
// same radius, write a cylinder
fprintf (out,"%s cylinder c%ld 0 0 7",color,++count);
for (i=0; in[0]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
for (i=0; in[1]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e\n",this->radius);
} else {
// different radius, write a cone
fprintf (out,"%s cone c%ld 0 0 8",color,++count);
for (i=0; in[0]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
for (i=0; in[1]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e %12.8e\n",this->radius,curr->r[0]->r);
}
} else {
// if we have no radii for either node...
// use the global radius and write a cylinder
fprintf (out,"%s cylinder c%ld 0 0 7",color,++count);
for (i=0; in[0]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
for (i=0; in[1]->x[i]);
for (i=thisdim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e\n",this->radius);
}
curr = curr->next;
if (++nlines%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
}
fprintf(stderr,"%d elements\n",nlines);
fflush(stderr);
// if things bomb, send back a nonzero
return (0);
}
/*
* Recursively write the Radiance spheres
*/
int write_radiance_nodes (FILE *out, node_group_ptr this, int dim,
double rad, int *count, int *nlines) {
int write_many_colors = FALSE;
int i;
double bigrad;
node_ptr curr;
char color[5];
strcpy(color,"def");
if (this->child[0]) {
(void) write_radiance_nodes (out, this->child[0], dim, rad, count, nlines);
(void) write_radiance_nodes (out, this->child[1], dim, rad, count, nlines);
} else {
curr = this->first;
while (curr) {
// choose the largest radius at this node
bigrad = -1.;
for (i=0; inumconn0; i++)
if (curr->conn0[i]->r[0])
if (curr->conn0[i]->r[0]->r > bigrad)
bigrad = curr->conn0[i]->r[0]->r;
for (i=0; inumconn1; i++)
if (curr->conn1[i]->r[1])
if (curr->conn1[i]->r[1]->r > bigrad)
bigrad = curr->conn1[i]->r[1]->r;
// print either the largest radius, or the default radius
if (bigrad < -EPSILON) bigrad = rad;
// but only print a sphere if it's large enough!
if (bigrad > EPSILON) {
if (write_many_colors) {
sprintf(color,"def%1d",(int)(-0.5*log(rand()/(double)(RAND_MAX))));
}
fprintf (out,"%s sphere s%ld 0 0 4",color,++(*count));
for (i=0; ix[i]);
for (i=dim; i<3; i++) fprintf (out," 0.");
fprintf (out," %12.8e\n",bigrad);
}
curr = curr->next;
if (++(*nlines)%DOTPER == 1) {
fprintf(stderr,".");
fflush(stderr);
}
}
}
return (0);
}
/*
* Calculate and write the appropriate dots for a given structure
*/
int write_dots (FILE *out, seg_group_ptr this, double res) {
int i;
int write_vtk_headers = FALSE;
double *thresh;
node_group_ptr dots;
// make a structure for the dots to print
dots = (NODE_GROUP*) malloc (sizeof(NODE_GROUP));
dots->num = 0;
dots->first = NULL;
dots->parent = NULL;
dots->child[0] = NULL;
dots->child[1] = NULL;
dots->axis = -1;
dots->min = NULL;
dots->max = NULL;
// set the dot proximity threshold
thresh = (double*) malloc (this->dim*sizeof(double));
for (i=0; idim; i++) thresh[i] = res;
// begin filling the dots structure
(void) write_dots_1 (this->nodes, dots, this->dim, thresh);
fprintf(stderr," made %d dots from nodes\n",dots->num);
// now, loop over all segments, adding any dots that will fit
(void) write_dots_2 (this, dots, this->dim, thresh);
fprintf(stderr," now %d dots after including segments\n",dots->num);
// finally, write the data
if (write_vtk_headers) {
fprintf (out,"# vtk DataFile Version 2.0\n");
fprintf (out,"%s\n",VERSION);
fprintf (out,"ASCII\n");
fprintf (out,"DATASET UNSTRUCTURED_GRID\n");
// note: we assume that the node count is correct
fprintf (out,"POINTS %d float\n",dots->num);
// finally, loop through the list of dots, writing what we've found
}
// write all of the dot locations
(void) write_dots_3 (out, dots, this->dim);
// then the rest of the data
if (write_vtk_headers) {
fprintf (out,"CELLS %d %d\n",dots->num,dots->num*2);
for (i=0; inum; i++) fprintf (out,"1 %d\n",i);
fprintf (out,"CELL_TYPES %d\n",dots->num);
for (i=0; inum; i++) fprintf (out,"1\n");
}
// free the tree structure and all of the dots
free_nodes (dots);
return (0);
}
/*
* Recursively include nodes as dots
*/
int write_dots_1 (node_group_ptr this, node_group_ptr dots,
int dim, double *thresh) {
int i;
double dist;
node_ptr curr, close_dot;
if (this->child[0] != NULL) {
(void) write_dots_1 (this->child[0], dots, dim, thresh);
(void) write_dots_1 (this->child[1], dots, dim, thresh);
} else {
curr = this->first;
while (curr != NULL) {
// set the threshold distance
dist = thresh[0];
// check to see if this node is too close to any other node
close_dot = find_closest_node (dots, dim, curr->x, NULL, &dist);
//close_dot = NULL;
//close_dot = curr;
// if it is not too close, add it to the dots list
if (close_dot == NULL) {
(void) add_node (dots, dim, curr->x, 0);
//fprintf(stderr," added dot, closest node was %g\n",dist);
//} else {
//fprintf(stderr," no dot, closest node was %g\n",dist);
}
curr = curr->next;
}
}
return (0);
}
/*
* Recursively include segments as dots
*/
int write_dots_2 (seg_group_ptr this, node_group_ptr dots,
int dim, double *thresh) {
int i,isubseg,nsubseg;
double dist,len;
double *testloc, *tangent;
seg_ptr curr;
node_ptr close_dot;
//if (this->child[0] != NULL) {
//(void) write_dots_2 (this->child[0], dots, dim, thresh);
//(void) write_dots_2 (this->child[1], dots, dim, thresh);
//} else {
testloc = (double*) malloc (dim*sizeof(double));
tangent = (double*) malloc (dim*sizeof(double));
curr = this->first;
while (curr != NULL) {
// how long is this segment?
len = seg_length (curr, dim);
nsubseg = (int)(len/thresh[0]);
// find the straight tangential vector along this edge
for (i=0; in[1]->x[i]-curr->n[0]->x[i];
// march along the segment, checking locations for particles
for (isubseg=1; isubsegn[0]->x[i] +
(double)isubseg * tangent[i] / (double)nsubseg;
// check to see if this dot is too close to any other dot
dist = thresh[0];
close_dot = find_closest_node (dots, dim, testloc, NULL, &dist);
// if it is not too close, add it to the dots list
if (close_dot == NULL) (void) add_node (dots, dim, testloc, 0);
}
// if the segment is between 1 and 2 thresholds long, try two points
//if (nsubseg == 1) {
if (FALSE) {
// set the test location
for (i=0; in[0]->x[i] +
(thresh[0]*1.01) * tangent[i] / len;
// check to see if this dot is too close to any other dot
dist = thresh[0];
close_dot = find_closest_node (dots, dim, testloc, NULL, &dist);
// if it is not too close, add it to the dots list
if (close_dot == NULL) (void) add_node (dots, dim, testloc, 0);
// now, try the dot thresh distance from the other node
for (i=0; in[1]->x[i] -
(thresh[0]*1.01) * tangent[i] / len;
dist = thresh[0];
close_dot = find_closest_node (dots, dim, testloc, NULL, &dist);
if (close_dot == NULL) (void) add_node (dots, dim, testloc, 0);
}
curr = curr->next;
}
//}
return (0);
}
/*
* Recursively write dots as ASCII file
*/
int write_dots_3 (FILE *out, node_group_ptr this, int dim) {
int i;
node_ptr curr;
if (this->child[0] != NULL) {
(void) write_dots_3 (out, this->child[0], dim);
(void) write_dots_3 (out, this->child[1], dim);
} else {
curr = this->first;
while (curr != NULL) {
// dump the data
fprintf(out,"%g",curr->x[0]);
for (i=1; ix[i]);
fprintf(out,"\n");
curr = curr->next;
}
}
return (0);
}
/*
* Create a Gaussian perturbation in as many dimensions as required,
* mean is 0, std dev is scale
*/
void perturb_gaussian (double *dx, double scale, int dim) {
int d;
double t1,t2;
// do two at a time
for (d=0; dn[0]->x[i]-this->n[1]->x[i],2);
return (sqrt(vds));
}
/* distance squared */
double seg_length_sqrd (seg_ptr this, int dim) {
int i;
double vds = 0.;
for (i=0; in[0]->x[i]-this->n[1]->x[i],2);
return (vds);
}
/* distance squared */
double vec_dist_sqrd (double *x, double *y, int dim) {
int i;
double vds = 0.;
for (i=0; ichild[0] != NULL) {
(void) free_nodes (this->child[0]);
(void) free_nodes (this->child[1]);
} else {
curr = this->first;
while (curr != NULL) {
free (curr->x);
free (curr->conn0);
free (curr->conn1);
next = curr->next;
free (curr);
curr = next;
}
}
free (this->min);
free (this->max);
free (this);
return (0);
}
/*
* Recursively free memory associated with tangents
*/
int free_tangents (tan_group_ptr this) {
tan_ptr curr, next;
if (this->child[0] != NULL) {
(void) free_tangents (this->child[0]);
(void) free_tangents (this->child[1]);
} else {
curr = this->first;
while (curr != NULL) {
free (curr->x);
next = curr->next;
free (curr);
curr = next;
}
}
free (this->min);
free (this->max);
free (this);
return (0);
}
/*
* Recursively free memory associated with radii
*/
int free_radii (rad_group_ptr this) {
rad_ptr curr, next;
if (this->child[0] != NULL) {
(void) free_radii (this->child[0]);
(void) free_radii (this->child[1]);
} else {
curr = this->first;
while (curr != NULL) {
next = curr->next;
free (curr);
curr = next;
}
}
free (this);
return (0);
}
/*
* Free memory associated with segments (everything, basically)
*/
int free_segs (seg_group_ptr this) {
seg_ptr curr, next;
// first, remove radii, tangents, and nodes
(void) free_radii (this->radii);
(void) free_tangents (this->tangents);
(void) free_nodes (this->nodes);
// then, delete all segments
curr = this->first;
while (curr != NULL) {
next = curr->next;
free (curr);
curr = next;
}
// finally, free the structure itself
free (this);
return (0);
}
/*
* This function writes basic usage information to stderr,
* and then quits. Too bad.
*/
int Usage(char progname[MAXSTR],int status) {
static char **cpp, *help_message[] =
{
"where [options] are one or more of the following:",
" ",
" -seg write an obj-like file of the segments (default)",
" ",
" -rad write a Radiance-readable file of the cylinders",
" ",
" -vtk write a VTK- and mayavi-readable file",
" ",
" -png write a PNG image file",
" ",
" -coarsen [l] coarsen structure so that no element is shorter than",
" length l; default l=1",
" ",
" -refine [l] refine structure so that no element is longer than",
" length l; default l=1",
" ",
" -roughen [s] roughen structure by perturbing nodes, s is scale",
" factor; default s=1",
" ",
" -treeradius [s [dim [val]]]",
" set all radii to mimic woody plant growth,",
" relative strength is s * 10^6 * E / (rho * g * l),",
" next two entries define which nodes are roots:",
" 'x 0' means whichever node is closest to the x=0 plane;",
" default= 1.0, x, 0.0",
" ",
" -translate [x [y [z]]]",
" translate structure by vector x,y,z; default= 0,0,0",
" ",
" -scale [x [y [z]]]",
" scale structure by magnitudes x,y,z; default= 1,1,1",
" mirroring can be done using negative values,",
" if fewer than (dimension) values are given, the",
" last value will be used for the rest of the dimensions",
" ",
" -info dump x,y,z,r min/max and other information",
" ",
" -split split into two .rad files along the longest dimension",
" ",
" -version returns version information",
" ",
" -help returns this help information",
" ",
"The input file should be in .seg or .rad format",
" ",
"Operations are done in the order that they appear on the command-",
"line, so make sure to list your input file first!",
" ",
"Options may be abbreviated to an unambiguous length (duh).",
"Output is always to stdout, so make sure to redirect (> or |)",
NULL
};
fprintf(stderr,"Usage:\n\n %s [infile | options]\n\n", progname);
for (cpp = help_message; *cpp; cpp++)
fprintf(stderr, "%s\n", *cpp);
fflush(stderr);
exit(status);
return(0);
}