0001 function eval_1d_app_model_obj_fn
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046 n_images = 100;
0047
0048 n_sets = 5;
0049
0050 image_width = 50;
0051
0052 bump_width = 0.2;
0053
0054 ref_shift = 0.2;
0055 max_shift = 0.2;
0056 step = 0.1;
0057
0058 los = zeros(n_images,floor((max_shift-ref_shift)/step));
0059 his = zeros(n_images,floor((max_shift-ref_shift)/step));
0060
0061 min_error = 0;
0062 max_error = 1;
0063 error_step = 0.1;
0064
0065
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
0074 spec_iters = 10;
0075
0076 gen_iters = 10;
0077
0078 build_unwarped_model = 0;
0079 build_perfect_warp_model = 0;
0080 find_gen_and_spec = 1;
0081
0082
0083
0084
0085
0086 for this_set = 1:n_sets
0087 tic;
0088 ['calculating set ', num2str(this_set)]
0089
0090
0091
0092 [imagelist images points his_no_norm los_no_norm] = make_1d_images(n_images, image_width, bump_width);
0093
0094 window = 9;
0095
0096
0097
0098 ref_hi_no_norm = his_no_norm(1);
0099
0100 ref_lo_no_norm = los_no_norm(1);
0101 ref_points_vec_no_norm = points(:,1);
0102
0103
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
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
0117
0118
0119
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
0124 intensity_model.sd = std(intensity_model.params);
0125 shape_model.sd = std(shape_model.params);
0126
0127 intensity_total_vars = intensity_model.total_var;
0128 shape_total_vars = shape_model.total_var;
0129
0130
0131 if (build_unwarped_model),
0132
0133 total_shape_variance = sum(shape_model.variances);
0134 total_intensity_variance = sum(intensity_model.variances);
0135
0136 combined_model.shape_weight = ones(1,size(shape_model.params,2));
0137
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
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
0146 [combined_model.mean_specificity, combined_model.std_specificity] = find_specificity(combined_model, images, spec_iters, ref_points_vec);
0147
0148 fig_title = 'No warp: ';
0149
0150
0151
0152
0153
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
0162
0163
0164 for i=1:n_images
0165
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
0168 warped_points_vec = warped_points(:,i);
0169
0170 if any(warped_points_vec > 1)
0171 disp('stop!');
0172
0173 end
0174 image_vec = images(:,i);
0175
0176
0177 warped_images(:,i) = interp1(ref_points_vec, image_vec, warped_points_vec, 'linear',0);
0178
0179 end
0180
0181 if(build_perfect_warp_model)
0182
0183
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
0189
0190
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
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
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
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
0212
0213
0214
0215
0216
0217
0218 error_ctr = 1;
0219 for error = min_error:error_step:max_error
0220 ['calculating error ', num2str(error)]
0221
0222 for im=1:n_images
0223 errors(error_ctr) = error;
0224
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
0229
0230 badly_warped_images(:,im) = interp1(ref_points_vec, image_vec, badly_warped_points_vec, 'linear', 0);
0231
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
0238
0239 factor_ctr = 1;
0240
0241 for i=5:10
0242 factor = 10^i * 1/(10^10);
0243
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
0246 range_score(this_set, error_ctr, factor_ctr) = measure_model(models(this_set, error_ctr, factor_ctr).model.variances, 50);
0247
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
0256
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
0272
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
0287 t = toc;
0288 disp(['Time for set ',num2str(this_set),': ',num2str(t)]);
0289 end
0290
0291
0292
0293
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
0331
0332
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
0343
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
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
0373
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
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