diff --git a/aberration_correction/CorrectByHyperspectralADMM.m b/aberration_correction/CorrectByHyperspectralADMM.m index 1d0767c..10c0a2b 100644 --- a/aberration_correction/CorrectByHyperspectralADMM.m +++ b/aberration_correction/CorrectByHyperspectralADMM.m @@ -58,6 +58,9 @@ % j-th colour channel or spectral band of the latent images. For example, % `sensor_map` is a matrix mapping discretized spectral power % distributions to RGB colours. +% - 'channel_mode': A Boolean value indicating whether the latent colour +% space is a set of colour channels (true) or a set of spectral bands +% (false). % The following variables are sometimes required: % - 'bands': A vector containing the wavelengths or colour channel indices % to use as the `lambda` input argument of 'polyfunToMatrix()' when @@ -198,12 +201,10 @@ 'downsampling_factor',... 'output_directory',... 'save_latent_image_files',... - 'add_border',... - 'full_GLambda',... + 'baek2017Algorithm2Options',... 'rho',... 'weights',... - 'tol',... - 'maxit'... + 'int_method'... }; %% Input data and parameters @@ -242,11 +243,11 @@ % Whether to expand the latent image relative to the input image to cover % all undistorted coordinates from the input image. This is the % `add_border` input argument. -add_border = false; +baek2017Algorithm2Options.add_border = false; % Whether to make the spectral gradient the same size as the image. This is % the `full_GLambda` input argument. -full_GLambda = false; +baek2017Algorithm2Options.full_GLambda = false; % Penalty parameters in ADMM, the `rho` input argument. % Sample values seem to be in the range 1-10 (see pages 89, 93, and 95 of @@ -260,10 +261,15 @@ % % Reasonable values for the third element are 10^-4 to 10^-3 (page 21 of % Boyd et al. 2011). -tol = [ 1e-3, 1e-2, 1e-3 ]; +baek2017Algorithm2Options.tol = [ 1e-3, 1e-2, 1e-3 ]; % Maximum number of inner and outer iterations, the `maxit` input argument -maxit = [ 20, 500 ]; +baek2017Algorithm2Options.maxit = [ 20, 500 ]; + +% If the latent space consists of wavelength bands, use this type of +% numerical integration in 'channelConversionMatrix()'. (Otherwise, a value +% of 'none' will automatically be used instead.) +int_method = 'trap'; % ## Debugging Flags baek2017Algorithm2Verbose = true; @@ -303,7 +309,7 @@ bands_polyfun = bands; bands = []; -model_variables_required = { 'sensor_map' }; +model_variables_required = { 'sensor_map', 'channel_mode' }; load(color_map_filename, model_variables_required{:}, optional_variable); if ~all(ismember(model_variables_required, who)) error('One or more of the required colour space conversion variables is not loaded.') @@ -311,6 +317,12 @@ bands_color = bands; +if channel_mode + baek2017Algorithm2Options.int_method = int_method; +else + baek2017Algorithm2Options.int_method = 'none'; +end + %% Preprocessing input data % Select the highest-priority value of `bands` @@ -376,9 +388,8 @@ [ I_latent{i}, image_bounds{i}, I_rgb, J_full, J_est ] = baek2017Algorithm2(... image_sampling, bayer_pattern, sensor_map_resampled,... - polyfun, bands, add_border,... - full_GLambda, I_raw, rho, weights,... - tol, maxit, baek2017Algorithm2Verbose... + polyfun, bands, I_raw, rho, weights,... + baek2017Algorithm2Options, baek2017Algorithm2Verbose... ); % Save the results diff --git a/aberration_correction/admm/baek2017Algorithm2.m b/aberration_correction/admm/baek2017Algorithm2.m index 5a93f3e..ff0d347 100644 --- a/aberration_correction/admm/baek2017Algorithm2.m +++ b/aberration_correction/admm/baek2017Algorithm2.m @@ -1,17 +1,15 @@ function [ I_3D, image_bounds, varargout ] = baek2017Algorithm2(... image_sampling, align, sensitivity,... - polyfun, lambda, add_border,... - full_GLambda, J_2D, rho, weights,... - tol, maxit, varargin... + polyfun, lambda, J_2D, rho, weights,... + options, varargin... ) % BAEK2017ALGORITHM2 Run ADMM as in Algorithm 2 of Baek et al. 2017 % % ## Syntax % I = baek2017Algorithm2(... % image_sampling, align, sensitivity,... -% polyfun, lambda, add_border,... -% full_GLambda, J, rho, weights,... -% tol, maxit [, verbose]... +% polyfun, lambda, J, rho, weights,... +% options [, verbose]... % ) % [ I, image_bounds ] = baek2017Algorithm2(___) % [ I, image_bounds, I_rgb ] = baek2017Algorithm2(___) @@ -21,9 +19,8 @@ % ## Description % I = baek2017Algorithm2(... % image_sampling, align, sensitivity,... -% polyfun, lambda, add_border,... -% full_GLambda, J, rho, weights,... -% tol, maxit [, verbose]... +% polyfun, lambda, J, rho, weights,... +% options [, verbose]... % ) % Estimate a latent RGB or hyperspectral image `I` from dispersion in % the input RAW image `J`. @@ -77,20 +74,6 @@ % encapsulated by `polyfun`. 'c' is the desired number of spectral bands % or colour channels in `I`. % -% add_border -- Padding flag -% A Boolean value indicating whether or not the `image_bounds` input -% argument of 'polyfunToMatrix()' should be empty. If `true`, the output -% image `I` will be large enough to contain the lateral chromatic -% aberration-corrected coordinates of all pixels in `J`. If `false`, -% the output image `I` will be clipped to the region occupied by -% `J`. -% -% full_GLambda -- Approximate spectral gradients at border bands -% A Boolean value used as the `replicate` input argument of -% 'spectralGradient()' when creating the spectral gradient matrix needed -% for regularizing `I`. Refer to the documentation of -% 'spectralGradient.m' for details. -% % J -- Input RAW image % A 2D array containing the raw colour-filter pattern data of an image. % @@ -105,25 +88,51 @@ % the `beta` weight on the regularization of the spectral gradient of the % spatial gradient of the image in the ADMM optimization problem. % -% If all elements of `weights` are zero, this function will find the -% minimum-norm least squares solution to the simpler problem: +% If all elements of `weights` are zero, this function will find a +% solution to the simpler problem: % argmin_i (||M * Omega * Phi * i - j||_2) ^ 2 % where `M` performs mosaicing, `Omega` converts colours to the RGB % colour space of the camera, and `Phi` warps the image according to the -% dispersion model. `j` is the input RAW image. -% -% tol -- Convergence tolerances -% The first element of `tol` is the tolerance value to use with MATLAB's -% 'pcg()' function when solving the I-minimization step of the ADMM -% algorithm. The second and third elements of `tol` are the absolute and -% relative tolerance values for the ADMM algorithm, as explained in -% Section 3.3.1 of Boyd et al. 2011. -% -% maxit -- Maximum number of iterations -% The first element of `maxit` contains the maximum number of iterations -% to use with MATLAB's 'pcg()' function when solving the I-minimization -% step of the ADMM algorithm. The second element of `maxit` contains the -% maximum number of ADMM iterations to perform. +% dispersion model. `j` is the input RAW image. If the linear system is +% underdetermined, the function will find the minimum-norm least squares +% solution. If the problem is sufficiently determined, as it may be in +% cases where `i` has fewer wavelength bands of colour channels than `j` +% has colour channels, then the function will find the exact or +% least-squares solution. +% +% This function will not terminate if `weights` has a mixture of zero and +% nonzero elements, because the ADMM algorithm will not converge. +% +% options -- Options and small parameters +% A structure with the following fields: +% - 'add_border': A Boolean value indicating whether or not the +% `image_bounds` input argument of 'polyfunToMatrix()' should be empty. +% If `true`, the output image `I` will be large enough to contain the +% lateral chromatic aberration-corrected coordinates of all pixels in +% `J`. If `false`, the output image `I` will be clipped to the region +% occupied by `J`. +% - 'full_GLambda': A Boolean value used as the `replicate` input +% argument of 'spectralGradient()' when creating the spectral gradient +% matrix needed for regularizing `I`. Refer to the documentation of +% 'spectralGradient.m' for details. +% - 'int_method': The numerical integration method used for spectral to +% colour space conversion. `int_method` is passed to +% 'channelConversionMatrix()' as its `int_method` argument. Refer to +% the documentation of 'channelConversionMatrix.m' for details. If +% 'int_method' is 'none', as should be the case when colour conversion +% is from a set of colour channels, not a spectral space, numerical +% integration will not be performed. +% - 'maxit': Maximum number of iterations: The first element of `maxit` +% contains the maximum number of iterations to use with MATLAB's +% 'pcg()' function when solving the I-minimization step of the ADMM +% algorithm. The second element of `maxit` contains the maximum number +% of ADMM iterations to perform. +% - 'tol': Convergence tolerances: The first element of `tol` is the +% tolerance value to use with MATLAB's 'pcg()' function when solving +% the I-minimization step of the ADMM algorithm. The second and third +% elements of `tol` are the absolute and relative tolerance values for +% the ADMM algorithm, as explained in Section 3.3.1 of Boyd et al. +% 2011. % % verbose -- Verbosity flag % If `true`, console output will be displayed to show the progress of the @@ -190,7 +199,7 @@ % File created May 27, 2018 nargoutchk(1, 5); -narginchk(12, 13); +narginchk(9, 10); if ~isempty(varargin) verbose = varargin{1}; @@ -202,14 +211,24 @@ error('The `image_sampling` input argument must contain an image height and width only.'); end +if isStringScalar(options.int_method) || ischar(options.int_method) + do_integration = ~strcmp(options.int_method, 'none'); +else + error('`options.int_method` must be a character vector or a string scalar.'); +end + J = J_2D(:); % Create constant matrices image_sampling_J = size(J_2D); n_bands = length(lambda); M = mosaicMatrix(image_sampling_J, align); -Omega = channelConversionMatrix(image_sampling_J, sensitivity); -if add_border +if do_integration + Omega = channelConversionMatrix(image_sampling_J, sensitivity, lambda, options.int_method); +else + Omega = channelConversionMatrix(image_sampling_J, sensitivity); +end +if options.add_border image_bounds = []; else image_bounds = [0, 0, image_sampling_J(2), image_sampling_J(1)]; @@ -219,13 +238,18 @@ ); if all(weights == 0) - % Minimum-norm least squares solution to the non-regularized problem A = M * Omega * Phi; - I = lsqminnorm(A, J); + if (size(A, 1) < size(A, 2)) || (rank(A) < size(A, 2)) + % Minimum-norm least squares solution to the non-regularized problem + I = lsqminnorm(A, J); + else + % Unique or least-squares solution + I = A \ J; + end else % Perform ADMM G_xy = spatialGradient([image_sampling, n_bands]); - G_lambda = spectralGradient([image_sampling, n_bands], full_GLambda); + G_lambda = spectralGradient([image_sampling, n_bands], options.full_GLambda); G_lambda_sz1 = size(G_lambda, 1); G_lambda_sz2 = size(G_lambda, 2); % The product `G_lambda * G_xy` must be defined, so `G_lambda` needs to be @@ -258,10 +282,10 @@ G_xy_T = G_xy.'; G_lambda_xy_T = G_lambda_xy.'; converged = false; - for iter = 1:maxit(2) + for iter = 1:options.maxit(2) % Optimization [ I, flag, relres, iter_pcg ] = pcg(... - A, b, tol(1), maxit(1), [], [], I... + A, b, options.tol(1), options.maxit(1), [], [], I... ); if(verbose) fprintf('%d: PCG (flag = %d, relres = %g, iter = %d)\n',... @@ -295,17 +319,17 @@ % Calculate stopping criteria % See Section 3.3.1 of Boyd et al. 2011. - epsilon_pri1 = sqrt(len_Z1) * tol(2) +... - tol(3) * max([norm(g_xy), norm(Z1)]); - epsilon_pri2 = sqrt(len_Z2) * tol(2) +... - tol(3) * max([norm(g_lambda_xy), norm(Z2)]); + epsilon_pri1 = sqrt(len_Z1) * options.tol(2) +... + options.tol(3) * max([norm(g_xy), norm(Z1)]); + epsilon_pri2 = sqrt(len_Z2) * options.tol(2) +... + options.tol(3) * max([norm(g_lambda_xy), norm(Z2)]); Y1 = rho(1) * U1; Y2 = rho(2) * U2; - epsilon_dual1 = sqrt(len_I) * tol(2) +... - tol(3) * norm(G_xy_T * Y1); - epsilon_dual2 = sqrt(len_I) * tol(2) +... - tol(3) * norm(G_lambda_xy_T * Y2); + epsilon_dual1 = sqrt(len_I) * options.tol(2) +... + options.tol(3) * norm(G_xy_T * Y1); + epsilon_dual2 = sqrt(len_I) * options.tol(2) +... + options.tol(3) * norm(G_lambda_xy_T * Y2); if(verbose) fprintf('Stop Crit. (e_p1 = %g, e_p2 = %g, e_d2 = %g, e_d2 = %g)\n',... @@ -335,7 +359,11 @@ I_3D = reshape(I, image_sampling(1), image_sampling(2), n_bands); if nargout > 2 - Omega_I = channelConversionMatrix(image_sampling, sensitivity); + if do_integration + Omega_I = channelConversionMatrix(image_sampling, sensitivity, lambda, options.int_method); + else + Omega_I = channelConversionMatrix(image_sampling, sensitivity); + end I_rgb = Omega_I * I; n_channels_rgb = 3; varargout{1} = reshape(I_rgb, image_sampling(1), image_sampling(2), n_channels_rgb); diff --git a/aberration_correction/channelConversionMatrix.m b/aberration_correction/channelConversionMatrix.m index 55540e0..e5fd28a 100644 --- a/aberration_correction/channelConversionMatrix.m +++ b/aberration_correction/channelConversionMatrix.m @@ -1,11 +1,15 @@ -function [ Omega ] = channelConversionMatrix(image_sampling, sensitivity) +function [ Omega ] = channelConversionMatrix(image_sampling, sensitivity, varargin) % CHANNELCONVERSIONMATRIX Create a sparse matrix to convert images between colour spaces % % ## Syntax -% Omega = channelConversionMatrix(image_sampling, sensitivity) +% Omega = channelConversionMatrix(... +% image_sampling, sensitivity [, lambda, int_method]... +% ) % % ## Description -% Omega = channelConversionMatrix(image_sampling, sensitivity) +% Omega = channelConversionMatrix(... +% image_sampling, sensitivity [, lambda, int_method]... +% ) % Returns a matrix for changing the spectral representation of an image. % % ## Input Arguments @@ -20,6 +24,26 @@ % hyperspectral representation to colours in the output colour space or % hyperspectral representation. % +% lambda -- Wavelength bands +% A vector, of length equal to the size of the second dimension of +% `sensitivity`, containing the wavelengths of the input discrete +% spectral space. `lambda(j)` is the wavelength or colour channel index +% corresponding to `sensitivity(:, j)`. `lambda` is used to scale the +% elements of the colour conversion matrix by the spacings between +% wavelengths, so that the matrix performs numerical integration of the +% products of spectra with the sensor's spectral sensitivities. +% +% `lambda` and `int_method` should not be passed if the input colour +% space consists of colour channels, not wavelength bands. +% +% int_method -- Numerical integration method +% The numerical integration method to use. `int_method` is passed to +% `integrationWeights()` as its `method` input argument. `int_method` +% should only be passed if the input colour space consists of wavelength +% bands, as discussed in the description of `lambda` above. +% +% Defaults to 'trap' if not passed. +% % ## Output Arguments % % Omega -- Colour channel conversion matrix @@ -35,7 +59,7 @@ % `J` is a vectorized form of the image in the output colour space, with `c2` % colour channels. % -% See also sonyQuantumEfficiency +% See also sonyQuantumEfficiency, integrationWeights % Bernard Llanos % Supervised by Dr. Y.H. Yang @@ -43,12 +67,23 @@ % File created May 24, 2018 nargoutchk(1, 1); -narginchk(2, 2); +narginchk(2, 4); if length(image_sampling) ~= 2 error('The `image_sampling` input argument must contain an image height and width only.'); end +channel_mode = true; +if ~isempty(varargin) + channel_mode = false; + lambda = varargin{1}; + if length(varargin) > 1 + int_method = varargin{2}; + else + int_method = 'trap'; + end +end + n_px = prod(image_sampling); c1 = size(sensitivity, 2); c2 = size(sensitivity, 1); @@ -63,10 +98,16 @@ % Column indices % Iterate over the input colour channels for each pixel, and repeat for each % output colour channel -columns = repmat(repelem((1:n_px).', c1), c2, 1) + repmat((0:(c1 - 1)).' * n_px, n_px_c2, 1); +columns = repmat(repelem((1:n_px).', c1), c2, 1) +... + repmat((0:(c1 - 1)).' * n_px, n_px_c2, 1); % Matrix values -elements = repmat(sensitivity.', n_px, 1); +if channel_mode + weights = ones(size(sensitivity)); +else + weights = repmat(integrationWeights(lambda, int_method), 1, c1); +end +elements = repmat((sensitivity.') .* weights, n_px, 1); elements = elements(:); % Assemble the sparse matrix diff --git a/aberration_data/BimaterialImages.m b/aberration_data/BimaterialImages.m index 8192620..fce593a 100644 --- a/aberration_data/BimaterialImages.m +++ b/aberration_data/BimaterialImages.m @@ -44,12 +44,6 @@ % - 'model_from_reference': A parameter of the above scripts, which % determines the frame of reference for the model of chromatic % aberration. It must be set to `false`. -% The following variables are sometimes required: -% - 'bands': A vector containing the wavelengths to use as the `lambda` -% input argument of 'polyfunToMatrix()' (to evaluate a dispersion model), -% and to form hyperspectral images. This variable is required only if not -% provided in the colour space conversion data file, or directly in this -% script (see below). % % The following variables are optional. If they are present, they are % assumed to define a conversion between a geometrical optics coordinate @@ -70,40 +64,24 @@ % - 'sensor_map': A 2D array, where `sensor_map(i, j)` is the sensitivity % of the i-th colour channel of the output sensor response images to the % j-th spectral band of the synthetic hyperspectral images. -% The following variables are sometimes required: % - 'bands': A vector containing the wavelengths to use as the `lambda` % input argument of 'polyfunToMatrix()' (to evaluate a dispersion model), -% and to form hyperspectral images. This variable takes precedence over -% the same variable provided with the dispersion model (see above), but -% can be overridden by a version of the variable provided directly in -% this script (see below). +% and to form hyperspectral images. This variable can be overridden by a +% version of the variable provided directly in this script (see below). +% However, it is still needed from this file to normalize spectral +% radiances, even when overridden in its other functions. % % ### Discrete spectral space % % The 'bands' variable defined in the parameters section of the code below -% can be empty (`[]`), in which case it is loaded from the dispersion model -% data or from the colour space conversion data (see above). -% -% The final value of 'bands' is determined according to the following list, -% in order by decreasing priority: -% - 'bands' defined in this script -% - 'bands' loaded from with colour space conversion data -% - 'bands' loaded from a polynomial model of dispersion +% can be empty (`[]`), in which case it is loaded from the colour space +% conversion data (see above). % % The final 'bands' vector is clipped to the intervals defined by the other % vectors of wavelengths, to avoid extrapolation when resampling data to % conform to the final value of 'bands'. Resampling is needed to express % all spectral quantities in the same discrete space of wavelengths. % -% If no variable 'bands' was found in the colour space conversion data, but -% the final value of 'bands' has a length equal to the size of the second -% dimension of 'sensor_map', 'bands' is assumed to be compatible with -% 'sensor_map'. Otherwise, 'sensor_map' is resampled along its second -% dimension, by assuming that its first and last columns correspond to the -% first and last elements of 'bands', and that its remaining columns -% correspond to equally-spaced wavelengths in-between the first and last -% elements of 'bands'. -% % ### Illuminant % The following two sources of data can be provided: % @@ -147,22 +125,27 @@ % the per-pixel weights in the soft segmentation of the chromaticity map, % and then by illuminating the spectra according to the illumination % weights in the illumination map. Image values are normalized radiances. -% - '*_3.tif': A colour image created by converting the hyperspectral -% image to the raw colour space of the camera. +% - '*_3.tif' and '*_3.mat': A colour image created by converting the +% hyperspectral image to the raw colour space of the camera. % - '*_hyperspectral_warped.mat': A warped version of the hyperspectral % image (stored in the variable 'I_hyper_warped') created by applying the % dispersion model to the image. -% - '*_3_warped.tif': A colour image created by converting the -% warped hyperspectral image to the raw colour space of the camera. -% - '*_raw_warped.tif': A colour-filter array image produced by mosaicing -% the warped sensor response image according to the colour-filter pattern -% of the camera. +% - '*_3_warped.tif' and '*_3_warped.mat': A colour image created by +% converting the warped hyperspectral image to the raw colour space of +% the camera. +% - '*_raw_warped.tif' and '*_raw_warped.mat': A colour-filter array image +% produced by mosaicing the warped sensor response image according to the +% colour-filter pattern of the camera. % % The raw colour space of the camera is determined by the colour space % conversion data provided as input to this script. A camera may apply % additional operations to convert sensor responses from the raw colour % space to, for example, sRGB colours. % +% Note that both '.mat' and '.tif' files are output for monochromatic or +% three-channel images to provide both easy display ('.tif' files) and +% lossless storage ('.mat' files). +% % ### Data file output % % A '.mat' file containing the following variables: @@ -170,13 +153,7 @@ % - 'bands': The final value of the 'bands' variable, determined as % discussed under the section "Discrete spectral space" above. % - 'bands_color': The 'bands' variable loaded from the colour space -% conversion data file, for reference. 'bands_color' is empty if no -% such variable was found in the data file, or if the value of the loaded -% variable was empty. -% - 'bands_polyfun': The 'bands' variable loaded from the dispersion model -% data file, for reference. 'bands_polyfun' is empty if no such variable -% was found in the data file, or if the value of the loaded variable was -% empty. +% conversion data file, for reference. % - 'chromaticity_filenames': A cell vector containing the input % chromaticity map filenames retrieved based on the wildcard provided in % the parameters section of the script. @@ -214,6 +191,8 @@ 'polynomial_model_filename',... 'dispersonFun_name',... 'color_map_filename',... + 'normalization_channel',... + 'int_method',... 'bands_interp_method',... 'use_cie_illuminant',... 'illuminant_filename',... @@ -221,7 +200,6 @@ 'illuminant_name',... 'illuminant_function_name',... 'xyzbar_filename',... - 'lambda_xyzbar',... 'reflectances_filename',... 'n_colors',... 'bayer_pattern',... @@ -248,7 +226,11 @@ end % Colour space conversion data -color_map_filename = '/home/llanos/GoogleDrive/ThesisResearch/Data and Results/20180615_TestingBimaterialImages/SonyColorMapData.mat'; +color_map_filename = '/home/llanos/GoogleDrive/ThesisResearch/Data and Results/20180629_TestingBimaterialImages/SonyColorMapData.mat'; +% Colour channel to use for radiance normalization +normalization_channel = 2; +% Integration method to use for colour calculations +int_method = 'trap'; % Override the wavelengths at which to evaluate the model of dispersion, if % desired. @@ -284,7 +266,7 @@ bayer_pattern = 'gbrg'; % Output directory for all images and saved data -output_directory = '/home/llanos/GoogleDrive/ThesisResearch/Data and Results/20180615_TestingBimaterialImages'; +output_directory = '/home/llanos/GoogleDrive/ThesisResearch/Data and Results/20180629_TestingBimaterialImages'; % ## Debugging Flags segmentColorsVerbose = false; @@ -320,7 +302,7 @@ lambda_illuminant, spd_illuminant,... lambda_colorChecker, reflectances,... lambda_xyzbar, xyzbar,... - illuminant_name... + illuminant_name, int_method... ); %% Load dispersion model @@ -359,8 +341,8 @@ %% Load colour space conversion data -model_variables_required = { 'sensor_map' }; -load(color_map_filename, model_variables_required{:}, optional_variable); +model_variables_required = { 'sensor_map', optional_variable }; +load(color_map_filename, model_variables_required{:}); if ~all(ismember(model_variables_required, who)) error('One or more of the required colour space conversion variables is not loaded.') end @@ -372,52 +354,33 @@ % Select the highest-priority value of `bands` if ~isempty(bands_script) bands = bands_script; -elseif ~isempty(bands_color) - bands = bands_color; -elseif ~isempty(bands_polyfun) - bands = bands_polyfun; else - error('The variable `bands` is not defined, or is empty'); + bands = bands_color; end % Use the intersection of all values of `bands` corresponding to arrays -if ~isempty(bands_color) - bands = bands(bands >= min(bands_color) & bands <= max(bands_color)); -end +bands = bands(bands >= min(bands_color) & bands <= max(bands_color)); bands = bands(bands >= min(lambda_colorChecker) & bands <= max(lambda_colorChecker)); -% Compare with colour space conversion data -n_bands = length(bands); -n_bands_sensor_map = size(sensor_map, 2); -resample_bands = false; -if ~isempty(bands_color) - if n_bands ~= length(bands_color) || any(bands ~= bands_color) - % Resampling is needed - resample_bands = true; - bands_for_interp = bands_color; - end -elseif n_bands_sensor_map ~= n_bands - % Resampling is needed, but will be "blind" - resample_bands = true; - bands_for_interp = linspace(bands(1), bands(end), n_bands_sensor_map); -end % Resample colour space conversion data if necessary -if resample_bands +if length(bands) ~= length(bands_color) || any(bands ~= bands_color) [sensor_map_resampled, bands] = resampleArrays(... - bands_for_interp, sensor_map.', bands,... + bands_color, sensor_map.', bands,... bands_interp_method... ); sensor_map_resampled = sensor_map_resampled.'; else sensor_map_resampled = sensor_map; end - n_bands = length(bands); %% Calculate spectral radiances [lambda_Rad, ~, Rad_normalized] = reflectanceToRadiance(... - lambda_illuminant, spd_illuminant, lambda_colorChecker, reflectances, lambda_xyzbar, xyzbar(:, 2)... + lambda_illuminant, spd_illuminant,... + lambda_colorChecker, reflectances,... + bands_color, sensor_map(normalization_channel, :).',... + int_method... ); % Resample radiances @@ -512,10 +475,10 @@ I_hyper = reshape(H, image_height, image_width, n_bands); % Compute the equivalent sensor response image - Omega = channelConversionMatrix(image_sampling, sensor_map_resampled); - raw_full = Omega * H; + Omega = channelConversionMatrix(image_sampling, sensor_map_resampled, bands, int_method); + raw_full = Omega * H; raw_full_3D = reshape(raw_full, image_height, image_width, n_channels_raw); - + % Simulate dispersion W = polyfunToMatrix(... dispersonFun, bands,... @@ -537,12 +500,22 @@ % Save the results output_filename = fullfile(output_directory, [chromaticity_names{i} '_hyperspectral.mat']); save(output_filename, 'I_hyper'); + + output_filename = fullfile(output_directory, [chromaticity_names{i} '_3.mat']); + save(output_filename, 'raw_full_3D'); output_filename = fullfile(output_directory, [chromaticity_names{i} '_3' ext]); imwrite(raw_full_3D, output_filename); + output_filename = fullfile(output_directory, [chromaticity_names{i} '_hyperspectral_warped.mat']); save(output_filename, 'I_hyper_warped'); + + output_filename = fullfile(output_directory, [chromaticity_names{i} '_3_warped.mat']); + save(output_filename, 'raw_warped_3D'); output_filename = fullfile(output_directory, [chromaticity_names{i} '_3_warped' ext]); imwrite(raw_warped_3D, output_filename); + + output_filename = fullfile(output_directory, [chromaticity_names{i} '_raw_warped.mat']); + save(output_filename, 'raw_2D'); output_filename = fullfile(output_directory, [chromaticity_names{i} '_raw_warped' ext]); imwrite(raw_2D, output_filename); end diff --git a/aberration_data/cieSpectralToColor.m b/aberration_data/cieSpectralToColor.m index 09ca3ec..0efe26d 100644 --- a/aberration_data/cieSpectralToColor.m +++ b/aberration_data/cieSpectralToColor.m @@ -2,14 +2,14 @@ % CIESPECTRALTOCOLOR Convert spectral radiance to sRGB using the CIE tristimulus functions % % ## Syntax -% rgb = cieSpectralToColor(lambda_C, C, lambda_R, R [, whitepoint]) -% [rgb, XYZ] = cieSpectralToColor(lambda_C, C, lambda_R, R [, whitepoint]) +% rgb = cieSpectralToColor(lambda_C, C, lambda_R, R [, whitepoint, int_method]) +% [rgb, XYZ] = cieSpectralToColor(lambda_C, C, lambda_R, R [, whitepoint, int_method]) % % ## Description -% rgb = cieSpectralToColor(lambda_C, C, lambda_R, R [, whitepoint]) +% rgb = cieSpectralToColor(lambda_C, C, lambda_R, R [, whitepoint, int_method]) % Returns sRGB values corresponding to the input spectral radiances % -% [rgb, XYZ] = cieSpectralToColor(lambda_C, C, lambda_R, R [, whitepoint]) +% [rgb, XYZ] = cieSpectralToColor(lambda_C, C, lambda_R, R [, whitepoint, int_method]) % Additionally returns the CIE 1931 tristimulus values corresponding to % the input spectral radiances % @@ -42,7 +42,15 @@ % whitepoint -- Illuminant whitepoint % A character vector describing the CIE standard illuminant with which % the radiances were obtained. `whitepoint` is passed to the MATLAB -% 'xyz2rgb()' function, and defaults to 'd65' if not passed. +% 'xyz2rgb()' function, and defaults to 'd65' if empty or if not passed. +% +% int_method -- Numerical integration method +% The numerical integration method to use when integrating over the +% tristimulus functions. `int_method` is passed to `integrationWeights()` +% as its `method` input argument. A value of 'rect' is required for this +% function to operate according to the ASTM E308 standard. +% +% Defaults to 'rect' if not passed. % % ## Output Arguments % @@ -73,9 +81,12 @@ % on June 5, 2018. % - Lindbloom, Bruce J. (2017). Computing XYZ From Spectral Data. Retrieved % from http://www.brucelindbloom.com on June 11, 2018. +% - ASTM E308-17 Standard Practice for Computing the Colors of Objects by +% Using the CIE System, ASTM International, West Conshohocken, PA, 2017, +% https://doi.org/10.1520/E0308-17 % % See also xyz2rgb, resampleArrays, reflectanceToColor, -% reflectanceToRadiance +% reflectanceToRadiance, integrationWeights % Bernard Llanos % Supervised by Dr. Y.H. Yang @@ -83,11 +94,17 @@ % File created June 7, 2018 nargoutchk(1, 2); -narginchk(4, 5); +narginchk(4, 6); +whitepoint = []; +int_method = 'rect'; if ~isempty(varargin) whitepoint = varargin{1}; -else + if length(varargin) > 1 + int_method = varargin{2}; + end +end +if isempty(whitepoint) whitepoint = 'd65'; end @@ -98,13 +115,9 @@ R_resampled = R_resampled.'; lambda = reshape(lambda, length(lambda), 1); -lambda_diff = diff(lambda); -lambda_diff = [ - lambda_diff; - lambda_diff(end) - ]; +weights = repmat(integrationWeights(lambda, int_method).', size(R, 2), 1); -XYZ = R_resampled * diag(lambda_diff) * C_resampled; +XYZ = (R_resampled .* weights) * C_resampled; XYZ(XYZ < 0) = 0; XYZ(XYZ > 1) = 1; diff --git a/aberration_data/reflectanceToColor.m b/aberration_data/reflectanceToColor.m index af8c67d..f2bf60a 100644 --- a/aberration_data/reflectanceToColor.m +++ b/aberration_data/reflectanceToColor.m @@ -3,21 +3,21 @@ % % ## Syntax % rgb = reflectanceToColor(... -% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, whitepoint]... +% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, whitepoint, int_method]... % ) % [rgb, XYZ] = reflectanceToColor(... -% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, whitepoint]... +% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, whitepoint, int_method]... % ) % % ## Description % rgb = reflectanceToColor(... -% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, whitepoint]... +% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, whitepoint, int_method]... % ) % Returns rgb values for the given reflectances seen under the given % illuminant % % [rgb, XYZ] = reflectanceToColor(... -% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, whitepoint]... +% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, whitepoint, int_method]... % ) % Additionally returns CIE 1931 tristimulus values for the given % reflectances seen under the given illuminant @@ -47,8 +47,9 @@ % % lambda_C -- Reference wavelength values % A vector of wavelengths at which the 'C(lambda)' distributions were -% sampled. An even sampling of 5 nm is required, spanning the range from -% 360 nm (inclusive) to 780 nm (inclusive). +% sampled. An even sampling of 5 nm, spanning the range from 360 nm +% (inclusive) to 780 nm (inclusive), is required for this function to +% conform to the ASTM E308 standard. % % C -- CIE 1931 color matching functions % A 3-column matrix, where `C(i, j)` contains the value of the j-th CIE @@ -59,7 +60,15 @@ % A character vector naming the CIE standard illuminant corresponding to % `L`. The 'cieSpectralToColor()' function will use a default value of % 'd65' if its `whitepoint` input argument is not passed by this -% function. +% function. `whitepoint` can be empty (`[]`). +% +% int_method -- Numerical integration method +% The numerical integration method to use when integrating over the +% tristimulus functions. `int_method` is passed to `integrationWeights()` +% as its `method` input argument. A value of 'rect' is required for this +% function to operate according to the ASTM E308 standard. +% +% Defaults to 'rect' if not passed. % % ## Output Arguments % @@ -83,7 +92,7 @@ % https://doi.org/10.1520/E0308-17 % % See also reflectanceToRadiance, ciedIlluminant, cieSpectralToColor, -% resampleArrays +% resampleArrays, integrationWeights % Bernard Llanos % Supervised by Dr. Y.H. Yang @@ -91,22 +100,28 @@ % File created June 12, 2018 nargoutchk(1, 2); -narginchk(6, 7); +narginchk(6, 8); + +if length(varargin) > 1 + int_method = varargin{2}; +else + int_method = 'rect'; +end [lambda_resampled, ~, Rad_normalized, lambda_C_resampled, C_resampled] = reflectanceToRadiance(... - lambda_L, L, lambda_Ref, Ref, lambda_C, C, 2 ... + lambda_L, L, lambda_Ref, Ref, lambda_C, C, 2, int_method... ); if isempty(varargin) [rgb, XYZ] = cieSpectralToColor(... lambda_resampled, C_resampled,... - lambda_resampled, Rad_normalized.'... + lambda_resampled, Rad_normalized.', [], int_method... ); else [rgb, XYZ] = cieSpectralToColor(... lambda_C_resampled, C_resampled,... lambda_resampled, Rad_normalized.',... - varargin{1}... + varargin{1}, int_method... ); end diff --git a/aberration_data/reflectanceToRadiance.m b/aberration_data/reflectanceToRadiance.m index 89fca39..8d2c3e3 100644 --- a/aberration_data/reflectanceToRadiance.m +++ b/aberration_data/reflectanceToRadiance.m @@ -6,12 +6,12 @@ % lambda_L, L, lambda_Ref, Ref... % ) % [lambda_Rad, Rad, Rad_normalized] = reflectanceToRadiance(... -% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, C_ind]... +% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, C_ind, int_method]... % ) % [... % lambda_Rad, Rad, Rad_normalized, lambda_C_resampled, C_resampled... % ] = reflectanceToRadiance(... -% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, C_ind]... +% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, C_ind, int_method]... % ) % % ## Description @@ -22,7 +22,7 @@ % illuminant % % [lambda_Rad, Rad, Rad_normalized] = reflectanceToRadiance(... -% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, C_ind]... +% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, C_ind, int_method]... % ) % Additionally returns radiances normalized by the sensor's response to % the illuminant. @@ -30,7 +30,7 @@ % [... % lambda_Rad, Rad, Rad_normalized, lambda_C_resampled, C_resampled... % ] = reflectanceToRadiance(... -% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, C_ind]... +% lambda_L, L, lambda_Ref, Ref, lambda_C, C [, C_ind, int_method]... % ) % Additionally returns the version of the sensor's response functions % which can be used for colour calculation, along with the wavelengths at @@ -80,6 +80,13 @@ % % Defaults to 1 if not passed. % +% int_method -- Numerical integration method +% The numerical integration method to use when computing the radiance +% normalization constant. `int_method` is passed to +% `integrationWeights()` as its `method` input argument. +% +% Defaults to 'rect' if not passed. +% % ## Output Arguments % % lambda_Rad -- Radiance wavelength values @@ -121,7 +128,7 @@ % https://doi.org/10.1520/E0308-17 % % See also ciedIlluminant, reflectanceToColor, cieSpectralToColor, -% resampleArrays +% resampleArrays, integrationWeights % Bernard Llanos % Supervised by Dr. Y.H. Yang @@ -164,22 +171,27 @@ nargoutchk(2, 2); narginchk(4, 4); normalize = false; -elseif length(varargin) == 2 || length(varargin) == 3 +elseif length(varargin) == 2 || length(varargin) == 3 || length(varargin) == 4 nargoutchk(3, 5); if nargout == 4 error('Either three or five output arguments must be requested.'); end - narginchk(6, 7); + narginchk(6, 8); normalize = true; lambda_C = varargin{1}; C = varargin{2}; + int_method = 'rect'; if length(varargin) > 2 C_ind = varargin{3}; + + if length(varargin) > 3 + int_method = varargin{4}; + end else C_ind = 1; end else - error('Either zero, two, or three optional input arguments must be provided, not %d.', length(varargin)); + error('Either zero, two, three, or four optional input arguments must be provided, not %d.', length(varargin)); end [... @@ -201,16 +213,11 @@ lambda_L, L, lambda_C, ybar,... 'spline'... ); - lambda_diff = diff(lambda_resampled2); - lambda_diff = [ - lambda_diff; - lambda_diff(end) - ]; - + weights = integrationWeights(lambda_resampled2, int_method); ybar_resampled = sumEnds(lambda_C, lambda_resampled2, ybar, ybar_resampled); % Compute the normalization constant - N = sum(L_resampled2 .* ybar_resampled .* lambda_diff); + N = dot(L_resampled2 .* ybar_resampled, weights); Rad_normalized = Rad ./ N; diff --git a/dispersion_model/makePolyfun.m b/dispersion_model/makePolyfun.m index 8728b85..4f65ddf 100644 --- a/dispersion_model/makePolyfun.m +++ b/dispersion_model/makePolyfun.m @@ -97,13 +97,18 @@ xy_normalized(:, 1:(end - 1)),... lambda_normalized(:, 1:(end - 1)),... ]; - xylambda_normalized_3d = repmat(permute(xylambda_normalized, [1 3 2]), 1, polyfun_data(d).n_powers, 1); - powers_rep = repmat(polyfun_data(d).powers, n_d, 1, 1); - vandermonde_matrix = prod(xylambda_normalized_3d .^ powers_rep, 3); - - disparity_x_normalized = vandermonde_matrix * polyfun_data(d).coeff_x; - disparity_y_normalized = vandermonde_matrix * polyfun_data(d).coeff_y; - disparity_normalized = [disparity_x_normalized disparity_y_normalized, ones(n_d, 1)]; + xylambda_normalized_3d = permute(xylambda_normalized, [1 3 2]); + + disparity_normalized = ones(n_d, 3); + for j = 1:n_d + vandermonde_vector = prod(... + repmat(xylambda_normalized_3d(j, :, :), 1, polyfun_data(d).n_powers, 1)... + .^ polyfun_data(d).powers, 3 ... + ); + disparity_normalized(j, 1) = dot(vandermonde_vector, polyfun_data(d).coeff_x); + disparity_normalized(j, 2) = dot(vandermonde_vector, polyfun_data(d).coeff_y); + end + disparity_d = (polyfun_data(d).T_disparity_inv * disparity_normalized.').'; disparity(filter_d, :) = disparity_d(:, 1:(end-1)); end diff --git a/integrationWeights.m b/integrationWeights.m new file mode 100644 index 0000000..e2dbe8b --- /dev/null +++ b/integrationWeights.m @@ -0,0 +1,64 @@ +function [w] = integrationWeights(x, method) +% INTEGRATIONWEIGHTS Calculate weights for numerical integration +% +% ## Syntax +% w = integrationWeights(x, method) +% +% ## Description +% w = integrationWeights(x, method) +% Returns weights to use in a weighted sum for numerical integration +% of a function, given the sampling locations in its 1D domain +% +% ## Input Arguments +% +% x -- Sampling locations +% A vector of values at which the function to be integrated is sampled. +% +% method -- Integration method +% A character vector or a string scalar indicating which numerical +% integration method to use. If `method` is: +% - 'rect', then weights will be produced according to a rectangle rule +% with the value of the function within each subinterval taken as its +% value on the left end of the subinterval. +% - 'trap', then weights will be produced according to the trapezoid rule +% +% ## Output Arguments +% +% w -- Integration weights +% A column vector with the same length as `x`, containing weights for +% numerical integration. A function, represented as a vector, `f`, of +% values sampled at the locations in `x`, has an integral in the interval +% [x(1), x(end)] approximated by `dot(f, w)`. +% +% See also trapz + +% Bernard Llanos +% Supervised by Dr. Y.H. Yang +% University of Alberta, Department of Computing Science +% File created June 28, 2018 + +nargoutchk(1, 1); +narginchk(2, 2); + +if ~isvector(x) + error('The domain of integration must be one-dimensional.'); +end + +if size(x, 1) < size(x, 2) + x = x.'; +end + +if isStringScalar(method) || ischar(method) + dx = diff(x); + if strcmp(method, 'rect') + w = [dx; 0]; + elseif strcmp(method, 'trap') + w = ([dx; 0] + [0; dx] ) / 2; + else + error('Unexpected value of `method`.'); + end +else + error('`method` must be a character vector or a string scalar.'); +end + +end \ No newline at end of file diff --git a/sensor/CIE1931ColorMap.m b/sensor/CIE1931ColorMap.m index b76b0cb..8d4d71e 100644 --- a/sensor/CIE1931ColorMap.m +++ b/sensor/CIE1931ColorMap.m @@ -19,6 +19,9 @@ % ### Sensor quantum efficiency data % % A '.mat' file containing the following variables: +% - 'channel_mode': A Boolean value set to `false` to indicate that the +% data in `sensor_map` represents spectral sensitivities, not colour +% channel mappings % - 'sensor_map': A 2D array, where `sensor_map(i, j)` is the sensitivity % of the i-th tristimulus function (x-bar, y-bar or z-bar) to light of % the j-th wavelength. @@ -37,12 +40,15 @@ % List of parameters to save with results parameters_list = { - 'data_source'... + 'data_source',... + 'channel_mode'... }; %% Input data and parameters data_source = '/home/llanos/GoogleDrive/ThesisResearch/Data and Results/20180614_ASTM_E308/Table1_CIE1931_2DegStandardObserver.csv'; +% Data represents spectral sensitivities, not colour channel mappings +channel_mode = false; %% Load the data @@ -64,6 +70,7 @@ legend('x-bar', 'y-bar', 'z-bar') %% Save to a file +sensor_map = sensor_map.'; save_variables_list = [ parameters_list, {... 'sensor_map',... 'bands'... diff --git a/sensor/RGBColorMap.m b/sensor/RGBColorMap.m index 9ca5d30..71b5f41 100644 --- a/sensor/RGBColorMap.m +++ b/sensor/RGBColorMap.m @@ -15,6 +15,9 @@ % ### Colour conversion data % % A '.mat' file containing the following variables: +% - 'channel_mode': A Boolean value set to `true` to indicate that the +% data in `sensor_map` represents colour channel mappings, not spectral +% sensitivities. % - 'sensor_map': A 2D array, where `sensor_map(i, j)` is the sensitivity % of the i-th colour channel (Red, Green, or Blue) to light of the j-th % colour channel (Red, Green, or Blue). @@ -28,9 +31,11 @@ bands = (1:3).'; sensor_map = eye(length(bands)); +channel_mode = true; save_variables_list = {... 'bands',... + 'channel_mode',... 'sensor_map'... }; uisave(save_variables_list,'RGBColorMapData'); \ No newline at end of file diff --git a/sensor/SonyColorMap.m b/sensor/SonyColorMap.m index 3e9d723..497b4ee 100644 --- a/sensor/SonyColorMap.m +++ b/sensor/SonyColorMap.m @@ -27,6 +27,9 @@ % ### Sensor quantum efficiency data % % A '.mat' file containing the following variables: +% - 'channel_mode': A Boolean value set to `false` to indicate that the +% data in `sensor_map` represents spectral sensitivities, not colour +% channel mappings % - 'sensor_map': A 2D array, where `sensor_map(i, j)` is the sensitivity % of the i-th colour channel (Red, Green, or Blue) to light of the j-th % wavelength. @@ -45,6 +48,7 @@ % List of parameters to save with results parameters_list = { + 'channel_mode',... 'data_source',... 'bands'... }; @@ -52,7 +56,8 @@ %% Input data and parameters fn_name = 'sonyQuantumEfficiency'; -bands = linspace(200, 1200, 1000); +bands = linspace(200, 1200, 1000).'; +channel_mode = false; %% Load the data