Demo data and documentation.

- Squash commit including some minor bug fixes.

Contents.m

@@ -0,0 +1,10 @@
+% Simultaneous Demosaicing and Chromatic Aberration Correction through Spectral Reconstruction
+% Version 2.0.0 18-Jul-2019
+%
+% Refer to README.md for setup instructions and high-level documentation.
+%
+% Subdirectories are documented with their own 'Contents.m' files, and are also
+% described in README.md
+%
+% Files
+%   SetFixedParameters - Set parameters used by multiple scripts in the codebase.

README.md

@@ -22,6 +22,12 @@
   on image datasets
 - Scripts for visualizing and analyzing the results
 
+## Citation
+- If you use this work, you may cite the following thesis, available at
+  https://era.library.ualberta.ca/
+
+B. Llanos. "Chromatic Aberration Correction and Spectral Reconstruction from Colour Images." M.Sc. thesis, University of Alberta, Canada, 2019.
+
 ## Environment
 - MATLAB R2018b was used for development, but other MATLAB releases will likely
   work.
@@ -46,7 +52,7 @@
 - A full list of dependencies for individual files can be determined as
   described at https://www.mathworks.com/help/matlab/matlab_prog/identify-dependencies.html
 
-### Third-party dependencies
+### Third-party code dependencies
 - All third-party dependencies are either optional, or are easily ommitted
   through minor modifications to the codebase.
 - Third-party dependencies should be added to the MATLAB path, after being
@@ -100,15 +106,57 @@
 - Optional dependency, used by [evaluation/RunOnDataset.m]() only when enabled
   in [evaluation/SetAlgorithms.m](). Also used by [data_analysis/TuneSunEtAl2017.m]().
 
+### Third-party data dependencies
+- All demos in [demo_data/README.md]() do not require third-party data.
+
+#### CIE 1931 Standard (2-degree) Observer spectral tristimulus functions
+- These functions were retrieved from Table 1 of the ASTM E308 standard
+  (cited below) and saved as a CSV file.
+- In the CSV file, the first column contained wavelength values, and the
+  following three columns contained the values of the 'x-bar', 'y-bar',
+  and 'z-bar' functions, respectively.
+- The data is loaded into the `xyzbar` variable in many scripts.
+
+#### Spectral Power Distribution of a CIE D-Illuminant
+- Data used for synthesizing the spectral radiance of a CIE D-illuminant
+  for an arbitrary correlated colour temperature was retrieved from Bruce Lindbloom's
+  website (cited below).
+- The original spreadsheet file published by Bruce Lindbloom was converted
+  to a CSV file, where the first column contained wavelength values, and the
+  following three columns contained the values of the 'S0', 'S1', and 'S2'
+  functions, respectively.
+- The data is loaded by several scripts as input for the [aberration_data/ciedIlluminant.m]()
+  function.
+
+#### ColorChecker spectral reflectances
+- Spectral reflectance measurements of the patches of the average X-Rite ColorChecker
+  chart from Danny Pascale's website (cited below) were used before in-house
+  spectral reflectance measurements were available
+  (saved as [demo_data/spectral_data/spectra_averaged.csv]()).
+- The reflectance data was stored in a CSV file, where the first column contained
+  wavelength values, and the following 24 columns contained the spectral reflectances
+  of the 24 ColorChecker CLASSIC patches.
+- It would be easy to create an equivalent file from
+  [demo_data/spectral_data/spectra_averaged.csv]() to use instead in most places.
+- Presently, only [aberration_data/BimaterialImages.m](),
+  [data_analysis/ColorCheckerClassicPreNov2014.m](), and
+  [data_analysis/SpectralDataJune262018.m]() rely on this data.
+
 ## Setup
 - Add all folders and subfolders to the MATLAB path, excluding [.git/]() and [deprecated/]().
 
 ## Basic usage instructions
-- Files starting with capital letters are MATLAB scripts, which can be used
-  as described in their documentation comments.
+- To run predefined demos, see [demo_data/README.md]().
+- Files with names starting with capital letters are MATLAB scripts, and can
+  be used as described in their documentation comments.
 - Remaining files are MATLAB functions called by the scripts.
   There are a few MATLAB functions used for ad-hoc analysis for which no scripts
   have yet been created to call them.
