Description of build_1d

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
