00001
#include <volume_io/internal_volume_io.h>
00002
#include <bicpl/deform.h>
00003
00004 public void get_model_shape_point(
00005 Point *origin,
00006 Vector *pos_model_dir,
00007 Vector *neg_model_dir,
00008 Real dist,
00009 Point *point )
00010 {
00011 Vector offset;
00012
00013
if( dist < 0.0 )
00014 {
00015 SCALE_VECTOR( offset, *neg_model_dir, dist );
00016 }
00017
else
00018 {
00019 SCALE_VECTOR( offset, *pos_model_dir, dist );
00020 }
00021
00022 ADD_POINT_VECTOR( *point, *origin, offset );
00023 }
00024
00025 public void compute_equilibrium_point(
00026
int point_index,
00027 BOOLEAN boundary_exists,
00028 Real boundary_dist,
00029 Real base_length,
00030 Real model_dist,
00031 Vector *pos_model_dir,
00032 Vector *neg_model_dir,
00033 Point *centroid,
00034
deformation_model_struct *deformation_model,
00035 Point *equilibrium_point )
00036 {
00037 Real curvature_offset, equil_dist;
00038
deform_model_struct *model;
00039
00040 model =
find_relevent_model( deformation_model, point_index );
00041
00042
00043
00044
if( !boundary_exists )
00045 boundary_dist = model_dist;
00046
00047 equil_dist = INTERPOLATE( model->
model_weight, boundary_dist, model_dist );
00048
00049
if( model->
min_curvature_offset <= model->
max_curvature_offset )
00050 {
00051 curvature_offset = (equil_dist - model_dist) / base_length;
00052
00053
if( curvature_offset < model->
min_curvature_offset ||
00054 curvature_offset > model->
max_curvature_offset )
00055 {
00056
if( curvature_offset < model->
min_curvature_offset )
00057 curvature_offset = model->
min_curvature_offset;
00058
else
00059 curvature_offset = model->
max_curvature_offset;
00060
00061 equil_dist = model_dist + curvature_offset * base_length;
00062 }
00063 }
00064
00065
get_model_shape_point( centroid, pos_model_dir, neg_model_dir,
00066 equil_dist, equilibrium_point );
00067 }
00068
00069 public void compute_model_dirs(
00070 Point *centroid,
00071 Vector *normal,
00072 Real base_length,
00073 Point *model_point,
00074 Real *model_dist,
00075 Point *search_origin,
00076 Vector *pos_model_dir,
00077 Vector *neg_model_dir )
00078 {
00079 Real to_model_dot_normal, normal_dot_normal, len;
00080 Real ratio, left_length, right_length, offset;
00081 Vector centroid_to_model, offset_vec, tmp, base_vector, perp;
00082
00083
if( EQUAL_POINTS( *centroid, *model_point ) )
00084 {
00085 *model_dist = 0.0;
00086 *search_origin = *centroid;
00087 *pos_model_dir = *normal;
00088 *neg_model_dir = *normal;
00089
return;
00090 }
00091
00092 SUB_POINTS( centroid_to_model, *model_point, *centroid );
00093 CROSS_VECTORS( perp, *normal, centroid_to_model );
00094 CROSS_VECTORS( base_vector, perp, *normal );
00095 len = MAGNITUDE( base_vector );
00096
00097
if( len == 0.0 )
00098 {
00099 *model_dist = DOT_VECTORS( centroid_to_model, *normal );
00100 *search_origin = *centroid;
00101 *pos_model_dir = *normal;
00102 *neg_model_dir = *normal;
00103
return;
00104 }
00105
00106 SCALE_VECTOR( base_vector, base_vector, base_length / 2.0 / len );
00107
00108 SUB_VECTORS( offset_vec, centroid_to_model, base_vector );
00109 left_length = MAGNITUDE( offset_vec );
00110 ADD_VECTORS( offset_vec, centroid_to_model, base_vector );
00111 right_length = MAGNITUDE( offset_vec );
00112
00113
if( left_length + right_length == 0.0 )
00114 {
00115 HANDLE_INTERNAL_ERROR(
"compute model dirs" );
00116 *model_dist = 0.0;
00117 *search_origin = *centroid;
00118 *pos_model_dir = *normal;
00119 *neg_model_dir = *normal;
00120
return;
00121 }
00122
00123 ratio = left_length / (left_length + right_length);
00124 ratio = 1.0 - 2.0 * ratio;
00125
00126
if( ratio < -0.0001 || ratio > 1.0001 )
00127 {
00128 HANDLE_INTERNAL_ERROR(
"compute model dirs ratio" );
00129 *model_dist = 0.0;
00130 *search_origin = *centroid;
00131 *pos_model_dir = *normal;
00132 *neg_model_dir = *normal;
00133
return;
00134 }
00135
00136 SCALE_VECTOR( offset_vec, base_vector, ratio );
00137 ADD_POINT_VECTOR( *search_origin, *centroid, offset_vec );
00138 SUB_POINT_VECTOR( *pos_model_dir, *model_point, *search_origin );
00139 *model_dist = MAGNITUDE( *pos_model_dir );
00140 NORMALIZE_VECTOR( *pos_model_dir, *pos_model_dir );
00141
00142 to_model_dot_normal = DOT_VECTORS( *pos_model_dir, *normal );
00143 normal_dot_normal = DOT_VECTORS( *normal, *normal );
00144
if( normal_dot_normal == 0.0 )
00145 normal_dot_normal = 1.0;
00146 offset = -2.0 * to_model_dot_normal / normal_dot_normal;
00147
00148 SCALE_VECTOR( offset_vec, *normal, offset );
00149 ADD_VECTORS( *neg_model_dir, *pos_model_dir, offset_vec );
00150 NORMALIZE_VECTOR( *neg_model_dir, *neg_model_dir );
00151
00152
if( DOT_VECTORS( *pos_model_dir, *normal ) < 0.0 )
00153 {
00154 tmp = *pos_model_dir;
00155 *pos_model_dir = *neg_model_dir;
00156 *neg_model_dir = tmp;
00157 *model_dist = -(*model_dist);
00158 }
00159
00160 SCALE_VECTOR( *neg_model_dir, *neg_model_dir, -1.0 );
00161 }
00162
00163 private void get_coords_in_nonortho_system(
00164 Vector *v,
00165 Vector *v1,
00166 Vector *v2,
00167 Vector *normalized_v3,
00168 Real *d1,
00169 Real *d2,
00170 Real *d3 )
00171 {
00172 Real mag_v1, mag_c, fact;
00173 Vector offset, v2d, c;
00174
00175 *d3 = DOT_VECTORS( *v, *normalized_v3 );
00176
00177 SCALE_VECTOR( offset, *normalized_v3, *d3 );
00178 SUB_VECTORS( v2d, *v, offset );
00179
00180 mag_v1 = MAGNITUDE( *v1 );
00181
if( mag_v1 == 0.0 )
00182 mag_v1 = 1.0;
00183
00184
if( MAGNITUDE( v2d ) / mag_v1 < 1.0e-5 )
00185 {
00186 *d2 = 0.0;
00187 *d1 = 0.0;
00188
return;
00189 }
00190
00191 fact = DOT_VECTORS( v2d, *v1 ) / mag_v1;
00192 SCALE_VECTOR( offset, *v1, fact/mag_v1 );
00193 SUB_VECTORS( c, v2d, offset );
00194 mag_c = MAGNITUDE( c );
00195
00196
if( mag_c / mag_v1 < 1.0e-5 )
00197 *d2 = 0.0;
00198
else
00199 {
00200 fact = DOT_VECTORS( *v2, c ) / mag_c;
00201
if( fact == 0.0 )
00202 *d2 = 0.0;
00203
else
00204 *d2 = mag_c / fact;
00205 }
00206
00207 SCALE_VECTOR( offset, *v2, *d2 );
00208 SUB_VECTORS( c, v2d, offset );
00209 mag_c = DOT_VECTORS( c, *v1 ) / mag_v1;
00210 *d1 = mag_c / mag_v1;
00211 }
00212
00213 #define MAX_MODEL_NEIGHS 2000
00214
00215 private void compute_model_point(
00216 Point points[],
00217
int point_index,
00218
int n_neighbours,
00219
int neighbours[],
00220 Point *centroid,
00221 Vector *normal,
00222 Point model_points[],
00223 Point model_centroids[],
00224 Vector model_normals[],
00225 Point *model_point )
00226 {
00227
int i;
00228 Real model_len_to_neighbour, len_to_neighbour, scale_factor;
00229 Real d1, d2, d3;
00230 Vector to_model_neighbour;
00231 Vector to_neighbour, model_offset;
00232 Vector to_model_point, third_model_vector, third_vector;
00233 Vector offset;
00234
00235 fill_Point( model_offset, 0.0, 0.0, 0.0 );
00236
00237 SUB_POINTS( to_model_point, model_points[point_index],
00238 model_centroids[point_index] );
00239
00240 for_less( i, 0, n_neighbours )
00241 {
00242 SUB_POINTS( to_model_neighbour,
00243 model_points[neighbours[i]], model_centroids[point_index] );
00244
00245 CROSS_VECTORS( third_model_vector,
00246 to_model_neighbour, model_normals[point_index] );
00247 NORMALIZE_VECTOR( third_model_vector, third_model_vector );
00248
00249
get_coords_in_nonortho_system( &to_model_point,
00250 &to_model_neighbour,
00251 &model_normals[point_index],
00252 &third_model_vector,
00253 &d1, &d2, &d3 );
00254
00255 model_len_to_neighbour = MAGNITUDE( to_model_neighbour );
00256
if( model_len_to_neighbour == 0.0 )
00257 model_len_to_neighbour = 1.0;
00258
00259 SUB_POINTS( to_neighbour, points[neighbours[i]], *centroid );
00260 CROSS_VECTORS( third_vector, to_neighbour, *normal );
00261 NORMALIZE_VECTOR( third_vector, third_vector );
00262
00263 len_to_neighbour = MAGNITUDE( to_neighbour );
00264 scale_factor = len_to_neighbour / model_len_to_neighbour;
00265 d2 *= scale_factor;
00266 d3 *= scale_factor;
00267
00268 SCALE_VECTOR( offset, to_neighbour, d1 );
00269 ADD_VECTORS( model_offset, model_offset, offset );
00270 SCALE_VECTOR( offset, *normal, d2 );
00271 ADD_VECTORS( model_offset, model_offset, offset );
00272 SCALE_VECTOR( offset, third_vector, d3 );
00273 ADD_VECTORS( model_offset, model_offset, offset );
00274 }
00275
00276 SCALE_VECTOR( model_offset, model_offset, 1.0 / (Real) n_neighbours );
00277
00278 ADD_POINT_VECTOR( *model_point, *centroid, model_offset );
00279
00280
#ifdef DEBUG
00281
if(
distance_between_points( model_point, &points[point_index] ) > 0.0001 )
00282 {
00283 print(
"%d: %g %g %g %g %g %g\n", point_index,
00284 Point_x(*model_point ),
00285 Point_y(*model_point ),
00286 Point_z(*model_point ),
00287 Point_x(points[point_index]),
00288 Point_y(points[point_index]),
00289 Point_z(points[point_index]) );
00290 HANDLE_INTERNAL_ERROR(
"dang" );
00291 }
00292
#endif
00293
}
00294
00295 public void get_model_point(
00296
deformation_model_struct *deformation_model,
00297 Point points[],
00298
int point_index,
00299
int n_neighbours,
00300
int neighbours[],
00301 Real curvatures[],
00302 Point *centroid,
00303 Vector *normal,
00304 Real base_length,
00305 Point *model_point )
00306 {
00307
int i;
00308 Real curvature, dist_from_centroid, normal_mag;
00309 Vector offset;
00310
deform_model_struct *model;
00311
Deformation_model_types model_type;
00312
00313 model =
find_relevent_model( deformation_model, point_index );
00314
00315 model_type = model->
model_type;
00316
00317
if( model_type ==
PARAMETRIC_MODEL && model->
n_model_points == 0 )
00318 model_type =
FLAT_MODEL;
00319
00320
switch( model_type )
00321 {
00322
case FLAT_MODEL:
00323 *model_point = *centroid;
00324
break;
00325
00326
case AVERAGE_MODEL:
00327 curvature = 0.0;
00328 for_less( i, 0, n_neighbours )
00329 curvature += curvatures[neighbours[i]];
00330 curvature /= (Real) n_neighbours;
00331
00332 dist_from_centroid = curvature * base_length;
00333 normal_mag = MAGNITUDE( *normal );
00334
if( normal_mag == 0.0 ) normal_mag = 1.0;
00335
00336 SCALE_VECTOR( offset, *normal, dist_from_centroid / normal_mag );
00337
00338 ADD_POINT_VECTOR( *model_point, *centroid, offset );
00339
break;
00340
00341
case PARAMETRIC_MODEL:
00342
case GENERAL_MODEL:
00343
compute_model_point( points, point_index, n_neighbours, neighbours,
00344 centroid, normal,
00345 model->
model_points,
00346 model->
model_centroids,
00347 model->
model_normals,
00348 model_point );
00349
break;
00350
00351
default:
00352 HANDLE_INTERNAL_ERROR(
"get_model_point" );
00353
break;
00354 }
00355 }
00356
00357 public void get_neighbours_of_line_vertex(
00358
lines_struct *lines,
00359
int vertex_index,
00360
int neighbours[2] )
00361 {
00362
int prev_index, next_index;
00363
00364
if( vertex_index == 0 )
00365 prev_index = lines->
n_points-1;
00366
else
00367 prev_index = (vertex_index - 1 + lines->
n_points) %
00368 lines->
n_points;
00369
00370 next_index = (vertex_index + 1 + lines->
n_points) % lines->
n_points;
00371
00372 neighbours[0] = lines->
indices[prev_index];
00373 neighbours[1] = lines->
indices[next_index];
00374 }
00375
00376 public BOOLEAN
deformation_model_includes_average(
00377
deformation_model_struct *model )
00378 {
00379
int i;
00380
00381 for_less( i, 0, model->
n_models )
00382
if( model->
models[i].
model_type ==
AVERAGE_MODEL )
00383
return(
TRUE );
00384
00385
return(
FALSE );
00386 }
00387
00388 public Real
compute_line_curvature(
00389
lines_struct *lines,
00390
int axis,
00391
int point_index,
00392
int prev_point_index,
00393
int next_point_index )
00394 {
00395 Real base_length, curvature;
00396 Point centroid;
00397 Vector normal, dir_to_point;
00398
00399
compute_line_centroid_and_normal( lines, axis, prev_point_index,
00400 next_point_index, ¢roid, &normal,
00401 &base_length );
00402
00403 SUB_POINTS( dir_to_point, lines->
points[point_index], centroid );
00404
00405 curvature = MAGNITUDE( dir_to_point ) / base_length;
00406
00407
if( DOT_VECTORS( dir_to_point, normal ) < 0.0 )
00408 curvature = -curvature;
00409
00410
return( curvature );
00411 }
00412
00413 public Real
deform_point(
00414
int point_index,
00415 Point points[],
00416 Point *equilibrium_point,
00417 Real fractional_step,
00418 Real max_step,
00419 BOOLEAN position_constrained,
00420 Real max_position_offset,
00421 Point original_positions[],
00422 Point *new_point )
00423 {
00424 Vector toward_equil;
00425 Real dist_to_equil, dist;
00426 Vector diff;
00427 Real factor;
00428
00429 SUB_POINTS( toward_equil, *equilibrium_point, points[point_index] );
00430 dist_to_equil = MAGNITUDE( toward_equil );
00431
00432
if( dist_to_equil * fractional_step > max_step )
00433 fractional_step = max_step / dist_to_equil;
00434
00435 SCALE_VECTOR( toward_equil, toward_equil, fractional_step );
00436
00437 ADD_POINT_VECTOR( *new_point, points[point_index], toward_equil );
00438
00439
if( position_constrained )
00440 {
00441 SUB_POINTS( diff, *new_point, original_positions[point_index] );
00442 dist = MAGNITUDE( diff );
00443
if( dist > max_position_offset )
00444 {
00445 factor = max_position_offset / dist;
00446 SCALE_VECTOR( diff, diff, factor );
00447 ADD_POINT_VECTOR( *new_point, original_positions[point_index],
00448 diff );
00449 }
00450 }
00451
00452
return( dist_to_equil );
00453 }
00454
00455 public void compute_line_centroid_and_normal(
00456
lines_struct *lines,
00457
int axis,
00458
int prev_point_index,
00459
int next_point_index,
00460 Point *centroid,
00461 Vector *normal,
00462 Real *base_length )
00463 {
00464
int a1, a2;
00465 Vector dir;
00466
00467 SUB_POINTS( dir, lines->
points[next_point_index],
00468 lines->
points[prev_point_index] );
00469
00470 *base_length = MAGNITUDE( dir );
00471
if( *base_length == 0.0 )
00472 *base_length = 1.0;
00473
00474 a1 = (axis + 1) % 3;
00475 a2 = (axis + 2) % 3;
00476
00477 Point_coord(*normal,axis) = Point_coord(dir,axis);
00478 Point_coord(*normal,a1) = Point_coord(dir,a2);
00479 Point_coord(*normal,a2) = -Point_coord(dir,a1);
00480 NORMALIZE_VECTOR( *normal, *normal );
00481
00482 INTERPOLATE_POINTS( *centroid,
00483 lines->
points[next_point_index],
00484 lines->
points[prev_point_index], 0.5 );
00485 }
00486
00487 public int get_subsampled_neighbours_of_point(
00488
deformation_model_struct *deformation_model,
00489
polygons_struct *polygons,
00490
int poly,
00491
int vertex_index,
00492
int neighbours[],
00493
int max_neighbours,
00494 BOOLEAN *interior_flag )
00495 {
00496
int point_index, model_vertex_index, model_poly;
00497
polygons_struct *model_polygons;
00498
deform_model_struct *model;
00499
00500 point_index = polygons->
indices[
00501
POINT_INDEX(polygons->
end_indices,poly,vertex_index)];
00502
00503 model =
find_relevent_model( deformation_model, point_index );
00504
00505
if( polygons->
n_points <= model->
up_to_n_points )
00506 {
00507 model_polygons = polygons;
00508 model_poly = poly;
00509 model_vertex_index = vertex_index;
00510 }
00511
else
00512 {
00513 model_polygons =
get_polygons_ptr( model->
model_object );
00514
00515
if( !
find_polygon_with_vertex( model_polygons, point_index,
00516 &model_poly, &model_vertex_index ) )
00517 {
00518 HANDLE_INTERNAL_ERROR(
"get_subsampled_neighbours_of_point" );
00519 }
00520 }
00521
00522
return(
get_neighbours_of_point( model_polygons, model_poly,
00523 model_vertex_index, neighbours,
00524 max_neighbours, interior_flag ) );
00525 }
