BICPL: deform_line.c Source File

00001 #include <volume_io/internal_volume_io.h> 00002 #include <bicpl/deform.h> 00003 00004 private void perturb_line_points( 00005 int axis, 00006 lines_struct *lines, 00007 Real curvature_factors[], 00008 Point new_points[], 00009 Real fractional_step, 00010 Real max_step, 00011 Real max_search_distance, 00012 int degrees_continuity, 00013 deform_data_struct *deform_data, 00014 boundary_definition_struct *boundary_def, 00015 deformation_model_struct *deformation_model, 00016 deform_stats *stats ); 00017 private Real one_iteration_lines( 00018 lines_struct *lines, 00019 deform_struct *deform_parms, 00020 int iteration ); 00021 00022 public void deform_lines( 00023 lines_struct *lines, 00024 deform_struct *deform_parms ) 00025 { 00026 int iteration; 00027 Real max_error; 00028 00029 iteration = 0; 00030 do 00031 { 00032 ++iteration; 00033 00034 max_error = one_iteration_lines( lines, deform_parms, iteration ); 00035 } 00036 while( max_error > deform_parms->stop_threshold && 00037 iteration < deform_parms->max_iterations ); 00038 } 00039 00040 public void deform_lines_one_iteration( 00041 lines_struct *lines, 00042 deform_struct *deform_parms, 00043 int iteration ) 00044 { 00045 (void) one_iteration_lines( lines, deform_parms, iteration ); 00046 } 00047 00048 private Real one_iteration_lines( 00049 lines_struct *lines, 00050 deform_struct *deform_parms, 00051 int iteration ) 00052 { 00053 int axis; 00054 Point *new_points, *tmp; 00055 Real *curvature_factors; 00056 deform_stats stats; 00057 00058 if( lines->n_items > 1 || 00059 lines->end_indices[0] != lines->n_points && 00060 lines->end_indices[0] != lines->n_points+1 ) 00061 { 00062 print_error( "Must use on single line.\n" ); 00063 return( 0.0 ); 00064 } 00065 00066 if( lines->n_points <= 0 ) 00067 { 00068 print_error( "Must use on nonempty line.\n" ); 00069 return( 0.0 ); 00070 } 00071 00072 if( !check_correct_deformation_lines( lines, 00073 &deform_parms->deformation_model ) ) 00074 { 00075 return( 0.0 ); 00076 } 00077 00078 ALLOC( new_points, lines->n_points ); 00079 ALLOC( curvature_factors, lines->n_points ); 00080 00081 axis = find_axial_plane( lines ); 00082 00083 initialize_deform_stats( &stats ); 00084 00085 perturb_line_points( axis, lines, curvature_factors, new_points, 00086 deform_parms->fractional_step, 00087 deform_parms->max_step, 00088 deform_parms->max_search_distance, 00089 deform_parms->degrees_continuity, 00090 &deform_parms->deform_data, 00091 &deform_parms->boundary_definition, 00092 &deform_parms->deformation_model, 00093 &stats ); 00094 00095 tmp = lines->points; 00096 lines->points = new_points; 00097 new_points = tmp; 00098 00099 print( "Iteration %d: ", iteration ); 00100 print_deform_stats( &stats, lines->n_points ); 00101 00102 FREE( new_points ); 00103 FREE( curvature_factors ); 00104 00105 return( stats.maximum ); 00106 } 00107 00108 public void get_line_equilibrium_point( 00109 lines_struct *lines, 00110 int axis, 00111 int point_index, 00112 int neighbours[], 00113 Real curvature_factors[], 00114 Real max_search_distance, 00115 int degrees_continuity, 00116 Volume volume, 00117 Volume label_volume, 00118 boundary_definition_struct *boundary_def, 00119 deformation_model_struct *deformation_model, 00120 Point *equilibrium_point, 00121 Point *boundary_point ) 00122 { 00123 BOOLEAN found_flag; 00124 Real base_length, model_distance, boundary_distance; 00125 Point centroid, model_point, search_origin; 00126 Vector normal, pos_model_dir, neg_model_dir; 00127 00128 compute_line_centroid_and_normal( lines, axis, neighbours[0], neighbours[1], 00129 &centroid, &normal, &base_length ); 00130 00131 get_model_point( deformation_model, lines->points, point_index, 00132 2, neighbours, curvature_factors, 00133 &centroid, &normal, base_length, &model_point ); 00134 00135 compute_model_dirs( &centroid, &normal, base_length, &model_point, 00136 &model_distance, &search_origin, 00137 &pos_model_dir, &neg_model_dir ); 00138 00139 found_flag = find_boundary_in_direction( volume, label_volume, NULL, 00140 NULL, NULL, model_distance, 00141 &search_origin, &pos_model_dir, &neg_model_dir, 00142 max_search_distance, max_search_distance, 00143 degrees_continuity, 00144 boundary_def, &boundary_distance ); 00145 00146 if( boundary_point != (Point *) NULL ) 00147 get_model_shape_point( &model_point, &pos_model_dir, &neg_model_dir, 00148 boundary_distance, boundary_point ); 00149 00150 compute_equilibrium_point( point_index, 00151 found_flag, boundary_distance, base_length, 00152 model_distance, &pos_model_dir, &neg_model_dir, 00153 &centroid, deformation_model, 00154 equilibrium_point ); 00155 } 00156 00157 private void perturb_line_points( 00158 int axis, 00159 lines_struct *lines, 00160 Real curvature_factors[], 00161 Point new_points[], 00162 Real fractional_step, 00163 Real max_step, 00164 Real max_search_distance, 00165 int degrees_continuity, 00166 deform_data_struct *deform_data, 00167 boundary_definition_struct *boundary_def, 00168 deformation_model_struct *deformation_model, 00169 deform_stats *stats ) 00170 { 00171 int point_index, vertex_index; 00172 int size, start_index, end_index; 00173 int neighbours[2], i; 00174 Point equilibrium_point; 00175 BOOLEAN closed; 00176 progress_struct progress; 00177 Real dist_from_equil; 00178 00179 for_less( point_index, 0, lines->n_points ) 00180 new_points[point_index] = lines->points[point_index]; 00181 00182 i = 0; 00183 size = GET_OBJECT_SIZE( *lines, i ); 00184 00185 closed = (size == lines->n_points+1); 00186 00187 if( closed ) 00188 { 00189 start_index = 0; 00190 end_index = size-2; 00191 } 00192 else 00193 { 00194 start_index = 1; 00195 end_index = size-2; 00196 } 00197 00198 if( deformation_model_includes_average(deformation_model) ) 00199 { 00200 for_inclusive( vertex_index, start_index, end_index ) 00201 { 00202 get_neighbours_of_line_vertex( lines, vertex_index, neighbours ); 00203 point_index = lines->indices[vertex_index]; 00204 curvature_factors[point_index] = compute_line_curvature( 00205 lines, axis, point_index, neighbours[0], neighbours[1] ); 00206 } 00207 } 00208 00209 initialize_progress_report( &progress, TRUE, end_index-start_index+1, 00210 "Deforming Line Points" ); 00211 00212 for_inclusive( vertex_index, start_index, end_index ) 00213 { 00214 get_neighbours_of_line_vertex( lines, vertex_index, neighbours ); 00215 point_index = lines->indices[vertex_index]; 00216 00217 get_line_equilibrium_point( lines, axis, point_index, neighbours, 00218 curvature_factors, 00219 max_search_distance, degrees_continuity, 00220 deform_data->volume, 00221 deform_data->label_volume, 00222 boundary_def, 00223 deformation_model, &equilibrium_point, 00224 (Point *) NULL ); 00225 00226 dist_from_equil = deform_point( point_index, lines->points, 00227 &equilibrium_point, 00228 fractional_step, max_step, 00229 deformation_model->position_constrained, 00230 deformation_model->max_position_offset, 00231 deformation_model->original_positions, 00232 &new_points[point_index] ); 00233 00234 record_error_in_deform_stats( stats, dist_from_equil ); 00235 00236 update_progress_report( &progress, vertex_index - start_index + 1 ); 00237 } 00238 00239 terminate_progress_report( &progress ); 00240 } 00241 00242 public int find_axial_plane( 00243 lines_struct *lines ) 00244 { 00245 int axis, p; 00246 BOOLEAN found_axis; 00247 00248 found_axis = FALSE; 00249 00250 for_less( axis, 0, 3 ) 00251 { 00252 found_axis = TRUE; 00253 for_less( p, 0, lines->n_points-1 ) 00254 { 00255 if( Point_coord(lines->points[p],axis) != 00256 Point_coord(lines->points[p+1],axis) ) 00257 { 00258 found_axis = FALSE; 00259 break; 00260 } 00261 } 00262 if( found_axis ) break; 00263 } 00264 00265 if( !found_axis ) 00266 { 00267 print_error( "No axis found.\n" ); 00268 axis = X; 00269 } 00270 00271 return( axis ); 00272 }