Description of eval_1d_app_model_obj

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
