Adjust reflectance. Comments and coordination in others

scene/sceneAdjustIlluminant.m

@@ -1,31 +1,35 @@
 function scene = sceneAdjustIlluminant(scene,illEnergy,preserveMean)
 %Adjust the current scene illuminant to the value in data
 %
-% Brief synopsis
+% Synopsis
 %  scene = sceneAdjustIlluminant(scene,illEnergy,preserveMean)
 %
 % Brief Description
-%  Change the scene illuminant.  The calculation scales the scene
-%  radiance by dividing out the current illuminant and then multiplying by
-%  the new illuminant. This preserves the reflectance. By default, we also
-%  preserve the scene mean luminance, which effectively scales the
-%  illuminant level.
+%  Change the scene illuminant.  
 %
-%  If you do not want the illuminant level to change, then set the
-%  preserveMean flag to false. It is true by default.
+% Inputs
+%  scene:      A scene structure, or the current scene will be assumed
+%  illEnergy:  Either a file name to spectral data or a vector (same length
+%              as scene wave) defining the illuminant in energy units or an
+%              ISETCam illuminant struct.
+%  preserveMean:  Scale result to preserve mean illuminant (default true)
 %
-%  If the current scene has no defined illuminant, we assume that the scene
-%  illuminant is D65.
+%  Description
+%    The calculation scales the scene radiance by dividing out the current
+%    illuminant and then multiplying by the new illuminant. This preserves
+%    the reflectance. By default, we also preserve the scene mean luminance,
+%    which effectively scales the illuminant level.
 %
-%  This appears to work if the scene is a spatial-spectral illuminant too. 
+%    If you do not want the illuminant level to change, then set the
+%    preserveMean flag to false. It is true by default.
 %
-% Parameters
-%  scene:      A scene structure, or the current scene will be assumed
-%  illuminant: Either a file name to spectral data or a vector (same length
-%              as scene wave) defining the illuminant in energy units
-%  preserveMean:  Scale result to preserve mean illuminant (default true)
+%    If the current scene has no defined illuminant, we assume that the scene
+%    illuminant is D65.
+%
+%    This appears to work if the scene is a spatial-spectral illuminant
+%    too.  
 %
-% Copyright ImagEval Consultants, LLC, 2010.
+% ieExamplesPrint('sceneAdjustIlluminant');
 %
 % See also
 %    sceneAdjustLuminance
@@ -49,7 +53,7 @@
 %}
 
 %%
-if ieNotDefined('scene'),        scene = vcGetObject('scene'); end
+if ieNotDefined('scene'),        scene = ieGetObject('scene'); end
 if ieNotDefined('preserveMean'), preserveMean = true; end
 
 % Make sure we have the illuminant data in the form of energy
@@ -64,6 +68,9 @@
     if ~exist(fullName,'file'), error('No file %s\n',fullName);
     else, illEnergy = ieReadSpectra(fullName,wave);
     end
+elseif isstruct(illEnergy) && isequal(illEnergy.type,'illuminant')
+    fullName = illuminantGet(illEnergy,'name');
+    illEnergy = illuminantGet(illEnergy,'energy');
 else
     % User sent numbers.  We check for numerical validity next.
     fullName = '';
@@ -73,10 +80,8 @@
 if max(illEnergy) > 10^5
     % Energy is not this big.
     warning('Illuminant energy values are high; may be photons, not energy.')
-elseif isequal(max(isnan(illEnergy(:))),1) || isequal(min(illEnergy(:)),0)
-    warndlg('NaNs or zeros present in proposed illuminant over this wavelength range. No transformation applied.');
-    pause(3);
-    return;
+elseif isequal(min(illEnergy(:)),0)
+    warning('Illuminant transformation cannot be inverted applied (zeroes).');
 end
 
 %% Start the conversion
@@ -86,7 +91,7 @@
     % sceneFromFile (or vcReadImage). Assume the illuminant is D65.  Lord
     % knows why.  Maybe we should do an illuminant estimation algorithm
     % here.
-    disp('Old scene.  Creating D65 illuminant')
+    disp('Old style scene.  Creating D65 illuminant')
     wave   = sceneGet(scene,'wave');
     curIll = ieReadSpectra('D65',wave);   % D65 in energy units
     scene  = sceneSet(scene,'illuminant energy',curIll);
@@ -149,5 +154,5 @@
 
 scene = sceneSet(scene,'illuminant comment',fullName);
 
-return;
+end
 

scene/sceneAdjustReflectance.m

@@ -1,35 +1,111 @@
 function scene = sceneAdjustReflectance(scene,newR)
-% Adjust the scene reflectances, keeping the illuminant unchanged
+% Adjust the reflectance, keeping the illuminant unchanged, for a scene
 %
-% Example:
+% Synopsis
+%    scene = sceneAdjustReflectance(scene,newR)
 %
