diff --git a/aberration_correction/CorrectByHyperspectralADMM.m b/aberration_correction/CorrectByHyperspectralADMM.m
index 8534991..7a38402 100644
--- a/aberration_correction/CorrectByHyperspectralADMM.m
+++ b/aberration_correction/CorrectByHyperspectralADMM.m
@@ -232,14 +232,7 @@
     ] = loadDispersionModel(reverse_dispersion_model_filename, false);
 end
 
-model_variables_required = { 'sensor_map', 'channel_mode', 'bands' };
-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
-
-bands_color = bands;
-
+[sensor_map, channel_mode, bands_color] = loadColorMap(color_map_filename);
 if channel_mode
     if has_dispersion && ...
        ((length(bands_color) ~= length(bands_dispersionfun)) ||...
diff --git a/aberration_correction/CorrectByWarping.m b/aberration_correction/CorrectByWarping.m
index d239fb5..8c83052 100644
--- a/aberration_correction/CorrectByWarping.m
+++ b/aberration_correction/CorrectByWarping.m
@@ -169,30 +169,14 @@
     dispersion_data, bands_dispersionfun, transform_data...
 ] = loadDispersionModel(forward_dispersion_model_filename, true, true);
 
-model_variables_required = { 'sensor_map', 'channel_mode', 'bands' };
-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
-
-bands_color = bands;
-
+[sensor_map, channel_mode, bands_color] = loadColorMap(color_map_filename);
 if channel_mode
     if ((length(bands) ~= length(bands_dispersionfun)) ||...
        any(bands(:) ~= bands_dispersionfun(:)))
         error('When warping colour images, the same colour channels must be used by the model of dispersion.');
     end
 else
-    findSamplingOptions = [];
-    bands = [];
-    variables_required = { 'findSamplingOptions', 'bands' };
-    load(sampling_filename, variables_required{:});
-    if isempty(bands)
-        error('`bands` was not loaded from "%s".', sampling_filename)
-    end
-    if isempty(findSamplingOptions)
-        error('`findSamplingOptions` was not loaded from "%s".', sampling_filename)
-    end
+    [findSamplingOptions, bands] = loadVariables(sampling_filename, {'findSamplingOptions', 'bands'});
     if disable_interpolation
         color_weights = colorWeights(...
             sensor_map, bands_color, bands, findSamplingOptions...
diff --git a/aberration_data/BimaterialImages.m b/aberration_data/BimaterialImages.m
index f1fcbf4..ba375fe 100644
--- a/aberration_data/BimaterialImages.m
+++ b/aberration_data/BimaterialImages.m
@@ -57,7 +57,7 @@
 % 'SonyColorMap.m', for example. The following variables are required:
 % - '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.
+%   j-th element of 'bands' (below).
 % - 'channel_mode': A Boolean value indicating whether the input colour
 %   space is a set of colour channels (true) or a set of spectral bands
 %   (false). A value of `false` is required.
@@ -271,7 +271,8 @@
 % ColorChecker spectral reflectances
 colorChecker_table = readtable(reflectances_filename);
 variable_names = colorChecker_table.Properties.VariableNames;
-lambda_colorChecker = colorChecker_table.(variable_names{1});
+bands = colorChecker_table.(variable_names{1});
+n_bands = length(bands);
 reflectances = colorChecker_table{:, 2:end};
 n_patches = length(variable_names) - 1;
 
@@ -284,7 +285,7 @@
         illuminant_temperature, lambda_illuminant, S_illuminant, lambda_illuminant...
     );
 else
-    lambda_illuminant = lambda_colorChecker;
+    lambda_illuminant = bands;
     spd_illuminant = feval(illuminant_function_name, lambda_illuminant);
 end
 
@@ -295,7 +296,7 @@
 
 colorChecker_rgb = reflectanceToColor(...
     lambda_illuminant, spd_illuminant,...
-    lambda_colorChecker, reflectances,...
+    bands, reflectances,...
     lambda_xyzbar, xyzbar, findSamplingOptions.int_method...
     );
 
