path: root/src/structure.c

/*
 * structure.c
 *
 * 3D analysis
 *
 * (c) 2007 Gordon Ball <gfb21@cam.ac.uk>
 *  dtr - Diffraction Tomography Reconstruction
 *
 */
 #include <string.h>
 #include <stdlib.h>
 #include <math.h>
 #include "reflections.h"
 #include "utils.h"
 #include "structure.h"
 
 static int reflect_count(ReflectionContext *rctx) {
 	Reflection *r;
 	r = rctx->reflections;
 	int count=0;
 	do {
 		count++;
 	} while ((r=r->next)!=NULL);	
 	return count;
 }
 
 
 /*
  * Find the largest single dimension variable in the world
  * Relevant for octtree
  */
 static double largest_dimension(ReflectionContext *rctx) {
 	Reflection *r;
 	double max=0;
 	double val;
 	r = rctx->reflections;
 	do {
 		if ((val = r->x) > max) max = val;
 		if ((val = r->y) > max) max = val;
 		if ((val = r->z) > max) max = val;
 	} while ((r = r->next) != NULL);
 	return max;
 }
 
 /*
  * Calculate the volume of a shell from r1->r2
  * Requires r2 > r1
  */
 static double get_shell_volume(double r2, double r1) {
 	return (4. * M_PI * (r2*r2*r2 - r1*r1*r1))/3.; 
 }
 
 /*
  * Attempts to calculate the optimum mask radius
  * We have a problem here with non-isotropicness since the points in question usually are close to planar, and not equally distributed about the origin
  * It may be more necessary to try and construct a bounding polygon - although that's going to get horribly messy
  * Will attempt an octtree based method, but for the moment this will return the radius for which 2/3 of the points are contained.
  */
 
 double get_mask_radius(ReflectionContext *rctx) {
 	Reflection *r;
 	double maxmod = 0.;
 	double modul;
 	double density;
 	double size;
 	double mask=0.;
 	int num_intervals=10;
 	int interval;
 	int count=0;
 	r = rctx->reflections;
 	do {
		if (r->type == REFLECTION_NORMAL) { 		
 			if ((modul = modulus(r->x,r->y,r->z)) > maxmod) maxmod = modul;
 			count++;
		}
 	} while ((r = r->next) != NULL);
 	
 	printf("g_m_r: count=%d, maxmod=%f\n",count,maxmod);
 	
 	int *bucket;
 	bucket = calloc(num_intervals,sizeof(int));
 	maxmod *= 1.001;
 	size = maxmod/num_intervals;
 	r = rctx->reflections;
 	do {
 		if (r->type == REFLECTION_NORMAL) {
 			modul = modulus(r->x,r->y,r->z);
 			interval = (int)((modul/maxmod)* (double)num_intervals);
 			bucket[interval] += 1;
 		}
 	} while ((r = r->next) != NULL);
 	
 	int i;
 	int sum=0;
 	for (i=0;i<num_intervals;i++) {
 		density = ( bucket[i] * 1e30 )/get_shell_volume(size*(i+1),size*i);
 		printf("interval %d - count=%d density=%f volume=%f r2=%f r1=%f\n",i,bucket[i],density,get_shell_volume(size*(i+1),size*i), size*(i+1), size*i);
 		if (( sum += bucket[i]) < (0.66*count)) mask = size*(i+1);
 	}
 	printf("g_m_r: returning %f\n",mask);
 	return mask;
 }
 
 /*
  * Returns a new ReflectionContext with all normal reflections outside mask_radius dumped
  * The new ReflectionContext will occupy a non-contigious block of memory
  */
 ReflectionContext *apply_mask_radius(double mask_radius, ReflectionContext *rctx) {
 	Reflection *r, *newr;
 	ReflectionContext *nrc;
 	double modul;
 	//printf("a_m_r: points starting %d\n",reflect_count(nrc));
 	
 	nrc = reflection_init();
 	r = rctx->reflections;
 	newr = nrc->reflections;
 	do {
 		if (r->type == REFLECTION_NORMAL) {
 			modul = modulus(r->x,r->y,r->z);
 			printf("a_m_r: modul/mask_radius %f\n",modul/mask_radius);
 			if (modul < mask_radius) {
 				newr = malloc(sizeof(Reflection));
 				memcpy(newr,r,sizeof(Reflection));
 				reflection_add_from_reflection(nrc,newr);
 			}
 		}
 	} while ((r=r->next) != NULL);
 	//printf("a_m_r: points remaining %d\n",reflect_count(nrc));
 	return nrc;
 }
 
 /*
  * checks if the supplied reflection falls in octtree volume vol
  * returns -1 if it falls outside it
  * returns 0-7 if it falls within depending which child volume it would occupy
  */
 static int in_octtree_volume(Reflection *r, OctTree *vol) {
 	int val=0;
 	double minx,miny,minz;
 	double maxx,maxy,maxz;
 	minx = vol->ox-vol->halfedge;
 	miny = vol->oy-vol->halfedge;
 	minz = vol->oz-vol->halfedge;
 	maxx = vol->ox+vol->halfedge;
 	maxy = vol->oy+vol->halfedge;
 	maxz = vol->oz+vol->halfedge;
 	
 	if ( r->x > maxx || r->x <= minx || r->y > maxy || r->y <= miny || r->z > maxz || r->z <= minz ) {
 		return -1;
 	} else {
 		if (r->x <= vol->ox) val += 1;
 		if (r->y <= vol->oy) val += 2;
 		if (r->z <= vol->oz) val += 4;
 	}
 	return val;
 }
 /*
  * set the x,y,z and halfedge params for a octtree child based on the parent and child#
  */
 static void set_octtree_origin(OctTree *parent, OctTree *child, int childnum) {
 	int val = childnum;
 	child->halfedge = parent->halfedge * 0.5;
 	if (val >= 4) {
 		child->oz = parent->oz - child->halfedge;
 		val %= 4;	
 	} else {
 		child->oz = parent->oz + child->halfedge;
 	}
 	if (val >= 2) {
 		child->oy = parent->oy - child->halfedge;
 		val %= 2;	
 	} else {
 		child->oy = parent->oy + child->halfedge;
 	}
 	if (val >= 1) {
 		child->ox = parent->ox - child->halfedge;	
 	} else {
 		child->ox = parent->ox + child->halfedge;
 	}
 }
 
 /*
  * checks to see if there are any reflections in this volume
  * if there are, create a new octtree and attach it to the appropriate branch of the parent
  * else attach null
  * if we haven't reached maximum depth, spawn 8 new requests until the desired resolution is reached
  */
 static void stack_octtree(Reflection **rl, int rcount, OctTree *parent, int childnum, int maxdepth) {
 	//if there are reflections here
 	//printf("s_o: starting depth=%d child=%d\n",parent->depth,childnum);
 	if (rcount > 0) {
 		//printf("s_o: reflections=%d\n",rcount);
 		OctTree *here = malloc(sizeof(OctTree)); //create a new OctTree node
 		here->child = calloc(8,sizeof(OctTree *));
 		here->parent = parent;
 		here->childnum = childnum;
 		here->list = NULL;
 		parent->child[childnum] = here; //attach it to the parent
 		set_octtree_origin(parent,here,childnum);
 		//printf("s_o: set origin (%f,%f,%f) halfedge %f\n",here->ox,here->oy,here->oz,here->halfedge);
 		
 		if ((here->depth = parent->depth+1) < maxdepth) { //only process children if we're not at maxdepth
 			//printf("s_o: allocating rcount=%d\n",rcount);
 			//int *dest = calloc(rcount,sizeof(int)); //list of the reflections and which child to route them to
 			//int *distrib = calloc(8,sizeof(int)); //count of reflections to route to each child
 			int dest[rcount];
 			int distrib[8] = {0,0,0,0,0,0,0,0};
 			Reflection **list;
	 		
 			int i,j;
 			for (i=0;i<rcount;i++) {
 				dest[i] = in_octtree_volume(rl[i],here);
 				//printf("s_o: reflection %d (%f,%f,%f) in volume %d\n",i,rl[i]->x,rl[i]->y,rl[i]->z,dest[i]);
		 		distrib[dest[i]] += 1;	
 			}
 			int n;
 			for (i=0;i<8;i++) {
 				if (distrib[i] > 0) {
	 				//printf("s_o: creating list for child %d, %d members\n",i,distrib[i]);
 					list = malloc(sizeof(Reflection *)*distrib[i]);
 					n=0;
 					for (j=0;j<rcount;j++) {
 						if (dest[j] == i) {
 							list[n++] = rl[j];
	 					}
 					}
 					stack_octtree(list,distrib[i],here,i,maxdepth);
	 				free(list);
 				} else {
 					//printf("s_o: no reflections for child %d\n",i);
 					here->child[i]=NULL;
 				}
 			}
 			//printf("s_o: [%d] ready to free distrib=%d dest=%d\n",here->depth,distrib,dest);
 			//free(distrib);
 			//printf("s_o: freed distrib\n");
 			//free(dest);
 			//printf("s_o: freed dest\n");
 		} else { //add a reflection list of the children
 			ReflectionList *l = malloc(sizeof(ReflectionList));
 			l->r = malloc(rcount*sizeof(Reflection *));
 			memcpy(l->r,rl,rcount*sizeof(Reflection *));
 			here->list = l;
 		}
 	} else { //if there are no reflections, just attach NULL and return
 		printf("s_o: no reflections, attaching null\n");
 		parent->child[childnum] = NULL;
 	}
 }
 
 
 /*
  * generate an octtree filling all space out to the largest dimension in the reflectionlist
  * then eliminate all volumes containing no reflections down to the desired accuracy
  * TODO: use the same basis as the reflection
  */
 
 OctTree *gen_octtree(ReflectionContext *rctx, int depth) {
 	printf("g_o: starting\n");
 	double max = largest_dimension(rctx)*1.01;
 	int count = reflect_count(rctx);
 	Reflection *r;
 	
 	OctTree *top;
 	top = malloc(sizeof(OctTree));
 	top->child = calloc(8,sizeof(OctTree *));
 	top->parent=NULL;
 	top->halfedge = max;
 	top->ox = 0;
 	top->oy = 0;
 	top->oz = 0;
 	top->depth=0;
 	top->childnum=-1;
 	
 	Reflection *rl[count];
 	r = rctx->reflections;
 	int n=0;
 	do {
 		rl[n++] = r;
 	} while ((r=r->next)!=NULL);
 	
 	int i;
 	for (i=0;i<8;i++) {
 		//printf("g_o: stack %d\n",i);
 		stack_octtree(rl,count,top,i,depth);
 	}
 	return top;
 }
 
 static void print_octtree_stack(OctTree *here, int* at_level) {
 		int i;
 		for (i=0;i<8;i++) {
 			if (here->child[i] != NULL) print_octtree_stack(here->child[i],at_level);	
 		}
 		at_level[here->depth]++;
 }
 
 void print_octtree(OctTree *tree) {
 		int* at_level = calloc(16,sizeof(int));
 		print_octtree_stack(tree,at_level);
 		int i;
 		for (i=0;i<16;i++) {
 			printf("level %d nodes %d\n",i,at_level[i]);
 		}
 }
 
 //return a list of pointers to the 26 surrounding nodes
 OctTreeLinkedList *get_adjacent_nodes(OctTree *o, int *count) {
 		
 }
 
 //return a list of all level x nodes
 void stack_get_depth(OctTree *o, int *maxdepth) {
 	int i;
 	if (o->depth > *maxdepth) *maxdepth = o->depth;
 	for (i=0;i<8;i++) {
 		if (o->child[i] != NULL) stack_get_depth(o->child[i],maxdepth);
 	}
 }
 
 OctTreeLinkedList *get_bottom_nodes(OctTree *o, int *count, int level) {
 		int depth=0;
 		stack_get_depth(o,&depth);
 		
 }
 
 void free_linked_list(OctTreeLinkedList *otll) {
 	OctTreeLinkedList *o,*next;
 	o = otll;
 	do {
 		next = o;
 		free(o);
 		o = next;
 	} while (o != NULL);	
 }
 
 void dump_histogram(ReflectionContext *rctx) {
 	Reflection *r;
 	int count = reflect_count(rctx);
 	int n=0;
 	double dist;
 	Reflection **rl = malloc(count*sizeof(Reflection *));
 	r = rctx->reflections;
 	do {
 		rl[n++] = r;
 	} while ((r = r->next) != NULL);
 	
 	FILE *f;
 	
 	f = fopen("histogram","w");
 	int i,j;
 	for (i=0;i<count;i++) {
 		for (j=i+1;j<count;j++) {
 			dist = modulus(rl[i]->x-rl[j]->x,rl[i]->y-rl[j]->y,rl[i]->z-rl[j]->z);
 			fprintf(f,"%f\n",dist); 
 		}	
 	}
 	fclose(f);	
 	
 }
 
 /*
  * look for sections of the tree with gaps of at least req_length between branches
  * add a node that branches after at least req_length or terminates
  */  
 OctTreeLinkedList *stack_find_sparse_trees(OctTree *o, OctTreeLinkedList *l, int *count, int req_length, int *cur_length, int allow_end) {
 	OctTreeLinkedList *nl;
 	int i;
 	int children=0;
 	for (i=0;i<8;i++) {
 		if (o->child[i] != NULL) children++;
 	}
 	if (children==0) { //end of a chain, add if sufficiently long
 		(*cur_length)++;
 		if (allow_end==1) {
 			if (*cur_length >= req_length) {
 				nl = malloc(sizeof(OctTreeLinkedList));
 				nl->o = o;
 				nl->next = NULL;
 				l->next = nl;
 				(*count)++;
 				printf("s_f_s_t: found end-of-chain length=%d depth=%d\n",*cur_length,o->depth);
 			} else {
 				nl = l;
 			}
 		} else {
 			nl = l;
 		}
 		(*cur_length)=0;
 		return nl; //return the current list pointer
 	} else {
 		if (children==1) { //middle of a singular chain, add to length counter
 			(*cur_length)++;
 			nl = l;
 		} else {
 			if (*cur_length >= req_length) { //branch point, see if the current chain is long enough to add
	 			nl = malloc(sizeof(OctTreeLinkedList));
 				nl->o = o;
 				nl->next = NULL;
 				l->next = nl;	
 				(*count)++;
 				printf("s_f_s_t: found branch point length=%d depth=%d\n",*cur_length,o->depth);
 			} else {
 				nl = l;
 			}
 			(*cur_length)=0; //regardless whether this section was long enough, zero the counter
 		}
 		
 		for (i=0;i<8;i++) {
 			if (o->child[i] != NULL) {
 				nl = stack_find_sparse_trees(o->child[i],nl,count,req_length,cur_length,allow_end);	
 			}	
 		}
 		return nl;	
 	}
 	
 	
 }
 
 OctTreeLinkedList *find_sparse_trees(OctTree *o, int req_length, int allow_end, int *count) {
 	printf("f_s_t: starting\n");
 	int cur_length=0;
 	OctTreeLinkedList *ol = malloc(sizeof(OctTreeLinkedList));
 	OctTreeLinkedList *ol2;
 	ol->o = NULL;
 	ol->next = NULL;
 	stack_find_sparse_trees(o,ol,count,req_length,&cur_length,allow_end);
 	ol2 = ol->next;
 	free(ol);
 	return ol2;
 }
 
 ReflectionContext *change_reflection_basis(ReflectionContext *rctx, Basis *basis) {
 	ReflectionContext *new = reflection_init();
 	Reflection *r;
 
 	//calculate a change-of-basis matrix, replace all the reflection coordinates thus
 	//hence we hopefully have a square-basis representation, which octtree statistics might work on. maybe.
 }