BUILD_1D_MODEL: Builds a 1-D model and returns some related statistics.
This is the main routine to be called in the code.
Code written by Katherine Smith, 2003
GENERAL
build_1d_model(spline_type, placement_type, objective_function)
INPUT/S
-spline_type:
The type of spline to be used. 'single_point' or
'multi_point' at present.
-placement_type:
Placement type for 'multi_point':
'grid' or 'edges' or 'random'
placement type for 'single_point':
'overlap_grid' or 'edges_and_scale'
or 'random_and_scale'.
-objective_function:
'msd_opt_separate' or 'model_opt_together'.
OUTPUT/S
-None. Figures produced and results displayed are main
data to be utilised.
PENDING WORK
-Improve readability of file.
-Add comments
KNOWN BUG/S
-The algorithm is known to be very slow. Some profiling
could possibly be used to hasten the process.
COMMENT/S
-This function needs some work. It is messy in my opinion.
-Smith: the function makes images and points
-Many comments added 27/11/03
-28/11/03: more comments added.
RELATED FUNCTION/S
See other functions in same directory, filenames
beginning with "build" in particular. There seem to be much
overlap due to cut-and-paste.
ABOUT
-Created: November 23rd, 2003
-Last update: November 28th, 2003
-Revision: 0.0.7
-Author: R. S. Schestowitz, University of Manchester
==============================================================0001 function build_1d_model(spline_type, placement_type, objective_function) 0002 % BUILD_1D_MODEL: Builds a 1-D model and returns some related statistics. 0003 % This is the main routine to be called in the code. 0004 % 0005 % Code written by Katherine Smith, 2003 0006 % 0007 % GENERAL 0008 % 0009 % build_1d_model(spline_type, placement_type, objective_function) 0010 % 0011 % INPUT/S 0012 % 0013 % -spline_type: 0014 % The type of spline to be used. 'single_point' or 0015 % 'multi_point' at present. 0016 % 0017 % -placement_type: 0018 % Placement type for 'multi_point': 0019 % 'grid' or 'edges' or 'random' 0020 % placement type for 'single_point': 0021 % 'overlap_grid' or 'edges_and_scale' 0022 % or 'random_and_scale'. 0023 % 0024 % -objective_function: 0025 % 'msd_opt_separate' or 'model_opt_together'. 0026 % 0027 % OUTPUT/S 0028 % 0029 % -None. Figures produced and results displayed are main 0030 % data to be utilised. 0031 % 0032 % PENDING WORK 0033 % 0034 % -Improve readability of file. 0035 % -Add comments 0036 % 0037 % KNOWN BUG/S 0038 % 0039 % -The algorithm is known to be very slow. Some profiling 0040 % could possibly be used to hasten the process. 0041 % 0042 % COMMENT/S 0043 % 0044 % -This function needs some work. It is messy in my opinion. 0045 % -Smith: the function makes images and points 0046 % -Many comments added 27/11/03 0047 % -28/11/03: more comments added. 0048 % 0049 % RELATED FUNCTION/S 0050 % 0051 % See other functions in same directory, filenames 0052 % beginning with "build" in particular. There seem to be much 0053 % overlap due to cut-and-paste. 0054 % 0055 % ABOUT 0056 % 0057 % -Created: November 23rd, 2003 0058 % -Last update: November 28th, 2003 0059 % -Revision: 0.0.7 0060 % -Author: R. S. Schestowitz, University of Manchester 0061 % ============================================================== 0062 0063 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0064 % variables (parameters) initialisation 0065 0066 frame_counter = 0; 0067 n_images = 10; 0068 % how many images to try to register in each set, default was 10 0069 n_iterations = 10; 0070 % how many iterations of the optimisation to run, default was 50. 0071 n_sets = 1; 0072 % [used to be 25] 0073 % how many sets of images to register 0074 % (in order to get better statistics). 0075 image_width = 50; 0076 % spline type 'single_point' or 'multi_point' 0077 spline_type = 'multi_point'; 0078 % placement type for multipoint: 'grid' or 'edges' or 'random' 0079 % placement type for singlepoint: 'overlap_grid' or 'edges_and_scale' or 'random_and_scale' 0080 placement_type = 'grid'; 0081 objective_function = 'msd_opt_together'; 0082 verbose = 'off'; 0083 % be verbose or not 0084 n_points = 5; 0085 % an argument that is passed to the model optimisation function. 0086 % meaning still unknown. 0087 do_plot = 1; 0088 % should model be plotted as it changes? (Boolean) 0089 0090 white_ctr = 0; 0091 % white counter? 0092 ref_shift = 0.2; 0093 % shift in reference image?? 0094 max_shift = 0.2; 0095 % maximum allowed shift? 0096 step = 0.1; 0097 % used below only 0098 los = zeros(n_images,floor((max_shift-ref_shift)/step)); 0099 his = zeros(n_images,floor((max_shift-ref_shift)/step)); 0100 % not understood yet 0101 blurred = 0; 0102 % image blurring enabled??? 0103 spec_iters = 25; 0104 % specificity iterations 0105 gen_iters = 25; 0106 % generalisability iterations 0107 min_error = 0; 0108 max_error = 1; 0109 % define error range?? 0110 error_step = 0.1; 0111 n_error_steps = (max_error-min_error)/error_step; 0112 0113 % END OF INITIALISATION 0114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0115 0116 tic; 0117 % start counting time 0118 for this_set = 1:n_sets 0119 % for all sets 0120 tic; 0121 % why is this called again?? 0122 % Smith: for each set of n_sets images, build a model 0123 % calculate obj fn values from these models, plot the mean with error 0124 % bars given by 1 stddev away 0125 ['calculating set ', num2str(this_set)] 0126 % user gets status 0127 [imagelist images points his los] = make_1d_images(n_images, image_width, 0.2); 0128 % create a set of images 0129 % window = 9; 0130 % images = average_smooth(images, window); blurred = 1; 0131 % images = gaussian_smooth(images, window); blurred = 1; 0132 % smooth all images - currently disabled 0133 0134 points = -1 + 2*(points-1)/(image_width-1); 0135 % Smith: normalise points from -1 to 1 0136 keep = 0.999999; 0137 % is this the precision required??? 0138 ref_points_vec = points(:,1); 0139 ref_image_vec = images(:,1); 0140 % set first image generated to be reference 0141 ref_hi = his(1); 0142 ref_lo = los(1); 0143 % and get its upper and lower boundaries 0144 0145 if(do_plot) 0146 subplot_fig = figure; 0147 H=figure(subplot_fig); 0148 for i=1:n_images, 0149 subplot(n_images,2,(2*i)-1); 0150 plot(images(:,i)),title('Unwarped images'); 0151 end 0152 end 0153 %% NOTE: Change above in traversal (RSS 26/11/03) 0154 0155 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0156 % Warping images, then modelling warped images 0157 % attempt to register images by warping. 0158 0159 warped_images = images; 0160 warped_points = points; 0161 % set this initial state to plot the examples before change 0162 0163 0164 % correctly_warped_points = points; 0165 % correctly_warped_images = images; 0166 % for i=1:n_images 0167 % % now bung in the right deformations to get them back 0168 % correctly_warped_points(:,i) = linear_warp(ref_points_vec, los(i,white_ctr), his(i,white_ctr), ref_lo, ref_hi); 0169 % correctly_warped_points_vec = correctly_warped_points(:,i); 0170 % image_vec = images(:,i); 0171 % % resample at the warped points - point on reg. grid has value of point at same index in warpy grid 0172 % % need to interpolate 0173 % correctly_warped_images(:,i) = interp1(ref_points_vec, image_vec, correctly_warped_points_vec); 0174 % correct_model = bufild_model(correctly_warped_images,correctly_warped_points,1,'','edge',0); 0175 % converging_score(i) = measure_model(correct_model.variances,50); 0176 % end 0177 % figure,plot(log(converging_score)),title('Log of score as model converges'); 0178 0179 0180 for n=1:n_iterations 0181 % iterations aim to get good statistical results by averaging 0182 disp(['iter ',num2str(n)]); 0183 % Smith: correct registration 0184 % attempt to register 0185 for i=1:n_images 0186 disp(['Warping image ',num2str(i),' of ',num2str(n_images)]); 0187 image_vec = warped_images(:,i); 0188 points_vec = warped_points(:,i); 0189 if(strcmp(objective_function,'model_opt_together')) 0190 if(do_plot) 0191 frame_counter = frame_counter + 1; 0192 M(frame_counter) = getframe(H); 0193 [param_list, warped_point, warped_image, score] = optimise_all_warps_model([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot, subplot_fig,i*2); 0194 else 0195 [param_list, warped_point, warped_image, score] = optimise_all_warps_model([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot); 0196 end 0197 elseif(strcmp(objective_function,'model_opt_separate')) 0198 if(do_plot) 0199 frame_counter = frame_counter + 1; 0200 M(frame_counter) = getframe(H); 0201 [param_list, warped_point, warped_image, score] = optimise_warps_model([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot, subplot_fig,i*2); 0202 else 0203 [param_list, warped_point, warped_image, score] = optimise_warps_model([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot); 0204 end 0205 elseif(strcmp(objective_function,'msd_opt_together')) 0206 if(do_plot) 0207 frame_counter = frame_counter + 1; 0208 M(frame_counter) = getframe(H); 0209 [param_list, warped_point, warped_image, score] = optimise_all_warps_msd(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig,i*2); 0210 else 0211 [param_list, warped_point, warped_image, score] = optimise_all_warps_msd(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot); 0212 end 0213 elseif(strcmp(objective_function,'msd_opt_separate')) 0214 if(do_plot) 0215 [param_list, warped_point, warped_image, score] = optimise_warps_msd(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig,i*2); 0216 frame_counter = frame_counter + 1; 0217 M(frame_counter) = getframe(H); 0218 else 0219 [param_list, warped_point, warped_image, score] = optimise_warps_msd(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot); 0220 end 0221 else 0222 error(['Unknown objective function: ', objective_function]); 0223 end 0224 0225 warped_images(:,i) = warped_image; 0226 warped_points(:,i) = warped_point; 0227 0228 % figure(subplot_fig),title(['Warped images, after image ',num2str(i),' iteration ',num2str(n)]); 0229 % dummy = waitforbuttonpress; 0230 if(n==n_iterations), 0231 final_score(i, this_set) = score; 0232 end 0233 end 0234 % end:images 0235 % filename = ['ex1num' num2str(n) '.png']; 0236 % set name of file 0237 %imwrite(warped_images(:,3), filename, 'png'); 0238 end 0239 % end:iteration 0240 0241 0242 % Smith: build model from warped images 0243 0244 w_c_model = build_model(warped_images, warped_points, keep,'Optimised warp', 'variance'); 0245 % This is where a model is being built for the warped images 0246 % [w_c_model.mean_specificity, w_c_model.std_specificity] = find_specificity(w_c_model, images, spec_iters, ref_points_vec); 0247 % show_combined_model(w_c_model,ref_points_vec, 2, 0.2, 'Combined model built from optimised warps'); 0248 0249 0250 0251 fig_title = 'Automatically aligned: '; 0252 % show_shape_model(w_shape_model, ref_points_vec, ref_image_vec, 2, white_width) 0253 % show_intensity_model(w_intensity_model, 2, white_width, fig_title) 0254 % show_combined_model(w_c_model, ref_points_vec, 2, white_width, fig_title) 0255 0256 w_intensity_total_vars = w_c_model.intensity_model.total_var; 0257 w_shape_total_vars = w_c_model.shape_model.total_var; 0258 % get these two variances to be used later to estimate ''goodness' 0259 % of the model, possibly using determinant of covariance matrix 0260 0261 model_score(this_set) = measure_model(w_c_model.variances, 20); 0262 msd_score(this_set) = measure_model_msd(warped_images); 0263 shape_modes(this_set) = w_c_model.n_shape_modes; 0264 intensity_modes(this_set) = size(w_c_model.intensity_model.pcs,2); 0265 shape_variance(this_set) = sum(w_c_model.shape_model.variances); 0266 intensity_variance(this_set) = sum(w_c_model.intensity_model.variances); 0267 0268 t = toc; 0269 disp(['Time for set ',num2str(this_set),': ',num2str(t)]); 0270 time(this_set)=t; 0271 % Get and show time statistics 0272 0273 end % end set 0274 0275 0276 % END OF MODEL-BUILDING 0277 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0278 0279 0280 0281 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0282 % now provide some statistics on the simulation/experiment 0283 0284 0285 0286 0287 disp('Mean match score:'); 0288 mean(final_score(:)) 0289 disp('Std match score:'); 0290 std(final_score(:)) 0291 0292 disp('Mean model score:'); 0293 mean(model_score(:)) 0294 disp('Std model score:'); 0295 std(model_score(:)) 0296 0297 disp('Mean msd score:'); 0298 mean(msd_score(:)) 0299 disp('Std msd score:'); 0300 std(msd_score(:)) 0301 0302 disp('Mean shape modes:'); 0303 mean(shape_modes(:)) 0304 disp('Std shape modes:'); 0305 std(shape_modes(:)) 0306 0307 disp('Mean intensity modes:'); 0308 mean(intensity_modes(:)) 0309 disp('Std intensity modes:'); 0310 std(intensity_modes(:)) 0311 0312 disp('Mean shape variance:'); 0313 mean(shape_variance(:)) 0314 disp('Std shape variance:'); 0315 std(shape_variance(:)) 0316 0317 disp('Mean intensity variance:'); 0318 mean(intensity_variance(:)) 0319 disp('Std intensity variance:'); 0320 std(intensity_variance(:)) 0321 0322 disp('Mean time: '); 0323 mean(time(:)) 0324 disp('Std time: '); 0325 std(time(:)) 0326 0327 t = toc; 0328 % get time from the counter initiated by <tic> 0329 disp(['Total time: ',num2str(t)]); 0330 % display total time 0331 0332 movie(M); 0333 movie2avi(M,'set10.avi'); 0334 0335 % 0336 % END OF STATISTICS 0337 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%