hdrvdp_visual_pathway.m 8.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290
  1. function [bands, L_adapt, band_freq, bb_padvalue] = hdrvdp_visual_pathway( img, name, metric_par, bb_padvalue )
  2. % HDRVDP_VISUAL_PATHWAY (internal) Process image along the visual pathway
  3. % to compute normalized perceptual response
  4. %
  5. % img - image data (can be multi-spectral)
  6. % name - string with the name of this map (shown in warnings and error
  7. % messages)
  8. % options - cell array with the 'option', value pairs
  9. % bands - CSF normalized freq. bands
  10. %
  11. % Copyright (c) 2011, Rafal Mantiuk <mantiuk@gmail.com>
  12. % Permission to use, copy, modify, and/or distribute this software for any
  13. % purpose with or without fee is hereby granted, provided that the above
  14. % copyright notice and this permission notice appear in all copies.
  15. %
  16. % THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  17. % WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  18. % MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  19. % ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  20. % WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  21. % ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  22. % OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  23. global hdrvdp_cache;
  24. if( any( isnan( img(:) ) ) )
  25. warning( 'hdrvdp:BadImageData', '%s image contains NaN pixel values', name );
  26. img(isnan(img)) = 1e-5;
  27. end
  28. % =================================
  29. % Precompute
  30. %
  31. % precompute common variables and structures or take them from cache
  32. % =================================
  33. width = size(img,2);
  34. height = size(img,1);
  35. img_sz = [height width]; % image size
  36. img_ch = size(img,3); % number of color channels
  37. rho2 = create_cycdeg_image( img_sz*2, metric_par.pix_per_deg ); % spatial frequency for each FFT coefficient, for 2x image size
  38. if( metric_par.do_mtf )
  39. mtf_filter = hdrvdp_mtf( rho2, metric_par ); % only the last param is adjusted
  40. else
  41. mtf_filter = ones( size(rho2) );
  42. end
  43. if( ~metric_par.debug )
  44. % memory savings for huge images
  45. clear rho2;
  46. end
  47. % Load spectral sensitivity curves
  48. [lambda, LMSR_S] = load_spectral_resp( 'log_cone_smith_pokorny_1975.csv' );
  49. LMSR_S(LMSR_S==0) = min(LMSR_S(:));
  50. LMSR_S = 10.^LMSR_S;
  51. [~, ROD_S] = load_spectral_resp( 'cie_scotopic_lum.txt' );
  52. LMSR_S(:,4) = ROD_S;
  53. IMG_E = metric_par.spectral_emission;
  54. % =================================
  55. % Precompute photoreceptor non-linearity
  56. % =================================
  57. pn = hdrvdp_get_from_cache( 'pn', [metric_par.rod_sensitivity metric_par.csf_sa], @() create_pn_jnd( metric_par ) );
  58. pn.jnd{1} = pn.jnd{1} * metric_par.sensitivity_correction;
  59. pn.jnd{2} = pn.jnd{2} * metric_par.sensitivity_correction;
  60. % =================================
  61. % Optical Transfer Function
  62. % =================================
  63. L_O = zeros( size(img) );
  64. for k=1:img_ch % for each color channel
  65. if( metric_par.do_mtf )
  66. % Use per-channel or one per all channels surround value
  67. pad_value = metric_par.surround_l( k );
  68. L_O(:,:,k) = clamp( fast_conv_fft( double(img(:,:,k)), mtf_filter, pad_value ), 1e-5, 1e10 );
  69. else
  70. % NO mtf
  71. L_O(:,:,k) = img(:,:,k);
  72. end
  73. end
  74. if( ~metric_par.debug )
  75. % memory savings for huge images
  76. clear mtf_filter;
  77. end
  78. %TODO - MTF reduces luminance values
  79. % =================================
  80. % Photoreceptor spectral sensitivity
  81. % =================================
  82. M_img_lmsr = zeros( img_ch, 4 ); % Color space transformation matrix
  83. for k=1:4
  84. for l=1:img_ch
  85. M_img_lmsr(l,k) = trapz( lambda, LMSR_S(:,k).*IMG_E(:,l) );
  86. end
  87. end
  88. % Color space conversion
  89. R_LMSR = reshape( reshape( L_O, width*height, img_ch )*M_img_lmsr, height, width, 4 );
  90. if( ~metric_par.debug )
  91. % memory savings for huge images
  92. clear L_O;
  93. end
  94. %surround_LMSR = metric_par.surround_l * M_img_lmsr;
  95. % =================================
  96. % Adapting luminance
  97. % =================================
  98. L_adapt = R_LMSR(:,:,1) + R_LMSR(:,:,2);
  99. % =================================
  100. % Photoreceptor non-linearity
  101. % =================================
  102. %La = mean( L_adapt(:) );
  103. P_LMR = zeros(height, width, 4);
  104. %surround_P_LMR = zeros(1,4);
  105. for k=[1:2 4] % ignore S - does not influence luminance
  106. if( k==4 )
  107. ph_type = 2; % rod
  108. ii = 3;
  109. else
  110. ph_type = 1; % cone
  111. ii = k;
  112. end
  113. P_LMR(:,:,ii) = pointOp( log10( clamp(R_LMSR(:,:,k), 10^pn.Y{ph_type}(1), 10^pn.Y{ph_type}(end)) ), ...
  114. pn.jnd{ph_type}, pn.Y{ph_type}(1), pn.Y{ph_type}(2)-pn.Y{ph_type}(1), 0 );
  115. % surround_P_LMR(ii) = interp1( pn_Y{ph_type}, pn_jnd{ph_type}, ...
  116. % log10( clamp(surround_LMSR(k), 10^pn_Y{ph_type}(1), 10^pn_Y{ph_type}(end)) ) );
  117. end
  118. if( ~metric_par.debug )
  119. % memory savings for huge images
  120. clear R_LMSR;
  121. end
  122. % =================================
  123. % Remove the DC component, from
  124. % cone and rod pathways separately
  125. % =================================
  126. % TODO - check if there is a better way to do it
  127. % cones
  128. P_C = P_LMR(:,:,1)+P_LMR(:,:,2);
  129. mm = mean(mean( P_C ));
  130. P_C = P_C - mm;
  131. % rods
  132. mm = mean(mean( P_LMR(:,:,3) ));
  133. P_R = P_LMR(:,:,3) - mm;
  134. % =================================
  135. % Achromatic response
  136. % =================================
  137. P = P_C + P_R;
  138. if( ~metric_par.debug )
  139. % memory savings for huge images
  140. clear P_LMR P_C P_R;
  141. end
  142. %surround_P = surround_P_LMR(1)+surround_P_LMR(2)+surround_P_LMR(3);
  143. % =================================
  144. % Multi-channel decomposition
  145. % =================================
  146. %[lo0filt,hi0filt,lofilt,bfilts,steermtx,harmonics] = eval(metric_par.steerpyr_filter);
  147. %max_ht = maxPyrHt(size(P), size(lofilt,1));
  148. %[bands.pyr,bands.pind] = buildSpyr( P, min(max_ht,6), metric_par.steerpyr_filter );
  149. [bands.pyr,bands.pind] = buildSpyr( P, 'auto', metric_par.steerpyr_filter );
  150. bands.sz = ones( spyrHt( bands.pind ) + 2, 1 );
  151. bands.sz(2:end-1) = spyrNumBands( bands.pind );
  152. band_freq = 2.^-(0:(spyrHt( bands.pind )+1)) * metric_par.pix_per_deg / 2;
  153. % CSF-filter the base band
  154. L_mean = mean( L_adapt(:) );
  155. BB = pyrBand( bands.pyr, bands.pind, sum(bands.sz) );
  156. % I wish to know why sqrt(2) works better, as it should be 2
  157. rho_bb = create_cycdeg_image( size(BB)*2, band_freq(end)*2*sqrt(2) );
  158. csf_bb = hdrvdp_ncsf( rho_bb, L_mean, metric_par );
  159. % pad value must be the same for test and reference images.
  160. if( bb_padvalue == -1 )
  161. bb_padvalue = mean(BB(:));
  162. end
  163. bands.pyr(pyrBandIndices(bands.pind,sum(bands.sz))) = fast_conv_fft( BB, csf_bb, bb_padvalue );
  164. if( 0 )
  165. for b=1:length(bands.sz)
  166. for o=1:bands.sz(b)
  167. b_ind = sum(bands.sz(1:(b-1)))+o;
  168. BB = pyrBand( bands.pyr, bands.pind, b_ind );
  169. rho_bb = create_cycdeg_image( size(BB)*2, band_freq(b)*4 );
  170. csf_bb = hdrvdp_csf( rho_bb, L_mean, metric_par );
  171. if( metric_par.surround_l == -1 )
  172. pad_value = mean(BB(:));
  173. else
  174. pad_value = metric_par.surround_l;
  175. end
  176. bands.pyr(pyrBandIndices(bands.pind,b_ind)) = fast_conv_fft( BB, csf_bb, pad_value );
  177. end
  178. end
  179. end
  180. function item = cache_get( item_name, item_signature, compute_func )
  181. sign_name = [ item_name '_sign' ];
  182. if( isfield( hdrvdp_cache, sign_name ) && all( hdrvdp_cache.( sign_name ) == item_signature) ) % caching
  183. item = hdrvdp_cache.( item_name );
  184. else
  185. item = compute_func();
  186. hdrvdp_cache.( sign_name ) = zeros(size(item_signature)); % in case of breaking at this point
  187. hdrvdp_cache.( item_name ) = item;
  188. hdrvdp_cache.( sign_name ) = item_signature;
  189. end
  190. end
  191. end
  192. function [Y jnd] = build_jndspace_from_S(l,S)
  193. L = 10.^l;
  194. dL = zeros(size(L));
  195. for k=1:length(L)
  196. thr = L(k)/S(k);
  197. % Different than in the paper because integration is done in the log
  198. % domain - requires substitution with a Jacobian determinant
  199. dL(k) = 1/thr * L(k) * log(10);
  200. end
  201. Y = l;
  202. jnd = cumtrapz( l, dL );
  203. end
  204. function pn = create_pn_jnd( metric_par )
  205. % Create lookup tables for intensity -> JND mapping
  206. c_l = logspace( -5, 5, 2048 );
  207. s_A = hdrvdp_joint_rod_cone_sens( c_l, metric_par );
  208. s_R = hdrvdp_rod_sens( c_l, metric_par ) * 10.^metric_par.rod_sensitivity;
  209. % s_C = s_L = S_M
  210. s_C = 0.5 * interp1( c_l, max(s_A-s_R, 1e-3), min( c_l*2, c_l(end) ) );
  211. pn = struct();
  212. [pn.Y{1} pn.jnd{1}] = build_jndspace_from_S( log10(c_l), s_C );
  213. [pn.Y{2} pn.jnd{2}] = build_jndspace_from_S( log10(c_l), s_R );
  214. end