function demosaicedImage = nhkDemosaic(vci,isa) 
% Color demosaicing interface routine.  
% This script is only modified in the sense that it calls our custom
% demosaicing routine under the 'nhkInter' method
% Hiroshi Shimamoto and Kunal Ghosh
%
%    demosaicedImage = Demosaic(vci,isa) 
%
% The input data are in the image sensory array structure. These data
% are a planar image (r,c) with color information.   The demosaic
% algorithms convert the sensor data  to an (r,c,3) using the
% information about the sensor in the ISA structure.
%    
% The returned DEMOSAICEDIMAGE is a (r,c,w) image, where w is usually 3.
%
% The available demosaic methods will grow over time. For now they are
%
%    {'bilinear'}
%    {'laplacian'}
%    {'adaptive laplacian'}
%    {'nearest neighbor'}
%        
%  These all function properly on GRBG Bayer images, but they do not function
%  properly on CMY. There are also some issues with RGGB that we are
%  examining.
%
% Examples:
%  d = Demosaic(vci,isa)
%  d = Demosaic(vci,isa)
%
% Copyright ImagEval Consultants, LLC, 2005.


% Programming:
%    The demosaic routines work  with GRBG Bayer format, but not other
%    sensor configurations.  This needs to be improved.

% Get the image input data.  Normally, these are just a copy of the ISA
% output data and are one-dimensional sensor plane. 
img = imageGet(vci,'input');   %figure; imagesc(img)

% If the data are in planar format, put them into RGB format.  Otherwise,
% if they are already in RGB format (i.e. multiplanar) copy them and move
% on. 
if ndims(img) == 2,  imgRGB = plane2mosaic(img,isa,0);
else                 imgRGB = img; 
end

% Figure out the ordering of the planes as a string of color filter letters
% (r,g,b...)
filterNames = sensorGet(isa,'filterNamesString');

m = imageGet(vci,'demosaicmethod');  
switch lower(m)
    case {'bilinear'}
        method = 'bilinear';
    case {'adaptive laplacian'}
        method = 'adaptivelaplacian';
    case {'laplacian'}
        method = 'laplacian';
    case {'lcc1'}
        method = 'lcc1';
    case {'nearest neighbor'}
        method = 'nearestneighbor';
    otherwise
        % warning('Assuming a customer selected demosaic method.');
        method = m;
end

% The data are now in RGB format.  We know which plane is which color
% because we have the filterNameString.
%
switch lower(method)
       
    case 'bilinear'
        demosaicedImage = Bilinear(imgRGB,sensorGet(isa,'pattern'));
        
    case 'adaptivelaplacian'
        clipToRange = 0;
        bPattern = sensorGet(isa,'colororder');
        if (strcmp(bPattern,'grbg') || ...
                strcmp(bPattern,'rggb') || ...
                strcmp(bPattern,'bggr'))
            demosaicedImage = AdaptiveLaplacian(imgRGB, bPattern, clipToRange);
        else
            hdl = ieSessionGet('vcimagehandles');
            str = 'Sensor color format incompatible with current Laplacian. Using Bilinear.';
            ieInWindowMessage(str,hdl)    
            demosaicedImage = Bilinear(imgRGB,sensorGet(isa,'pattern'));
        end
        
    case 'laplacian'
        clipToRange = 0;
        bPattern = sensorGet(isa,'colororder');
        if (strcmp(bPattern,'grbg') || ...
                strcmp(bPattern,'rggb') || ...
                strcmp(bPattern,'bggr'))
            demosaicedImage = Laplacian(imgRGB,bPattern,clipToRange);
        else
            hdl = ieSessionGet('vcimagehandles');
            str = 'Sensor color format incompatible with current Laplacian. Using Bilinear.';
            ieInWindowMessage(str,hdl)
            demosaicedImage = Bilinear(imgRGB,sensorGet(isa,'pattern'));
        end
        
    case 'lcc1'
        % Not properly implemented yet.
        demosaicedImage = lcc1(imgRGB);
        
    case 'nearestneighbor'
        % Should work on CMY, too, no?
        demosaicedImage = NearestNeighbor(imgRGB,sensorGet(isa,'colororder'));
        
    case 'nhkinter'
        % This is our custom method. We try out various demosaicing
        % algorithms by uncommenting appropriate lines.

%         demosaicedImage = Bilinear(imgRGB,sensorGet(isa,'pattern'));
%%          demosaicedImage = nhkInter(imgRGB,sensorGet(isa,'pattern'));
%        demosaicedImage = nhkIP_4x2_nearestNeighbor(imgRGB,sensorGet(isa,'pattern'));
%          demosaicedImage = nhkIP_4x2_bilinear(imgRGB,sensorGet(isa,'pattern'));
%%         demosaicedImage = nhkIP_3x3_nearestNeighbor(imgRGB,sensorGet(isa,'pattern'));
%%         demosaicedImage = nhkIP_3x3_bilinear(imgRGB,sensorGet(isa,'pattern'));
%         demosaicedImage = nhkIP_3x3_bilinear2D(imgRGB,sensorGet(isa,'pattern'));

        % new demosaicing??
          demosaicedImage = nhkNewDemosaic_nearestNeighbor(imgRGB,sensorGet(isa,'pattern'));
%          demosaicedImage = nhkNewDemosaic_bilinear(imgRGB,sensorGet(isa,'pattern'));
%          demosaicedImage = nhkNewDemosaic_bilinear2D(imgRGB,sensorGet(isa,'pattern'));

        if ~isempty(which(method))
            % demosaicedImage = method(imgRGB);
            %demosaicedImage = feval(method,imgRGB,isa);
        else
            error(sprintf('Demosaic method (%s) not found.',method));
        end
end

return;