@@ -318,24 +319,13 @@
 
 %% Load colour space conversion data
 
-model_variables_required = { 'sensor_map', 'channel_mode', 'bands' };
-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
-if channel_mode
-    error('The input space of the colour conversion data must be a spectral space, not a space of colour channels.')
-end
-
-bands_color = bands;
-bands = lambda_colorChecker;
-n_bands = length(bands);
+[sensor_map, ~, bands_color] = loadColorMap(color_map_filename, false);
 
 %% Calculate spectral radiances
 
 [lambda_Rad, ~, Rad_normalized] = reflectanceToRadiance(...
     lambda_illuminant, spd_illuminant,...
-    lambda_colorChecker, reflectances,...
+    bands, reflectances,...
     bands_color, sensor_map.',...
     normalization_channel, findSamplingOptions.int_method...
 );
diff --git a/data_analysis/ColorCorrection.m b/data_analysis/ColorCorrection.m
index 0a3dc30..dd44e4b 100644
--- a/data_analysis/ColorCorrection.m
+++ b/data_analysis/ColorCorrection.m
@@ -73,19 +73,8 @@
 n_images = n_spectral + n_color;
 
 if n_spectral ~= 0
-    variables_required = { 'sensor_map', 'bands' };
-    load(color_map_filename, variables_required{:});
-    if ~all(ismember(variables_required, who))
-        error('One or more of the required colour space conversion variables is not loaded.')
-    end
-    bands_color = bands;
-    clear bands
-
-    variables_required = { 'findSamplingOptions', 'bands' };
-    load(sampling_filename, variables_required{:});
-    if ~all(ismember(variables_required, who))
-        error('One or more of the required sampling variables is not loaded.')
-    end
+    [sensor_map, ~, bands_color] = loadColorMap(color_map_filename);
+    [findSamplingOptions, bands] = loadVariables(sampling_filename, {'findSamplingOptions', 'bands'});
     color_weights = colorWeights(...
         sensor_map, bands_color, bands, findSamplingOptions...
     );
@@ -93,20 +82,10 @@
 
 if n_images ~= 0
     if use_chromadapt
-        variables_required = {'wb_illum', 'wb_scale'};
-        load(wb_filename, variables_required{:});
-        if ~all(ismember(variables_required, who))
-            error('One or more of the required white balancing variables is not loaded.')
-        end
+        [wb_illum, wb_scale] = loadVariables(wb_filename, {'wb_illum', 'wb_scale'});
         postfix = '_wb';
     else
-        load(xyz_weights_filename, xyz_weights_variable);
-        if exist(xyz_weights_variable, 'var')
-            xyz_weights = eval(xyz_weights_variable);
-        end
-        if ~exist(xyz_weights_variable, 'var') || isempty(xyz_weights)
-            error('No raw colour to XYZ conversion data loaded.')
-        end
+        xyz_weights = loadVariables(xyz_weights_filename, xyz_weights_variable);
         postfix = ['_' xyz_weights_variable];
     end
 end
diff --git a/data_analysis/DispersionPSFs.m b/data_analysis/DispersionPSFs.m
index 767e16b..bbc5981 100644
--- a/data_analysis/DispersionPSFs.m
+++ b/data_analysis/DispersionPSFs.m
@@ -107,15 +107,7 @@
 );
 df_color = makeDispersionfun(dispersion_data, T_roi);
 
-model_variables_required = { 'sensor_map', 'channel_mode', 'bands' };
-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
-if channel_mode
-    error('The input space of the colour conversion data must be a spectral space, not a space of colour channels.')
-end
-bands_color = bands;
+[sensor_map, ~, bands_color] = loadColorMap(color_map_filename, false);
 
 %% Load reflectance data
 
