From caa749f94e580ec6042b42fabd455b9289ee1eed Mon Sep 17 00:00:00 2001 From: Wade Schwartzkopf Date: Tue, 7 Apr 2020 07:01:03 -0400 Subject: [PATCH] Add RGIQE --- Processing/RGIQE.m | 389 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 389 insertions(+) create mode 100644 Processing/RGIQE.m diff --git a/Processing/RGIQE.m b/Processing/RGIQE.m new file mode 100644 index 0000000..535d1b0 --- /dev/null +++ b/Processing/RGIQE.m @@ -0,0 +1,389 @@ +function [image_information_metric, rniirs, diagnostic] = RGIQE(filename, varargin) +% RGIQE Radar generalized image quality equation +% [image_information_metric, rniirs] = RGIQE(filename, 'PropertyName', PropertyValue, ...) +% +% An implementation of a generalized image quality equation to provide a +% quantitative image quality metric and an RNIIRS prediction. +% +% RNIIRS (Radar National Imagery Interpretation Rating Scale) is an image +% quality scale from 1 to 9. This "quality" refers to local image fidelity +% and a users ability to perform certain detection and recognition task +% with an image. The "quality" that RNIIRS attempts to capture does NOT +% refer to image size, geopositioning accuracy, radiometric accuracy, or +% any other host of characterists that may improve an image's utility. +% Although RNIIRS is subjective, generally assessed by a trained analyst, +% an image quality equation (IQE) may be used to predict RNIIRS for a given +% SAR collect. +% +% The General Image Quality Equation (GIQE) provided here uses generic, +% fundamental, sensor-indepenedent parameters of radar collection to +% measure local image fidelity by quantifying the information density in +% the ground in terms of bits\m^2, using Shannon's channel capacity +% theorem. Unlike RNIIRS, this value is quantitative and not subjective or +% based on analysts' assessments. We can however compute an RNIIRS +% prediction based off of it by fitting the information metric to analysts +% ratings. Both the information theory metric and the RNIIRS prediction +% are computed in this function. +% +% INPUTS: +% +% filename Must either be a complex image in a format +% recognized by the MATLAB SAR Toolbox or an XML +% file with SICD-formatted XML. If empty or not +% passed, the user will be prompted to select a +% file with a dialog box. +% +% Property name Description +% multi_noise Multiplicative noise. This value is equal to: +% +% MNR * sigma_0_background +% +% where MNR is the multiplicative noise ratio +% and sigma_0_background is the assumed power +% level of the background and ambiguous +% areas. +% MNR can further be broken down into: +% +% MNR = ISLR + QNR + AMBR +% +% where ISLR is the integrated sidelobe +% ratio, QNR is the quantization noise ratio, +% and AMBR is the ambiguity ratio. +% (Default = 0) +% interactive "AUTO", "NOISE", "SIGNAL", "BOTH", or "NONE". +% Enable interactive measurement of noise and +% signal levels. "AUTO" only runs +% interactive measurment if metadata from +% source file does not provide it. (Default +% is "AUTO"). +% frames If a file contains more than one complex +% dataset, this value specifies for which of +% those datasets to compute metrics. (Default +% is all datasets in the file.) +% signal_sigma_override Allows caller to override signal level used +% for computation (and INTERACTIVE input +% parameter). In general, this is not +% recommended, as RNIIRS estimation function +% was trained using the default values. If +% the dataset provided is radiometrically +% calibrated, this value will be treated as +% calibrated signal level in the ground plane +% (sigma-0). If the dataset is not +% calibrated, this value will be assume to be +% average pixel power. (Default is 1 square +% meter of radar return per square meter on +% the ground for co-pol, 0.25 for cross-pol.) +% +% Note that the selection of signal level and noise level are fundamentally +% different in this model. Noise is measureable, quantifiable, and likely +% fairy accuractely predictable (at least for additive noise). There is a +% right answer for what the noise level is for any given dataset. On the +% other hand, signal level is somewhat arbitrarily chosen, selected to be +% the power level of the scene content that one may be most interested in +% examining. It is for this reason that this function allows for a signal +% level override, but not a direct noise level override (although noise can +% be adjusted with the MULTI_NOISE input parameter.) +% +% Author: Wade Schwartzkopf, NGA/Research + +% Of the lines in this file, >95% is file and input parameter parsing and +% documention, and <5% actually computing the image quality metrics... + +%% Open input file +if ((nargin<1)||isempty(filename)) % If no filename was give, use dialog box to ask for one + % Recall last interactively selected path used + if ispref('matlab_sar_toolbox','last_used_directory') + pathname = getpref('matlab_sar_toolbox','last_used_directory'); + if ~ischar(pathname)||~exist(pathname,'dir') + pathname = pwd; + end + else + pathname = pwd; + end + [filename,pathname]=uigetfile(sar_file_extensions(),... + 'Open SAR Data File',pathname); + setpref('matlab_sar_toolbox','last_used_directory',pathname); +else % Path was already passed in with filenames + pathname=''; +end +if(filename) + fullfilename=[pathname filename]; +else % filename=0. Cancel was pressed, instead of a file being chosen. + return; +end +if ~isempty(pathname) + setpref('matlab_sar_toolbox','last_used_directory',pathname); %store path +end +try + ro = open_reader(fullfilename); + % Check for multiple images + if ~iscell(ro) + ro = {ro}; + end + meta = cell(numel(ro),1); + for i = 1:numel(ro) + meta{i} = ro{i}.get_meta(); + end +catch + try + fid = fopen(fullfilename); + meta = {sicdxml2struct(xmlread( java.io.StringBufferInputStream( char(fread(fid)'))))}; + fclose(fid); + catch + error('RGIQE:InvalidFile','Invalid File.'); + end +end + +%% Parse other input parameters +p = inputParser; +p.addParameter('multi_noise', 0, @(x) isnumeric(x) && isscalar(x)); +p.addParameter('interactive', 'AUTO', @(x) ismember(upper(x), ... + {'AUTO', 'NOISE', 'SIGNAL', 'BOTH', 'NONE'})); +p.addParameter('frames', 1:numel(meta), @(x) all(ismember(x, 1:numel(meta)))); +p.addParameter('signal_sigma_override', [], @(x) isnumeric(x) && isscalar(x)); +parse(p,varargin{:}) + +%% Iterate through datasets provided +meta = meta(p.Results.frames); +image_information_metric = zeros(numel(meta),1); +rniirs = zeros(numel(meta),1); +for i = 1:numel(meta) + meta{i} = derived_sicd_fields(meta{i}); % Assure all radiometric fields are derived if not there + + %% Get noise level + if isfield(meta{i}, 'Radiometric') && ... + isfield(meta{i}.Radiometric,'NoiseLevel') && ... + isfield(meta{i}.Radiometric.NoiseLevel,'NoiseLevelType') && ... + strcmpi(meta{i}.Radiometric.NoiseLevel.NoiseLevelType,'ABSOLUTE') && ... + isfield(meta{i}.Radiometric.NoiseLevel,'NoisePoly') && ... + ~ismember(upper(p.Results.interactive), {'NOISE', 'BOTH'}) + % SICD metadata can only tell us additive noise. Estimation of + % multiplicative noise must come from elsewhere. + noise_sigma = 10^(meta{i}.Radiometric.NoiseLevel.NoisePoly(1)/10) + p.Results.multi_noise; + elseif ismember(upper(p.Results.interactive), {'AUTO', 'NOISE', 'BOTH'}) + % Measured noise contains both additive and multiplicative noise + uiwait(msgbox('Please select a noise region for measurement.','Noise measurement')); + aoi = mitm_viewer(fullfilename,'mode','aoi','closeAfterSelect',true,'initialFrame',i); + if isempty(aoi), return; end % Cancel was selected + cdata = double(ro{i}.read_chip(... + [aoi(1) aoi(1)+aoi(3)-1],[aoi(2) aoi(2)+aoi(4)-1])); + noise_sigma = mean(abs(cdata(:)).^2); + else + noise_sigma = NaN; + end + if isfield(meta{i}, 'Radiometric') && ... + isfield(meta{i}.Radiometric, 'SigmaZeroSFPoly') + % Scale to be equal to noise equivalent sigma zero + noise_sigma = meta{i}.Radiometric.SigmaZeroSFPoly * noise_sigma; + end + + %% Get signal level + % One issue with using the Shannon-Hartley channel capacity formulation + % is that it requires a signal level. Certainly the amount of + % information a SAR image conveys is related to its signal level, so + % this is how we hope it would behave. A target with brighter returns + % is always easier to interpret than a target that is in the noise. + % However, it presents a problem to image quality prediction, since the + % signal level can't be known prior to measurement. This constraint is + % implicitly true for traditional RNIIRS assessments as well, which + % assumes specific types of targets. + % One approach for setting the signal level used in the channel + % capacity equation would be to assume a fixed radiometric signal level + % that is representative of the sorts of scene content we are generally + % interested in and use that same signal level across all images for + % which we compute this metric. This allows for an even + % apples-to-apples comparison across images of varying content and is + % consistent with the traditional RNIIRS viewpoint, which assumes + % specific types of tasks and scene content. + % Another approach would be to use the signal level in the actual scene + % by measuring the average pixel power across a full image or AOI. This + % approach may be more accurate at measuring the true information in + % any given image area. However, it may be less consistent with the + % traditional RNIIRS viewpoint, and this approach can only be used + % after the collection has been made, not prior to it for planning + % collection and predicting image quality. + % We also allow the caller of the function to pass their own signal + % level. Passing a non-default signal level is not recommended for + % RNIIRS estimation, since the RNIIRS function was trained on the + % default values. However, we provide this flexibility in case it + % might be educational to see how it affects the information theory + % metric. + if ~ismember('signal_sigma_override', p.UsingDefaults) + signal_sigma = p.Results.signal_sigma_override; + elseif isfield(meta{i}, 'Radiometric') && ... + isfield(meta{i}.Radiometric, 'SigmaZeroSFPoly') && ... + ~ismember(upper(p.Results.interactive), {'SIGNAL', 'BOTH'}) + % Default value: 1 square meter of return per square meter on + % the ground (roughly similar to a large vehicle over its area). + % Ulaby plots show cross-pol generally 5-10dB lower than + % co-pol. We pick 1/4 as simple round number. + if isfield(meta{i},'ImageFormation') && ... + isfield(meta{i}.ImageFormation,'TxRcvPolarizationProc') && ... + numel(unique(split(meta{i}.ImageFormation.TxRcvPolarizationProc,':')))>1 + signal_sigma = 0.25; + else + signal_sigma = 1; + end + elseif ismember(upper(p.Results.interactive), {'AUTO', 'SIGNAL', 'BOTH'}) + ButtonName = questdlg(['Do you want to use average signal '... + 'power across entire image or select an AOI?'], ... + 'Signal level estimation', ... + 'Image', 'AOI', 'Image'); + if isempty(ButtonName), return; end %window was closed CANCEL + if strcmpi(ButtonName,'AOI') + aoi = mitm_viewer(fullfilename,'mode','aoi','closeAfterSelect',true,'initialFrame',i); + if isempty(aoi), return; end % Cancel was selected + cdata = double(ro{i}.read_chip(... + [aoi(1) aoi(1)+aoi(3)-1],[aoi(2) aoi(2)+aoi(4)-1])); + signal_sigma = mean(abs(cdata(:)).^2) - noise_sigma; + else + samplesize=[1000 1000]; % Exract 1000x1000 array of samples to estimate mean + subsample=ceil(double([meta{i}.ImageData.NumCols meta{i}.ImageData.NumRows])./samplesize); + data = abs(single(ro{i}.read_chip([1 meta{i}.ImageData.NumCols],... + [1 meta{i}.ImageData.NumRows],subsample))); + signal_sigma = mean(abs(data(:)).^2) - noise_sigma; + end + if isfield(meta{i}, 'Radiometric') && ... + isfield(meta{i}.Radiometric, 'SigmaZeroSFPoly') + % Scale to be equal to noise equivalent sigma zero + signal_sigma = meta{i}.Radiometric.SigmaZeroSFPoly * signal_sigma; + end + warning('RGIQE:SignalLevelMeasured',['More uncertainty with ' ... + 'RNIIRS estimates exists when signal level is sampled, ' ... + 'rather than when using a radiometric constant.']); + else + signal_sigma = NaN; + end + + %% Save some diagnostic statistics + if isfield(meta{i}, 'Radiometric') && ... + isfield(meta{i}.Radiometric, 'SigmaZeroSFPoly') + diagnostic.nesz(i) = 10*log10(noise_sigma); + end + diagnostic.snr(i) = signal_sigma/noise_sigma; + diagnostic.sp_resolution(i) = sqrt(meta{i}.Grid.Row.ImpRespWid * ... + meta{i}.Grid.Col.ImpRespWid); + diagnostic.ellipicity(i) = max(meta{i}.Grid.Row.ImpRespWid, ... + meta{i}.Grid.Col.ImpRespWid)/min(meta{i}.Grid.Row.ImpRespWid, ... + meta{i}.Grid.Col.ImpRespWid); + if numel(meta)>1 + diagnostic.meta{i} = meta{i}; + else + diagnostic.meta = meta{i}; + end + + %% Compute bandwidth area + % The cosine of the slope angle is proper scale factor to project the + % bandwidth area into the ground plane. + diagnostic.bandwidth_area(i) = ... + meta{i}.Grid.Col.ImpRespBW * meta{i}.Grid.Row.ImpRespBW * ... + cosd(meta{i}.SCPCOA.SlopeAng); + % Computing area from vertices is much more work, but it may help + % visualize. + % Project slant plane spatial frequency bounds to ground plane + % coords_slant_x = meta{i}.Grid.Col.KCtr + ... + % meta{i}.Grid.Col.ImpRespBW * [-1 1 1 -1 -1]/2; + % coords_slant_y = meta{i}.Grid.Row.KCtr + ... + % meta{i}.Grid.Row.ImpRespBW * [1 1 -1 -1 1]/2; + % ruvect = [meta{i}.Grid.Row.UVectECF.X meta{i}.Grid.Row.UVectECF.Y meta{i}.Grid.Row.UVectECF.Z]; + % cuvect = [meta{i}.Grid.Col.UVectECF.X meta{i}.Grid.Col.UVectECF.Y meta{i}.Grid.Col.UVectECF.Z]; + % coords_slant_3d = ruvect'*coords_slant_y + cuvect'*coords_slant_x; + % Compute basis vectors for ground plane + % gpn = wgs_84_norm([meta{i}.GeoData.SCP.ECF.X meta{i}.GeoData.SCP.ECF.Y meta{i}.GeoData.SCP.ECF.Z])'; + % gruvect = project_ground(ruvect, gpn); % Project range vector to ground plane + % gruvect = gruvect/norm(gruvect); + % gcuvect = cross(gruvect,gpn); + % gcuvect = gcuvect/norm(gcuvect); + % Project spatial frequencies to ground plane + % coords_ground_3d = coords_slant_3d; + % coords_ground_x = zeros(size(coords_slant_x)); + % coords_ground_y = zeros(size(coords_slant_y)); + % for j=1:numel(coords_slant_x) + % % Planes in ECF + % coords_ground_3d(:,j) = project_ground(coords_slant_3d(:,j)',gpn)'; + % % Convert to 2D basis vectors in ground plane + % coords_ground_x(j) = dot(coords_ground_3d(:,j),gcuvect); + % coords_ground_y(j) = dot(coords_ground_3d(:,j),gruvect); + % end + % This is a way to compute the area of a any polygon from its vertices. + % diagnostic.bandwidth_area = sum((coords_ground_x(1:(end-1)).*coords_ground_y(2:end)) - ... + % (coords_ground_y(1:(end-1)).*coords_ground_x(2:end)))/2; + % Visualize spatial frequency projections in slant and ground planes + % figure; + % title('Spatial frequency bounds 3D'); + % g_plane = gruvect'*[100 100 -100 -100 100] + gcuvect'*[-100 100 100 -100 -100]; + % plot3(coords_slant_3d(1,:),coords_slant_3d(2,:),coords_slant_3d(3,:),... + % coords_ground_3d(1,:),coords_ground_3d(2,:),coords_ground_3d(3,:),... + % g_plane(1,:),g_plane(2,:),g_plane(3,:)); + % legend('Slant plane bounds', 'Ground plane bounds', 'Ground plane') + % figure; + % plot(coords_slant_x, coords_slant_y, coords_ground_x, coords_ground_y); + % title('Spatial frequency bounds'); + % legend('Slant plane bounds', 'Ground plane bounds') + + %% After 300+ line of setup code and comments, we finally compute quality metric... + % Shannon-Hartley theorem of channel capacity + % C = B * log2(1 + (S/N)) + % We propose a 2-dimension information rate based on a 2-dimensional + % bandwidth. This 2D bandwidth is an area in a space where each axis + % is cycles/m, as opposed to 1D distance in Hz (cycles per second). + % Therefore the units of this metric are bits/m^2. + image_information_metric(i) = diagnostic.bandwidth_area(i) * log2(1 + diagnostic.snr(i)); + + %% Map bits/m^2 into RNIIRS + if ismember('signal_sigma_override', p.UsingDefaults) + % The coefficients here were derived by fitting hundreds of SAR + % datasets with varying parameters to analysts ratings. This + % author prefers the information metric but also provides this as + % well for comparison to this legacy scale. + coeffs = [.3960 3.7555]; % Probably more digits of precision than we really have + rniirs(i) = polyval(coeffs,log2(image_information_metric(i))); + % For completeness, we handle the very low RNIIRS case where the + % logarithmic equation would have been negative. Negative is not + % allowed in the RNIIRS scale. We compute where a line tangent to + % the bits/m^2-to-RNIIRS function intersects with the origin. We + % will use this tangent line as the mapping for very low RNIIRS, + % since it approaches zero as information approaches zero and + % remains non-negative. This way the entire mapping function + % between bits/m^2 and RNIIRS is continuous and has a continuous + % derivative. The function will be logarithmic with respect to the + % information metric above this point (presumably nearly all data + % will fit in here) and linear with respect to the metric below + % this point (extremely low RNIIRS). Note that since NGA has no + % precedent or analysts ratings for RNIIRS<1, we are really free to + % define this however we want. + linlog_transition = exp(1 - coeffs(2)*log(2)/coeffs(1)); + % Derivative of linear portion (as well as intersection with log portion) + d_rniirs = coeffs(1)/(log(2)*linlog_transition); + % This is what the two lines would look line at very low RNIIRS: + % x=linspace(0,.02,1000); + % figure; plot(x, polyval(coeffs,log2(x)), x, x*d_rniirs); + % That was a lot of explanation for a case that will likely rarely + % be used, wasn't it? + if image_information_metric(i)