0001 function y = convert_img_units(img,arg1,arg2)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 if ischar(img) && strcmp(img,'UNIT_TEST'); y = do_unit_test; return; end
0020
0021 if nargin < 2
0022 error('EIDORS:ConvertUnits:Input', ...
0023 'Not enough input arguments.');
0024 end
0025
0026
0027 if isnumeric(img)
0028
0029 if nargin < 3
0030 error('EIDORS:ConvertUnits:Input', ...
0031 'Two units are required to convert numeric values.');
0032 else
0033 x = img;
0034 cur_unit = arg1;
0035 new_unit = arg2;
0036 y = map_data(x, cur_unit, new_unit);
0037 return
0038 end
0039 end
0040
0041
0042 if any(isfield(img,{'elem_data','node_data'}))
0043 if nargin == 3
0044 cur_unit = arg1;
0045 new_unit = arg2;
0046 else
0047 try
0048 cur_unit = img.current_params;
0049 new_unit = arg1;
0050 catch
0051 error(['Don''t know what to convert. '...
0052 'Need either 3 arguments or img.current_params.']);
0053 end
0054 end
0055 else
0056 if nargin ==3
0057 cur_unit = arg1;
0058 img.data_mapper = cur_unit;
0059 new_unit = arg2;
0060 else
0061 new_unit = arg1;
0062 end
0063 img = data_mapper(img);
0064 cur_unit = img.current_params;
0065
0066
0067
0068 end
0069
0070
0071 [x do_nodes] = get_data(img);
0072
0073 y = map_data(x, cur_unit, new_unit);
0074
0075
0076 if do_nodes
0077 img.node_data = y;
0078 else
0079 img.elem_data = y;
0080 end
0081 img.current_params = new_unit;
0082 y = img;
0083
0084
0085 img = eidors_obj('image',img);
0086
0087 end
0088
0089 function x = map_data(x, cur_unit, new_unit)
0090
0091 if isempty(cur_unit) || isempty(new_unit)
0092 error('Unit string must not be empty.');
0093 end
0094
0095 if strcmp(cur_unit, new_unit)
0096 return
0097 end
0098
0099 try
0100 f = str2func(sprintf('%s2%s',cur_unit,new_unit));
0101 x = feval(f,x);
0102 x = fix_and_test(x);
0103 catch err
0104 if isempty(err.identifier)
0105 if strcmp(err.message(end-(30:-1:0)), ...
0106 'no function and no method found')
0107 err.identifier='MATLAB:UndefinedFunction';
0108 end
0109 end
0110
0111 if strcmp(err.identifier,'MATLAB:UndefinedFunction') || ( ...
0112 (length(err.message)>=24) && ...
0113 strcmp(err.message(1:23),'invalid function handle') );
0114 cur_pre = regexp(cur_unit,'^(.*?)_','match');
0115 if ~isempty(cur_pre) && ismember(cur_pre{1}, {'log_', 'log10_'});
0116 l = length(cur_pre{1});
0117 x = reverse(cur_pre{1}, x);
0118 x = map_data(x, cur_unit(l+1:end), new_unit);
0119 x = fix_and_test(x);
0120 return
0121 end
0122 new_pre = regexp(new_unit,'^(.*?)_','match');
0123 if ~isempty(new_pre) && ismember(new_pre{1}, {'log_', 'log10_', 'abs_'});
0124 l = length(new_pre{1});
0125 x = map_data(x, cur_unit, new_unit(l+1:end));
0126 x = apply(new_pre{1}, x);
0127 x = fix_and_test(x);
0128 return
0129 end
0130
0131
0132 error('EIDORS:ConversionNotSupported', errorstr(cur_unit, new_unit));
0133 else
0134 rethrow(err);
0135 end
0136 end
0137
0138 end
0139
0140 function x = fix_and_test(x)
0141 x(x == +inf) = +realmax;
0142 x(x == -inf) = -realmax;
0143 err_if_inf_or_nan(x, 'map_data-post');
0144 end
0145
0146 function err_if_inf_or_nan(x, str)
0147 if any(isnan(x) | isinf(x))
0148 error(sprintf('bad %s (%d NaN, %d Inf of %d)', ...
0149 str, ...
0150 length(find(isnan(x))), ...
0151 length(find(isinf(x))), ...
0152 length(x)));
0153 end
0154
0155 end
0156
0157 function x = conductivity2resistivity(x)
0158 x = 1./x;
0159 end
0160
0161 function x = resistivity2conductivity(x)
0162 x = 1./x;
0163 end
0164
0165 function erstr = errorstr(cur_unit, new_unit)
0166 erstr = sprintf( ['Unit conversion from %s to %s is not supported.\n'...
0167 'Please write a function %s2%s.'], cur_unit, new_unit, cur_unit, new_unit);
0168 end
0169
0170 function x = reverse(str, x);
0171 switch str
0172 case 'log10_'
0173 if any(x >= log10(realmax)-eps) warning('loss of precision -> inf'); end
0174 x = 10.^x;
0175 case 'log_'
0176 if any(x >= log(realmax)-eps) warning('loss of precision -> inf'); end
0177 x = exp(x);
0178 otherwise
0179 error('Cannot reverse %s', str);
0180 end
0181 end
0182
0183
0184 function x = apply(str, x)
0185 switch str
0186 case 'log10_'
0187 if any(x <= 0 + eps) warning('loss of precision -> -inf'); end
0188 x = log10(x);
0189 case 'log_'
0190 if any(x <= 0 + eps) warning('loss of precision -> -inf'); end
0191 x = log(x);
0192 case 'abs_'
0193 x = abs(x);
0194 otherwise
0195 error('Cannot apply %s', str);
0196 end
0197 end
0198
0199 function [x, do_nodes] = get_data(img)
0200 do_nodes = 0;
0201 try
0202 x = img.elem_data;
0203 catch
0204 try
0205 x = img.node_data;
0206 do_nodes = 1;
0207 catch
0208 error('No elem_data or node_data found');
0209 end
0210 end
0211 end
0212
0213
0214 function pass = do_unit_test
0215 pass = 1;
0216 d = 1;
0217 while (d ~= 1) && (d ~= 0)
0218 d = rand(1);
0219 end
0220 disp('TEST: convert_img_units()');
0221 elem_types = {'conductivity', 'log_conductivity', 'log10_conductivity', ...
0222 'resistivity', 'log_resistivity', 'log10_resistivity'};
0223 expected = [d log(d) log10(d) 1./d log(1./d) log10(1./d); ...
0224 exp(d) d log10(exp(d)) 1./exp(d) log(1./exp(d)) log10(1./exp(d)); ...
0225 10.^d log(10.^d ) d 1./10.^d log(1./10.^d ) log10(1./10.^d ); ...
0226 1./d log(1./d ) log10(1./d) d log(d) log10(d); ...
0227 1./exp(d) log(1./exp(d)) log10(1./exp(d)) exp(d) d log10(exp(d)); ...
0228 1./10.^d log(1./10.^d) log10(1./10.^d) 10.^d log(10.^d) d ];
0229 for i = 1:length(elem_types)
0230 for j = 1:length(elem_types)
0231
0232 in = elem_types{i};
0233 out= elem_types{j};
0234 unit_test_cmp(sprintf('convert_img_units(%s -> %s)\n', in, out), ...
0235 convert_img_units(d, in, out), ...
0236 expected(i,j) );
0237 end
0238 end
0239 end
0240
0241
0242 function pass = test_map_data(data, in, out, expected, pass)
0243 fprintf('TEST: convert_img_units(%s -> %s)\n', in, out);
0244 if convert_img_units(data, in, out) ~= expected
0245 pass = 0;
0246 fprintf('TEST: FAIL convert_img_units(%s -> %s)\n', in, out);
0247 end
0248 end