diff --git a/data_analysis/ImageEvaluation.m b/data_analysis/ImageEvaluation.m
index 005b6a5..96a7f47 100644
--- a/data_analysis/ImageEvaluation.m
+++ b/data_analysis/ImageEvaluation.m
@@ -79,13 +79,7 @@
 end
 
 if ~is_rgb
-    load(bands_filename, bands_variable);
-    if exist(bands_variable, 'var')
-        bands = eval(bands_variable);
-    end
-    if ~exist(bands_variable, 'var') || isempty(bands)
-        error('No wavelength band or colour channel index information loaded.')
-    end
+    bands = loadVariables(bands_filename, bands_variable);
 end
 
 if is_rgb
@@ -124,14 +118,7 @@
         if isempty(spectral_weights_filenames{i})
             spectral_weights = eye(length(bands));
         else
-            spectral_weights = [];
-            load(spectral_weights_filenames{i}, spectral_weights_variables{i});
-            if exist(spectral_weights_variables{i}, 'var')
-                spectral_weights = eval(spectral_weights_variables{i});
-            end
-            if ~exist(spectral_weights_variables{i}, 'var') || isempty(spectral_weights)
-                error('No colour space conversion information loaded for image %d.', i)
-            end
+            spectral_weights = loadVariables(spectral_weights_filenames{i}, spectral_weights_variables{i});
         end
         
         dp.evaluation.global_spectral.plot_color =...
diff --git a/data_analysis/IntraChannelDispersion.m b/data_analysis/IntraChannelDispersion.m
index 5d9b3d5..bdf22d5 100644
--- a/data_analysis/IntraChannelDispersion.m
+++ b/data_analysis/IntraChannelDispersion.m
@@ -34,14 +34,7 @@
 
 % Find the wavelengths bounding the high-sensitivity regions of the camera's
 % colour channels
-model_variables_required = { 'sensor_map', 'channel_mode', 'bands' };
-load(color_map_filename, model_variables_required{:});
-if ~all(ismember(model_variables_required, who))
-    error('One or more of the required camera sensitivity variables is not loaded.')
-end
-if channel_mode
-    error('The input space of the camera sensitivity data must be a spectral space, not a space of colour channels.')
-end
+[sensor_map, ~, bands] = loadColorMap(color_map_filename, false);
 
 cutoff_wavelength_sensitivites = max(sensor_map, [], 2);
 sensor_map_relative = sensor_map ./ repmat(cutoff_wavelength_sensitivites, 1, size(sensor_map, 2));
diff --git a/data_analysis/ShowPatches.m b/data_analysis/ShowPatches.m
index cf476fa..c815a1f 100644
--- a/data_analysis/ShowPatches.m
+++ b/data_analysis/ShowPatches.m
@@ -151,13 +151,7 @@
 n_images = length(image_prefixes);
 n_groups = length(algorithm_groups);
 
-load(xyz_weights_filename, xyz_weights_variable);
-if exist(xyz_weights_variable, 'var')
-    xyz_weights = eval(xyz_weights_variable);
-end
-if ~exist(xyz_weights_variable, 'var') || isempty(xyz_weights)
-    error('No raw colour to XYZ conversion matrix loaded.')
-end
+xyz_weights = loadVariables(xyz_weights_filename, xyz_weights_variable);
 
 % Find the class of the image data
 I_filename = fullfile(input_directory, sprintf('%s%s%s', image_prefixes{1}, algorithm_groups{1}{1}, input_postfix));
diff --git a/dispersion_model/RegistrationDispersion.m b/dispersion_model/RegistrationDispersion.m
index f1e6712..4dd953a 100644
--- a/dispersion_model/RegistrationDispersion.m
+++ b/dispersion_model/RegistrationDispersion.m
@@ -216,13 +216,7 @@
 %% Load calibration data
 
 if ~rgb_mode