+- The tokens `${DIRPATH}` and `${FILEPATH}` that appear in character vector
+  literals need to be replaced with directory paths and file paths, respectively,
+  that are specific to your data file locations. Some directory and file paths
+  are presently set to point to the demo data in [./demo_data/](), but can be
+  changed to point to your data.
 
 ## Tips
 - Raw colour-filter array images as used in this codebase are truly raw. No gamma
@@ -119,6 +167,11 @@
 - Running more parallel workers will increase memory consumption, but not
   excessively. The code has been designed to avoid having MATLAB send
   redundant data to parallel workers.
+- If you want to run the code on your own dataset of images, adjust parameters
+  in [SetFixedParameters.m](). The following parameters are most relevant:
+  - `bayer_pattern`
+  - `findSamplingOptions`
+  - `dispersionfunToMatrixOptions.resolution`
 
 ### Troubleshooting
 
@@ -228,6 +281,8 @@
     of the output of other code in the repository. In particular, there are
     scripts for creating colour-corrected images. There are also scripts created
     for tangential or one-time experiments.
+  - 'demo_data': Sample data for testing the codebase, and sample results on the
+    demo data.
   - 'deprecated': Unmaintained code, that is no longer called by the programs
     in other folders. Code in this folder is present for archival purposes only.
   - 'disk_fitting': Disk keypoint-based calibration of lateral chromatic aberration,
@@ -328,6 +383,10 @@ and applications,” IEEE Transactions on Pattern Analysis and Ma-
 chine Intelligence, vol. 41, no. 1, pp. 20–33, 2019.
 doi: 10.1109/TPAMI.2017.2760833.
 
+FLIR. (Jan. 27, 2017). FLIR FLEA3 GigE Vision Imaging Performance
+Specification. version 1.1, [Online]. Available:
+https://www.ptgrey.com/support/downloads/10109 (visited on 05/08/2017).
+
 D.H. Foster. (2018). Tutorial on Transforming Hyperspectral Images to
 RGB Colour Images, [Online].
 Available: http://personalpages.manchester.ac.uk/staff/d.h.foster/Tutorial_HSI2RGB/Tutorial_HSI2RGB.html

SetFixedParameters.m

@@ -1,5 +1,5 @@
 %% Set Fixed Parameters
-% Set values of parameters for image correction that seldomly need to be
+% Set values of parameters, common to multiple scripts, that seldomly need to be
 % changed.
 %
 % ## Usage
@@ -144,13 +144,13 @@
 % 'findSampling()' should not even be called.)
 findSamplingOptions.int_method = 'trap';
 
-findSamplingOptions.power_threshold = 0.95;
+findSamplingOptions.power_threshold = 1;
 % As an alternative to automatically determining the number of spectral
 % bands, according to `findSamplingOptions.power_threshold`, set it
 % explicitly (if the following option is an integer greater than zero).
-findSamplingOptions.n_bands = 8;
+findSamplingOptions.n_bands = 0;
 
-findSamplingOptions.support_threshold = 0.05;
+findSamplingOptions.support_threshold = 0;
 
 findSamplingOptions.bands_padding = 1000;
 
@@ -161,7 +161,7 @@
 % between samples in the sequence of samples being interpolated. The
 % interpolation function `f(x)` returns the weight for a sample at location
 % `x` relative to the current  interpolation location.
-findSamplingOptions.interpolant = @gaussian;
+findSamplingOptions.interpolant = @triangle;
 
 % Interpolation function for other spectral data, such as sensor spectral
 % sensitivities or ground truth spectral radiances. In contrast to
@@ -176,7 +176,7 @@
 % details.
 
 % Resolution at which to sample spectral dispersion
-dispersionfunToMatrixOptions.resolution = 0.1; % pixels
+dispersionfunToMatrixOptions.resolution = 0; % pixels
 
 dispersionfunToMatrixOptions.int_method = findSamplingOptions.int_method;
 dispersionfunToMatrixOptions.support_threshold = findSamplingOptions.support_threshold;
