unpack_reconst_matrix

 UNPACK_RECONST_MATRIX: unpack a compressed, stored reconstruction matrix
 Reconstruction matrices, especially on a circular body, have many
  symmetries (left-right, up-down, flip, and recoprocity = 16x).
 To save space, a matrix can be stored in this format, and then
  used to save space. This function is used to rebuild the matrix

 RM = unpack_reconst_matrix(packed_matrix, Nelec, Ngrid);
   Nelec = number of electrodes
   Ngrid = number of grid points

 options = {'circ'} or {'octagon'} => region shape

 Example:
  load GREIT_v10_Circ_Matrix.mat
  RM = unpack_reconst_matrix(GREIT_v10_Circ_Matrix, 16, 32);

0001 function RM = unpack_reconst_matrix(packed_matrix, Nelec, Ngrid, options);
0002 % UNPACK_RECONST_MATRIX: unpack a compressed, stored reconstruction matrix
0003 % Reconstruction matrices, especially on a circular body, have many
0004 %  symmetries (left-right, up-down, flip, and recoprocity = 16x).
0005 % To save space, a matrix can be stored in this format, and then
0006 %  used to save space. This function is used to rebuild the matrix
0007 %
0008 % RM = unpack_reconst_matrix(packed_matrix, Nelec, Ngrid);
0009 %   Nelec = number of electrodes
0010 %   Ngrid = number of grid points
0011 %
0012 % options = {'circ'} or {'octagon'} => region shape
0013 %
0014 % Example:
0015 %  load GREIT_v10_Circ_Matrix.mat
0016 %  RM = unpack_reconst_matrix(GREIT_v10_Circ_Matrix, 16, 32);
0017 
0018 
0019 % (C) 2008-2010 Andy Adler. License: GPL version 2 or version 3
0020 % $Id: unpack_reconst_matrix.m 3017 2012-06-06 13:56:06Z aadler $
0021 
0022 [RM,sel1,sel2] = unpack_matrix(packed_matrix, Nelec, Ngrid);
0023 if nargin>=4; for i=1:length(options)
0024    if strcmp(options{i},'octagon')
0025 size(RM)
0026 size(sel1)
0027 size(sel2)
0028       RM= RM(sel2,:);
0029    elseif strcmp(options{i},'flipud')
0030       RM = reshape(RM, Ngrid,Ngrid, []);
0031       RM = RM(:,Ngrid:-1:1,:);
0032       RM = reshape(RM, Ngrid*Ngrid, []);
0033    else
0034       error(['Huh? Option=(',options{i},') not understood']);
0035    end
0036 end; end
0037 
0038 function [RM,sel1,sel2] = unpack_matrix(PM, Nelec, Ngrid);
0039    elec_idx = 1:Nelec;
0040    grid_idx = 1:Ngrid; grid_idx = grid_idx - mean(grid_idx);
0041    Nelec2 = Nelec/2;
0042    Nelec34= Nelec*3/4;
0043 
0044    % Take a slice
0045    [x,y]= meshgrid( elec_idx, elec_idx);
0046    ss1 = (y-x)>1 & (y-x)<15;
0047    sel1 = abs(x-y)>1 & abs(x-y)<15;
0048    
0049    [x,y]= meshgrid( grid_idx, grid_idx);
0050    ss2 = abs(x-y)<25 & abs(x+y)<25 ...
0051        & x<0 & y<0 & x>=y ;
0052    sel2 = abs(x-y)<25 & abs(x+y)<25;
0053  
0054    % Build up
0055    BP  = zeros(Nelec^2, Ngrid^2);
0056    BP(ss1,ss2) = PM;
0057    BP  = reshape(BP, Nelec,Nelec,Ngrid,Ngrid);
0058 
0059    % Reciprocity
0060    BP  = BP + permute(BP, [2,1,3,4]);
0061 
0062    % FLIP LR
0063    el= Nelec:-1:1;
0064    BP= BP + BP(el,el,[Ngrid:-1:1],:);
0065    % FLIP UD
0066    el= [Nelec2:-1:1,Nelec:-1:Nelec2+1];
0067    BP= BP + BP(el,el,:,[Ngrid:-1:1]);
0068    % Transpose
0069    el= [Nelec34:-1:1,Nelec:-1:Nelec34+1];
0070    BP= BP + permute(BP(el,el,:,:), [1,2,4,3]);
0071 
0072    % Final UD flip to match radiological view (upward toward patient)
0073    % Here electrodes are connected CW starting from TDC
0074    BP= BP(:,:,:,[Ngrid:-1:1]);
0075 
0076    RM= reshape(BP, Nelec^2, [])';
0077    RM= RM(:,sel1);
0078 % This creates the diamond shape, but we want to leave shape choice later
0079 %  RM= RM(sel2,sel1);