-    load(bands_filename, bands_variable);
-    if exist(bands_variable, 'var')
-        bands = eval(bands_variable);
-    end
-    if ~exist(bands_variable, 'var') || isempty(bands)
-        error('No wavelengths loaded.')
-    end
+    bands = loadVariables(bands_filename, bands_variable);
     n_bands = length(bands);
     
     [~, reference_index] = min(abs(bands - reference_wavelength));
diff --git a/evaluation/Choi2017Input.m b/evaluation/Choi2017Input.m
index 0111f18..9299d3f 100644
--- a/evaluation/Choi2017Input.m
+++ b/evaluation/Choi2017Input.m
@@ -78,34 +78,12 @@
 
 %% Load wavelengths
 
-load(bands_filename, bands_variable);
-if exist(bands_variable, 'var')
-    bands_spectral = eval(bands_variable);
-end
-if ~exist(bands_variable, 'var') || isempty(bands_spectral)
-    error('No wavelength band information loaded.')
-end
-
-load(desired_bands_filename, desired_bands_variable);
-if exist(desired_bands_variable, 'var')
-    desired_bands_spectral = eval(desired_bands_variable);
-end
-if ~exist(desired_bands_variable, 'var') || isempty(desired_bands_spectral)
-    error('No desired wavelength band information loaded.')
-end
+bands_spectral = loadVariables(bands_filename, bands_variable);
+desired_bands_spectral = loadVariables(desired_bands_filename, desired_bands_variable);
 
 %% Load sensor data
 
-model_variables_required = { 'sensor_map', 'channel_mode', 'bands' };
-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
-if channel_mode
-    error('The input space of the colour conversion data must be a spectral space, not a space of colour channels.')
-end
-
-bands_color = bands;
+[sensor_map, ~, bands_color] = loadColorMap(color_map_filename, false);
 
 %% Load the illuminant
 
diff --git a/evaluation/Choi2017Output.m b/evaluation/Choi2017Output.m
index d5da9eb..c99ef77 100644
--- a/evaluation/Choi2017Output.m
+++ b/evaluation/Choi2017Output.m
@@ -36,11 +36,9 @@
 
 %% Load conversion data
 
-variables_required = { 'spectral_weights_output', 'color_weights' };
-load(conversion_filename, variables_required{:});
-if ~all(ismember(variables_required, who))
-    error('One or more of the required colour/spectral space conversion variables is not loaded.')
-end
+[spectral_weights_output, color_weights] = loadVariables(...
+    conversion_filename, {'spectral_weights_output', 'color_weights'}...
+);
 
 %% Output images
 
diff --git a/evaluation/CorrectChirpImage.m b/evaluation/CorrectChirpImage.m
index 0f8621b..8156332 100644
--- a/evaluation/CorrectChirpImage.m
+++ b/evaluation/CorrectChirpImage.m
@@ -318,15 +318,7 @@
 
 %% Load calibration data
 
-model_variables_required = { 'sensor_map', 'channel_mode', 'bands' };
-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
-if channel_mode
-    error('The input space of the colour conversion data must be a spectral space, not a space of colour channels.')
-end
-bands_color = bands;
+[sensor_map, ~, bands_color] = loadColorMap(color_map_filename, false);
 n_channels = size(sensor_map, 1);
 
 % Normalize spectral sensitivities to avoid out-of-range colour values.
diff --git a/evaluation/EvaluateColorChecker.m b/evaluation/EvaluateColorChecker.m
index fdbcba7..8b345c0 100644
--- a/evaluation/EvaluateColorChecker.m
+++ b/evaluation/EvaluateColorChecker.m
@@ -52,7 +52,7 @@
 % 'SonyColorMap.m', for example. The following variables are required:
 % - '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 spectral images.
+%   j-th element of 'bands' (below).
 % - 'channel_mode': A Boolean value indicating whether the input colour
 %   space is a set of colour channels (true) or a set of spectral bands
 %   (false). A value of `false` is required.
@@ -294,33 +294,13 @@
 bands_qhyper = [];
 color_weights = [];
 if has_spectral
