hdrvdp.m 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526
  1. function res = hdrvdp( test, reference, color_encoding, pixels_per_degree, options )
  2. % HDR-VDP-2 compute visually significant differences between an image pair.
  3. %
  4. % diff = HDRVDP( test, reference, color_encoding, pixels_per_degree, options )
  5. %
  6. % Parameters:
  7. % test - image to be tested (e.g. with distortions)
  8. % reference - reference image (e.g. without distortions)
  9. % color_encoding - color representation for both input images. See below.
  10. % pixels_per_degree - visual resolution of the image. See below.
  11. % options - cell array with { 'option', value } pairs. See the list of options
  12. % below. Note if unknown or misspelled option is passed, no warning
  13. % message is issued.
  14. %
  15. % The function returns a structure with the following fields:
  16. % P_map - probability of detection per pixel (matrix 0-1)
  17. % P_det - a single valued probability of detection (scalar 0-1)
  18. % C_map - threshold normalized contrast map, so that C_max=1
  19. % corresponds to the detection threshold (P_det=0.5).
  20. % C_max - maximum threshold normalized contrast, so that C_max=1
  21. % corresponds to the detection threshold (P_det=0.5).
  22. % Q - Quality correlate, which is 100 for the best quality and gets
  23. % lower for lower quality. Q can be negative in case of very large
  24. % differences.
  25. %
  26. % Test and references images are matrices of size (height, width,
  27. % channel_count). If there is only one channel, it is assumed to be
  28. % achromatic (luminance) with D65 spectra. The values must be given in
  29. % absolute luminance units (cd/m^2). If there are three color channels,
  30. % they are assumed to be RGB of an LCD display with red, green and blue LED
  31. % backlight. If different number of color channels is passed, their spectral
  32. % emission curve should be stored in a comma separated text file and its
  33. % name should be passed as 'spectral_emission' option.
  34. %
  35. % Note that the current version of HDR-VDP does NOT take color differences
  36. % into account. Spectral channels are used to properly compute luminance
  37. % sensed by rods and cones.
  38. %
  39. % COLOR ENCODING:
  40. %
  41. % HDR-VDP-2 requires that the color encoding is specified explicitly to avoid
  42. % mistakes when passing images. HDR-VDP operates on absolue physical units,
  43. % not pixel values that can be found in images. Therefore, it is necessary
  44. % to specify how the metric should interpret input images. The available
  45. % options are:
  46. %
  47. % 'luminance' - images contain absolute luminance values provided in
  48. % photopic cd/m^2. The images must contain exactly one color channel.
  49. %
  50. % 'luma-display' - images contain grayscale pixel values, sometimes known
  51. % as gamma-corrected luminance. The images must contain exactly one color
  52. % channel and the maximum pixel value must be 1 (not 256).
  53. % It corresponds to a gray-scale channel in
  54. % YCrCb color spaces used for video encoding. Because 'luma' alone does not
  55. % specify luminance, HDR-VDP-2 assumes the following display model:
  56. %
  57. % L = 99 * V^2.2 + 1,
  58. %
  59. % where L is luminance and V is luma. This corresponds to a display with
  60. % 2.2 gamma, 100 cd/m^2 maximum luminance and 1 cd/m^2 black level.
  61. %
  62. % 'sRGB-display' - use this color encoding for standard (LDR) color images.
  63. % This encoding assumes an sRGB display, with 100 cd/m^2 peak luminance and
  64. % 1 cd/m^2 black level. Note that this is different from sRGB->XYZ
  65. % transform, which assumes the peak luminance of 80 cd/m^2 and 1 cd/m^2
  66. % black level. The maximum pixel value must be 1 (not 256).
  67. %
  68. % 'rgb-bt.709' - use this color space for HDR images AFTER they are at
  69. % least roughly calibrated in absolute photometric units. The encoding
  70. % assumes the ITU-R BT.709 RGB color primaries (the same as for sRGB), but also
  71. % non-gamma corrected RGB color space.
  72. %
  73. % 'XYZ' - input image is provided as ABSOLUTE trichromatic color values in
  74. % CIE XYZ (1931) color space. The Y channel must be equal luminance in
  75. % cd/m^2.
  76. %
  77. % PIXELS_PER_DEGREE:
  78. %
  79. % This argument specifies the angular resolution of the image in terms of the
  80. % number of pixels per one visual degree. It
  81. % will change with a viewing distance and display resolution. A typical
  82. % value for a stanard resolution computer display is around 30. A
  83. % convenience function hdrvdp_pix_per_deg can be used to compute pixels per
  84. % degree parameter.
  85. %
  86. % OPTIONS:
  87. % The options must be given as name-value pairs in a cell array.
  88. % Default values are given in square brackets.
  89. %
  90. % 'surround_l' - [mean value of each color channel] luminance/intensity of the
  91. % surround, which is assumed to be uniform outside the input images.
  92. % The default value is 1e-5 (10^-5), which corresponds to almost
  93. % absolute darkness. Note that surround_l=0 should be avoided. It is
  94. % unrealistic to expect absolutely no light in physically plausible
  95. % environment. The special value -1 makes the surround surround_l
  96. % equal to the geometric mean of the image.
  97. % 'spectral_emission' - name of the comma separated file that contains
  98. % spectral emission curves of color channels of reference and test
  99. % images.
  100. % 'rgb_display' - [ccfl-lcd] use this option to specify one of the
  101. % predefined emission spectra for typical displays. Availaeble options
  102. % are:
  103. % crt - a typical CRT display
  104. % ccfl-lcd - an LCD display with CCFL backlight
  105. % led-lcd - an LCD display with LED backlight
  106. %
  107. % The following are the most important options used to fine-tune and calibrate
  108. % HDR-VDP:
  109. % 'peak_sensitivity' - absolute sensitivity of the HDR-VDP
  110. % 'mask_p' - excitation of the visual contrast masking model
  111. % 'mask_q' - inhibition of the visual contrast masking model
  112. %
  113. % EXAMPLE:
  114. % The following example creates a luminance ramp (gradient), distorts it
  115. % with a random noise and computes detection probabilities using HDR-VDP.
  116. %
  117. % reference = logspace( log10(0.1), log10(1000), 512 )' * ones( 1, 512 );
  118. % test = reference .* (1 + (rand( 512, 512 )-0.5)*0.1);
  119. % res = hdrvdp( test, reference, 'luminance', 30, { 'surround_l', 13 } );
  120. % clf;
  121. % imshow( hdrvdp_visualize( res.P_map, test ) );
  122. %
  123. % BUGS and LIMITATIONS:
  124. %
  125. % If you suspect the predictions are wrong, check first the Frequently
  126. % Asked Question section on the HDR-VDP-2 web-site
  127. % (http://hdrvdp.sourceforge.net). If it does not help, post your problem to
  128. % the HDR-VDP discussion group: http://groups.google.com/group/hdrvdp
  129. % (preffered) or send an e-mail directly to the author.
  130. %
  131. % REFERENCES
  132. %
  133. % The metric is described in the paper:
  134. % Mantiuk, Rafat, Kil Joong Kim, Allan G. Rempel, and Wolfgang Heidrich.
  135. % ?HDR-VDP-2: A Calibrated Visual Metric for Visibility and
  136. % Quality Predictions in All Luminance Conditions.?
  137. % ACM Trans. Graph (Proc. SIGGRAPH) 30, no. 4 (July 01, 2011): 1.
  138. % doi:10.1145/2010324.1964935.
  139. %
  140. % When refering to the metric, please cite the paper above and include the
  141. % version number, for example "HDR-VDP 2.2.0". To check the version number,
  142. % see the ChangeLog.txt. To check the version in the code, call
  143. % hdrvdp_version.txt.
  144. %
  145. % Copyright (c) 2011, Rafal Mantiuk <mantiuk@gmail.com>
  146. % Permission to use, copy, modify, and/or distribute this software for any
  147. % purpose with or without fee is hereby granted, provided that the above
  148. % copyright notice and this permission notice appear in all copies.
  149. %
  150. % THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  151. % WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  152. % MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  153. % ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  154. % WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  155. % ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  156. % OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  157. if( any( size(test) ~= size(reference) ) )
  158. error( 'reference and test images must be the same size' );
  159. end
  160. if( ~exist( 'options', 'var' ) )
  161. options = {};
  162. end
  163. if( ~exist( 'reconSpyr', 'file' ) )
  164. % If matlabPyrTools not in the path, add them now
  165. % Get the path to the hdrvdp directory
  166. pathstr = fileparts(mfilename( 'fullpath' ));
  167. addpath( fullfile( pathstr, 'matlabPyrTools_1.4_fixed' ) );
  168. % Re-check if everything went OK
  169. if( ~exist( 'reconSpyr', 'file' ) )
  170. error( 'Failed to add matlabPyrTools to the path.' );
  171. end
  172. end
  173. metric_par = hdrvdp_parse_options( options );
  174. % The parameters overwrite the options
  175. if( ~isempty( pixels_per_degree ) )
  176. metric_par.pix_per_deg = pixels_per_degree;
  177. end
  178. if( ~isempty( color_encoding ) )
  179. metric_par.color_encoding = color_encoding;
  180. end
  181. % Load spectral emission curves
  182. img_channels = size( test, 3 );
  183. switch lower( metric_par.color_encoding )
  184. case 'luminance'
  185. if( img_channels ~= 1 )
  186. error( 'Only one channel must be provided for "luminance" color encoding' );
  187. end
  188. check_if_values_plausible( test, metric_par );
  189. case 'luma-display'
  190. if( img_channels ~= 1 )
  191. error( 'Only one channel must be provided for "luma-display" color encoding' );
  192. end
  193. test = display_model( test, 2.2, 99, 1 );
  194. reference = display_model( reference, 2.2, 99, 1 );
  195. case 'srgb-display'
  196. if( img_channels ~= 3 )
  197. error( 'Exactly 3 color channels must be provided for "sRGB-display" color encoding' );
  198. end
  199. test = display_model_srgb( test );
  200. reference = display_model_srgb( reference );
  201. case 'rgb-bt.709'
  202. if( img_channels ~= 3 )
  203. error( 'Exactly 3 color channels must be provided for "rgb-bt.709" color encoding' );
  204. end
  205. check_if_values_plausible( test(:,:,2), metric_par );
  206. case 'xyz'
  207. if( img_channels ~= 3 )
  208. error( 'Exactly 3 color channels must be provided for "XYZ" color encoding' );
  209. end
  210. test = xyz2rgb( test );
  211. reference = xyz2rgb( reference );
  212. check_if_values_plausible( test(:,:,2), metric_par );
  213. case 'generic'
  214. if( isempty( metric_par.spectral_emission ) )
  215. error( '"spectral_emission" option must be specified when using the "generic" color encoding' );
  216. end
  217. end
  218. if( metric_par.surround_l == -1 )
  219. % If surround_l is set to -1, use the geometric (log) mean of each color
  220. % channel of the reference image
  221. metric_par.surround_l = geomean( reshape( reference, [size(reference,1)*size(reference,2) size(reference,3)] ) );
  222. elseif( length(metric_par.surround_l) == 1 )
  223. metric_par.surround_l = repmat( metric_par.surround_l, [1 img_channels] );
  224. end
  225. if( img_channels == 3 && ~isfield( metric_par, 'rgb_display' ) )
  226. metric_par.rgb_display = 'ccfl-lcd';
  227. end
  228. if( ~isempty( metric_par.spectral_emission ) )
  229. [tmp, IMG_E] = load_spectral_resp( [metric_par.base_dir '/' metric_par.spectral_emission] );
  230. elseif( isfield( metric_par, 'rgb_display' ) )
  231. [tmp, IMG_E] = load_spectral_resp( sprintf( 'emission_spectra_%s.csv', metric_par.rgb_display ) );
  232. elseif( img_channels == 1 )
  233. [tmp, IMG_E] = load_spectral_resp( 'd65.csv' );
  234. else
  235. error( '"spectral_emissiom" option needs to be specified' );
  236. end
  237. if( img_channels == 1 && size(IMG_E,2)>1 )
  238. % sum-up spectral responses of all channels for luminance-only data
  239. IMG_E = sum( IMG_E, 2 );
  240. end
  241. if( img_channels ~= size( IMG_E, 2 ) )
  242. error( 'Spectral response data is either missing or is specified for different number of color channels than the input image' );
  243. end
  244. metric_par.spectral_emission = IMG_E;
  245. % Compute spatially- and orientation-selective bands
  246. % Process reference first to reuse bb_padvalue
  247. [B_R L_adapt_reference band_freq bb_padvalue] = hdrvdp_visual_pathway( reference, 'reference', metric_par, -1 );
  248. [B_T L_adapt_test] = hdrvdp_visual_pathway( test, 'test', metric_par, bb_padvalue );
  249. L_adapt = (L_adapt_test + L_adapt_reference)./2;
  250. % precompute CSF
  251. csf_la = logspace( -5, 5, 256 );
  252. csf_log_la = log10( csf_la );
  253. b_count = length(B_T.sz);
  254. CSF = zeros( length(csf_la), b_count );
  255. for b=1:b_count
  256. CSF(:,b) = hdrvdp_ncsf( band_freq(b), csf_la', metric_par );
  257. end
  258. log_La = log10(clamp( L_adapt, csf_la(1), csf_la(end) ));
  259. D_bands = B_T;
  260. % Pixels that are actually different
  261. diff_mask = any((abs( test-reference ) ./ reference) > 0.001, 3);
  262. if( metric_par.do_quality_raw_data )
  263. qres = quality_pred_init();
  264. end
  265. Q = 0; % quality correlate
  266. % For each band
  267. for b = 1:b_count
  268. %masking params
  269. p = 10.^metric_par.mask_p;
  270. q = 10.^metric_par.mask_q;
  271. pf = (10.^metric_par.psych_func_slope)/p;
  272. % accumulate masking activity across orientations (cross-orientation
  273. % masking)
  274. mask_xo = zeros( get_band_size(B_T,b,1) );
  275. for o=1:B_T.sz(b)
  276. mask_xo = mask_xo + mutual_masking( b, o ); %min( abs(B_T{b,o}), abs(B_R{b,o}) );
  277. end
  278. log_La_rs = clamp( imresize(log_La,get_band_size(B_T,b,1)), csf_log_la(1), csf_log_la(end) );
  279. % per-pixel contrast sensitivity
  280. CSF_b = interp1( csf_log_la, CSF(:,b), log_La_rs );
  281. % REMOVED: Transform CSF linear sensitivity to the non-linear space of the
  282. % photoreceptor response
  283. % CSF_b = CSF_b .* 1./hdrvdp_joint_rod_cone_sens_diff( 10.^log_La_rs, metric_par );
  284. band_norm = 2^(b-1); % = 1/n_f from the paper
  285. band_mult = 1; %2^-(b-1);
  286. for o=1:B_T.sz(b)
  287. %excitation difference
  288. band_diff = (get_band(B_T,b,o) - get_band(B_R,b,o))*band_mult;
  289. if( metric_par.do_si_gauss )
  290. band_scale = size(band_diff,1)/size(test,1);
  291. % Sigma grows with lower frequencies to subtend a similar number of
  292. % cycles. Note that the scale differs between the bands.
  293. sigma = 10^metric_par.si_size / band_freq(b) * metric_par.pix_per_deg * band_scale;
  294. local_sum = fast_gauss( abs(band_diff), sigma, false );
  295. ex_diff = (sign_pow(band_diff/band_norm, p-1) * band_norm) .* local_sum;
  296. else
  297. ex_diff = sign_pow(band_diff/band_norm, p) * band_norm;
  298. end
  299. if( b == b_count )
  300. % base band
  301. N_nCSF = 1;
  302. else
  303. N_nCSF = (1./CSF_b);
  304. end
  305. if( metric_par.do_masking )
  306. k_mask_self = 10.^metric_par.mask_self; % self-masking
  307. k_mask_xo = 10.^metric_par.mask_xo; % masking across orientations
  308. k_mask_xn = 10.^metric_par.mask_xn; % masking across neighboring bands
  309. self_mask = mutual_masking( b, o );
  310. mask_xn = zeros( size( self_mask ) );
  311. if( b > 1 )
  312. mask_xn = max( imresize( mutual_masking( b-1, o ), size( self_mask ) ), 0 )/(band_norm/2);
  313. end
  314. if( b < (b_count-1) )
  315. mask_xn = mask_xn + max( imresize( mutual_masking( b+1, o ), size( self_mask ) ), 0 )/(band_norm*2);
  316. end
  317. % correct activity for this particular channel
  318. band_mask_xo = max( mask_xo - self_mask, 0 );
  319. N_mask = band_norm * (k_mask_self*(self_mask./N_nCSF/band_norm).^q + ...
  320. k_mask_xo*(band_mask_xo./N_nCSF/band_norm).^q + ...
  321. k_mask_xn*(mask_xn./N_nCSF).^q);
  322. D = ex_diff./sqrt( N_nCSF.^(2*p) + N_mask.^2 );
  323. % if( b == b_count )
  324. % D_bands = set_band( D_bands, b, o, D );
  325. % else
  326. D_bands = set_band( D_bands, b, o, sign_pow(D/band_norm, pf)*band_norm );
  327. % end
  328. else
  329. % NO masking
  330. D = ex_diff./N_nCSF.^p;
  331. D_bands = set_band( D_bands, b, o, sign_pow(D/band_norm, pf)*band_norm );
  332. end
  333. % Quality prediction
  334. % test
  335. w_f = interp1( metric_par.quality_band_freq, metric_par.quality_band_w, ...
  336. clamp( band_freq(b), metric_par.quality_band_freq(end), metric_par.quality_band_freq(1) ) );
  337. epsilon = 1e-12;
  338. % Mask the pixels that are almost identical in test and
  339. % reference images. Improves predictions for small localized
  340. % differences.
  341. diff_mask_b = imresize( double(diff_mask), size( D ) );
  342. D_p = D .* diff_mask_b;
  343. Q = Q + (log( msre( D_p ) + epsilon )-log(epsilon)) * w_f;
  344. if( metric_par.do_quality_raw_data )
  345. qres = quality_pred_band( qres, D_p, b, o );
  346. end
  347. end
  348. end
  349. S_map = abs(reconSpyr( D_bands.pyr, D_bands.pind, metric_par.steerpyr_filter ));
  350. % TODO: localized distortions may cause prediction of visibilble differences
  351. % in other parts of an image because they affect low frequencies. This is
  352. % especially apparent for super-threshold differences. A mechanism to
  353. % restrict location of such changes is needed.
  354. %
  355. %S_map = S_map .* double(diff_mask);
  356. if( metric_par.do_spatial_pooling )
  357. S_map = sum(S_map(:))/(max(S_map(:))+1e-12)*S_map;
  358. end
  359. res.P_map = 1 - exp( log(0.5)*S_map );
  360. res.P_det = max( res.P_map(:) );
  361. res.C_map = S_map;
  362. res.C_max = max( S_map(:) );
  363. res.Q = 100-Q;
  364. % This MOS fitting did not work very well, disabled for now
  365. %res.Q_MOS = 100 / (1+exp(metric_par.quality_logistic_q1*(Q+metric_par.quality_logistic_q2)));
  366. if( metric_par.do_quality_raw_data )
  367. res.qres = qres;
  368. end
  369. function m = mutual_masking( b, o )
  370. m = min( abs(get_band(B_T,b,o)), abs(get_band(B_R,b,o)) );
  371. % simplistic phase-uncertainty mechanism
  372. % TODO - improve it
  373. if( metric_par.do_si_gauss )
  374. m = blur_gaussian( m, 10^metric_par.si_size );
  375. else
  376. F = ones( 3, 3 );
  377. m = conv2( m, F/numel(F), 'same');
  378. end
  379. end
  380. end
  381. function y = sign_pow( x, e )
  382. y = sign(x) .* abs(x).^e;
  383. end
  384. function B = get_band( bands, b, o )
  385. oc = min( o, bands.sz(b) );
  386. B = pyrBand( bands.pyr, bands.pind, sum(bands.sz(1:(b-1)))+oc );
  387. end
  388. function sz = get_band_size( bands, b, o )
  389. sz = bands.pind(sum(bands.sz(1:(b-1)))+o,:);
  390. end
  391. function bands = set_band( bands, b, o, B )
  392. bands.pyr(pyrBandIndices(bands.pind,sum(bands.sz(1:(b-1)))+o)) = B;
  393. end
  394. function d = msre( X )
  395. d = sqrt(sum(X(:).^2))/numel(X);
  396. end
  397. function L = display_model( V, gamma, peak, black_level )
  398. L = peak * V.^gamma + black_level;
  399. end
  400. function RGB = display_model_srgb( sRGB )
  401. a = 0.055;
  402. thr = 0.04045;
  403. RGB = zeros(size(sRGB));
  404. RGB(sRGB<=thr) = sRGB(sRGB<=thr)/12.92;
  405. RGB(sRGB>thr) = ((sRGB(sRGB>thr)+a)/(1+a)).^2.4;
  406. RGB = 99*RGB + 1;
  407. end
  408. function RGB = xyz2rgb( XYZ )
  409. % Transform image color space using matrix M
  410. % dest = M * src
  411. M = [ 3.240708 -1.537259 -0.498570;
  412. -0.969257 1.875995 0.041555;
  413. 0.055636 -0.203996 1.057069 ];
  414. RGB = reshape( (M * reshape( XYZ, [size(XYZ,1)*size(XYZ,2) 3] )')', ...
  415. [size(XYZ,1) size(XYZ,2) 3] );
  416. end
  417. function Y = spatial_summation( X, sigma )
  418. % Essentilally a non-normalized Gaussian filter
  419. %
  420. ksize = round(sigma*6);
  421. h = fspecial( 'gaussian', ksize, sigma );
  422. h = h / max(h(:));
  423. Y = imfilter( X, h, 'replicate' );
  424. end
  425. function check_if_values_plausible( img, metric_par )
  426. % Check if the image is in plausible range and report a warning if not.
  427. % This is because the metric is often misused and used for with
  428. % non-absolute luminace data.
  429. if( ~metric_par.disable_lowvals_warning )
  430. if( max(img(:)) <= 1 )
  431. warning( 'hdrvdp:lowvals', [ 'The images contain very low physical values, below 1 cd/m^2. ' ...
  432. 'The passed values are most probably not scaled in absolute units, as requied for this color encoding. ' ...
  433. 'See ''doc hdrvdp'' for details. To disable this wanrning message, add option { ''disable_lowvals_warning'', ''true'' }.' ] );
  434. end
  435. end
  436. end