vmquantc.m

有关matlab的电子书籍有一定的帮助希望有用

function [im, map] = vmquantc( r, g, b, k, Q, dith, Qe )
%VMQUANTC Variance Minimization Color Quantization (MEX file).
%   VMQUANTC implements the color quantization algorithm for
%   the function VMQUANT.  It has the same syntax as VMQUANT.
%
%   Don't use this function directly, use VMQUANT or RGB2IND
%   instead.

%   Joseph M. Winograd 6-93
%   Copyright 1993-1998 The MathWorks, Inc.  All Rights Reserved.
%   $Revision: 5.6 $  $Date: 1997/11/24 15:56:50 $

%       Reference: Xiaolin Wu, "Efficient Statistical Computation for
%       Optimal Color QUantization," Graphics Gems II, (ed. James
%       Arvo).  Academic Press: Boston. 1991.

%        Data Structures:
%       P, M0, M1r, M1g, M1b, M2 are 3d matrices packed into 2d
%       form.  They are all three maps of the r, g, b color cube
%       where the i=1, j=1, k=1 planes are all outside the color
%       space (color population=0) to allow direct calculation
%       of statistics for r=0, g=0, b=0 planes.  Then, r,g,b=0:1
%       are mapped into i,j,k=2:Qx+1.

%       
%        Inefficiencies:
%       Mean of each box is calculated during scan as well as for color mapping

%    Issues:
%       All colors in image are quantized to [Qr Qg Qb] bins and
%       color distortion is thereby introduced.  Turn the below
%       flag realmap on, and an unquantized color map will be generated
%       by calculating the mean of the population of each box,
%       at appreciable computational cost.
%

%%%     realmap = 1;
%
        
error('Missing VMQUANTC.MEX');        