-    load(bands_filename, bands_variable);
-    if exist(bands_variable, 'var')
-        bands_estimated = eval(bands_variable);
-    end
-    if ~exist(bands_variable, 'var') || isempty(bands_estimated)
-        error('No wavelengths loaded.')
-    end
+    bands_estimated = loadVariables(bands_filename, bands_variable);
     n_bands_estimated = length(bands_estimated);
 
-    load(qhyper_bands_filename, qhyper_bands_variable);
-    if exist(qhyper_bands_variable, 'var')
-        bands_qhyper = eval(qhyper_bands_variable);
-    end
-    if ~exist(qhyper_bands_variable, 'var') || isempty(bands_qhyper)
-        error('No wavelengths loaded.')
-    end
+    bands_qhyper = loadVariables(qhyper_bands_filename, qhyper_bands_variable);
     n_bands_qhyper = length(bands_qhyper);
 
-    model_variables_required = { 'sensor_map', 'channel_mode', 'bands' };
-    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
-    if channel_mode
-        error('The input space of the colour conversion data must be a spectral space, not a space of colour channels.')
-    end
-    bands_color = bands;
+    [sensor_map, ~, bands_color] = loadColorMap(color_map_filename, false);
     color_weights = colorWeights(...
         sensor_map, bands_color, bands_estimated, findSamplingOptions...
     );
diff --git a/evaluation/PointSpectralEvaluation.m b/evaluation/PointSpectralEvaluation.m
index 3759b83..cb0f3bb 100644
--- a/evaluation/PointSpectralEvaluation.m
+++ b/evaluation/PointSpectralEvaluation.m
@@ -296,24 +296,12 @@
     
 bands_cell = cell(n_all, 1);
 for i = 1:n_images
-    load(bands_filenames_rep{i}, bands_variable);
-    if exist(bands_variable, 'var')
-        bands_cell{i + n_measurement_sets} = eval(bands_variable);
-    end
-    if ~exist(bands_variable, 'var') || isempty(bands_cell{i + n_measurement_sets})
-        error('No wavelengths loaded from file "%s".', bands_filenames_rep{i})
-    end
+    bands_cell{i + n_measurement_sets} = loadVariables(bands_filenames_rep{i}, bands_variable);
 end
 bands = bands_cell;
 
 if has_images
-    load(location_filename, location_variable);
-    if exist(location_variable, 'var')
-        patch_centers = eval(location_variable);
-    end
-    if ~exist(location_variable, 'var') || isempty(patch_centers)
-        error('No image patch center locations loaded from file "%s".', location_filename)
-    end
+    patch_centers = loadVariables(location_filename, location_variable);
 end
 
 %% Load all spectral measurements and spectral images
diff --git a/evaluation/RunOnDataset.m b/evaluation/RunOnDataset.m
index 749270d..a2a0a21 100644
--- a/evaluation/RunOnDataset.m
+++ b/evaluation/RunOnDataset.m
@@ -256,7 +256,7 @@
 %% Finalize the set of algorithms to run
 
 if ~use_automatic_weights
-    load(admm_algorithms_filename, 'admm_algorithms');
+    admm_algorithms = loadVariables(admm_algorithms_filename, 'admm_algorithms');
 end
 admm_algorithm_fields = fieldnames(admm_algorithms);
 n_admm_algorithms = length(admm_algorithm_fields);
diff --git a/evaluation/helper_script/PreprocessDataset.m b/evaluation/helper_script/PreprocessDataset.m
index 924290d..37eb134 100644
--- a/evaluation/helper_script/PreprocessDataset.m
+++ b/evaluation/helper_script/PreprocessDataset.m
@@ -48,19 +48,13 @@
     error('Spectral images are provided, but not the camera''s spectral sensitivity data.');
 end
 
