Description of eval_1d_app_model_obj

0001 function eval_1d_app_model_obj_fn
0002 % EVAL_1D_APP_MODEL_OBJ_FN: Evaluate 1-D appearance model objective function. In
0003 %                           practice it is a procedure that carries out many
0004 %                           experiments that evaluate the objective function.
0005 %
0006 % Code written by Katherine Smith, 2003
0007 %
0008 %    GENERAL
0009 %
0010 %      eval_1d_app_model_obj_fn
0011 %
0012 %    INPUT/S
0013 %
0014 %      -
0015 %
0016 %    OUTPUT/S
0017 %
0018 %      -
0019 %
0020 %    PENDING WORK
0021 %
0022 %      -
0023 %
0024 %    KNOWN BUG/S
0025 %
0026 %      -
0027 %
0028 %    COMMENT/S
0029 %
0030 %
0031 %
0032 %    RELATED FUNCTION/S
0033 %
0034 %      There are some near duplicates in the same directory.
0035 %
0036 %    ABOUT
0037 %
0038 %      -Created:     November 23rd, 2003
0039 %      -Last update: Novermber 30th, 2003
0040 %      -Revision:    0.0.1
0041 %      -Author:      R. S. Schestowitz, University of Manchester
0042 % ==============================================================
0043 
0044        % Smith: make images and points
0045 n_images = 100;
0046        % Number of images to use for evaluation
0047 n_sets = 5;
0048        % Number of sets of images (for averaging results); deafult was 50
0049 image_width = 50;
0050        % Width of the images
0051 bump_width = 0.2;
0052        % Width of the bump in the images
0053 ref_shift = 0.2;
0054 max_shift = 0.2;
0055 step = 0.1;
0056        % Warping parameters??
0057 los = zeros(n_images,floor((max_shift-ref_shift)/step));
0058 his = zeros(n_images,floor((max_shift-ref_shift)/step));
0059        % Will hold the number of positions where bump is high and low
0060 min_error = 0;
0061 max_error = 1;
0062 error_step = 0.1;
0063        % See section on <Warping images badly,
0064        % then modelling poorly warped images> where these are used
0065 
0066 n_error_steps = (max_error-min_error)/error_step;
0067 range_score = zeros(n_sets, n_error_steps, 11);
0068 factor_step = 0.1;
0069 factor_min = 0;
0070 factor_max = 0.4;
0071 blurred = 0;
0072        % blurred images?? (Boolean)
0073 spec_iters = 10;
0074        % iterations for specificity
0075 gen_iters = 10;
0076        % iterations for generalisability
0077 build_unwarped_model = 0;
0078 build_perfect_warp_model = 0;
0079 find_gen_and_spec = 1;
0080        % Booleans that are self-explanatory. They contol of the flow of this
0081        % routine as it performs many unrelated tests
0082 
0083         
0084        % The following is similar to build_1d_model. REORGANISE?
0085 for this_set = 1:n_sets
0086   tic;
0087   ['calculating set ', num2str(this_set)]
0088   % for each set of n_sets images, build a model
0089   % calculate obj fn values from these models, plot the mean with error
0090   % bars given by 1 stddev away
0091   [imagelist images points his_no_norm los_no_norm] = make_1d_images(n_images, image_width, bump_width);
0092   % smooth images
0093   window = 9;
0094 %    images = average_smooth(images, window); blurred = 1;
0095 %    images = gaussian_smooth(images, window); blurred = 1;
0096 
0097   ref_hi_no_norm = his_no_norm(1);
0098                                       % Set reference to be image 1??
0099   ref_lo_no_norm = los_no_norm(1);
0100   ref_points_vec_no_norm = points(:,1);
0101 
0102   % Smith: normalise points from -1 to 1
0103   points = -1 + 2*(points-1)/(image_width-1);
0104   los = -1 + 2*(los_no_norm-1)/(image_width-1);
0105   his = -1 + 2*(his_no_norm-1)/(image_width-1);
0106          % Is this some normalisation that is required for reasonable warps
0107 
0108 
0109   keep_pts = 0.999999;
0110   keep_int = 0.8;
0111   ref_points_vec = points(:,1);
0112   ref_images_vec = images(:,1);
0113   ref_hi = his(1);
0114   ref_lo = los(1);
0115           % again use image 1 as reference?? Or is it just initialisation?
0116 
0117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0118 % Smith: Modelling unwarped images
0119 
0120   [intensity_model.pcs,intensity_model.variances,intensity_model.params, intensity_model.mean, intensity_model.var_r, intensity_model.total_var, intensity_model.impdata] = st_pca(images', keep_int);
0121   [shape_model.pcs,shape_model.variances,shape_model.params, shape_model.mean, shape_model.var_r, shape_model.total_var, shape_model.impdata] = st_pca(points', keep_pts);
0122           % perform PCA on the images
0123   intensity_model.sd = std(intensity_model.params);
0124   shape_model.sd = std(shape_model.params);
0125           % get standard deviation from the models (of shape and intensity)
0126   intensity_total_vars = intensity_model.total_var;
0127   shape_total_vars = shape_model.total_var;
0128           % get total variance for each model
0129 
0130   if (build_unwarped_model),
0131     % Smith: combine params - initially normalise to variance described in model
0132     total_shape_variance = sum(shape_model.variances);
0133     total_intensity_variance = sum(intensity_model.variances);
0134           % get total variances
0135     combined_model.shape_weight = ones(1,size(shape_model.params,2));
0136           % set weight Ws to initially be 1??
0137     [combined_model.pcs, combined_model.variances, combined_model.params, combined_model.mean, combined_model.total_var, combined_model.impdata] = make_combined_model(shape_model, intensity_model, combined_model.shape_weight);
0138           % get combined model
0139     combined_model.sd = std(combined_model.params);
0140     combined_model.intensity_model = intensity_model;
0141     combined_model.shape_model = shape_model;
0142     combined_model.n_shape_modes = size(shape_model.pcs,2);
0143     combined_model.label = 'No warp';
0144              % set some parameters of the model
0145     [combined_model.mean_specificity, combined_model.std_specificity] = find_specificity(combined_model, images, spec_iters, ref_points_vec);
0146              % find the specificity of the combined model
0147     fig_title = 'No warp: ';
0148 
0149 %  show_intensity_model(intensity_model, 2, fig_title)
0150 
0151       %This isn't very interesting so don't show it!
0152 %  show_shape_model(shape_model, ref_points_vec, ref_images_vec, 2, white_width)
0153 
0154     show_combined_model(combined_model, ref_points_vec, 2, 5, fig_title)
0155     disp(['Score of unwarped model: ', num2str(measure_model(combined_model.variances, 50))]);
0156   end
0157 
0158 
0159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0160 % Smith: Warping images, then modelling warped images
0161 % now warp them using the known correspondences and build another model.
0162 
0163   for i=1:n_images
0164     % now bung in the right deformations to get them back
0165     warped_points(:,i) = -1 + 2*(linear_warp(ref_points_vec_no_norm, los_no_norm(i), his_no_norm(i), ref_lo_no_norm, ref_hi_no_norm)-1)/(image_width-1);
0166             % warping images? Why -1 + 2*warp()???
0167     warped_points_vec = warped_points(:,i);
0168             % make a simple copy (duplicate)
0169     if any(warped_points_vec > 1)
0170       disp('stop!');
0171         % RSS: Why should warp have this effect anyway??
0172     end
0173     image_vec = images(:,i);
0174     % resample at the warped points - point on reg. grid has value of point at same index in warpy grid
0175     % need to interpolate
0176     warped_images(:,i) = interp1(ref_points_vec, image_vec, warped_points_vec, 'linear',0);
0177          % apply the warp to the images.
0178   end
0179   
0180   if(build_perfect_warp_model)
0181       % This flaf appears at the initialisation section of this function
0182       % This is similar to the section above? --> MERGER??
0183     [w_intensity_model.pcs,w_intensity_model.variances,w_intensity_model.params, w_intensity_model.mean, w_intensity_model.var_r, w_intensity_model.total_var, w_intensity_model.impdata] = st_pca(warped_images', keep_int);
0184     [w_shape_model.pcs,w_shape_model.variances,w_shape_model.params, w_shape_model.mean, w_shape_model.var_r, w_shape_model.total_var, w_shape_model.impdata] = st_pca(warped_points', keep_pts);
0185     w_intensity_model.sd = std(w_intensity_model.params);
0186     w_shape_model.sd = std(w_shape_model.params);
0187          % As above
0188     
0189     % Smith: combine params - initially normalise to variance described in model
0190     total_shape_variance = sum(w_shape_model.variances);
0191     total_intensity_variance = sum(w_intensity_model.variances);
0192     w_c_model.shape_weight = ones(1,size(w_shape_model.pcs,2));
0193     if(total_shape_variance ~= 0)
0194       w_c_model.shape_weight(:) = sqrt(total_intensity_variance/total_shape_variance);
0195     end
0196            % RSS: What is this about??
0197     [w_c_model.pcs, w_c_model.variances, w_c_model.params, w_c_model.mean, w_c_model.total_var, w_c_model.impdata] = make_combined_model(w_shape_model, w_intensity_model, w_c_model.shape_weight);
0198       % combined model
0199     w_c_model.sd = sqrt(w_c_model.variances);
0200     w_c_model.intensity_model = w_intensity_model;
0201     w_c_model.shape_model = w_shape_model;
0202     w_c_model.n_shape_modes = size(w_shape_model.pcs,2);
0203     w_c_model.label = 'Perfect warp';
0204 %    [w_c_model.mean_specificity, w_c_model.std_specificity] = find_specificity(w_c_model, images, spec_iters, ref_points_vec);
0205 
0206     fig_title = 'Aligned: ';
0207     show_combined_model(w_c_model, ref_points_vec, 2, 5, fig_title)
0208     disp(['Score of perfectly warped model: ', num2str(measure_model(w_c_model.variances, 50))]);
0209   end
0210            % RSS: What is it that makes it 'perfect'?
0211 
0212 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0213 % Warping images badly, then modelling poorly warped images
0214 
0215 % Again this repeats some of the above in an inefficient way
0216 
0217   error_ctr = 1;
0218   for error = min_error:error_step:max_error
0219     ['calculating error ', num2str(error)]
0220 %    for error = -1*(image_width*0.2):image_width*0.2
0221     for im=1:n_images
0222       errors(error_ctr) = error;
0223       % Smith: now bung in deformations to get them back(ish)
0224       badly_warped_points(:,im) = perturb_points_cps(warped_points(:,im), error);
0225       badly_warped_points_vec = badly_warped_points(:,im);
0226       image_vec = images(:,im);
0227       % Smith: resample at the warped points - point on reg. grid has value of point at same index in warpy grid
0228       % need to interpolate
0229       badly_warped_images(:,im) = interp1(ref_points_vec, image_vec, badly_warped_points_vec, 'linear', 0);
0230       % RSS: apply these bad warps by interpolating.
0231     end
0232     [b_w_intensity_model.pcs,b_w_intensity_model.variances,b_w_intensity_model.params, b_w_intensity_model.mean, b_w_intensity_model.var_r, b_w_intensity_model.total_var, b_w_intensity_model.impdata] = st_pca(badly_warped_images', keep_int);
0233     [b_w_shape_model.pcs,b_w_shape_model.variances,b_w_shape_model.params, b_w_shape_model.mean, b_w_shape_model.var_r, b_w_shape_model.total_var, b_w_shape_model.impdata] = st_pca(badly_warped_points', keep_pts);
0234     b_w_intensity_model.sd = std(b_w_intensity_model.params);
0235     b_w_shape_model.sd = std(b_w_shape_model.params);
0236            % Create the models and get the standard deviation of these models
0237     % Smith: try using some different normalisation factors
0238     factor_ctr = 1;
0239           % This <factor> appears to be the coefficient Ws
0240     for i=5:10
0241       factor = 10^i * 1/(10^10); 
0242           % Silly numbers... 'i' could range from 1 to 5 with same effect
0243       models(this_set, error_ctr, factor_ctr).model = build_model_from_models(b_w_shape_model, b_w_intensity_model, ['constant: ', num2str(factor)],'constant', factor);
0244         % create a model from the models available
0245       range_score(this_set, error_ctr, factor_ctr) = measure_model(models(this_set, error_ctr, factor_ctr).model.variances, 50);
0246         % and score it
0247       if(find_gen_and_spec)
0248         [models(this_set, error_ctr, factor_ctr).mean_specificity, models(this_set, error_ctr, factor_ctr).std_specificity] = find_specificity(models(this_set, error_ctr, factor_ctr).model, images, spec_iters, ref_points_vec);
0249         [models(this_set, error_ctr, factor_ctr).mean_generalisability, models(this_set, error_ctr, factor_ctr).std_generalisability] = find_generalisability(models(this_set, error_ctr, factor_ctr).model, gen_iters);
0250       end
0251       factor_ctr = factor_ctr+1;
0252     end
0253 
0254     % Smith: variance
0255     % RSS: Isn't the following the same as above?
0256     
0257     total_shape_variance = sum(b_w_shape_model.variances);
0258     total_intensity_variance = sum(b_w_intensity_model.variances);
0259     if(total_shape_variance ~= 0)
0260       factor = sqrt(total_intensity_variance/total_shape_variance);
0261     end
0262     models(this_set, error_ctr, factor_ctr).model = build_model_from_models(b_w_shape_model, b_w_intensity_model, ['variance: ', num2str(factor)],'constant', factor);
0263     if(find_gen_and_spec)
0264       [models(this_set, error_ctr, factor_ctr).mean_specificity, models(this_set, error_ctr, factor_ctr).std_specificity] = find_specificity(models(this_set, error_ctr, factor_ctr).model, images, spec_iters, ref_points_vec);
0265       [models(this_set, error_ctr, factor_ctr).mean_generalisability, models(this_set, error_ctr, factor_ctr).std_generalisability] = find_generalisability(models(this_set, error_ctr, factor_ctr).model, gen_iters);
0266     end
0267     range_score(this_set, error_ctr, factor_ctr) = measure_model(models(this_set, error_ctr, factor_ctr).model.variances, 50);
0268     factor_ctr = factor_ctr+1;
0269 
0270     % Smith: edges
0271     % RSS: Not sure what this is yet...
0272     edges = abs(diff(badly_warped_images,1,2));
0273     factor = sum(edges(:))/size(edges(:),1);
0274     models(this_set, error_ctr, factor_ctr).model = build_model_from_models(b_w_shape_model, b_w_intensity_model, ['edges: ', num2str(factor)],'constant', factor);
0275     range_score(this_set, error_ctr, factor_ctr) = measure_model(models(this_set, error_ctr, factor_ctr).model.variances, 50);
0276     if(find_gen_and_spec)
0277       [models(this_set, error_ctr, factor_ctr).mean_specificity, models(this_set, error_ctr, factor_ctr).std_specificity] = find_specificity(models(this_set, error_ctr, factor_ctr).model, images, spec_iters, ref_points_vec);
0278       [models(this_set, error_ctr, factor_ctr).mean_generalisability, models(this_set, error_ctr, factor_ctr).std_generalisability] = find_generalisability(models(this_set, error_ctr, factor_ctr).model, gen_iters);
0279     end
0280     factor_ctr = factor_ctr+1;
0281     
0282     shape_modes_ctr(error_ctr,this_set) = size(b_w_shape_model.pcs,2);
0283     intensity_modes_ctr(error_ctr,this_set) = size(b_w_intensity_model.pcs,2);
0284     error_ctr = error_ctr+1;
0285   end % error
0286   t = toc;
0287   disp(['Time for set ',num2str(this_set),': ',num2str(t)]);
0288 end % set
0289 
0290 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0291 % Plot graphs showing mean generalisability and specificity
0292 % against degree of misregistration
0293 
0294 if(find_gen_and_spec)
0295    for j=1:error_ctr-1
0296       error_labels(j).label = num2str(errors(j));
0297       for k=1:factor_ctr-1
0298         spec_labels(k).label = models(1,1,k).model.label;
0299         for stupid = 1:n_sets
0300           temp_spec_means(stupid,k) = models(stupid,j,k).mean_specificity;
0301           temp_spec_stds(stupid,k) = models(stupid,j,k).std_specificity;
0302           temp_gen_means(stupid,k) = models(stupid,j,k).mean_generalisability;
0303           temp_gen_stds(stupid,k) = models(stupid,j,k).std_generalisability;
0304         end
0305         spec_means(k,j) = mean(mean(temp_spec_means));
0306         spec_stds(k,j) = mean(mean(temp_spec_stds));
0307         gen_means(k,j) = mean(mean(temp_gen_means));
0308         gen_stds(k,j) = mean(mean(temp_gen_stds));
0309       end
0310       error_spec_means(j,:) = mean(temp_spec_means);
0311       error_gen_means(j,:) = mean(temp_gen_means);
0312    end
0313    spec_v_gen_fig = figure;, figure(spec_v_gen_fig), plot(error_spec_means', error_gen_means'), title('Specificity vs generalisability'),...
0314    xlabel('Specificity'),ylabel('Generalisability');
0315    cmap = colormap(jet(factor_ctr-1));
0316    for i=1:factor_ctr-1
0317      colour = cmap(i,:);
0318      hold on, plot(error_spec_means(:,i),error_gen_means(:,i),'Marker','o','MarkerFaceColor',colour,'LineStyle','none');
0319    end
0320     
0321    legend(error_labels(:).label,-1);
0322    cb = colorbar;
0323    for i=1:factor_ctr-1
0324       cb_ticks(i) = i+0.5;
0325       cb_labels{i} = models(1,1,i).model.label;
0326    end
0327    set(cb,'YTick',cb_ticks);
0328    set(cb,'YTickLabel',cb_labels);
0329     %spec_v_gen_fig = figure;, figure(spec_v_gen_fig), plot(mean(temp_spec_means), mean(temp_gen_means)), title(['Specificity vs generalisability, error ',num2str(j)]),...
0330     %  xlabel('Specificity'),ylabel('Generalisability');
0331     %legend(spec_labels(:).label,-1);
0332    spec_fig = figure;,figure(spec_fig),errorbar(log(spec_means)', spec_stds'./spec_means'), ...
0333    title('Specificity'),xlabel('Error'),ylabel('Log Specificity (msd of rand example from training example)');
0334    legend(spec_labels(:).label,-1);
0335    gen_fig = figure;,figure(gen_fig),errorbar(log(gen_means)', gen_stds'./gen_means'), ...
0336    title('Generalisability'),xlabel('Error'),ylabel('Log Generalisability (msd of example constructed from model to example');
0337    legend(spec_labels(:).label,-1);
0338 end
0339   
0340 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0341 % Plot graphs showing log of determinant of c params
0342 % against degree of misregistration
0343 
0344 if(blurred == 1)
0345   fig = figure, cmap = jet(256);, title('Blurred');
0346 else
0347   fig = figure, cmap = jet(256);, title('Unblurred');
0348 end
0349 line_ctr = 1;
0350 mode_fig = figure;,figure(mode_fig),errorbar([mean(shape_modes_ctr,2) mean(intensity_modes_ctr,2)], [std(shape_modes_ctr,1,2) std(intensity_modes_ctr,1,2)]),title('number of modes'),xlabel('Max error'),ylabel('number of modes');
0351 legend('Shape', 'Intensity');
0352 
0353 %for(factor = 0:factor_step:factor_max)
0354 mean_score = zeros(size(models,2), size(models,3));
0355 for(i = 1:size(models,3)) 
0356   mean_score(:,i) = mean(range_score(:,:,i))';
0357   stddev(:,i) = std(range_score(:,:,i))';
0358   legend_labels(line_ctr).name = models(1,1,line_ctr).model.label;
0359   line_ctr = line_ctr+1;
0360 end
0361 
0362 figure(fig),hold on,lines(line_ctr).line = errorbar(log(mean_score), stddev);
0363 for i=1:size(models,3)
0364   shape_weights(i) = models(1,1,i).model.shape_weight(1);
0365 end
0366 legend(legend_labels(:).name,-1);
0367 xlabel('Max error');
0368 ylabel('log of determinant of covariance matrix');
0369 
0370 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0371 % Plot graphs showing generalisability and specificity
0372 % against number of modes retained for each degree of misregistration
0373 
0374 if(find_gen_and_spec)
0375   for i=1:error_ctr-1
0376       for(j=1:factor_ctr-1)
0377          for(k=1:n_sets)
0378            num_modes(j,k) = size(models(k,i,j).model.pcs,2);
0379          end
0380       end
0381       most_modes = max(max(num_modes));
0382       mean_specs = zeros(factor_ctr-1,most_modes);
0383       std_specs = zeros(factor_ctr-1,most_modes);
0384       for j=1:factor_ctr-1
0385          [a,b] = show_spec_v_modes(models(:,i,j), images, spec_iters, ref_points_vec);
0386          mean_specs(j,1:size(a,2)) = a;
0387          std_specs(j,1:size(b,2)) = b;
0388       %  show_gen_v_modes(models, gen_iters);
0389       end
0390       figure,errorbar(mean(mean_specs),mean(std_specs)),xlabel('Number of modes'),ylabel('Specificity'), title(['Error ',num2str(i)]);
0391       clear mean_specs;
0392       clear std_specs;
0393   end
0394 end
eval_1d_app_model_obj_fn

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