@@ -265,9 +265,9 @@
 % Only use even integers for the patch and padding sizes, to ensure that patches
 % are valid colour filter array images.
 patch_sizes = [ % Each row contains a (number of rows, number of columns) pair
-   128, 128;
+   64, 64;
 ]; 
-paddings = 16;
+paddings = 8;
 
 solvePatchesColorOptions.patch_options = struct;
 solvePatchesColorOptions.patch_options.patch_size = patch_sizes(1, :);

aberration_correction/Contents.m

@@ -0,0 +1,19 @@
+% ABERRATION_CORRECTION
+% Version 2.0.0 18-Jul-2019
+%
+% Spectral or colour image reconstruction and chromatic aberration correction.
+%
+% Image reconstruction or correction scripts
+%   CorrectByHyperspectralADMM - Demosaicing and hyperspectral ADMM-based correction of chromatic aberration
+%   CorrectByWarping           - Correction of chromatic aberration by image warping
+%
+% Image transformation
+%   channelConversionMatrix    - Create a sparse matrix to convert images between colour spaces
+%   mosaicMatrix               - Create a sparse matrix to mosaic an image
+%   warpImageSpectral          - Patch-wise warping and conversion to colour of a spectral image
+%
+% Regularization operators for optimization problems
+%   spatialGradient            - Create sparse matrices acting as image spatial gradient operators
+%   spatialLaplacian           - Create a sparse matrix acting as an image spatial Laplacian operator
+%   spectralGradient           - Create a sparse matrix acting as an image spectral gradient operator
+%   TestMatrices               - Test script for matrix generation functions

aberration_correction/CorrectByHyperspectralADMM.m

@@ -183,33 +183,33 @@
 
 % Wildcard for 'ls()' to find the images to process.
 % '.mat' or image files can be loaded
-input_images_wildcard = '/home/llanos/GoogleDrive/ThesisResearch/Results/20190208_ComputarLens/dataset/exposure_blending/d2_colorChecker30cm_unfiltered.mat';
+input_images_wildcard = fullfile('.', 'demo_data', 'multispectral_images', 'd1_disks32cmV2_raw.mat');
 input_images_variable_name = 'I_raw'; % Used only when loading '.mat' files
 
 % Model of dispersion
 % Can be empty
-reverse_dispersion_model_filename = '/home/llanos/GoogleDrive/ThesisResearch/Results/20190208_ComputarLens/dispersion/spectral/full_image/RAWDiskDispersionResults_spectral_polynomial_fromNonReference.mat';
+reverse_dispersion_model_filename = fullfile('.', 'demo_data', 'dispersion_models', 'disk_fitting', 'RAWDiskDispersionResults_spectral_polynomial_fromNonReference.mat');
 
 % Colour space conversion data
-color_map_filename = '/home/llanos/GoogleDrive/ThesisResearch/Results/20190208_ComputarLens/dataset/SonyColorMapData.mat';
+color_map_filename = fullfile('.', 'demo_data', 'multispectral_images', 'sensor.mat');
 
 % Output directory for all images and saved parameters
-output_directory = '/home/llanos/Downloads';
+output_directory = fullfile('.', 'demo_data', 'aberration_correction');
 
 % ## Options for patch-wise image estimation
 
 % Only estimate a single patch, with its top-left corner at the given (row,
 % column) location. If empty (`[]`), the entire image will be estimated. The
 % patch corner indices must be odd integers to avoid creating a patch with a
 % different colour filter array pattern from the whole image.
-target_patch = [1405, 271];
+target_patch = [];
 
 % Only select regularization weights for a single patch, with its top-left
 % corner at the given (row, column) location. If empty (`[]`), regularization
 % weights will be selected for each patch separately. (THIS IS SLOW) Again, the
 % patch corner indices must be odd integers to avoid creating a patch with a
 % different colour filter array pattern from the whole image.
-target_patch_weights = [1405, 271];
+target_patch_weights = [99, 99];
 
 % Also compare with (or only run) whole image estimation, meaning that the image
 % is treated as a single patch. Only enable this for small images.

aberration_correction/CorrectByWarping.m

@@ -132,16 +132,16 @@
 
 % Wildcard for 'ls()' to find the images to process.
 % '.mat' or image files can be loaded