-bands = [];
-bands_variable = 'bands';
 if has_spectral
     spectral_filenames = listFiles(dp.spectral_images_wildcard);
     n_images_spectral = length(spectral_filenames);
     n_images = n_images_spectral;
     names = trimCommon(spectral_filenames, [false, true]);
     
-    load(dp.wavelengths, bands_variable);
-    if isempty(bands)
-        error('No wavelength band information is associated with the spectral images.')
-    end
-    bands_spectral = bands;
+    bands_spectral = loadVariables(dp.wavelengths, 'bands');
 end
 if has_rgb
     rgb_filenames = listFiles(dp.rgb_images_wildcard);
@@ -98,19 +92,9 @@
     end
 end
 
-bands = [];
 if has_color_map
     % Load colour conversion data
-    
-    model_variables_required = { 'sensor_map', 'channel_mode' };
-    load(dp.color_map, model_variables_required{:}, bands_variable);
-    if ~all(ismember(model_variables_required, who))
-        error('One or more of the required colour space conversion variables is not loaded.')
-    end
-    if isempty(bands)
-        error('No wavelength band information is associated with the colour conversion data.')
-    end
-    bands_color = bands;
+    [sensor_map, channel_mode, bands_color] = loadColorMap(dp.color_map);
     
     if channel_mode
         if (length(bands_spectral) ~= length(bands_color)) ||...
diff --git a/evaluation/loadDispersionModel.m b/evaluation/loadDispersionModel.m
index a9787c9..6e8a810 100644
--- a/evaluation/loadDispersionModel.m
+++ b/evaluation/loadDispersionModel.m
@@ -82,21 +82,11 @@
     bands_req = varargin{1};
 end
 
-dispersion_data = [];
-model_from_reference = [];
-bands = [];
-model_space = [];
-fill = [];
-
-model_variables = { 'dispersion_data', 'model_from_reference', 'bands' };
-variables_transform = { 'model_space', 'fill' };
-load(filename, model_variables{:}, variables_transform{:});
-if isempty(dispersion_data) || isempty(model_from_reference)
-    error('One or more of the dispersion model variables is not loaded.')
-end
-if bands_req && isempty(bands)
-    error('The `bands` variable is empty or is not loaded.')
-end
+varnames = { 'dispersion_data', 'model_from_reference', 'model_space', 'fill', 'bands' };
+not_required = [false(1, 2), true(1, 2), bands_req];
+[...
+    dispersion_data, model_from_reference, model_space, fill, bands...
+] = loadVariables(filename, varnames, not_required, not_required);
 if mfr_expected ~= model_from_reference
     error('Dispersion model is in the wrong frame of reference.')
 end
