123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526 |
- function res = hdrvdp( test, reference, color_encoding, pixels_per_degree, options )
- % HDR-VDP-2 compute visually significant differences between an image pair.
- %
- % diff = HDRVDP( test, reference, color_encoding, pixels_per_degree, options )
- %
- % Parameters:
- % test - image to be tested (e.g. with distortions)
- % reference - reference image (e.g. without distortions)
- % color_encoding - color representation for both input images. See below.
- % pixels_per_degree - visual resolution of the image. See below.
- % options - cell array with { 'option', value } pairs. See the list of options
- % below. Note if unknown or misspelled option is passed, no warning
- % message is issued.
- %
- % The function returns a structure with the following fields:
- % P_map - probability of detection per pixel (matrix 0-1)
- % P_det - a single valued probability of detection (scalar 0-1)
- % C_map - threshold normalized contrast map, so that C_max=1
- % corresponds to the detection threshold (P_det=0.5).
- % C_max - maximum threshold normalized contrast, so that C_max=1
- % corresponds to the detection threshold (P_det=0.5).
- % Q - Quality correlate, which is 100 for the best quality and gets
- % lower for lower quality. Q can be negative in case of very large
- % differences.
- %
- % Test and references images are matrices of size (height, width,
- % channel_count). If there is only one channel, it is assumed to be
- % achromatic (luminance) with D65 spectra. The values must be given in
- % absolute luminance units (cd/m^2). If there are three color channels,
- % they are assumed to be RGB of an LCD display with red, green and blue LED
- % backlight. If different number of color channels is passed, their spectral
- % emission curve should be stored in a comma separated text file and its
- % name should be passed as 'spectral_emission' option.
- %
- % Note that the current version of HDR-VDP does NOT take color differences
- % into account. Spectral channels are used to properly compute luminance
- % sensed by rods and cones.
- %
- % COLOR ENCODING:
- %
- % HDR-VDP-2 requires that the color encoding is specified explicitly to avoid
- % mistakes when passing images. HDR-VDP operates on absolue physical units,
- % not pixel values that can be found in images. Therefore, it is necessary
- % to specify how the metric should interpret input images. The available
- % options are:
- %
- % 'luminance' - images contain absolute luminance values provided in
- % photopic cd/m^2. The images must contain exactly one color channel.
- %
- % 'luma-display' - images contain grayscale pixel values, sometimes known
- % as gamma-corrected luminance. The images must contain exactly one color
- % channel and the maximum pixel value must be 1 (not 256).
- % It corresponds to a gray-scale channel in
- % YCrCb color spaces used for video encoding. Because 'luma' alone does not
- % specify luminance, HDR-VDP-2 assumes the following display model:
- %
- % L = 99 * V^2.2 + 1,
- %
- % where L is luminance and V is luma. This corresponds to a display with
- % 2.2 gamma, 100 cd/m^2 maximum luminance and 1 cd/m^2 black level.
- %
- % 'sRGB-display' - use this color encoding for standard (LDR) color images.
- % This encoding assumes an sRGB display, with 100 cd/m^2 peak luminance and
- % 1 cd/m^2 black level. Note that this is different from sRGB->XYZ
- % transform, which assumes the peak luminance of 80 cd/m^2 and 1 cd/m^2
- % black level. The maximum pixel value must be 1 (not 256).
- %
- % 'rgb-bt.709' - use this color space for HDR images AFTER they are at
- % least roughly calibrated in absolute photometric units. The encoding
- % assumes the ITU-R BT.709 RGB color primaries (the same as for sRGB), but also
- % non-gamma corrected RGB color space.
- %
- % 'XYZ' - input image is provided as ABSOLUTE trichromatic color values in
- % CIE XYZ (1931) color space. The Y channel must be equal luminance in
- % cd/m^2.
- %
- % PIXELS_PER_DEGREE:
- %
- % This argument specifies the angular resolution of the image in terms of the
- % number of pixels per one visual degree. It
- % will change with a viewing distance and display resolution. A typical
- % value for a stanard resolution computer display is around 30. A
- % convenience function hdrvdp_pix_per_deg can be used to compute pixels per
- % degree parameter.
- %
- % OPTIONS:
- % The options must be given as name-value pairs in a cell array.
- % Default values are given in square brackets.
- %
- % 'surround_l' - [mean value of each color channel] luminance/intensity of the
- % surround, which is assumed to be uniform outside the input images.
- % The default value is 1e-5 (10^-5), which corresponds to almost
- % absolute darkness. Note that surround_l=0 should be avoided. It is
- % unrealistic to expect absolutely no light in physically plausible
- % environment. The special value -1 makes the surround surround_l
- % equal to the geometric mean of the image.
- % 'spectral_emission' - name of the comma separated file that contains
- % spectral emission curves of color channels of reference and test
- % images.
- % 'rgb_display' - [ccfl-lcd] use this option to specify one of the
- % predefined emission spectra for typical displays. Availaeble options
- % are:
- % crt - a typical CRT display
- % ccfl-lcd - an LCD display with CCFL backlight
- % led-lcd - an LCD display with LED backlight
- %
- % The following are the most important options used to fine-tune and calibrate
- % HDR-VDP:
- % 'peak_sensitivity' - absolute sensitivity of the HDR-VDP
- % 'mask_p' - excitation of the visual contrast masking model
- % 'mask_q' - inhibition of the visual contrast masking model
- %
- % EXAMPLE:
- % The following example creates a luminance ramp (gradient), distorts it
- % with a random noise and computes detection probabilities using HDR-VDP.
- %
- % reference = logspace( log10(0.1), log10(1000), 512 )' * ones( 1, 512 );
- % test = reference .* (1 + (rand( 512, 512 )-0.5)*0.1);
- % res = hdrvdp( test, reference, 'luminance', 30, { 'surround_l', 13 } );
- % clf;
- % imshow( hdrvdp_visualize( res.P_map, test ) );
- %
- % BUGS and LIMITATIONS:
- %
- % If you suspect the predictions are wrong, check first the Frequently
- % Asked Question section on the HDR-VDP-2 web-site
- % (http://hdrvdp.sourceforge.net). If it does not help, post your problem to
- % the HDR-VDP discussion group: http://groups.google.com/group/hdrvdp
- % (preffered) or send an e-mail directly to the author.
- %
- % REFERENCES
- %
- % The metric is described in the paper:
- % Mantiuk, Rafat, Kil Joong Kim, Allan G. Rempel, and Wolfgang Heidrich.
- % ?HDR-VDP-2: A Calibrated Visual Metric for Visibility and
- % Quality Predictions in All Luminance Conditions.?
- % ACM Trans. Graph (Proc. SIGGRAPH) 30, no. 4 (July 01, 2011): 1.
- % doi:10.1145/2010324.1964935.
- %
- % When refering to the metric, please cite the paper above and include the
- % version number, for example "HDR-VDP 2.2.0". To check the version number,
- % see the ChangeLog.txt. To check the version in the code, call
- % hdrvdp_version.txt.
- %
- % Copyright (c) 2011, Rafal Mantiuk <mantiuk@gmail.com>
- % Permission to use, copy, modify, and/or distribute this software for any
- % purpose with or without fee is hereby granted, provided that the above
- % copyright notice and this permission notice appear in all copies.
- %
- % THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
- % WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
- % MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
- % ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
- % WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
- % ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
- % OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
- if( any( size(test) ~= size(reference) ) )
- error( 'reference and test images must be the same size' );
- end
- if( ~exist( 'options', 'var' ) )
- options = {};
- end
- if( ~exist( 'reconSpyr', 'file' ) )
- % If matlabPyrTools not in the path, add them now
-
- % Get the path to the hdrvdp directory
- pathstr = fileparts(mfilename( 'fullpath' ));
-
- addpath( fullfile( pathstr, 'matlabPyrTools_1.4_fixed' ) );
-
- % Re-check if everything went OK
- if( ~exist( 'reconSpyr', 'file' ) )
- error( 'Failed to add matlabPyrTools to the path.' );
- end
- end
- metric_par = hdrvdp_parse_options( options );
- % The parameters overwrite the options
- if( ~isempty( pixels_per_degree ) )
- metric_par.pix_per_deg = pixels_per_degree;
- end
- if( ~isempty( color_encoding ) )
- metric_par.color_encoding = color_encoding;
- end
- % Load spectral emission curves
- img_channels = size( test, 3 );
- switch lower( metric_par.color_encoding )
- case 'luminance'
- if( img_channels ~= 1 )
- error( 'Only one channel must be provided for "luminance" color encoding' );
- end
- check_if_values_plausible( test, metric_par );
- case 'luma-display'
- if( img_channels ~= 1 )
- error( 'Only one channel must be provided for "luma-display" color encoding' );
- end
- test = display_model( test, 2.2, 99, 1 );
- reference = display_model( reference, 2.2, 99, 1 );
- case 'srgb-display'
- if( img_channels ~= 3 )
- error( 'Exactly 3 color channels must be provided for "sRGB-display" color encoding' );
- end
- test = display_model_srgb( test );
- reference = display_model_srgb( reference );
- case 'rgb-bt.709'
- if( img_channels ~= 3 )
- error( 'Exactly 3 color channels must be provided for "rgb-bt.709" color encoding' );
- end
- check_if_values_plausible( test(:,:,2), metric_par );
- case 'xyz'
- if( img_channels ~= 3 )
- error( 'Exactly 3 color channels must be provided for "XYZ" color encoding' );
- end
- test = xyz2rgb( test );
- reference = xyz2rgb( reference );
- check_if_values_plausible( test(:,:,2), metric_par );
- case 'generic'
- if( isempty( metric_par.spectral_emission ) )
- error( '"spectral_emission" option must be specified when using the "generic" color encoding' );
- end
- end
- if( metric_par.surround_l == -1 )
- % If surround_l is set to -1, use the geometric (log) mean of each color
- % channel of the reference image
- metric_par.surround_l = geomean( reshape( reference, [size(reference,1)*size(reference,2) size(reference,3)] ) );
- elseif( length(metric_par.surround_l) == 1 )
- metric_par.surround_l = repmat( metric_par.surround_l, [1 img_channels] );
- end
- if( img_channels == 3 && ~isfield( metric_par, 'rgb_display' ) )
- metric_par.rgb_display = 'ccfl-lcd';
- end
- if( ~isempty( metric_par.spectral_emission ) )
- [tmp, IMG_E] = load_spectral_resp( [metric_par.base_dir '/' metric_par.spectral_emission] );
- elseif( isfield( metric_par, 'rgb_display' ) )
- [tmp, IMG_E] = load_spectral_resp( sprintf( 'emission_spectra_%s.csv', metric_par.rgb_display ) );
- elseif( img_channels == 1 )
- [tmp, IMG_E] = load_spectral_resp( 'd65.csv' );
- else
- error( '"spectral_emissiom" option needs to be specified' );
- end
- if( img_channels == 1 && size(IMG_E,2)>1 )
- % sum-up spectral responses of all channels for luminance-only data
- IMG_E = sum( IMG_E, 2 );
- end
- if( img_channels ~= size( IMG_E, 2 ) )
- error( 'Spectral response data is either missing or is specified for different number of color channels than the input image' );
- end
- metric_par.spectral_emission = IMG_E;
- % Compute spatially- and orientation-selective bands
- % Process reference first to reuse bb_padvalue
- [B_R L_adapt_reference band_freq bb_padvalue] = hdrvdp_visual_pathway( reference, 'reference', metric_par, -1 );
- [B_T L_adapt_test] = hdrvdp_visual_pathway( test, 'test', metric_par, bb_padvalue );
- L_adapt = (L_adapt_test + L_adapt_reference)./2;
- % precompute CSF
- csf_la = logspace( -5, 5, 256 );
- csf_log_la = log10( csf_la );
- b_count = length(B_T.sz);
- CSF = zeros( length(csf_la), b_count );
- for b=1:b_count
- CSF(:,b) = hdrvdp_ncsf( band_freq(b), csf_la', metric_par );
- end
- log_La = log10(clamp( L_adapt, csf_la(1), csf_la(end) ));
- D_bands = B_T;
- % Pixels that are actually different
- diff_mask = any((abs( test-reference ) ./ reference) > 0.001, 3);
- if( metric_par.do_quality_raw_data )
- qres = quality_pred_init();
- end
- Q = 0; % quality correlate
- % For each band
- for b = 1:b_count
-
- %masking params
- p = 10.^metric_par.mask_p;
- q = 10.^metric_par.mask_q;
- pf = (10.^metric_par.psych_func_slope)/p;
-
- % accumulate masking activity across orientations (cross-orientation
- % masking)
- mask_xo = zeros( get_band_size(B_T,b,1) );
- for o=1:B_T.sz(b)
- mask_xo = mask_xo + mutual_masking( b, o ); %min( abs(B_T{b,o}), abs(B_R{b,o}) );
- end
-
- log_La_rs = clamp( imresize(log_La,get_band_size(B_T,b,1)), csf_log_la(1), csf_log_la(end) );
- % per-pixel contrast sensitivity
- CSF_b = interp1( csf_log_la, CSF(:,b), log_La_rs );
-
- % REMOVED: Transform CSF linear sensitivity to the non-linear space of the
- % photoreceptor response
- % CSF_b = CSF_b .* 1./hdrvdp_joint_rod_cone_sens_diff( 10.^log_La_rs, metric_par );
-
- band_norm = 2^(b-1); % = 1/n_f from the paper
- band_mult = 1; %2^-(b-1);
-
- for o=1:B_T.sz(b)
-
- %excitation difference
- band_diff = (get_band(B_T,b,o) - get_band(B_R,b,o))*band_mult;
-
- if( metric_par.do_si_gauss )
- band_scale = size(band_diff,1)/size(test,1);
- % Sigma grows with lower frequencies to subtend a similar number of
- % cycles. Note that the scale differs between the bands.
- sigma = 10^metric_par.si_size / band_freq(b) * metric_par.pix_per_deg * band_scale;
- local_sum = fast_gauss( abs(band_diff), sigma, false );
- ex_diff = (sign_pow(band_diff/band_norm, p-1) * band_norm) .* local_sum;
- else
- ex_diff = sign_pow(band_diff/band_norm, p) * band_norm;
- end
-
-
- if( b == b_count )
- % base band
- N_nCSF = 1;
- else
- N_nCSF = (1./CSF_b);
- end
-
- if( metric_par.do_masking )
-
- k_mask_self = 10.^metric_par.mask_self; % self-masking
- k_mask_xo = 10.^metric_par.mask_xo; % masking across orientations
- k_mask_xn = 10.^metric_par.mask_xn; % masking across neighboring bands
- self_mask = mutual_masking( b, o );
-
- mask_xn = zeros( size( self_mask ) );
- if( b > 1 )
- mask_xn = max( imresize( mutual_masking( b-1, o ), size( self_mask ) ), 0 )/(band_norm/2);
- end
- if( b < (b_count-1) )
- mask_xn = mask_xn + max( imresize( mutual_masking( b+1, o ), size( self_mask ) ), 0 )/(band_norm*2);
- end
-
- % correct activity for this particular channel
- band_mask_xo = max( mask_xo - self_mask, 0 );
-
- N_mask = band_norm * (k_mask_self*(self_mask./N_nCSF/band_norm).^q + ...
- k_mask_xo*(band_mask_xo./N_nCSF/band_norm).^q + ...
- k_mask_xn*(mask_xn./N_nCSF).^q);
-
-
- D = ex_diff./sqrt( N_nCSF.^(2*p) + N_mask.^2 );
- % if( b == b_count )
- % D_bands = set_band( D_bands, b, o, D );
- % else
- D_bands = set_band( D_bands, b, o, sign_pow(D/band_norm, pf)*band_norm );
- % end
- else
- % NO masking
- D = ex_diff./N_nCSF.^p;
- D_bands = set_band( D_bands, b, o, sign_pow(D/band_norm, pf)*band_norm );
- end
-
- % Quality prediction
- % test
- w_f = interp1( metric_par.quality_band_freq, metric_par.quality_band_w, ...
- clamp( band_freq(b), metric_par.quality_band_freq(end), metric_par.quality_band_freq(1) ) );
- epsilon = 1e-12;
-
- % Mask the pixels that are almost identical in test and
- % reference images. Improves predictions for small localized
- % differences.
- diff_mask_b = imresize( double(diff_mask), size( D ) );
- D_p = D .* diff_mask_b;
-
- Q = Q + (log( msre( D_p ) + epsilon )-log(epsilon)) * w_f;
-
- if( metric_par.do_quality_raw_data )
- qres = quality_pred_band( qres, D_p, b, o );
- end
-
- end
-
- end
- S_map = abs(reconSpyr( D_bands.pyr, D_bands.pind, metric_par.steerpyr_filter ));
- % TODO: localized distortions may cause prediction of visibilble differences
- % in other parts of an image because they affect low frequencies. This is
- % especially apparent for super-threshold differences. A mechanism to
- % restrict location of such changes is needed.
- %
- %S_map = S_map .* double(diff_mask);
- if( metric_par.do_spatial_pooling )
- S_map = sum(S_map(:))/(max(S_map(:))+1e-12)*S_map;
- end
- res.P_map = 1 - exp( log(0.5)*S_map );
- res.P_det = max( res.P_map(:) );
- res.C_map = S_map;
- res.C_max = max( S_map(:) );
- res.Q = 100-Q;
- % This MOS fitting did not work very well, disabled for now
- %res.Q_MOS = 100 / (1+exp(metric_par.quality_logistic_q1*(Q+metric_par.quality_logistic_q2)));
- if( metric_par.do_quality_raw_data )
- res.qres = qres;
- end
- function m = mutual_masking( b, o )
- m = min( abs(get_band(B_T,b,o)), abs(get_band(B_R,b,o)) );
- % simplistic phase-uncertainty mechanism
- % TODO - improve it
-
- if( metric_par.do_si_gauss )
- m = blur_gaussian( m, 10^metric_par.si_size );
-
- else
- F = ones( 3, 3 );
- m = conv2( m, F/numel(F), 'same');
- end
- end
- end
- function y = sign_pow( x, e )
- y = sign(x) .* abs(x).^e;
- end
- function B = get_band( bands, b, o )
- oc = min( o, bands.sz(b) );
- B = pyrBand( bands.pyr, bands.pind, sum(bands.sz(1:(b-1)))+oc );
- end
- function sz = get_band_size( bands, b, o )
- sz = bands.pind(sum(bands.sz(1:(b-1)))+o,:);
- end
- function bands = set_band( bands, b, o, B )
- bands.pyr(pyrBandIndices(bands.pind,sum(bands.sz(1:(b-1)))+o)) = B;
- end
- function d = msre( X )
- d = sqrt(sum(X(:).^2))/numel(X);
- end
- function L = display_model( V, gamma, peak, black_level )
- L = peak * V.^gamma + black_level;
- end
- function RGB = display_model_srgb( sRGB )
- a = 0.055;
- thr = 0.04045;
- RGB = zeros(size(sRGB));
- RGB(sRGB<=thr) = sRGB(sRGB<=thr)/12.92;
- RGB(sRGB>thr) = ((sRGB(sRGB>thr)+a)/(1+a)).^2.4;
- RGB = 99*RGB + 1;
- end
- function RGB = xyz2rgb( XYZ )
- % Transform image color space using matrix M
- % dest = M * src
- M = [ 3.240708 -1.537259 -0.498570;
- -0.969257 1.875995 0.041555;
- 0.055636 -0.203996 1.057069 ];
- RGB = reshape( (M * reshape( XYZ, [size(XYZ,1)*size(XYZ,2) 3] )')', ...
- [size(XYZ,1) size(XYZ,2) 3] );
- end
- function Y = spatial_summation( X, sigma )
- % Essentilally a non-normalized Gaussian filter
- %
- ksize = round(sigma*6);
- h = fspecial( 'gaussian', ksize, sigma );
- h = h / max(h(:));
- Y = imfilter( X, h, 'replicate' );
-
- end
- function check_if_values_plausible( img, metric_par )
- % Check if the image is in plausible range and report a warning if not.
- % This is because the metric is often misused and used for with
- % non-absolute luminace data.
- if( ~metric_par.disable_lowvals_warning )
- if( max(img(:)) <= 1 )
- warning( 'hdrvdp:lowvals', [ 'The images contain very low physical values, below 1 cd/m^2. ' ...
- 'The passed values are most probably not scaled in absolute units, as requied for this color encoding. ' ...
- 'See ''doc hdrvdp'' for details. To disable this wanrning message, add option { ''disable_lowvals_warning'', ''true'' }.' ] );
- end
- end
- end
|