-input_images_wildcard = '/home/llanos/GoogleDrive/ThesisResearch/Results/20190421_ComputarLens_revisedAlgorithms/channel_scaling/*_dHyper.mat';
-input_images_variable_name = 'I_hyper'; % Used only when loading '.mat' files
+input_images_wildcard = fullfile('.', 'demo_data', 'multispectral_images', 'd1_disks32cmV2_d3.mat');
+input_images_variable_name = 'I_3'; % Used only when loading '.mat' files
 
 % Whether the input images are spectral images or colour images
-is_spectral = true;
+is_spectral = false;
 
 % Colour space conversion data
 convert_to_color = false; % Whether or not to convert images to a different colour space
 if convert_to_color
-    color_map_filename = '/home/llanos/GoogleDrive/ThesisResearch/Results/20190421_ComputarLens_revisedAlgorithms/channel_scaling/sensor.mat';
+    color_map_filename = '${FILEPATH}';
 else
     color_map_filename = []; % Unused
 end
@@ -153,7 +153,7 @@
     % The file must contain a vector, `bands`, of wavelengths corresponding to the
     % spectral bands of the spectral images. It must also contain a variable,
     % `findSamplingOptions`, such as defined in 'SetFixedParameters.m'.
-    sampling_filename = '/home/llanos/GoogleDrive/ThesisResearch/Results/20190421_ComputarLens_revisedAlgorithms/run_on_dataset_dHyper_dispersion/RunOnDataset_20190421_ComputarLens_dHyper_dispersion.mat';
+    sampling_filename = fullfile('.', 'demo_data', 'multispectral_images', 'sensor.mat');
 
     % Whether or not to prevent spectral resampling and interpolation during
     % warping. When `disable_interpolation` is `true`, the `options` input
@@ -167,10 +167,10 @@
 end
 
 % Model of dispersion
-forward_dispersion_model_filename = '/home/llanos/GoogleDrive/ThesisResearch/Results/20190208_ComputarLens/dispersion/spectral/full_image/RAWDiskDispersionResults_spectral_polynomial_fromReference.mat';
+forward_dispersion_model_filename = fullfile('.', 'demo_data', 'dispersion_models', 'disk_fitting', 'RAWDiskDispersionResults_RGB_polynomial_fromReference.mat');
 
 % Output directory for all images and saved parameters
-output_directory = '/home/llanos/GoogleDrive/ThesisResearch/Results/20190421_ComputarLens_revisedAlgorithms/warping_filtered_images/bandpassFiltered_warpCorrected_asSpectralImages';
+output_directory = fullfile('.', 'demo_data', 'aberration_correction');
 
 % Parameters which do not usually need to be changed. Some of these
 % parameters will be overridden by the input data.
@@ -234,8 +234,12 @@
         'interpolant_ref', findSamplingOptions.interpolant_ref,...
         'bands_in', bands...
     );
-else
+elseif is_spectral
     dispersion_options = struct('bands_in', bands);
+else
+    n_channels_rgb = 3;
+    bands_rgb = 1:n_channels_rgb;
+    dispersion_options = struct('bands_in', bands_rgb);
 end
 
 %% Process the images

aberration_correction/admm/Contents.m

@@ -0,0 +1,16 @@
+% ADMM
+% Version 2.0.0 18-Jul-2019
+%
+% Image reconstruction optimization problem setup and resolution.
+%
+% Alternating Direction Method of Multipliers
+%   baek2017Algorithm2LowMemory     - Run ADMM (loosely) as in Algorithm 2 of Baek et al. 2017
+%   initBaek2017Algorithm2LowMemory - Allocate memory for 'baek2017Algorithm2LowMemory()'
+%   solvePatchesColor               - Run ADMM (loosely) as in Algorithm 2 of Baek et al. 2017, with weight selection and patch-wise decomposition
+%   solvePatchesSpectral            - Run ADMM for spectral image estimation
+%
+% Regularization weight selection for the optimization problem
+%   initPenalties                   - Allocate memory for 'penalties()'
+%   initWeightsLowMemory            - Allocate memory for 'weightsLowMemory()'
+%   penalties                       - Calculate data fitting and regularization errors
+%   weightsLowMemory                - Select regularization weights using a grid search

