show_current

 SHOW_CURRENT: show current or other quantity defined
  on nodes onto the image

 show_current( img, volt )
   img -> img object 
   volt-> voltage on nodes (if not specified, img is solved
      via fwd_solve, or the value on node_data is used)

 Without specified values, show_current will
   create one current vector for each element
 The points are specified as either
   img.fwd_model.mdl_slice_mapper.npx   - number of points in horizontal direction
   img.fwd_model.mdl_slice_mapper.npy   - number of points in vertical 
    or
   img.fwd_model.mdl_slice_mapper.x_pts - vector of points in horizontal direction
   img.fwd_model.mdl_slice_mapper.y_pts - vector of points in vertical

 For 3D models, the slice may be specified as (see mdl_slice_mapper for information)
   img.fwd_model.mdl_slice_mapper.level 
 
 If an output is specified, then no image is draws
 
 Examples:
   img = mk_image( mk_common_model('b2c2',8));
   show_current(img); 
 OR
   img = mk_image( mk_common_model('b2c2',8));
   img.fwd_model.mdl_slice_mapper.npx  = 64;
   img.fwd_model.mdl_slice_mapper.npy  = 32;
   show_current(img); 
 OR
   img = mk_image( mk_common_model('b2c2',8));
   q = show_current(img); 
   quiver(q.xp,q.yp, q.xc,q.yc, 2,'b','LineWidth',2);