diff --git a/sensor/loadColorMap.m b/sensor/loadColorMap.m
new file mode 100644
index 0000000..7e152e5
--- /dev/null
+++ b/sensor/loadColorMap.m
@@ -0,0 +1,81 @@
+function [sensor_map, channel_mode, bands] = loadColorMap(filename, channel_mode_expected)
+% LOADCOLORMAP Load colour space conversion data
+%
+% ## Syntax
+% sensor_map = loadColorMap(filename [, channel_mode])
+% [sensor_map, channel_mode] = loadColorMap(filename [, channel_mode])
+% [sensor_map, channel_mode, bands] = loadColorMap(filename [, channel_mode])
+%
+% ## Description
+% sensor_map = loadColorMap(filename [, channel_mode])
+%   Returns the `sensor_map` variable loaded from the file.
+%
+% [sensor_map, channel_mode] = loadColorMap(filename [, channel_mode])
+%   Additionally returns the `channel_mode` variable loaded from the file.
+%
+% [sensor_map, channel_mode, bands] = loadColorMap(filename [, channel_mode])
+%   Additionally returns the `bands` variable loaded from the file.
+%
+% ## Input Arguments
+%
+% filename -- Filename
+%   A character vector containing the filename and path of the '.mat' file from
+%   which to load the data. `filename` must include the file extension.
+%
+% channel_mode -- Expected value of `channel_mode`
+%   A logical scalar indicating the expected value of `channel_mode` in the
+%   file. An error will be thrown if the value of `channel_mode` in the file
+%   does not match this value. If the expected value is empty, or `channel_mode`
+%   is not passed as an input argument, then the value of `channel_mode` in the
+%   file will not be validated.
+%
+% ## Output Arguments
+%
+% sensor_map -- Channel sensitivity data
+%   A 2D array, where `sensor_map(i, j)` is the sensitivity of the i-th colour
+%   channel or spectral band of an output colour space to the j-th colour
+%   channel or spectral band of an input colour space. For example, `sensor_map`
+%   is a matrix mapping discretized spectral power distributions to RGB colours.
+%
+% channel_mode -- Type of sensitivity data
+%   A Boolean value indicating whether the input colour space is a set of
+%   colour channels (true) or a set of spectral bands (false).
+%
+% bands -- Input colour space
+%   A vector containing the wavelengths or colour channel indices corresponding
+%   to the second dimension of 'sensor_map'. The file need not contain this
+%   variable if `bands` is not requested as an output argument.
+%
+% See also loadVariables
+
+% Bernard Llanos
+% Supervised by Dr. Y.H. Yang
+% University of Alberta, Department of Computing Science
+% File created May 23, 2019
+
+nargoutchk(1, 3);
+narginchk(1, 2);
+
+if nargin < 2
+    channel_mode_expected = [];
+end
+
+if ~isempty(channel_mode_expected) || nargout > 1
+    channel_mode = loadVariables(filename, 'channel_mode');
+    if ~isempty(channel_mode_expected)
+        if channel_mode_expected ~= channel_mode
+            error(...
+                '`channel_mode` in "%s" is expected to be `logical(%d)`.',...
+                filename, channel_mode_expected...
+            );
+        end
+    end
+end
+
+if nargout > 2
+    [sensor_map, bands] = loadVariables(filename, {'sensor_map', 'bands'});
+else
+    sensor_map = loadVariables(filename, 'sensor_map');
+end
+
+end
\ No newline at end of file
diff --git a/utilities/loadImage.m b/utilities/loadImage.m
index c3d1b68..bfdc496 100644
--- a/utilities/loadImage.m
+++ b/utilities/loadImage.m
@@ -17,7 +17,7 @@
 %
 % ## Input Arguments
 %
-% filename -- Dispersion model filename
+% filename -- Image filename
 %   A character vector containing the filename and path of the '.mat' or
 %   image format file containing the image, including the file extension.
 %
@@ -59,7 +59,7 @@
 %     ACM Transactions on Graphics (Proc. SIGGRAPH Asia 2017), 36(6),
 %     218:1–13. doi:10.1145/3130800.3130810
 %
-% See also saveImages, imread, load, im2double
+% See also saveImages, imread, loadVariables, im2double
 
 % Bernard Llanos
 % Supervised by Dr. Y.H. Yang
@@ -79,24 +79,16 @@
     else
         variable_name = varargin{1};
     end
-    load(filename, variable_name);
-    if exist(variable_name, 'var')
-        I = eval(variable_name);
-        if ~isnumeric(I)
-            error(...
-                'The input image variable `%s` in %s does not have a numeric class.',...
-                variable_name, filename...
-            );
-        elseif isinteger(I)
-            I = im2double(I);
-        elseif isa(I, 'single')
-            I = double(I);
-        end
-    else
+    I = loadVariables(filename, variable_name);
+    if ~isnumeric(I)
         error(...
-            'The input image variable %s was not loaded from %s.',...
+            'The input image variable `%s` in "%s" does not have a numeric class.',...
             variable_name, filename...
-            );
+        );
+    elseif isinteger(I)
+        I = im2double(I);
+    elseif isa(I, 'single')
+        I = double(I);
     end
 elseif strcmp(ext, exr_ext)
     disp('Loading an OpenEXR image...');
