123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290 |
- function [bands, L_adapt, band_freq, bb_padvalue] = hdrvdp_visual_pathway( img, name, metric_par, bb_padvalue )
- % HDRVDP_VISUAL_PATHWAY (internal) Process image along the visual pathway
- % to compute normalized perceptual response
- %
- % img - image data (can be multi-spectral)
- % name - string with the name of this map (shown in warnings and error
- % messages)
- % options - cell array with the 'option', value pairs
- % bands - CSF normalized freq. bands
- %
- % 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.
- global hdrvdp_cache;
- if( any( isnan( img(:) ) ) )
- warning( 'hdrvdp:BadImageData', '%s image contains NaN pixel values', name );
- img(isnan(img)) = 1e-5;
- end
- % =================================
- % Precompute
- %
- % precompute common variables and structures or take them from cache
- % =================================
- width = size(img,2);
- height = size(img,1);
- img_sz = [height width]; % image size
- img_ch = size(img,3); % number of color channels
- rho2 = create_cycdeg_image( img_sz*2, metric_par.pix_per_deg ); % spatial frequency for each FFT coefficient, for 2x image size
- if( metric_par.do_mtf )
- mtf_filter = hdrvdp_mtf( rho2, metric_par ); % only the last param is adjusted
- else
- mtf_filter = ones( size(rho2) );
- end
- if( ~metric_par.debug )
- % memory savings for huge images
- clear rho2;
- end
- % Load spectral sensitivity curves
- [lambda, LMSR_S] = load_spectral_resp( 'log_cone_smith_pokorny_1975.csv' );
- LMSR_S(LMSR_S==0) = min(LMSR_S(:));
- LMSR_S = 10.^LMSR_S;
- [~, ROD_S] = load_spectral_resp( 'cie_scotopic_lum.txt' );
- LMSR_S(:,4) = ROD_S;
- IMG_E = metric_par.spectral_emission;
- % =================================
- % Precompute photoreceptor non-linearity
- % =================================
- pn = hdrvdp_get_from_cache( 'pn', [metric_par.rod_sensitivity metric_par.csf_sa], @() create_pn_jnd( metric_par ) );
- pn.jnd{1} = pn.jnd{1} * metric_par.sensitivity_correction;
- pn.jnd{2} = pn.jnd{2} * metric_par.sensitivity_correction;
- % =================================
- % Optical Transfer Function
- % =================================
- L_O = zeros( size(img) );
- for k=1:img_ch % for each color channel
- if( metric_par.do_mtf )
- % Use per-channel or one per all channels surround value
- pad_value = metric_par.surround_l( k );
- L_O(:,:,k) = clamp( fast_conv_fft( double(img(:,:,k)), mtf_filter, pad_value ), 1e-5, 1e10 );
- else
- % NO mtf
- L_O(:,:,k) = img(:,:,k);
- end
- end
- if( ~metric_par.debug )
- % memory savings for huge images
- clear mtf_filter;
- end
- %TODO - MTF reduces luminance values
- % =================================
- % Photoreceptor spectral sensitivity
- % =================================
- M_img_lmsr = zeros( img_ch, 4 ); % Color space transformation matrix
- for k=1:4
- for l=1:img_ch
- M_img_lmsr(l,k) = trapz( lambda, LMSR_S(:,k).*IMG_E(:,l) );
- end
- end
- % Color space conversion
- R_LMSR = reshape( reshape( L_O, width*height, img_ch )*M_img_lmsr, height, width, 4 );
- if( ~metric_par.debug )
- % memory savings for huge images
- clear L_O;
- end
- %surround_LMSR = metric_par.surround_l * M_img_lmsr;
- % =================================
- % Adapting luminance
- % =================================
- L_adapt = R_LMSR(:,:,1) + R_LMSR(:,:,2);
- % =================================
- % Photoreceptor non-linearity
- % =================================
- %La = mean( L_adapt(:) );
- P_LMR = zeros(height, width, 4);
- %surround_P_LMR = zeros(1,4);
- for k=[1:2 4] % ignore S - does not influence luminance
- if( k==4 )
- ph_type = 2; % rod
- ii = 3;
- else
- ph_type = 1; % cone
- ii = k;
- end
-
- P_LMR(:,:,ii) = pointOp( log10( clamp(R_LMSR(:,:,k), 10^pn.Y{ph_type}(1), 10^pn.Y{ph_type}(end)) ), ...
- pn.jnd{ph_type}, pn.Y{ph_type}(1), pn.Y{ph_type}(2)-pn.Y{ph_type}(1), 0 );
-
- % surround_P_LMR(ii) = interp1( pn_Y{ph_type}, pn_jnd{ph_type}, ...
- % log10( clamp(surround_LMSR(k), 10^pn_Y{ph_type}(1), 10^pn_Y{ph_type}(end)) ) );
- end
- if( ~metric_par.debug )
- % memory savings for huge images
- clear R_LMSR;
- end
- % =================================
- % Remove the DC component, from
- % cone and rod pathways separately
- % =================================
- % TODO - check if there is a better way to do it
- % cones
- P_C = P_LMR(:,:,1)+P_LMR(:,:,2);
- mm = mean(mean( P_C ));
- P_C = P_C - mm;
- % rods
- mm = mean(mean( P_LMR(:,:,3) ));
- P_R = P_LMR(:,:,3) - mm;
- % =================================
- % Achromatic response
- % =================================
- P = P_C + P_R;
- if( ~metric_par.debug )
- % memory savings for huge images
- clear P_LMR P_C P_R;
- end
- %surround_P = surround_P_LMR(1)+surround_P_LMR(2)+surround_P_LMR(3);
- % =================================
- % Multi-channel decomposition
- % =================================
- %[lo0filt,hi0filt,lofilt,bfilts,steermtx,harmonics] = eval(metric_par.steerpyr_filter);
- %max_ht = maxPyrHt(size(P), size(lofilt,1));
- %[bands.pyr,bands.pind] = buildSpyr( P, min(max_ht,6), metric_par.steerpyr_filter );
- [bands.pyr,bands.pind] = buildSpyr( P, 'auto', metric_par.steerpyr_filter );
- bands.sz = ones( spyrHt( bands.pind ) + 2, 1 );
- bands.sz(2:end-1) = spyrNumBands( bands.pind );
- band_freq = 2.^-(0:(spyrHt( bands.pind )+1)) * metric_par.pix_per_deg / 2;
- % CSF-filter the base band
- L_mean = mean( L_adapt(:) );
-
- BB = pyrBand( bands.pyr, bands.pind, sum(bands.sz) );
- % I wish to know why sqrt(2) works better, as it should be 2
- rho_bb = create_cycdeg_image( size(BB)*2, band_freq(end)*2*sqrt(2) );
- csf_bb = hdrvdp_ncsf( rho_bb, L_mean, metric_par );
-
- % pad value must be the same for test and reference images.
- if( bb_padvalue == -1 )
- bb_padvalue = mean(BB(:));
- end
-
- bands.pyr(pyrBandIndices(bands.pind,sum(bands.sz))) = fast_conv_fft( BB, csf_bb, bb_padvalue );
- if( 0 )
- for b=1:length(bands.sz)
- for o=1:bands.sz(b)
-
- b_ind = sum(bands.sz(1:(b-1)))+o;
- BB = pyrBand( bands.pyr, bands.pind, b_ind );
- rho_bb = create_cycdeg_image( size(BB)*2, band_freq(b)*4 );
- csf_bb = hdrvdp_csf( rho_bb, L_mean, metric_par );
-
- if( metric_par.surround_l == -1 )
- pad_value = mean(BB(:));
- else
- pad_value = metric_par.surround_l;
- end
-
- bands.pyr(pyrBandIndices(bands.pind,b_ind)) = fast_conv_fft( BB, csf_bb, pad_value );
- end
-
- end
- end
- function item = cache_get( item_name, item_signature, compute_func )
- sign_name = [ item_name '_sign' ];
- if( isfield( hdrvdp_cache, sign_name ) && all( hdrvdp_cache.( sign_name ) == item_signature) ) % caching
- item = hdrvdp_cache.( item_name );
- else
- item = compute_func();
- hdrvdp_cache.( sign_name ) = zeros(size(item_signature)); % in case of breaking at this point
- hdrvdp_cache.( item_name ) = item;
- hdrvdp_cache.( sign_name ) = item_signature;
- end
- end
- end
- function [Y jnd] = build_jndspace_from_S(l,S)
- L = 10.^l;
- dL = zeros(size(L));
- for k=1:length(L)
- thr = L(k)/S(k);
- % Different than in the paper because integration is done in the log
- % domain - requires substitution with a Jacobian determinant
- dL(k) = 1/thr * L(k) * log(10);
- end
- Y = l;
- jnd = cumtrapz( l, dL );
- end
- function pn = create_pn_jnd( metric_par )
- % Create lookup tables for intensity -> JND mapping
- c_l = logspace( -5, 5, 2048 );
- s_A = hdrvdp_joint_rod_cone_sens( c_l, metric_par );
- s_R = hdrvdp_rod_sens( c_l, metric_par ) * 10.^metric_par.rod_sensitivity;
- % s_C = s_L = S_M
- s_C = 0.5 * interp1( c_l, max(s_A-s_R, 1e-3), min( c_l*2, c_l(end) ) );
- pn = struct();
- [pn.Y{1} pn.jnd{1}] = build_jndspace_from_S( log10(c_l), s_C );
- [pn.Y{2} pn.jnd{2}] = build_jndspace_from_S( log10(c_l), s_R );
- end
|