BRL-CAD
voxels.c
Go to the documentation of this file.
1 /* V O X E L S . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2009-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 
21 #include "common.h"
22 
23 #include <stdlib.h>
24 #include <math.h>
25 #include <string.h>
26 #include <stdio.h>
27 
28 #include "vmath.h" /* vector math macros */
29 #include "raytrace.h" /* librt interface definitions */
30 
31 #include "analyze.h"
32 
33 
34 /**
35  * Function to get the corresponding region entry to a region name.
36  */
37 HIDDEN struct voxelRegion *
38 getRegionByName(struct voxelRegion *head, const char *regionName) {
39  struct voxelRegion *ret = NULL;
40 
41  BU_ASSERT(regionName != NULL);
42 
43  if (head->regionName == NULL) { /* the first region on this voxel */
44  head->regionName = bu_strdup(regionName);
45  ret = head;
46  }
47  else {
48  while (head->nextRegion != NULL) {
49  if (bu_strcmp(head->regionName, regionName) == 0) {
50  ret = head;
51  break;
52  }
53 
54  head = head->nextRegion;
55  }
56 
57  if (ret == NULL) { /* not found until here */
58  if (bu_strcmp(head->regionName ,regionName) == 0) /* is it the last one on the list? */
59  ret = head;
60  else {
61  BU_ALLOC(ret, struct voxelRegion);
62  head->nextRegion = ret;
63  ret->regionName = bu_strdup(regionName);
64  }
65  }
66  }
67 
68  return ret;
69 }
70 
71 
72 /**
73  * rt_shootray() was told to call this on a hit.
74  *
75  * This callback routine utilizes the application structure which
76  * describes the current state of the raytrace.
77  *
78  * This callback routine is provided a circular linked list of
79  * partitions, each one describing one in and out segment of one
80  * region for each region encountered.
81  *
82  * The 'segs' segment list is unused in this example.
83  */
84 int
85 hit_voxelize(struct application *ap, struct partition *PartHeadp, struct seg *UNUSED(segs))
86 {
87  struct partition *pp = PartHeadp->pt_forw;
88  struct rayInfo *voxelHits = (struct rayInfo*) ap->a_uptr;
89  fastf_t sizeVoxel = voxelHits->sizeVoxel;
90  fastf_t *fillDistances = voxelHits->fillDistances;
91 
92  while (pp != PartHeadp) {
93  /**
94  * hitInp, hitOutp are hit structures to save distances where
95  * ray entered and exited the present partition. hitDistIn,
96  * hitDistOut are the respective distances from the origin of
97  * ray. voxelNumIn, voxelNumOut are the voxel numbers where
98  * ray entered and exited the present partition.
99  */
100 
101  struct hit *hitInp = pp->pt_inhit;
102  struct hit *hitOutp = pp->pt_outhit;
103  fastf_t hitDistIn = hitInp->hit_dist - 1.;
104  fastf_t hitDistOut = hitOutp->hit_dist - 1.;
105  int voxelNumIn = (int)(hitDistIn / sizeVoxel);
106  int voxelNumOut = (int)(hitDistOut / sizeVoxel);
107 
108  if (EQUAL((hitDistOut / sizeVoxel), floor(hitDistOut / sizeVoxel)))
109  voxelNumOut = FMAX(voxelNumIn, voxelNumOut - 1);
110 
111  /**
112  * If voxel entered and voxel exited are same then nothing can
113  * be evaluated till we see the next partition too. If not,
114  * evaluate entry voxel. Also, all the intermediate voxels are
115  * in.
116  */
117  if (voxelNumIn == voxelNumOut) {
118  fillDistances[voxelNumIn] += hitDistOut - hitDistIn;
119  getRegionByName(voxelHits->regionList + voxelNumIn, pp->pt_regionp->reg_name)->regionDistance += hitDistOut - hitDistIn;
120  } else {
121  int j;
122 
123  fillDistances[voxelNumIn] += (voxelNumIn + 1) * sizeVoxel - hitDistIn;
124  getRegionByName(voxelHits->regionList + voxelNumIn, pp->pt_regionp->reg_name)->regionDistance += (voxelNumIn + 1) * sizeVoxel - hitDistIn;
125 
126  for (j = voxelNumIn + 1; j < voxelNumOut; ++j) {
127  fillDistances[j] += sizeVoxel;
128  getRegionByName(voxelHits->regionList + j, pp->pt_regionp->reg_name)->regionDistance += sizeVoxel;
129  }
130 
131  fillDistances[voxelNumOut] += hitDistOut - (voxelNumOut * sizeVoxel);
132  getRegionByName(voxelHits->regionList + voxelNumOut, pp->pt_regionp->reg_name)->regionDistance += hitDistOut - (voxelNumOut * sizeVoxel);
133  }
134 
135  pp = pp->pt_forw;
136  }
137 
138  return 0;
139 
140 }
141 
142 
143 /**
144  * voxelize function takes raytrace instance and user parameters as inputs
145  */
146 void
147 voxelize(struct rt_i *rtip, fastf_t sizeVoxel[3], int levelOfDetail, void (*create_boxes)(void *callBackData, int x, int y, int z, const char *regionName, fastf_t percentageFill), void *callBackData)
148 {
149  struct rayInfo voxelHits;
150  int numVoxel[3];
151  int yMin;
152  int zMin;
153  int i, j, k, rayNum;
154  fastf_t *voxelArray;
155  fastf_t rayTraceDistance;
156  fastf_t effectiveDistance;
157 
158  /* get bounding box values etc. */
159  rt_prep_parallel(rtip, 1);
160 
161  /* calculate number of voxels in each dimension */
162  numVoxel[0] = (int)(((rtip->mdl_max)[0] - (rtip->mdl_min)[0])/sizeVoxel[0]) + 1;
163  numVoxel[1] = (int)(((rtip->mdl_max)[1] - (rtip->mdl_min)[1])/sizeVoxel[1]) + 1;
164  numVoxel[2] = (int)(((rtip->mdl_max)[2] - (rtip->mdl_min)[2])/sizeVoxel[2]) + 1;
165 
166  if (EQUAL(numVoxel[0] - 1, (((rtip->mdl_max)[0] - (rtip->mdl_min)[0])/sizeVoxel[0])))
167  numVoxel[0] -=1;
168 
169  if (EQUAL(numVoxel[1] - 1, (((rtip->mdl_max)[1] - (rtip->mdl_min)[1])/sizeVoxel[1])))
170  numVoxel[1] -=1;
171 
172  if (EQUAL(numVoxel[2] - 1, (((rtip->mdl_max)[2] - (rtip->mdl_min)[2])/sizeVoxel[2])))
173  numVoxel[2] -=1;
174 
175  voxelHits.sizeVoxel = sizeVoxel[0];
176 
177  /* voxelArray stores the distance in path of ray inside a voxel which is filled
178  * initialize with 0s */
179  voxelArray = (fastf_t *)bu_calloc(numVoxel[0], sizeof(fastf_t), "voxelize:voxelArray");
180 
181  /* regionList holds the names of voxels inside the voxels
182  * initialize with NULLs */
183  voxelHits.regionList = (struct voxelRegion *)bu_calloc(numVoxel[0], sizeof(struct voxelRegion), "voxelize:regionList");
184 
185  /* minimum value of bounding box in Y and Z directions */
186  yMin = (int)((rtip->mdl_min)[1]);
187  zMin = (int)((rtip->mdl_min)[2]);
188 
189  BU_ASSERT_LONG(levelOfDetail, >, 0);
190  /* 1.0 / (levelOfDetail + 1) and effectiveDistance have to be used multiple times in the following loops */
191  rayTraceDistance = 1. / levelOfDetail;
192  effectiveDistance = levelOfDetail * levelOfDetail * sizeVoxel[0];
193 
194  /* start shooting */
195  for (i = 0; i < numVoxel[2]; ++i) {
196  for (j = 0; j < numVoxel[1]; ++j) {
197  struct application ap;
198 
199  RT_APPLICATION_INIT(&ap);
200  ap.a_rt_i = rtip;
201  ap.a_onehit = 0;
202  VSET(ap.a_ray.r_dir, 1., 0., 0.);
203 
204  ap.a_hit = hit_voxelize;
205  ap.a_miss = NULL;
206  ap.a_uptr = &voxelHits;
207 
208  voxelHits.fillDistances = voxelArray;
209 
210  for (rayNum = 0; rayNum < levelOfDetail; ++rayNum) {
211  for (k = 0; k < levelOfDetail; ++k) {
212 
213  /* ray is hit through evenly spaced points of the unit sized voxels */
214  VSET(ap.a_ray.r_pt, (rtip->mdl_min)[0] - 1.,
215  yMin + (j + (k + 0.5) * rayTraceDistance) * sizeVoxel[1],
216  zMin + (i + (rayNum + 0.5) * rayTraceDistance) * sizeVoxel[2]);
217  rt_shootray(&ap);
218  }
219  }
220 
221  /* output results via a call-back supplied by user*/
222  for (k = 0; k < numVoxel[0]; ++k) {
223  if (voxelHits.regionList[k].regionName == NULL)
224  /* an air voxel */
225  create_boxes(callBackData, k, j, i, NULL, 0.);
226  else {
227  struct voxelRegion *tmp = voxelHits.regionList + k;
228  struct voxelRegion *old = tmp->nextRegion;
229 
230  create_boxes(callBackData, k, j, i, tmp->regionName, tmp->regionDistance / effectiveDistance);
231 
232  if (tmp->regionName != 0)
233  bu_free(tmp->regionName, "voxelize:voxelRegion:regionName");
234 
235  while (old != NULL) {
236  tmp = old;
237  create_boxes(callBackData, k, j, i, tmp->regionName, tmp->regionDistance / effectiveDistance);
238  old = tmp->nextRegion;
239 
240  /* free the space allocated for new regions */
241  if (tmp->regionName != 0)
242  bu_free(tmp->regionName, "voxelize:voxelRegion:regionName");
243 
244  BU_FREE(tmp, struct voxelRegion);
245  }
246  }
247 
248  voxelArray[k] = 0.;
249  voxelHits.regionList[k].regionName = NULL;
250  voxelHits.regionList[k].nextRegion = NULL;
251  voxelHits.regionList[k].regionDistance = 0.;
252  }
253  }
254  }
255 
256  bu_free(voxelArray, "voxelize:voxelArray");
257  bu_free(voxelHits.regionList, "voxelize:regionList");
258 }
259 
260 
261 /*
262  * Local Variables:
263  * tab-width: 8
264  * mode: C
265  * indent-tabs-mode: t
266  * c-file-style: "stroustrup"
267  * End:
268  * ex: shiftwidth=4 tabstop=8
269  */
struct xray a_ray
Actual ray to be shot.
Definition: raytrace.h:1583
struct region * pt_regionp
ptr to containing region
Definition: raytrace.h:580
struct hit * pt_outhit
OUT hit ptr.
Definition: raytrace.h:579
point_t mdl_min
min corner of model bounding RPP
Definition: raytrace.h:1769
#define VSET(a, b, c, d)
Definition: color.c:53
Definition: raytrace.h:368
#define BU_ASSERT_LONG(_lhs, _relation, _rhs)
Definition: defines.h:240
Definition: raytrace.h:248
fastf_t * fillDistances
Definition: analyze.h:100
Header file for the BRL-CAD common definitions.
const char * reg_name
Identifying string.
Definition: raytrace.h:539
void voxelize(struct rt_i *rtip, fastf_t sizeVoxel[3], int levelOfDetail, void(*create_boxes)(void *callBackData, int x, int y, int z, const char *regionName, fastf_t percentageFill), void *callBackData)
Definition: voxels.c:147
#define BU_ASSERT(_equation)
Definition: defines.h:216
int a_onehit
flag to stop on first hit
Definition: raytrace.h:1586
int(* a_hit)(struct application *, struct partition *, struct seg *)
called when shot hits model
Definition: raytrace.h:1584
struct hit * pt_inhit
IN hit pointer.
Definition: raytrace.h:577
#define HIDDEN
Definition: common.h:86
struct rt_i * a_rt_i
this librt instance
Definition: raytrace.h:1588
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
fastf_t sizeVoxel
Definition: analyze.h:99
#define UNUSED(parameter)
Definition: common.h:239
char * regionName
Definition: analyze.h:89
struct voxelRegion * nextRegion
Definition: analyze.h:91
#define BU_FREE(_ptr, _type)
Definition: malloc.h:229
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
int bu_strcmp(const char *string1, const char *string2)
Definition: str.c:171
point_t r_pt
Point at which ray starts.
Definition: raytrace.h:218
point_t mdl_max
max corner of model bounding RPP
Definition: raytrace.h:1770
#define FMAX(a, b)
Definition: common.h:102
struct voxelRegion * regionList
Definition: analyze.h:101
HIDDEN struct voxelRegion * getRegionByName(struct voxelRegion *head, const char *regionName)
Definition: voxels.c:38
void rt_prep_parallel(struct rt_i *rtip, int ncpu)
int rt_shootray(struct application *ap)
void * a_uptr
application-specific pointer
Definition: raytrace.h:1618
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
int(* a_miss)(struct application *)
called when shot misses
Definition: raytrace.h:1585
int hit_voxelize(struct application *ap, struct partition *PartHeadp, struct seg *segs)
Definition: voxels.c:85
HIDDEN void create_boxes(void *callBackData, int x, int y, int z, const char *a, fastf_t fill)
Definition: voxelize.c:51
struct partition * pt_forw
forwards link
Definition: raytrace.h:574
#define RT_APPLICATION_INIT(_p)
Definition: raytrace.h:1676
fastf_t hit_dist
dist from r_pt to hit_point
Definition: raytrace.h:250
double fastf_t
Definition: defines.h:300
fastf_t regionDistance
Definition: analyze.h:90
#define bu_strdup(s)
Definition: str.h:71