xy= [270.6959  335.0508  375.7458  400.3521  422.1192  391.8346  353.9787  298.1414 ...
     230.9473  182.6811  131.5758  104.1303   83.3096  108.8623  140.0933  187.4131;
      32.9571   51.2429  127.1286  192.0429  245.9857  303.5857  347.4714  370.3286 ...
     370.3286  369.4143  341.0714  292.6143  227.7000  178.3286  110.6714   46.6714]';
load ~/docs/eidors/htdocs/tutorial/netgen/extrusion/CT2.mat
clf; subplot(231);
%A= double(A); for i=1:3; I(:,:,i) = A'/max(A(:)); end; imagesc(I);
hold on;
plot(trunk(:,1),trunk(:,2),'m*-');
plot(lung(:,1),lung(:,2),'r*-');
hh=plot(xy(:,1),xy(:,2), 'b.'); set(hh,'MarkerSize',20);
hold off
axis equal; axis([25,450,0,400]); axis off; title('image and electrod posns');

trunk = trunk*.01;
lung  = lung*.01;
xy    = xy*.01;

centroid= [2.5174 2.1323]; % centroid = mean(trunk)
elecx = xy(:,1) - centroid(1);
elecy = xy(:,2) - centroid(2);
electh= atan2(elecy,elecx)*180/pi; % electh=electh(1:3);
%[fmdl, mat_idx] = ng_mk_extruded_model({1,trunk*[.01]*[1,0;0,1],[4,50]},[ 6,0,1.0],[0.1]);
%[fmdl, mat_idx] = ng_mk_extruded_model({2,trunk*[.01]*[1,0;0,1],[4,50]},[electh,1+0*electh],[0.1]);

 lung = flipud(lung(1:3:end,:));
 [fmdl, mat_idx] = ng_mk_extruded_model({2,{trunk, lung} ,[4,50],.1},[electh,1+0*electh],[0.1]);
[stim,meas_sel] = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current'}, 1);
fmdl.stimulation = stim;

flipy= diag([1,-1,1]);
fmdl.nodes = fmdl.nodes*flipy;
img=mk_image(fmdl,1);
img.elem_data(mat_idx{2})= 0.5;
jac_backgnd = img.elem_data;
subplot(232);
show_fem(img); view(0,70);
hold on; 
hh=plot3(xy(:,1),xy(:,2),1+xy*[0;0], 'b.'); set(hh,'MarkerSize',20);
hold off;
xlim( mean(xlim) + .55*diff(xlim)*[-1,1]);
ylim( mean(ylim) + .55*diff(ylim)*[-1,1]);
title('3D FEM model and elecs');

fmdlr = ng_mk_extruded_model({0,trunk,[4,50],.1},[0,0],[0.1]);
flipy= diag([1,-1]);
fmdlr.nodes = fmdlr.nodes*flipy;
subplot(233);
show_fem(fmdlr);
title('2D FEM model and elecs');

stim = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current'}, 1);
img.fwd_model.stimulation = stim;

vh= fwd_solve(img);

targ= mk_c2f_circ_mapping(fmdl, [2.8;-2.8;1;0.3]);
img.elem_data = img.elem_data + targ*.5;
subplot(234); show_fem(img); view(2);
title('Simulated change region');

vi = fwd_solve(img);
vi = add_noise( 3, vi, vh );
subplot(235); plot([vh.meas, 100*(vi.meas - vh.meas)]);
title('simulated data');

imdl = mk_common_model('c2c2',16);
subplot(235);
imdl.fwd_model.normalize_measurements= 1;
imr= inv_solve(imdl, vh, vi);
show_fem(imr); ax= [min(imr.fwd_model.nodes), max(imr.fwd_model.nodes)]; axis(ax([1,3,2,4]));
title('2D circ reconstruction (SNR=3)');

% COARSE FINE MAPPING
cmdl.mk_coarse_fine_mapping.f2c_offset = [0,0,1];
cmdl.mk_coarse_fine_mapping.f2c_project = speye(3); % Scaling not required
cmdl.mk_coarse_fine_mapping.z_depth = inf;
cmdl.mk_coarse_fine_mapping.z_depth = 0.2;
c2f= mk_coarse_fine_mapping( fmdl, fmdlr);


imdl.name = 'CT pig 3D model';
imdl.fwd_model = fmdl;
imdl.rec_model = fmdlr;
imdl.fwd_model.coarse2fine = c2f;


subplot(236);
imdl.fwd_model.normalize_measurements= 1;
imdl.hyperparameter.value = 1e-3;
imdl.jacobian_bkgnd.value = jac_backgnd;
imr= inv_solve(imdl, vh, vi);

% NOW DIG IT OUT OF THE CACHE
imds = imdl;
   J = calc_jacobian( mk_image( imdl ));
   RtR = calc_RtR_prior( imdl );
   W   = calc_meas_icov( imdl );
   hp  = calc_hyperparameter( imdl );
   RM= (J'*W*J +  hp^2*RtR)\J'*W;
imds.solve = @solve_use_matrix; 
imds.solve_use_matrix.RM  = RM;

imds = rmfield(imds,{'hyperparameter','RtR_prior','jacobian_bkgnd'});
imds.fwd_model = rmfield(imds.fwd_model,{'nodes','elems','boundary','coarse2fine'});
imds.fwd_model.meas_select = meas_sel;
save pig_mdl imds vi vh
imr= inv_solve(imds, vh, vi);
show_fem(imr); ax= [min(imr.fwd_model.nodes), max(imr.fwd_model.nodes)]; axis(ax([1,3,2,4]));
title('3D mdl reconstruction (SNR=3)');



