aboutsummaryrefslogtreecommitdiff
path: root/src/structure.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/structure.c')
-rw-r--r--src/structure.c444
1 files changed, 0 insertions, 444 deletions
diff --git a/src/structure.c b/src/structure.c
deleted file mode 100644
index 8ce0888..0000000
--- a/src/structure.c
+++ /dev/null
@@ -1,444 +0,0 @@
-/*
- * 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.
- }
-
-
-
-
-
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-03 04:41:24 +0000