convert_img_units

PURPOSE ^

CONVERT_IMG_UNITS change image data units

SYNOPSIS ^

function y = convert_img_units(img,arg1,arg2)

DESCRIPTION ^

CONVERT_IMG_UNITS change image data units 
  img = convert_img_units(img,new_unit) converts img.elem_data or img.node_data
  expressed in img.current_params into different units for the same 
  physical property
 
 Examples: 
  img = mk_image(mdl, 2, 'resisitivity');
  img = convert_img_units(img, 'conductivity');

  img = mk_image(mdl, 1);
  img.current_params = 'log_resistivity';
  img = convert_img_units(img, 'conductivity');

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function y = convert_img_units(img,arg1,arg2)
0002 %CONVERT_IMG_UNITS change image data units
0003 %  img = convert_img_units(img,new_unit) converts img.elem_data or img.node_data
0004 %  expressed in img.current_params into different units for the same
0005 %  physical property
0006 %
0007 % Examples:
0008 %  img = mk_image(mdl, 2, 'resisitivity');
0009 %  img = convert_img_units(img, 'conductivity');
0010 %
0011 %  img = mk_image(mdl, 1);
0012 %  img.current_params = 'log_resistivity';
0013 %  img = convert_img_units(img, 'conductivity');
0014 
0015 % (C) 2012-2014 Alistair Boyle and Bartlomiej Grychtol
0016 % License: GPL version 2 or 3
0017 % $Id: convert_img_units.m 6926 2024-05-31 15:34:13Z bgrychtol $
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 % convert_img_units(x, in, out);
0027 if isnumeric(img)
0028   % require 3 arguments
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 % convert_img_units(img, out);
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   % what to do about old parametrization??
0067 %   img = rmfield(img, 'cur_unit');
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 % standard field order
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 %nothing to do
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) % Octave error off str2func
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     %no direct function for conversion, see if we can help
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 % test functions
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 %     pass = test_map_data(d, elem_types{i}, elem_types{j}, expected(i,j), pass);
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

Generated on Sun 29-Dec-2024 11:41:59 by m2html © 2005