Description of sigmatome2

0001 function [Filter, stim_pattern]= sigmatome2_filter(test);
0002 % SIGMATOME2_FILTER:  Hardware filter and stim_patterns for Sigmatome II device
0003 %
0004 % Usage:
0005 %  [Filter, stim_pat]= sigmatome2_filter;
0006 %    Filter   is a 416x208 Hardware filter
0007 %    stim_pat is a EIDORS stim_pattern structure
0008 %
0009 % Example:
0010 %   imdl= mk_common_model('c2t2',16);
0011 %   [Filter,stim_pat] = sigmatome2_filter;
0012 %   imdl.fwd_model.stimulation = stim_pat;
0013 %   imdl.fwd_model.jacobian = @jacobian_filtered;
0014 %   imdl.fwd_model.jacobian_filtered.jacobian = @jacobian_adjoint;
0015 %   imdl.fwd_model.jacobian_filtered.filter   = Filter;
0016 %   imdl.meas_icov = speye( size(Filter,1) );
0017 %   imdl.fwd_model = rmfield(imdl.fwd_model, 'meas_select');
0018 %
0019 %   img= mk_image(imdl);
0020 %   vh= fwd_solve(img);  vh = Filter * vh.meas;
0021 %   img.elem_data(50)=1.1;
0022 %   vi= fwd_solve(img);  vi = Filter * vi.meas;
0023 %
0024 %   imdl.hyperparameter.value = 0.03;
0025 %   imgr= inv_solve(imdl, vh, vi);
0026 %   show_slices(imgr);
0027 
0028 % (C) 2015 Andy Adler. License: GPL version 2 or version 3
0029 %   Based on information from Robert Guardo and Herve Gagnon
0030 % $Id: sigmatome2_filter.m 6365 2022-05-05 13:33:14Z aadler $
0031 
0032 if nargin>=1 && strcmp(test,'UNIT_TEST'); do_unit_test; return; end
0033 
0034 Filter   = calc_sigmatome2_filter;
0035 protocol = ComputeAllConfigs([0 1 2 3], 16) + 1;
0036 stim_pattern = stim_meas_list( protocol );
0037 
0038 
0039 % Generate the filter function of the sigmatome system
0040 function filter = calc_sigmatome2_filter
0041    Filter3= [0,         0.002992185737274, 0.002885352918496,-0.001992209356618, ...
0042      0.000247868161165, 0.001320969614758,-0.002994506380400, 0.000217287704644, ...
0043      0.000820760030266,-0.002388560206626, 0.005467093174149,-0.000465034175608, ...
0044     -0.014049400948852, 0.011968526958771, 0.017269336581277,-0.033327865383464, ...
0045     -0.000687996077078, 0.065884070897957,-0.055509942269519,-0.084217637095161, ...
0046      0.414790173689917, 0.947607669916644, 0.723749489957398, 0.088544295521622, ...
0047     -0.118281839343952, 0.032044694797133, 0.043642225555493,-0.035476390241121, ...
0048     -0.007655680300629, 0.022891001199580,-0.005076656322158,-0.007794874431280, ...
0049      0.006724761261718,-0.000064773581282,-0.002313251865051, 0.001234620138549, ...
0050     -0.002848224928827,-0.000772726894395, 0.002559973165365,-0.001297698809235, ...
0051      0.000569993353282, 0.002828133181959, 0.000426156082041];
0052    filter = sparse(416,208);
0053    filter(1:length(Filter3),3) = Filter3';
0054    for i=1:208; 
0055       filter(:,i) = circshift(filter(:,3),+2*(i-3));
0056    end
0057 
0058 function [Configs Mesptr Curptr Index] = ComputeAllConfigs(InitConfig, N);
0059 %
0060 %ComputeAllConfigs: Cette fonction calcule toutes les configurations
0061 %                  [Source, Puit, Suiveur, Inverseur] a partir de la
0062 %                  configuration ititiale et du nombre d'electrodes.
0063 %
0064 %SYNTAXE:  Configs = ComputeAllConfigs(InitConfig, NbreElectrodes);
0065 %
0066 %INPUT:    InitConfig(1 x 4):  Configuration initiale de mesure exprimee dans
0067 %                              le format suivant : [Source, Puit, Suiveur,
0068 %                              Inverseur]. Les electrodes sont numerotees de
0069 %                              0 a (NbreElectrodes-1).
0070 %
0071 %          N:                  Nombre d'electrodes.
0072 %
0073 %OUTPUT:   Configs(m x 4):     Configurations d'electrodes : [Source, Puit,
0074 %                              Suiveur, Inverseur] pour les "m" mesures prises
0075 %                              par le scanhead.
0076 %
0077 %          Mesptr(2 x N):      Chaque colonne correspond aux positions
0078 %                              d'electrodes qui forment une des "N" paires
0079 %                              d'electrodes utilises pour effectuer les mesures
0080 %                              de tension.
0081 %
0082 %          Curptr(2 x N):      Chaque colonne correspond aux numeros des noeuds
0083 %                              qui forment une des "N" paires d'electrodes
0084 %                              utilises pour injecter les courants.
0085 %
0086 %          Index(m x 2):       Pour les "m" mesures prises par le scanhead,
0087 %                              doublet indiquant le numero de la paire
0088 %                              d'electrodes qui sert a applique le courant et
0089 %                              le numero de la paire qui sert a effectuer la
0090 %                              mesure.
0091 
0092 % Copyright (c) 2009 Hervé Gagnon, Ecole Polytechnique de Montréal.
0093 
0094    % Validation du parametre InitConfig
0095    if (size(InitConfig,1) ~= 1 || size(InitConfig,2) ~= 4)
0096        error('Le parametre "InitConfig" doit etre de dimension (1 x 4)!');
0097    end
0098 
0099    % Validation du parametre NbreElectrodes
0100    if (length(N) ~= 1)
0101        error('Le parametre "NbreElectrodes" doit etre un scalaire!');
0102    end
0103 
0104    SrcePosInit = InitConfig(1);
0105    SinkPosInit = InitConfig(2);
0106    FollPosInit = InitConfig(3);
0107    InvtPosInit = InitConfig(4);
0108 
0109    GapCourant = abs(SrcePosInit - SinkPosInit);
0110    GapTension = abs(FollPosInit - InvtPosInit);
0111 
0112    if ((GapCourant == N/2) && (GapTension == N/2))
0113        m = N*(N - 2);
0114    elseif((GapCourant == GapTension) || (GapCourant + GapTension == N))
0115        m = N*(N - 3);
0116    else
0117        m = N*(N - 4);
0118    end
0119       
0120    SrcePos = SrcePosInit;
0121    SinkPos = SinkPosInit;
0122    FollPos = FollPosInit;
0123    InvtPos = InvtPosInit;
0124 
0125    Configs = zeros(m, 4);
0126 
0127    for i = 1:m
0128        Configs(i,:) = [SrcePos SinkPos FollPos InvtPos];
0129        
0130        SrcePos = IncrementVarModuloN(SrcePos, N);
0131        SinkPos = IncrementVarModuloN(SinkPos, N);
0132        FollPos = IncrementVarModuloN(FollPos, N);
0133        InvtPos = IncrementVarModuloN(InvtPos, N);
0134        
0135        if ((SrcePos == SrcePosInit) && (SinkPos == SinkPosInit))
0136           FollPos = IncrementVarModuloN(FollPos, N);
0137           InvtPos = IncrementVarModuloN(InvtPos, N);
0138           while (FollPos == SrcePos || FollPos == SinkPos || InvtPos == SrcePos || InvtPos == SinkPos)
0139               FollPos = IncrementVarModuloN(FollPos, N);
0140               InvtPos = IncrementVarModuloN(InvtPos, N);
0141           end
0142        end
0143    end
0144 
0145    %%%%%%%%%%%%%%%%%%%%%%%%%%
0146    Curptr = zeros(2,N);
0147    Mesptr = zeros(2,N);
0148 
0149    Curptr(:,1) = [SrcePosInit; SinkPosInit];
0150    Mesptr(:,1) = [FollPosInit; InvtPosInit];
0151 
0152    for i = 2:N
0153        Curptr(1,i) = IncrementVarModuloN(Curptr(1,i-1), N);
0154        Curptr(2,i) = IncrementVarModuloN(Curptr(2,i-1), N);
0155        Mesptr(1,i) = IncrementVarModuloN(Mesptr(1,i-1), N);
0156        Mesptr(2,i) = IncrementVarModuloN(Mesptr(2,i-1), N);
0157    end
0158 
0159    Index = zeros(m, 2);
0160    for i = 1:m
0161        Index(i,1) = find((Configs(i,1) == Curptr(1,:)) & (Configs(i,2) == Curptr(2,:)));
0162        Index(i,2) = find((Configs(i,3) == Mesptr(1,:)) & (Configs(i,4) == Mesptr(2,:)));
0163    end
0164 
0165 function VarOut = IncrementVarModuloN(VarIn, N);
0166    VarOut = VarIn + 1;
0167    if (VarOut == N)
0168        VarOut = 0;
0169    end
0170 
0171 function do_unit_test
0172    imdl= mk_common_model('c2t2',16);
0173    [Filter,stim_pat] = sigmatome2_filter;
0174    imdl.fwd_model.stimulation = stim_pat;
0175    imdl.fwd_model.jacobian = @jacobian_filtered;
0176    imdl.fwd_model.jacobian_filtered.jacobian = @jacobian_adjoint;
0177    imdl.fwd_model.jacobian_filtered.filter   = Filter;
0178    imdl.meas_icov = speye( size(Filter,1) );
0179    imdl.fwd_model = rmfield(imdl.fwd_model, 'meas_select');
0180 
0181    img= mk_image(imdl);
0182    vh= fwd_solve(img);  vh = Filter * vh.meas;
0183    img.elem_data(50)=1.1;
0184    vi= fwd_solve(img);  vi = Filter * vi.meas;
0185 
0186    imdl.hyperparameter.value = 0.03;
0187    imgr= inv_solve(imdl, vh, vi);
0188    show_slices(imgr);
0189 
0190    max_imgr = max(imgr.elem_data(:));
0191    unit_test_cmp('Max', max_imgr, 6.683411990674756e-04, 1e-8);
0192    tst= [17;25;26;37;38;50;66;82];
0193    unit_test_cmp('Loc', find(imgr.elem_data>0.95*max_imgr), tst);
sigmatome2_filter

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