%%%error( nargchk( 3, 7, nargin ) );
%%%
%%%if nargin < 6
%%%    dith = 0;
%%%end
%%%
%%%if nargin < 5
%%%    Qr = 33;
%%%    Qg = 33;
%%%    Qb = 33;
%%%else
%%%    Qr = 2^Q(1)+1;
%%%    Qg = 2^Q(2)+1;
%%%    Qb = 2^Q(3)+1;    
%%%end
%%%
%%%if nargin < 4
%%%    k = 256;
%%%else
%%%    k = fix(k);
%%%end
%%%
%%%%    Quantize [R G B] to [Qr Qg Qb] levels.
%%%
%%%bar = waitbar( 0, 'Performing color conversion...' );
%%%rgbsize = [Qr Qg Qb];
%%%rq = round( r * (Qr-2) + 2 );                % Map 0:1 onto 2:Qx+1
%%%gq = round( g * (Qg-2) + 2 );
%%%bq = round( b * (Qb-2) + 2 );
%%%waitbar( 0.01 );
%%%
%%%%    Calculate P(c).
%%%
%%%[row col] = elem3d( rgbsize, rq, gq, bq );
%%%P = sparse( row, col, 1, Qr*Qg, Qb );                
%%%waitbar( .025 );
%%%
%%%%    Calculate zeroth moment of P(c).
%%%
%%%M0 = cumsum3d( rgbsize, P );
%%%
%%%%    Calculate first moment of P(c). 
%%%%    (Use three 3D vectors r, g, b to simulate a 4 dimensional space).
%%%
%%%M1r = cumsum3d( rgbsize, reshape( cumsum( ones( Qr, Qg*Qb ) ), ...
%%%     Qr*Qg, Qb ) .* P );
%%%
%%%M1g = cumsum3d( rgbsize, reshape( cumsum( ones( Qg, Qr*Qb ) )', ...
%%%     Qb, Qr*Qg )' .* P );
%%%
%%%M1b = cumsum3d( rgbsize, cumsum( ones( Qb, Qr*Qg ) )' .* P );
%%%waitbar( .05 );
%%%
%%%%    Calculate second moment of P(c).
%%%
%%%M2 = cumsum3d( rgbsize, (reshape( cumsum( ones( Qr, Qg*Qb ) ), ...
%%%        Qr*Qg, Qb ) .^ 2 + reshape( cumsum( ones( Qg, ...
%%%     Qr*Qb ) )', Qb, Qr*Qg )' .^ 2 + cumsum( ones( ...
%%%     Qb, Qr*Qg ) )' .^ 2) .* P );
%%%waitbar( .1 );
%%%
%%%%    Find k-1 partition locations.
%%%
%%%%    Initialize variables for first pass.
%%%
%%%box = 1;
%%%boxll = [1 1 1];
%%%boxur = rgbsize;
%%%E = zeros(k,1);
%%%E(box) = sum( M2( elem3d( rgbsize, ...
%%%            [boxll(box,1) boxll(box,1) boxur(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxur(box,2) boxll(box,2) boxur(box,2)], ...
%%%            [boxur(box,3) boxll(box,3) boxll(box,3) boxur(box,3)] ) ) - ...
%%%            M2( elem3d( rgbsize, ...
%%%            [boxll(box,1) boxll(box,1) boxur(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxur(box,2) boxll(box,2) boxur(box,2)], ...
%%%            [boxll(box,3) boxur(box,3) boxur(box,3) boxll(box,3)] ) ) );
%%%wbox = sum( M0( elem3d( rgbsize, ...
%%%            [boxll(box,1) boxll(box,1) boxur(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxur(box,2) boxll(box,2) boxur(box,2)], ...
%%%            [boxur(box,3) boxll(box,3) boxll(box,3) boxur(box,3)] ) ) - ...
%%%            M0( elem3d( rgbsize, ...
%%%            [boxll(box,1) boxur(box,1) boxll(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxll(box,2) boxur(box,2) boxur(box,2)], ...
%%%            [boxll(box,3) boxur(box,3) boxur(box,3) boxll(box,3)] ) ) );
%%%M1box = [sum( M1r( elem3d( rgbsize, ...
%%%                [boxll(box,1) boxur(box,1) boxll(box,1) boxur(box,1)], ...
%%%                [boxur(box,2) boxll(box,2) boxll(box,2) boxur(box,2)], ...
%%%                [boxll(box,3) boxll(box,3) boxur(box,3) boxur(box,3)] ) ) - ...
%%%                M1r( elem3d( rgbsize, ...
%%%                [boxll(box,1) boxur(box,1) boxll(box,1) boxur(box,1)], ...
%%%                [boxll(box,2) boxur(box,2) boxur(box,2) boxll(box,2)], ...
%%%                [boxll(box,3) boxll(box,3) boxur(box,3) boxur(box,3)] ) ) ) ...
%%%                sum( M1g( elem3d( rgbsize, ...
%%%                [boxll(box,1) boxur(box,1) boxll(box,1) boxur(box,1)], ...
%%%                [boxur(box,2) boxll(box,2) boxll(box,2) boxur(box,2)], ...
%%%                [boxll(box,3) boxll(box,3) boxur(box,3) boxur(box,3)] ) ) - ...
%%%                M1g( elem3d( rgbsize, ...
%%%                [boxll(box,1) boxur(box,1) boxll(box,1) boxur(box,1)], ...
%%%                [boxll(box,2) boxur(box,2) boxur(box,2) boxll(box,2)], ...
%%%                [boxll(box,3) boxll(box,3) boxur(box,3) boxur(box,3)] ) ) ) ...
%%%                sum( M1b( elem3d( rgbsize, ...
%%%                [boxll(box,1) boxur(box,1) boxll(box,1) boxur(box,1)], ...
%%%                [boxur(box,2) boxll(box,2) boxll(box,2) boxur(box,2)], ...
%%%                [boxll(box,3) boxll(box,3) boxur(box,3) boxur(box,3)] ) ) - ...
%%%                M1b( elem3d( rgbsize, ...
%%%                [boxll(box,1) boxur(box,1) boxll(box,1) boxur(box,1)], ...
%%%                [boxll(box,2) boxur(box,2) boxur(box,2) boxll(box,2)], ...
%%%                [boxll(box,3) boxll(box,3) boxur(box,3) boxur(box,3)] ) ) )];
%%%E(box) = E(box) - (M1box * M1box') / wbox;
%%%
%%%for color=1:k-1
%%%
%%%%    Find the box with the greatest variance.
%%%
%%%    [maxvar box] = max( E );
%%%
%%%    if ~maxvar               %       If no boxes have any variance than we're
%%%     waitbar( .7 );  %       done.  (Thanks, Clay!)
%%%     break;
%%%    end
%%%
%%%%    Find point of minimum variance in r dimension of box.
%%%
%%%    M0base = sum( M0( elem3d( rgbsize, [boxll(box,1) boxll(box,1)], ...
%%%     [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) - ...
%%%     M0( elem3d( rgbsize, [boxll(box,1) boxll(box,1)], ...
%%%     [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) );
%%%    M1base = [sum( M1r( elem3d( rgbsize, [boxll(box,1) boxll(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) - ...
%%%         M1r( elem3d( rgbsize, [boxll(box,1) boxll(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) ) ...
%%%     sum( M1g( elem3d( rgbsize, [boxll(box,1) boxll(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) - ...
%%%         M1g( elem3d( rgbsize, [boxll(box,1) boxll(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) ) ...
%%%     sum( M1b( elem3d( rgbsize, [boxll(box,1) boxll(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) - ...
%%%         M1b( elem3d( rgbsize, [boxll(box,1) boxll(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) )];
%%%    wbox = M0base + sum( M0( elem3d( rgbsize, [boxur(box,1) boxur(box,1)], ...
%%%     [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%     M0( elem3d( rgbsize, [boxur(box,1) boxur(box,1)], ...
%%%     [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) );
%%%    M1box = M1base + ...
%%%     [sum( M1r( elem3d( rgbsize, [boxur(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%         M1r( elem3d( rgbsize, [boxur(box,1) boxur(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) ) ...
%%%     sum( M1g( elem3d( rgbsize, [boxur(box,1) boxur(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%         M1g( elem3d( rgbsize, [boxur(box,1) boxur(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) ) ...
%%%     sum( M1b( elem3d( rgbsize, [boxur(box,1) boxur(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%         M1b( elem3d( rgbsize, [boxur(box,1) boxur(box,1)], ...
%%%         [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) )];
%%%
%%%    varr = zeros(1,boxur(box,1)-boxll(box,1)-1); % Allocate space for the array.
%%%    for i=boxll(box,1)+1:boxur(box,1)-1
%%%     wi = M0base + sum( M0( elem3d( rgbsize, [i i], ...
%%%         [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%         M0( elem3d( rgbsize, [i i], ...
%%%         [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) );
%%%     M1i = M1base + [sum( M1r( elem3d( rgbsize, [i i], ...
%%%            [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%            M1r( elem3d( rgbsize, [i i], ...
%%%            [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) ) ...
%%%            sum( M1g( elem3d( rgbsize, [i i], ...
%%%            [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%            M1g( elem3d( rgbsize, [i i], ...
%%%            [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) ) ...
%%%            sum( M1b( elem3d( rgbsize, [i i], ...
%%%            [boxll(box,2) boxur(box,2)], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%            M1b( elem3d( rgbsize, [i i], ...
%%%            [boxll(box,2) boxur(box,2)], [boxur(box,3) boxll(box,3)] ) ) )];
%%%     vind = i - boxll(box,1);        % Inefficient to calc each time.
%%%     if wi
%%%         varr(vind) = (M1i * M1i') / wi;
%%%     end
%%%     M1diff = M1box - M1i;
%%%     wdiff = wbox - wi;
%%%     if wdiff
%%%         varr(vind) = varr(vind) + (M1diff * M1diff') / wdiff;
%%%     end
%%%    end      % for i (red)
%%%    if isempty(varr)
%%%     varr = -1;
%%%    end
%%%
%%%%       Find point of minimum variance in g dimension of box.
%%% 
%%%    M0base = sum( M0( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%        [boxll(box,2) boxll(box,2)], [boxur(box,3) boxll(box,3)] ) ) - ...
%%%        M0( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%        [boxll(box,2) boxll(box,2)], [boxll(box,3) boxur(box,3)] ) ) );
%%%    M1base = [sum( M1r( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxll(box,2)], [boxur(box,3) boxll(box,3)] ) ) - ...
%%%            M1r( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxll(box,2)], [boxll(box,3) boxur(box,3)] ) ) ) ...
%%%        sum( M1g( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxll(box,2)], [boxur(box,3) boxll(box,3)] ) ) - ...
%%%            M1g( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxll(box,2)], [boxll(box,3) boxur(box,3)] ) ) ) ...
%%%        sum( M1b( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxll(box,2)], [boxur(box,3) boxll(box,3)] ) ) - ...
%%%            M1b( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [boxll(box,2) boxll(box,2)], [boxll(box,3) boxur(box,3)] ) ) )];
%%% 
%%%    varg = zeros(1,boxur(box,2)-boxll(box,2)-1); % Allocate space for the array.
%%%    for i=boxll(box,2)+1:boxur(box,2)-1
%%%        wi = M0base + ...
%%%         sum( M0( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [i i], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%            M0( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [i i], [boxur(box,3) boxll(box,3)] ) ) );
%%%        M1i = M1base + ...
%%%         [sum( M1r( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [i i], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%            M1r( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [i i], [boxur(box,3) boxll(box,3)] ) ) ) ...
%%%            sum( M1g( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [i i], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%            M1g( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [i i], [boxur(box,3) boxll(box,3)] ) ) ) ...
%%%            sum( M1b( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [i i], [boxll(box,3) boxur(box,3)] ) ) - ...
%%%            M1b( elem3d( rgbsize, [boxll(box,1) boxur(box,1)], ...
%%%            [i i], [boxur(box,3) boxll(box,3)] ) ) )];
%%%     vind = i - boxll(box,2);        % Could just be updated...