aberration_correction/patch/Contents.m

@@ -0,0 +1,7 @@
+% PATCH
+% Version 2.0.0 18-Jul-2019
+%
+% Helper code for image splitting into patches.
+%
+% Files
+%   patchBoundaries - Find the inner and outer boundaries of an image patch

aberration_data/BimaterialImages.m

@@ -212,56 +212,56 @@
 %% Input data and parameters
 
 % Wildcard for 'ls()' to find the chromaticity maps.
-input_chromaticity_maps_wildcard = '/home/llanos/GoogleDrive/ThesisResearch/Results/20181127_ColorCheckerReflectanceImages/colorIDs.png';
+input_chromaticity_maps_wildcard = fullfile('.', 'demo_data', 'hdr_averaged_images', 'd1_disks32cmV2_550nm_rgb01.tif');
 
 % Threshold for assuming that an RGB image is actually a greyscale image
 grey_difference_threshold = 1;
 
 % Treat greyscale images as colour class IDs rather than pictures. An error
 % will be thrown if the images are not greyscale, and this option is
 % `true`.
-is_classes = true;
+is_classes = false;
 
 % Wildcard for 'ls()' to find the shading maps.
 % Set to an empty array to use constant shading
-input_shading_maps_wildcard = [];
+input_shading_maps_wildcard = fullfile('.', 'demo_data', 'hdr_averaged_images', 'd1_disks32cmV2_550nm_rgb01.tif');
 
 % Model of dispersion (can be empty, for no dispersion)
-reverse_dispersion_model_filename = []; %'/home/llanos/GoogleDrive/ThesisResearch/Results/20181020_DoubleConvexThickLensDispersion_Final/Models/DoubleConvexThickLensDispersionResults_spectral_spline_fromNonReference.mat';
+reverse_dispersion_model_filename = fullfile('.', 'demo_data', 'dispersion_models', 'disk_fitting', 'RAWDiskDispersionResults_spectral_polynomial_fromNonReference.mat');
 
 % Colour space conversion data
-color_map_filename = '/home/llanos/GoogleDrive/ThesisResearch/Results/20181127_TestingChoiEtAl2017/NikonD5100ColorMapData.mat';
+color_map_filename = fullfile('.', 'demo_data', 'multispectral_images', 'sensor.mat');
 % Colour channel to use for radiance normalization
 normalization_channel = 2;
 
 use_cie_illuminant = true;
 if use_cie_illuminant
     % CIE D-illuminant
-    illuminant_filename = '/home/llanos/GoogleDrive/ThesisResearch/Data/20180604_Spectral power distributions_BruceLindbloom/DIlluminants.csv';
+    illuminant_filename = '${FILEPATH}';
     illuminant_temperature = 6504; % From https://en.wikipedia.org/wiki/Standard_illuminant#Illuminant_series_D
     illuminant_function_name = []; % Unused
 else
     % Arbitrary illuminant function
-    illuminant_function_name = 'one';
+    illuminant_function_name = [];
     illuminant_filename = []; % Unused
     illuminant_temperature = []; % Unused
 end
 
 % CIE tristimulus functions
-xyzbar_filename = '/home/llanos/GoogleDrive/ThesisResearch/Data/20180614_ASTM_E308/Table1_CIE1931_2DegStandardObserver.csv';
+xyzbar_filename = '${FILEPATH}';
 
 % ColorChecker spectral reflectances
-reflectances_filename = '/home/llanos/GoogleDrive/ThesisResearch/Data/20180604_ColorCheckerSpectralData_BabelColor/ColorChecker_spectra_reformatted_llanos.csv';
+reflectances_filename = '${FILEPATH}';
 
 % Number of colours to blend in each output image (overridden by the actual
 % number of colour classes in the image, if `is_classes` is `true`)
 n_colors = 2;
 
 % Output directory for all images and saved data
-output_directory = '/home/llanos/Downloads';
+output_directory = '${DIRPATH}';
 
 % ## Debugging Flags
-segmentColorsVerbose = false;
+segmentColorsVerbose = true;
 
 % Parameters which do not usually need to be changed
 run('SetFixedParameters.m')