hdrvdp_visualize.m 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322
  1. function map = hdrvdp_visualize( p1, varargin )
  2. % HDRVDP_VISUALIZE produces color or grayscale visualization of the metric
  3. % predictions
  4. %
  5. % map = HDRVDP_VISUALIZE( P )
  6. % map = HDRVDP_VISUALIZE( P, context_image )
  7. % map = HDRVDP_VISUALIZE( P, context_image, target )
  8. % map = HDRVDP_VISUALIZE( P, context_image, target, colormap )
  9. % map = HDRVDP_VISUALIZE( 'pmap', P, options )
  10. % map = HDRVDP_VISUALIZE( 'diff', P, test, reference, options )
  11. %
  12. % Creates a colorful visualization from the HDR-VDP-2 results. Depending on
  13. % whether 'diff' or 'pmap' type is selected, the function visualizes a
  14. % different aspect of the prediction. The older syntax with the first
  15. % parameter 'P' is equivalent to 'pmap'.
  16. %
  17. % 'pmap' produces the probability of detection map. Such a map shows where
  18. % and how likely a difference will be noticed. However, It does not
  19. % show what this differece is.
  20. %
  21. % The required argument 'P' is the probability of detection map. P
  22. % must be within the range 0-1
  23. %
  24. % 'diff' shows the contrast-normalized per-pixel difference weighted by the
  25. % probability of detection. The resulting images do NOT show
  26. % probabilities. However, they better correspond to the perceived
  27. % differences and thus are easier to interpret.
  28. %
  29. % The required argument 'P' is the probability of detection map. P
  30. % must be within the range 0-1
  31. %
  32. % The two required argumants are the test and references images, the
  33. % same as passed to the hdrvdp.
  34. %
  35. % 'options' a cell array with { 'name', value } pairs. Refer to the examples
  36. % below. Available options:
  37. %
  38. % 'context_image' (default: []) - display context image on top of a
  39. % colour map. This is usually a 'test' image passed to the hdrvdp.
  40. %
  41. %
  42. % 'target' (default 'screen') - specifies desired output device and can
  43. % be one of:
  44. %
  45. % 'screen' - (default) - the map be shown on a color screen.
  46. % The map will contain good reproduction of the context
  47. % image 'img'.
  48. %
  49. % 'print' - the map can be printed on a gray-scale printer, so color
  50. % information will be lost. If this target is selected,
  51. % the color map will contain luma differences in addition
  52. % to color differences. To ensure that the context image
  53. % does not interfere with errors, only low-contrast and
  54. % high frequency content of the image will be preserved.
  55. %
  56. % 'colormap' (default 'trichromatic') parameter can be one of:
  57. %
  58. % 'trichromatic' - Errors are represented as multiple colors: blue,
  59. % cyan, green, yellow and red, which correspond to P equal
  60. % 0, 0.25, 0.5, 0.75 and 1.
  61. %
  62. % 'dichromatic' - more appropriate if observes may be color
  63. % deficient. The hue changes from cyan (0.0) to gray (0.5)
  64. % and then yellow (1.0). The look of images for
  65. % color-deficient observers can be tested at:
  66. % http://www.colblindor.com/coblis-color-blindness-simulator/
  67. %
  68. % 'monochromatic' - use only grayscale. Makes sense only with
  69. % target='print' or when no context image is specified
  70. %
  71. %
  72. % Tone-mapping is applied to the context image to reduce the dynamic range
  73. % so that highly saturated colors can be used also in bright image regions.
  74. % Tone-mapping will enhance contrast if it is lower than the available
  75. % contrast. This behavior is useful for images that have contrast that is
  76. % near detection threshold (e.g. ModelFest stimuli).
  77. %
  78. % The function returns a gamma-corrected sRGB image.
  79. %
  80. % Example:
  81. %
  82. % reference = logspace( log10(0.1), log10(1000), 512 )' * ones( 1, 512 );
  83. % test = reference .* (1 + (rand( 512, 512 )-0.5)*0.02);
  84. % res = hdrvdp( test, reference, 'luminance', 30 );
  85. %
  86. % imshow( hdrvdp_visualize( 'pmap', res.P_map, { 'context_image', test } ) )
  87. %
  88. % Legend for the color-scales can be found in the color_scales
  89. % directory.
  90. %
  91. % Copyright (c) 2011, Rafal Mantiuk <mantiuk@gmail.com>
  92. % Permission to use, copy, modify, and/or distribute this software for any
  93. % purpose with or without fee is hereby granted, provided that the above
  94. % copyright notice and this permission notice appear in all copies.
  95. %
  96. % THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  97. % WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  98. % MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  99. % ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  100. % WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  101. % ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  102. % OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  103. if( ischar( p1 ) ) % new syntax
  104. switch( p1 )
  105. case 'diff'
  106. [P, test_img, reference_img, options] = set_params( varargin, 'required', 'required', 'required', {} );
  107. opt = parse_options( options );
  108. P = norm_diff_img( test_img, reference_img ) .* P;
  109. case 'pmap'
  110. [P, options] = set_params( varargin, 'required', {} );
  111. opt = parse_options( options );
  112. otherwise
  113. error( 'Urecognizaed visuzalization type' );
  114. end
  115. else % old syntax
  116. P = p1;
  117. opt = struct();
  118. [opt.context_image, opt.target, opt.colormap] = set_params( varargin, [], 'screen', 'trichromatic' );
  119. end
  120. if( isempty(opt.context_image) )
  121. tmo_img = ones(size(P))*0.5;
  122. elseif( strcmp( opt.target, 'print' ) )
  123. l = log_luminance( opt.context_image );
  124. hp_img = (l - blur_gaussian( l, 2 ) + mean(l(:)));
  125. tmo_img = vis_tonemap(hp_img, 0.1) + 0.5;
  126. elseif( strcmp( opt.target, 'screen' ) )
  127. tmo_img = vis_tonemap( log_luminance(opt.context_image), 0.6 );
  128. else
  129. error( 'Unknown target: %s', opt.target );
  130. end
  131. P(P<0) = 0;
  132. P(P>1) = 1;
  133. if( strcmp( opt.colormap, 'trichromatic' ) || strcmp( opt.colormap, 'print' ) )
  134. color_map = [0.2 0.2 1.0;
  135. 0.2 1.0 1.0;
  136. 0.2 1.0 0.2;
  137. 1.0 1.0 0.2;
  138. 1.0 0.2 0.2];
  139. color_map_in = [0 0.25 0.5 0.75 1];
  140. elseif( strcmp( opt.colormap, 'dichromatic' ) )
  141. color_map = [0.2 1.0 1.0;
  142. 1.0 1.0 1.0
  143. 1.0 1.0 0.2];
  144. color_map_in = [0 0.5 1];
  145. elseif( strcmp( opt.colormap, 'monochromatic' ) )
  146. color_map = [1.0 1.0 1.0;
  147. 1.0 1.0 1.0];
  148. color_map_in = [0 1];
  149. else
  150. error( 'Unknown colormap: %s', opt.colormap );
  151. end
  152. if( strcmp( opt.target, 'screen' ) )
  153. color_map_l = color_map * [0.2126 0.7152 0.0722]'; %sum(color_map,2);
  154. color_map_ch = color_map ./ repmat( color_map_l, [1 3] );
  155. else
  156. if( strcmp( opt.colormap, 'monochromatic' ) )
  157. color_map_l = (color_map * [0.2126 0.7152 0.0722]') ./ color_map_in';
  158. else
  159. % luminance map start at 0.3, so that colors are visible
  160. color_map_l = (color_map * [0.2126 0.7152 0.0722]') ./ (color_map_in'*0.8+0.2);
  161. end
  162. color_map_ch = color_map ./ repmat( color_map_l, [1 3] );
  163. end
  164. %The line belows display pixel values
  165. %round(min( color_map_ch*255, 255 ))
  166. map = zeros( size(P,1), size(P,2), 3 );
  167. map(:,:,1) = interp1( color_map_in, color_map_ch(:,1), P );
  168. map(:,:,2) = interp1( color_map_in, color_map_ch(:,2), P );
  169. map(:,:,3) = interp1( color_map_in, color_map_ch(:,3), P );
  170. %map(:,:,3) = 1 - map(:,:,2) - map(:,:,1);
  171. %map = repmat( tmo_img, [1 1 3] );
  172. map = map .* repmat( tmo_img, [1 1 3] );
  173. end
  174. function l = log_luminance( X )
  175. if( size(X,3) == 3 )
  176. Y = X(:,:,1) * 0.212656 + X(:,:,2) * 0.715158 + X(:,:,3) * 0.072186;
  177. else
  178. Y = X;
  179. end
  180. Y(Y<=0) = min(Y(Y>0));
  181. l = log(Y);
  182. end
  183. function tmo_img = vis_tonemap( b, dr )
  184. t = 3;
  185. b_min = min(b(:));
  186. b_max = max(b(:));
  187. b_scale = linspace( b_min, b_max, 1024 );
  188. b_p = hist( b(:), b_scale );
  189. b_p = b_p / sum( b_p(:) );
  190. sum_b_p = sum( b_p.^(1/t) );
  191. dy = b_p.^(1/t) / sum_b_p;
  192. v = cumsum( dy )*dr + (1-dr)/2;
  193. tmo_img = interp1( b_scale, v, b );
  194. end
  195. function Y = blur_gaussian( X, sigma )
  196. ksize2 = round(sigma*3);
  197. gauss_1d = exp( -(-ksize2:ksize2).^2/(2*sigma^2) );
  198. gauss_1d = gauss_1d/sum(gauss_1d);
  199. Y = conv2( X, gauss_1d, 'same' );
  200. Y = conv2( Y, gauss_1d', 'same' );
  201. end
  202. function varargout = set_params( vals, varargin )
  203. for kk=1:length(varargin)
  204. if( length(vals) >= kk )
  205. varargout(kk) = vals(kk);
  206. else
  207. if( ischar( varargin{kk} ) && strcmp( 'required', varargin{kk} ) )
  208. error( 'Required parameter missing' );
  209. else
  210. varargout(kk) = varargin(kk);
  211. end
  212. end
  213. end
  214. end
  215. function I = norm_diff_img( test, reference )
  216. % Computer contrast-normalized difference image
  217. D = get_luminance(test)-get_luminance(reference);
  218. sigma = 5;
  219. window = fspecial('gaussian', round(sigma*4), sigma);
  220. img_mu = filter2(window, D, 'same');
  221. sigma_sq = max(0, filter2(window, D.^2, 'same' ) - img_mu.^2);
  222. v = ( sigma_sq ).^(1/2);
  223. I = min(D./(v+1),1);
  224. end
  225. function Y = get_luminance( img )
  226. % Return 2D matrix of luminance values for 3D matrix with an RGB image
  227. %dims = sum(nnz( size(img)>1 ));
  228. dims = find(size(img)>1,1,'last');
  229. if( dims == 3 )
  230. Y = img(:,:,1) * 0.212656 + img(:,:,2) * 0.715158 + img(:,:,3) * 0.072186;
  231. elseif( dims == 1 || dims == 2 )
  232. Y = img;
  233. else
  234. error( 'get_luminance: wrong matrix dimension' );
  235. end
  236. end
  237. function opt = parse_options( options )
  238. opt = struct();
  239. valid_opts = { 'colormap', 'trichromatic', 'target', 'screen', 'context_image', [] };
  240. for kk=1:2:length(options)
  241. opt_found = false;
  242. for ll=1:2:length(valid_opts)
  243. if( strcmp( options{kk}, valid_opts{ll} ) )
  244. opt_found = true;
  245. valid_opts{ll+1} = 'done';
  246. break;
  247. end
  248. end
  249. if( ~opt_found )
  250. error( 'Unrecognized option: %s', options{kk} );
  251. end
  252. opt.(options{kk}) = options{kk+1};
  253. end
  254. %Insert default options
  255. for kk=1:2:length(valid_opts)
  256. if( ~ischar( valid_opts{kk+1} ) || ~strcmp( valid_opts{kk+1}, 'done' ) )
  257. opt.(valid_opts{kk}) = valid_opts{kk+1};
  258. end
  259. end
  260. end