%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                      MAIN SIMULATION SCRIPT                         %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This is the main simulation script that we use to create the scene,
% apply the optical low-pass filter implemented with a 2D Gaussian and
% create the optical image, create the sensor image using the default
% sensor, and finally call our custom demosaicing routines through the
% processor. The demosaicing routines are called using nhkRender. 
% We intend to also use this script to call our custom sensor.

% Created by Kunal Ghosh and Hiroshi Shimamoto with the help of Brian
% Wandell and Joyce Farrell for EE 362 / Pych 221

iset

% Create the scene

 scene = sceneCreate('slantedbar');
% scene = sceneCreate;
sensor = sensorCreate;
sensor = sensorSet(sensor,'size',[256 256]); % when frequencyOrientation?
% sensor = sensorSet(sensor,'size',[385 385]); % when slantedbar?
sfov = sensorGet(sensor,'fov');
scene = sceneSet(scene, 'fov', sfov);
%      {'macbethd65'}  - Create a Macbeth D65 image.  Optional
%      {'frequencyOrientation'} - Demosaicking test pattern, equal photon spectrum
%      {'slantedbar'}  - Used for ISO spatial resolution, equal photon spectrum
%      {'zonePlate'}   - Circular zone plot, equal photon spectrum
vcReplaceObject(scene); sceneWindow;

% Create an optical image of the scene with optional optical low-pass filter
oi = oiCreate;
oi = oiSet(oi,'diffuserMethod','none'); % 'blur' to add optical LPF 
oi = oiSet(oi,'diffuserBlur',[2.8 0.0000001]*10^-6); % adjust 2D Gaussian filter
oi = oiCompute(scene,oi);
vcReplaceObject(oi); oiWindow;

% Create the sensor and generate the sensor image with custom color filters

pattern_4x2 = ...
     [2, 1; 
     3, 2; 
     2, 3; 
    1, 2];
pattern_3x3 = ...
    [2, 1, 3; 
     3, 2, 1; 
     1, 3, 2];
 pattern_bayer = ...
     [2, 1;
     3, 2];
sensor = sensorSet(sensor,'pattern',pattern_bayer); % set color filter

% fov = sensorGet(sensor,'fov');
sensor = sensorSetSizeToFOV(sensor,sfov);
sensor = sensorCompute(sensor,vcGetObject('oi'),0);
vcReplaceAndSelectObject(sensor); sensorImageWindow;

% Check the pattern and cfa appearance with:
%   sensorGet(sensor,'pattern')
%   [cfa,cfaN,mp] = sensorDetermineCFA; image(cfaN); colormap(mp);

% Call the custom demosaicing routine and generate the demosaiced image
vci = vcImageCreate;
vci = imageSet(vci,'demosaicMethod','nhkInter');
% vci = imageSet(vci,'colorbalancemethod','grayworld');
vci = nhkRender(vci,sensor); % call the demosaicing script
vcReplaceObject(vci); vcimageWindow;