-% See also:  sceneAdjustIlluminant, sceneAdjustLuminance
+% Inputs:
+%   scene:  A scene (ISETCam scene struct)
+%   newR:   The new reflectance (1D, or row x col x wave) of the scene
+%           reflectance
+%
+% Returns
+%  scene:  The modified scene
+%
+% Description:
+%   Each scene radiance can be decomposed into the illumination data and
+%   the reflectance data.  We do that decomposition first, we replace the
+%   reflectance data with the new reflectance data (newR) and recompute and
+%   return the new scene radiance.
+%
+%   In the general case, we replace the scene reflectances, point by point,
+%   with the reflectances in the row x col x wave newR.
+%
+%   As a special case, when the scene is uniform, newR can be a reflectance
+%   vector that is used at every point in the scene.
+%
+% See also:  
+%  sceneAdjustIlluminant, sceneAdjustLuminance
 %
-% (c) Imageval, 2012
 
+% Examples:
+%{
+  scene = sceneCreate('uniform ee');
+  wave = sceneGet(scene,'wave');
+  fName  = fullfile(oreyeRootPath,'data','tongue','meanReflectance');
+  tongue = ieReadSpectra(fName,wave);
+  scene = sceneAdjustReflectance(scene,tongue);
+%}
+%{
+  % We should test sceneAdjustReflectance more
+%}
+
+%% Parameters
 if ieNotDefined('scene'), error('Scene required.'); end
 if ieNotDefined('newR'), error('RGB format reflectance required'); end
 
-% Check that the new reflectance size matches the current scene
-sz = sceneGet(scene,'size');
-[row,col,w] = size(newR);
-if sz ~= [row,col], error('Reflectance size mis-match'); end
-
-nWave = sceneGet(scene,'nWave');
-if nWave ~=w, error('Wavelength dimension mis-match'); end
+% Check that newR is a reflectance
+assert(max(newR(:)) <= 1.0);
+assert(min(newR(:)) >= 0.0);
 
-%% Format the reflectance for fast multiply
+wave = sceneGet(scene,'wave');
 
-% We multiply the scene illuminant by the reflectance to get the new
-% photons
+illFormat  = sceneGet(scene,'illuminant format');
 illPhotons = sceneGet(scene,'illuminant photons');
+if isvector(newR)
+    % This is the uniform scene case
+    if ~(length(wave) == length(newR))
+        error('newR does not match the scene wavelength');
+    end
+
+    switch illFormat
+        case 'spectral'
+            sz = sceneGet(scene,'size');
+            photons = newR .* illPhotons;
+            photons = repmat(photons(:),1,sz(1),sz(2));
+            photons = permute(photons,[2,3,1]);
+        case 'spatial spectral'
+            warning('Illuminant is spatial-spectral; reflectance is a vector.\n');
+            % Get spatial-spectral illuminant and reformat
+            [photons,row,col] = RGB2XWFormat(illPhotons);
+
+            % Multiply by the one reflectance.  The result is a map of the
+            % illuminant, really.
+            photons = photons*diag(newR);
+
+            % Put the photons back in shape
+            photons = XW2RGBFormat(photons,row,col);
+
+        otherwise
+            error('Unknown illuminant format %s\n',illFormat);
+    end
+
+elseif ndims(newR) == 3
+    if ~(length(wave) == size(newR,3))
+        error('newR does not match the scene wavelength');
+    end
 
-% Format anbd multiply and reformat
-[newR,row,col] = RGB2XWFormat(newR);
-photons = newR*diag(illPhotons);
-photons = XW2RGBFormat(photons,row,col);
+        switch illFormat
+        case 'spectral'
+            % Format, multiply and reformat
+            [newR,row,col] = RGB2XWFormat(newR);
+
+            % Need to allow for the spatial-spectral case
+            photons        = newR*diag(illPhotons);
+
+            photons        = XW2RGBFormat(photons,row,col);
+
+        case 'spatial spectral'
+            photons = illPhotons .* newR;
+
+        otherwise
+            error('Unknown illuminant format %s\n',illFormat);
+        end
+
+
+end
 
-% Set
+% Replace the photons.
 scene = sceneSet(scene,'photons',photons);
 
 end

scripts/image/s_ipGamutReflectance.m

@@ -45,7 +45,7 @@
 % Store the result
 s2 = sceneAdjustReflectance(s,r3);
 s2 = sceneSet(s2,'name','Gaussian');
-ieAddObject(s2); sceneWindow;
+sceneWindow(s2);
 
 %% Show the the modified scene
 
@@ -73,6 +73,7 @@
     r3 = conv2(r3,g,'same');
 end
 r3 = XW2RGBFormat(r3,row,col);
+r3 = ieClip(r3,0,1);
 s3 = sceneAdjustReflectance(s,r3);
 s3 = sceneSet(s3,'name','DoG');
 ieAddObject(s3); sceneWindow;