diff --git a/utilities/loadVariables.m b/utilities/loadVariables.m
new file mode 100644
index 0000000..c05e156
--- /dev/null
+++ b/utilities/loadVariables.m
@@ -0,0 +1,109 @@
+function varargout = loadVariables(filename, varnames, empty, none)
+% LOADVARIABLES Load variables from a '.mat' file
+%
+% ## Syntax
+% [val1, val2, ... ] = loadVariables(filename, varnames [, empty, none])
+%
+% ## Description
+% [val1, val2, ... ] = loadVariables(filename, varnames [, empty, none])
+%   Returns the variables loaded from the file.
+%
+% ## Input Arguments
+%
+% filename -- Filename
+%   A character vector or string containing the filename and path of the '.mat'
+%   file, including the file extension.
+%
+% varnames -- Variable names
+%   The name of the single variable to be loaded from the file (a character
+%   vector or string scalar), or a cell vector or string vector of variable
+%   names.
+%
+% empty -- Allow empty variables
+%   A flag indicating whether an empty value of the variable is allowed
+%   (`true`), or whether an error should be raised if the variable has an empty
+%   value (`false`). Defaults to `false` if not passed. If there are multiple
+%   variables to load, `empty` can optionally be a logical vector of
+%   corresponding per-variable flags.
+%
+% none -- Allow missing variables
+%   A flag indicating whether the variable can be missing from the file
+%   (`true`), or whether an error should be raised if the variable is missing
+%   (`false`). Defaults to `false` if not passed. If `none` is `true`, and the
+%   variable is not found in the file, `val` will be empty (`[]`). If there are
+%   multiple variables to load, `empty` can optionally be a logical vector of
+%   corresponding per-variable values.
+%
+% ## Output Arguments
+%
+% val1, val2, ... -- Variable values
+%   The results of loading variables named in `varnames` from the '.mat' file.
+%
+% ## Notes
+% - An error will be thrown if any variables are allowed to be absent from the
+%   '.mat' file, but are expected to be non-empty.
+%
+% See also load, loadImage
+
+% Bernard Llanos
+% Supervised by Dr. Y.H. Yang
+% University of Alberta, Department of Computing Science
+% File created May 23, 2019
+
+
+narginchk(2, 4);
+
+if ischar(varnames)
+    varnames = {varnames};
+elseif ~(iscell(varnames) || isstring(varnames))
+    error('`varnames` must be a character vector, string scalar, string array, or cell vector.');
+end
+n_vars = length(varnames);
+nargoutchk(n_vars, n_vars);
+
+if nargin < 3
+    empty = false(n_vars, 1);
+elseif length(empty) < n_vars
+    if ~isscalar(empty)
+        error('`empty` must have the same length as `varnames`, or must be a scalar.');
+    end
+    empty = repmat(empty, n_vars, 1);
+end
+
+if nargin < 4
+    none = false(n_vars, 1);
+elseif length(none) < n_vars
+    if ~isscalar(none)
+        error('`none` must have the same length as `varnames`, or must be a scalar.');
+    end
+    none = repmat(none, n_vars, 1);
+end
+
+if any(none & ~empty)
+    error('Variables cannot be allowed to be absent, yet expected to be non-empty.');
+end
+
+file = matfile(filename, 'Writable', false);
+varargout = cell(n_vars, 1);
+for i = 1:n_vars
+    try
+        val = file.(varnames{i});
+    catch ME
+        switch ME.identifier
+            case 'MATLAB:MatFile:VariableNotInFile'
+                if none(i)
+                    val = [];
+                else
+                    rethrow(ME)
+                end
+            otherwise
+                rethrow(ME)
+        end
+    end
+    if isempty(val) && ~empty(i)
+        error('Empty value of `%s` loaded from "%s".', varnames{i}, filename);
+    end
+    varargout{i} = val;
+end
+
+end
\ No newline at end of file