0001 function quiv = show_current( img, vv )
0002 % SHOW_CURRENT: show current or other quantity defined
0003 %  on nodes onto the image
0004 %
0005 % show_current( img, volt )
0006 %   img -> img object
0007 %   volt-> voltage on nodes (if not specified, img is solved
0008 %      via fwd_solve, or the value on node_data is used)
0009 %
0010 % Without specified values, show_current will
0011 %   create one current vector for each element
0012 % The points are specified as either
0013 %   img.fwd_model.mdl_slice_mapper.npx   - number of points in horizontal direction
0014 %   img.fwd_model.mdl_slice_mapper.npy   - number of points in vertical
0015 %    or
0016 %   img.fwd_model.mdl_slice_mapper.x_pts - vector of points in horizontal direction
0017 %   img.fwd_model.mdl_slice_mapper.y_pts - vector of points in vertical
0018 %
0019 % For 3D models, the slice may be specified as (see mdl_slice_mapper for information)
0020 %   img.fwd_model.mdl_slice_mapper.level
0021 %
0022 % If an output is specified, then no image is draws
0023 %
0024 % Examples:
0025 %   img = mk_image( mk_common_model('b2c2',8));
0026 %   show_current(img);
0027 % OR
0028 %   img = mk_image( mk_common_model('b2c2',8));
0029 %   img.fwd_model.mdl_slice_mapper.npx  = 64;
0030 %   img.fwd_model.mdl_slice_mapper.npy  = 32;
0031 %   show_current(img);
0032 % OR
0033 %   img = mk_image( mk_common_model('b2c2',8));
0034 %   q = show_current(img);
0035 %   quiver(q.xp,q.yp, q.xc,q.yc, 2,'b','LineWidth',2);
0036 
0037 
0038 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0039 % $Id: show_current.html 2819 2011-09-07 16:43:11Z aadler $
0040 
0041 if isstr(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0042 
0043 dims = size(img.fwd_model.nodes,2);
0044 
0045 if nargin==1; % We need to calculate
0046    if isfield(img,'elem_data')
0047       img.fwd_solve.get_all_meas = 1;
0048       vh = fwd_solve(img);
0049       vv = vh.volt(:,1);
0050    elseif isfield(img,'node_data');
0051       vv = img.node_data(:,1);
0052       error('show_current: cannot interpolate conductivity onto elements (yet)');
0053    else
0054       error('show_current: one parameter provided, and cannot solve for node voltages');
0055    end
0056 end 
0057 
0058 Nel = size(img.fwd_model.elems,1);
0059 elemcur = zeros(Nel+1,dims);
0060 % Calc field as I = sigma E
0061 %V1 = V0 + Ex*x1 + Ey*y1   [ 1 x1 y1 ] [ V0 ]
0062 %V2 = V0 + Ex*x2 + Ey*y2 = [ 1 x2 y2 ]*[ Ex ]
0063 %V3 = V0 + Ex*x3 + Ey*y    [ 1 x3 y3 ] [ Ey ]
0064 oo = ones(dims+1,1);
0065 for i=1:Nel
0066   idx = img.fwd_model.elems(i,:);
0067   nod = img.fwd_model.nodes(idx,:);
0068   VE  = ([oo, nod])\vv(idx);
0069    
0070   elemcur(i+1,:) = img.elem_data(i)*VE(2:end)';
0071 end
0072 
0073 if isfield(img.fwd_model, 'mdl_slice_mapper');
0074    if dims == 2;
0075       img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0076    end
0077    elem_ptr = mdl_slice_mapper( img.fwd_model, 'elem' );
0078    szep = size(elem_ptr);
0079 
0080    [xp,yp] = grid_the_space( img.fwd_model);
0081 
0082    xc = reshape( elemcur(elem_ptr+1,1), szep);
0083    yc = reshape( elemcur(elem_ptr+1,2), szep);
0084 else 
0085    pp = interp_mesh( img.fwd_model);
0086    xp = pp(:,1);
0087    yp= pp(:,2);
0088    xc = elemcur(2:end,1);
0089    yc = elemcur(2:end,2);
0090 end
0091 
0092 quiv.xp = xp; 
0093 quiv.yp = yp; 
0094 quiv.xc = xc; 
0095 quiv.yc = yc; 
0096 quiver(xp,yp,xc,yc,2,'k');
0097 
0098 function  [x,y] = grid_the_space( fmdl);
0099 
0100   xspace = []; yspace = [];
0101   try 
0102      xspace = fmdl.mdl_slice_mapper.x_pts;
0103      yspace = fmdl.mdl_slice_mapper.y_pts;
0104   end
0105 
0106   if isempty(xspace)
0107      npx  = fmdl.mdl_slice_mapper.npx;
0108      npy  = fmdl.mdl_slice_mapper.npy;
0109 
0110      xmin = min(fmdl.nodes(:,1));    xmax = max(fmdl.nodes(:,1));
0111      xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0112 
0113      ymin = min(fmdl.nodes(:,2));    ymax = max(fmdl.nodes(:,2));
0114      ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0115 
0116      range= max([xrange, yrange]);
0117      xspace = linspace( xmean - range*0.5, xmean + range*0.5, npx );
0118      yspace = linspace( ymean + range*0.5, ymean - range*0.5, npy );
0119   end
0120 
0121   [x,y]=meshgrid( xspace, yspace );
0122 
0123 function do_unit_test
0124    fmdl.nodes = [0,0;0,1;1,0;1,1];
0125    fmdl.elems = [1,2,3;2,3,4];
0126    fmdl.electrode(1).nodes = [1,2]; fmdl.electrode(1).z_contact = 0.01;
0127    fmdl.electrode(2).nodes = [3,4]; fmdl.electrode(2).z_contact = 0.01;
0128    fmdl.gnd_node = 1;
0129    fmdl.stimulation(1).stim_pattern = [1;-1];
0130    fmdl.stimulation(1).meas_pattern = [1,-1];
0131    fmdl.solve = @aa_fwd_solve;
0132    fmdl.system_mat = @aa_calc_system_mat;
0133    fmdl.type = 'fwd_model'
0134    img = mk_image(fmdl,[1,1]); 
0135    img.fwd_solve.get_all_meas = 1;
0136 
0137    img.fwd_model.mdl_slice_mapper.npx = 6;
0138    img.fwd_model.mdl_slice_mapper.npy = 6;
0139    show_current(img);
0140 
0141    imdl= mk_common_model('d2d2c',8);
0142    img = calc_jacobian_bkgnd( imdl );
0143    img.fwd_model.mdl_slice_mapper.npx = 64;
0144    img.fwd_model.mdl_slice_mapper.npy = 64;
0145    show_current(img);
0146 
0147    img.fwd_model.mdl_slice_mapper.x_pts = linspace(0,1,62);
0148    img.fwd_model.mdl_slice_mapper.y_pts = linspace(0,1,56);
0149    img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0150 
0151    img.fwd_solve.get_all_meas = 1;
0152    vh = fwd_solve(img);
0153    show_current(img, vh.volt(:,1));