BICPL: plane_polygon_intersect.c Source File

00001 #include <volume_io/internal_volume_io.h> 00002 #include <bicpl/geom.h> 00003 00004 private void coalesce_lines( 00005 lines_struct *lines ) 00006 { 00007 int p, n, l, *n_neighbours, **neighbours, total_neighbours; 00008 int size, p1, p2, n_items, n_indices, *indices, *end_indices, edge; 00009 00010 if( lines->n_points <= 0 || 00011 lines->n_items <= 0 ) 00012 { 00013 return; 00014 } 00015 00016 ALLOC( n_neighbours, lines->n_points ); 00017 00018 for_less( p, 0, lines->n_points ) 00019 n_neighbours[p] = 0; 00020 00021 for_less( l, 0, lines->n_items ) 00022 { 00023 size = GET_OBJECT_SIZE( *lines, l ); 00024 00025 for_less( edge, 0, size-1 ) 00026 { 00027 p1 = lines->indices[POINT_INDEX(lines->end_indices,l,edge)]; 00028 p2 = lines->indices[POINT_INDEX(lines->end_indices,l,edge+1)]; 00029 00030 ++n_neighbours[p1]; 00031 ++n_neighbours[p2]; 00032 } 00033 } 00034 00035 total_neighbours = 0; 00036 for_less( p, 0, lines->n_points ) 00037 { 00038 if( n_neighbours[p] != 0 && n_neighbours[p] != 2 ) 00039 { 00040 return; 00041 } 00042 total_neighbours += n_neighbours[p]; 00043 } 00044 00045 ALLOC( neighbours, lines->n_points ); 00046 ALLOC( neighbours[0], total_neighbours ); 00047 for_less( p, 1, lines->n_points ) 00048 { 00049 neighbours[p] = &neighbours[p-1][n_neighbours[p-1]]; 00050 } 00051 00052 for_less( p, 0, lines->n_points ) 00053 n_neighbours[p] = 0; 00054 00055 for_less( l, 0, lines->n_items ) 00056 { 00057 size = GET_OBJECT_SIZE( *lines, l ); 00058 00059 for_less( edge, 0, size-1 ) 00060 { 00061 p1 = lines->indices[POINT_INDEX(lines->end_indices,l,edge)]; 00062 p2 = lines->indices[POINT_INDEX(lines->end_indices,l,edge+1)]; 00063 00064 neighbours[p1][n_neighbours[p1]] = p2; 00065 ++n_neighbours[p1]; 00066 neighbours[p2][n_neighbours[p2]] = p1; 00067 ++n_neighbours[p2]; 00068 } 00069 } 00070 00071 n_items = 0; 00072 end_indices = NULL; 00073 n_indices = 0; 00074 indices = NULL; 00075 00076 for_less( l, 0, lines->n_items ) 00077 { 00078 size = GET_OBJECT_SIZE( *lines, l ); 00079 00080 for_less( edge, 0, size-1 ) 00081 { 00082 p1 = lines->indices[POINT_INDEX(lines->end_indices,l,edge)]; 00083 p2 = lines->indices[POINT_INDEX(lines->end_indices,l,edge+1)]; 00084 00085 for_less( n, 0, n_neighbours[p1] ) 00086 { 00087 if( neighbours[p1][n] == p2 ) 00088 break; 00089 } 00090 00091 if( n >= n_neighbours[p1] ) 00092 continue; 00093 00094 ADD_ELEMENT_TO_ARRAY( indices, n_indices, p1, DEFAULT_CHUNK_SIZE ); 00095 00096 while( n < n_neighbours[p1] ) 00097 { 00098 ADD_ELEMENT_TO_ARRAY( indices, n_indices, p2, 00099 DEFAULT_CHUNK_SIZE ); 00100 00101 neighbours[p1][n] = -1; 00102 00103 for_less( n, 0, n_neighbours[p2] ) 00104 { 00105 if( neighbours[p2][n] == p1 ) 00106 break; 00107 } 00108 00109 if( n >= n_neighbours[p1] ) 00110 { 00111 handle_internal_error( "coalesce_lines" ); 00112 return; 00113 } 00114 00115 neighbours[p2][n] = -1; 00116 00117 p1 = p2; 00118 00119 for_less( n, 0, n_neighbours[p1] ) 00120 { 00121 if( neighbours[p1][n] >= 0 ) 00122 break; 00123 } 00124 00125 p2 = neighbours[p1][n]; 00126 } 00127 00128 ADD_ELEMENT_TO_ARRAY( end_indices, n_items, n_indices, 00129 DEFAULT_CHUNK_SIZE ); 00130 } 00131 } 00132 00133 FREE( lines->end_indices ); 00134 FREE( lines->indices ); 00135 00136 lines->n_items = n_items; 00137 lines->indices = indices; 00138 lines->end_indices = end_indices; 00139 } 00140 00141 public void intersect_planes_with_polygons( 00142 polygons_struct *polygons, 00143 Point *plane_origin, 00144 Vector *plane_normal, 00145 lines_struct *lines ) 00146 { 00147 int n_points, p, index; 00148 int poly, edge, size; 00149 int point_index1, point_index2; 00150 Vector v1, v2; 00151 Point point; 00152 Real t1, t2, ratios[2]; 00153 int p1s[2], p2s[2]; 00154 hash2_table_struct hash; 00155 00156 initialize_lines( lines, WHITE ); 00157 00158 initialize_hash2_table( &hash, polygons->n_items, (int) sizeof(int), 00159 0.5, .25 ); 00160 00161 index = 0; 00162 00163 for_less( poly, 0, polygons->n_items ) 00164 { 00165 n_points = 0; 00166 00167 size = GET_OBJECT_SIZE( *polygons, poly ); 00168 00169 for_less( edge, 0, size ) 00170 { 00171 point_index1 = polygons->indices[ 00172 POINT_INDEX(polygons->end_indices,poly,edge)]; 00173 point_index2 = polygons->indices[ 00174 POINT_INDEX(polygons->end_indices,poly,(edge+1)%size)]; 00175 00176 SUB_POINTS( v1, polygons->points[point_index1], *plane_origin ); 00177 SUB_POINTS( v2, polygons->points[point_index2], *plane_origin ); 00178 00179 t1 = DOT_VECTORS( v1, *plane_normal ); 00180 t2 = DOT_VECTORS( v2, *plane_normal ); 00181 00182 /*--- check if it intersects the line segment [p1 .. p2) */ 00183 00184 if( t1 == 0.0 || (t1 > 0.0 && t2 < 0.0) || (t1 < 0.0 && t2 > 0.0) ) 00185 { 00186 if( n_points < 2 ) 00187 { 00188 if( t1 == 0.0 ) 00189 ratios[n_points] = 0.0; 00190 else 00191 ratios[n_points] = t1 / (t1 - t2); 00192 00193 p1s[n_points] = MIN( point_index1, point_index2 ); 00194 p2s[n_points] = MAX( point_index1, point_index2 ); 00195 if( p1s[n_points] != point_index1 ) 00196 ratios[n_points] = 1.0 - ratios[n_points]; 00197 00198 if( ratios[n_points] == 0.0 ) 00199 { 00200 p2s[n_points] = n_points; 00201 ratios[n_points] = 0.0; 00202 } 00203 else if( ratios[n_points] == 1.0 ) 00204 { 00205 p1s[n_points] = p2s[n_points]; 00206 p2s[n_points] = n_points; 00207 ratios[n_points] = 0.0; 00208 } 00209 } 00210 00211 ++n_points; 00212 } 00213 } 00214 00215 if( n_points == 2 && 00216 (ratios[0] != 0.0 || ratios[1] != 0.0 || p1s[0] != p1s[1]) ) 00217 { 00218 start_new_line( lines ); 00219 00220 for_less( p, 0, 2 ) 00221 { 00222 if( !lookup_in_hash2_table( &hash, p1s[p], p2s[p], 00223 (void *) &index ) ) 00224 { 00225 if( ratios[p] == 0.0 ) 00226 point = polygons->points[p1s[p]]; 00227 else 00228 { 00229 INTERPOLATE_POINTS( point, 00230 polygons->points[p1s[p]], 00231 polygons->points[p2s[p]], 00232 ratios[p] ); 00233 } 00234 00235 index = lines->n_points; 00236 insert_in_hash2_table( &hash, p1s[p], p2s[p], 00237 (void *) &index ); 00238 00239 ADD_ELEMENT_TO_ARRAY( lines->points, lines->n_points, 00240 point, DEFAULT_CHUNK_SIZE ); 00241 } 00242 00243 ADD_ELEMENT_TO_ARRAY( lines->indices, 00244 lines->end_indices[lines->n_items-1], 00245 index, DEFAULT_CHUNK_SIZE ); 00246 } 00247 } 00248 } 00249 00250 delete_hash2_table( &hash ); 00251 00252 coalesce_lines( lines ); 00253 } 00254 00255 public void intersect_planes_with_quadmesh( 00256 quadmesh_struct *quadmesh, 00257 Point *plane_origin, 00258 Vector *plane_normal, 00259 lines_struct *lines ) 00260 { 00261 int n_points, p, index, m, n, x_index, y_index; 00262 int edge; 00263 int point_index1, point_index2, indices[4]; 00264 Vector v1, v2; 00265 Point point, *points; 00266 Real t1, t2, ratios[2]; 00267 int p1s[2], p2s[2]; 00268 hash2_table_struct hash; 00269 00270 initialize_lines( lines, WHITE ); 00271 00272 get_quadmesh_n_objects( quadmesh, &m, &n ); 00273 00274 initialize_hash2_table( &hash, m * n, (int) sizeof(int), 00275 0.5, .25 ); 00276 00277 points = quadmesh->points; 00278 00279 index = 0; 00280 00281 for_less( x_index, 0, m ) 00282 { 00283 for_less( y_index, 0, n ) 00284 { 00285 n_points = 0; 00286 00287 indices[0] = IJ(x_index,y_index,quadmesh->n); 00288 indices[1] = IJ((x_index+1)%quadmesh->m,y_index,quadmesh->n); 00289 indices[2] = IJ((x_index+1)%quadmesh->m,(y_index+1)%quadmesh->n, 00290 quadmesh->n); 00291 indices[3] = IJ(x_index,(y_index+1)%quadmesh->n, 00292 quadmesh->n); 00293 00294 for_less( edge, 0, 4 ) 00295 { 00296 point_index1 = indices[edge]; 00297 point_index2 = indices[(edge+1) % 4]; 00298 00299 SUB_POINTS( v1, points[point_index1], *plane_origin ); 00300 SUB_POINTS( v2, points[point_index2], *plane_origin ); 00301 00302 t1 = DOT_VECTORS( v1, *plane_normal ); 00303 t2 = DOT_VECTORS( v2, *plane_normal ); 00304 00305 /*--- check if it intersects the line segment [p1 .. p2) */ 00306 00307 if( t1 == 0.0 || (t1 > 0.0 && t2 < 0.0) || (t1 < 0.0 && t2 > 0.0) ) 00308 { 00309 if( n_points < 2 ) 00310 { 00311 if( t1 == 0.0 ) 00312 ratios[n_points] = 0.0; 00313 else 00314 ratios[n_points] = t1 / (t1 - t2); 00315 00316 p1s[n_points] = MIN( point_index1, point_index2 ); 00317 p2s[n_points] = MAX( point_index1, point_index2 ); 00318 if( p1s[n_points] != point_index1 ) 00319 ratios[n_points] = 1.0 - ratios[n_points]; 00320 00321 if( ratios[n_points] == 0.0 ) 00322 { 00323 p2s[n_points] = n_points; 00324 ratios[n_points] = 0.0; 00325 } 00326 else if( ratios[n_points] == 1.0 ) 00327 { 00328 p1s[n_points] = p2s[n_points]; 00329 p2s[n_points] = n_points; 00330 ratios[n_points] = 0.0; 00331 } 00332 } 00333 00334 ++n_points; 00335 } 00336 } 00337 00338 if( n_points == 2 && 00339 (ratios[0] != 0.0 || ratios[1] != 0.0 || p1s[0] != p1s[1]) ) 00340 { 00341 start_new_line( lines ); 00342 00343 for_less( p, 0, 2 ) 00344 { 00345 if( !lookup_in_hash2_table( &hash, p1s[p], p2s[p], 00346 (void *) &index ) ) 00347 { 00348 if( ratios[p] == 0.0 ) 00349 point = points[p1s[p]]; 00350 else 00351 { 00352 INTERPOLATE_POINTS( point, 00353 points[p1s[p]], points[p2s[p]], 00354 ratios[p] ); 00355 } 00356 00357 index = lines->n_points; 00358 insert_in_hash2_table( &hash, p1s[p], p2s[p], 00359 (void *) &index ); 00360 00361 ADD_ELEMENT_TO_ARRAY( lines->points, lines->n_points, 00362 point, DEFAULT_CHUNK_SIZE ); 00363 } 00364 00365 ADD_ELEMENT_TO_ARRAY( lines->indices, 00366 lines->end_indices[lines->n_items-1], 00367 index, DEFAULT_CHUNK_SIZE ); 00368 } 00369 } 00370 } 00371 } 00372 00373 delete_hash2_table( &hash ); 00374 00375 coalesce_lines( lines ); 00376 }