images_private_dct.cpp

一篇有关数字水印的程序

//
// MATLAB Compiler: 3.0
// Date: Wed Apr 25 16:41:03 2007
// Arguments: "-B" "macro_default" "-O" "all" "-O" "fold_scalar_mxarrays:on"
// "-O" "fold_non_scalar_mxarrays:on" "-O" "optimize_integer_for_loops:on" "-O"
// "array_indexing:on" "-O" "optimize_conditionals:on" "-B" "sglcpp" "-p" "-W"
// "main" "-L" "Cpp" "-t" "-T" "link:exe" "-h" "libmmfile.mlib" "-W" "mainhg"
// "libmwsglm.mlib" "-l" "-p" "-W" "main" "-L" "Cpp" "-t" "-T" "link:exe" "-h"
// "libmmfile.mlib" "-h" "dct2.m" 
//
#include "images_private_dct.hpp"
#include "images_private_checknargin.hpp"
#include "libmatlbm.hpp"
#include "libmmfile.hpp"
static mwArray _mxarray0_ = mclInitializeDouble(1.0);
static mwArray _mxarray1_ = mclInitializeDouble(2.0);

static mxChar _array3_[6] = { 'd', 'o', 'u', 'b', 'l', 'e' };
static mwArray _mxarray2_ = mclInitializeString(6, _array3_);
static mwArray _mxarray4_ = mclInitializeDouble(0.0);
static mwArray _mxarray5_ = mclInitializeComplex(-0.0, -1.0);
static mwArray _mxarray6_ = mclInitializeDouble(3.141592653589793);
static mwArray _mxarray7_ = mclInitializeDouble(-2.0);

void InitializeModule_images_private_dct() {
}

void TerminateModule_images_private_dct() {
}

static mwArray Mimages_private_dct(int nargout_, mwArray a, mwArray n);

_mexLocalFunctionTable _local_function_table_images_private_dct
  = { 0, (mexFunctionTableEntry *)NULL };

//
// The function "images_private_dct" contains the normal interface for the
// "images/private/dct" M-function from file
// "d:\matlab6p5\toolbox\images\images\private\dct.m" (lines 1-85). This
// function processes any input arguments and passes them to the implementation
// version of the function, appearing above.
//
mwArray images_private_dct(mwArray a, mwArray n) {
    int nargout = 1;
    mwArray b = mwArray::UNDEFINED;
    b = Mimages_private_dct(nargout, a, n);
    return b;
}

//
// The function "mlxImages_private_dct" contains the feval interface for the
// "images/private/dct" M-function from file
// "d:\matlab6p5\toolbox\images\images\private\dct.m" (lines 1-85). The feval
// function calls the implementation version of images/private/dct through this
// function. This function processes any input arguments and passes them to the
// implementation version of the function, appearing above.
//
void mlxImages_private_dct(int nlhs,
                           mxArray * plhs[],
                           int nrhs,
                           mxArray * prhs[]) {
    MW_BEGIN_MLX();
    {
        mwArray mprhs[2];
        mwArray mplhs[1];
        int i;
        mclCppUndefineArrays(1, mplhs);
        if (nlhs > 1) {
            error(
              mwVarargin(
                mwArray(
                  "Run-time Error: File: images/private/dct Line: 1 Colum"
                  "n: 1 The function \"images/private/dct\" was called wi"
                  "th more than the declared number of outputs (1).")));
        }
        if (nrhs > 2) {
            error(
              mwVarargin(
                mwArray(
                  "Run-time Error: File: images/private/dct Line: 1 Colu"
                  "mn: 1 The function \"images/private/dct\" was called "
                  "with more than the declared number of inputs (2).")));
        }
        for (i = 0; i < 2 && i < nrhs; ++i) {
            mprhs[i] = mwArray(prhs[i], 0);
        }
        for (; i < 2; ++i) {
            mprhs[i].MakeDIN();
        }
        mplhs[0] = Mimages_private_dct(nlhs, mprhs[0], mprhs[1]);
        plhs[0] = mplhs[0].FreezeData();
    }
    MW_END_MLX();
}

//
// The function "Mimages_private_dct" is the implementation version of the
// "images/private/dct" M-function from file
// "d:\matlab6p5\toolbox\images\images\private\dct.m" (lines 1-85). It contains
// the actual compiled code for that M-function. It is a static function and
// must only be called from one of the interface functions, appearing below.
//
//
// function b=dct(a,n)
//
static mwArray Mimages_private_dct(int nargout_, mwArray a, mwArray n) {
    mclMlineEnterFunction(
      "d:\\matlab6p5\\toolbox\\images\\images\\private\\dct.m",
      "images/private/dct")
    mwLocalFunctionTable save_local_function_table_
      = &_local_function_table_images_private_dct;
    int nargin_ = nargin(2, mwVarargin(a, n));
    mwArray b = mwArray::UNDEFINED;
    mwArray W = mwArray::UNDEFINED;
    mwArray ww = mwArray::UNDEFINED;
    mwArray yy = mwArray::UNDEFINED;
    mwArray y = mwArray::UNDEFINED;
    mwArray aa = mwArray::UNDEFINED;
    mwArray m = mwArray::UNDEFINED;
    mwArray do_trans = mwArray::UNDEFINED;
    mwArray ans = mwArray::UNDEFINED;
    //
    // %DCT Discrete cosine transform.
    // %
    // %   Y = DCT(X) returns the discrete cosine transform of X. The
    // %   vector Y is the same size as X and contains the discrete
    // %   cosine transform coefficients.  
    // %
    // %   Y = DCT(X,N) pads or truncates the vector X to length N
    // %   before transforming.  
    // %
    // %   If X is a matrix, the DCT operation is applied to each
    // %   column.  This transform can be inverted using IDCT.  
    // %
    // %   See also FFT,IFFT, and IDCT.
    // 
    // %   Copyright 1993-2002 The MathWorks, Inc.  
    // %   $Revision: 5.12 $  $Date: 2002/03/15 15:57:11 $
    // 
    // %   References: 
    // %       1) A. K. Jain, "Fundamentals of Digital Image
    // %          Processing", pp. 150-153.
    // %       2) Wallace, "The JPEG Still Picture Compression Standard",
    // %          Communications of the ACM, April 1991.
    // 
    // checknargin(1,2,nargin,mfilename);
    //
    mclMline(25);
    images_private_checknargin(
      _mxarray0_, _mxarray1_, nargin_, "images/private/dct");
    //
    // 
    // if ~isa(a, 'double')
    //
    mclMline(27);
    if (mclNotBool(isa(mwVa(a, "a"), _mxarray2_))) {
        //
        // a = double(a);
        //
        mclMline(28);
        a = double_func(mwVa(a, "a"));
    //
    // end
    //
    mclMline(29);
    }
    //
    // 
    // if min(size(a))==1 
    //
    mclMline(31);
    if (mclEqBool(min(size(mwValueVarargout(), mwVa(a, "a"))), _mxarray0_)) {
        //
        // if size(a,2)>1
        //
        mclMline(32);
        if (mclGtBool(
              size(mwValueVarargout(), mwVa(a, "a"), _mxarray1_), _mxarray0_)) {
            //
            // do_trans = 1;
            //
            mclMline(33);
            do_trans = _mxarray0_;
        //
        // else
        //
        mclMline(34);
        } else {
            //
            // do_trans = 0;
            //
            mclMline(35);
            do_trans = _mxarray4_;
        //
        // end
        //
        mclMline(36);
        }
        //
        // a = a(:); 
        //
        mclMline(37);
        a = mclArrayRef(mwVa(a, "a"), colon());
    //
    // else  
    //
    mclMline(38);
    } else {
        //
        // do_trans = 0; 
        //
        mclMline(39);
        do_trans = _mxarray4_;
    //
    // end
    //
    mclMline(40);
    }
    //
    // if nargin==1,
    //
    mclMline(41);
    if (nargin_ == 1) {
        //
        // n = size(a,1);
        //
        mclMline(42);
        n = size(mwValueVarargout(), mwVa(a, "a"), _mxarray0_);
    //
    // end
    //
    mclMline(43);
    }
    //
    // m = size(a,2);
    //
    mclMline(44);
    m = size(mwValueVarargout(), mwVa(a, "a"), _mxarray1_);
    //
    // 
    // % Pad or truncate a if necessary
    // if size(a,1)<n,
    //
    mclMline(47);
    if (mclLtBool(
          size(mwValueVarargout(), mwVa(a, "a"), _mxarray0_), mwVa(n, "n"))) {
        //
        // aa = zeros(n,m);
        //
        mclMline(48);
        aa = zeros(mwVarargin(mwVa(n, "n"), mwVv(m, "m")));
        //
        // aa(1:size(a,1),:) = a;
        //
        mclMline(49);
        mclArrayAssign(
          &aa,
          mwVa(a, "a"),
          colon(_mxarray0_, size(mwValueVarargout(), mwVa(a, "a"), _mxarray0_)),
          colon());
    //
    // else
    //
    mclMline(50);
    } else {
        //
        // aa = a(1:n,:);
        //
        mclMline(51);
        aa
          = mclArrayRef(
              mwVa(a, "a"), colon(_mxarray0_, mwVa(n, "n")), colon());
    //
    // end
    //
    mclMline(52);
    }
    //
    // 
    // if rem(n,2)==1 | ~isreal(a), % odd case
    //
    mclMline(54);
    {
        mwArray a_ = rem(mwVa(n, "n"), _mxarray1_) == _mxarray0_;
        if (tobool(a_) || tobool(a_ | ~ isreal(mwVa(a, "a")))) {
            //
            // % Form intermediate even-symmetric matrix.
            // y = zeros(2*n,m);
            //
            mclMline(56);
            y = zeros(mwVarargin(_mxarray1_ * mwVa(n, "n"), mwVv(m, "m")));
            //
            // y(1:n,:) = aa;
            //
            mclMline(57);
            mclArrayAssign(
              &y, mwVv(aa, "aa"), colon(_mxarray0_, mwVa(n, "n")), colon());
            //
            // y(n+1:n+n,:) = flipud(aa);
            //
            mclMline(58);
            mclArrayAssign(
              &y,
              flipud(mwVv(aa, "aa")),
              colon(mwVa(n, "n") + _mxarray0_, mwVa(n, "n") + mwVa(n, "n")),
              colon());
            //
            // 
            // % Perform FFT
            // yy = fft(y);
            //
            mclMline(61);
            yy = fft(mwVv(y, "y"));
            //
            // 
            // % Compute DCT coefficients
            // ww = (exp(-i*(0:n-1)*pi/(2*n))/sqrt(2*n)).';
            //
            mclMline(64);
            ww
              = transpose(
                  exp(
                    _mxarray5_ * colon(_mxarray4_, mwVa(n, "n") - _mxarray0_)
                    * _mxarray6_
                    / (_mxarray1_ * mwVa(n, "n")))
                  / sqrt(_mxarray1_ * mwVa(n, "n")));
            //
            // ww(1) = ww(1) / sqrt(2);
            //
            mclMline(65);
            mclIntArrayAssign(
              &ww, mclIntArrayRef(mwVv(ww, "ww"), 1) / sqrt(_mxarray1_), 1);
            //
            // b = ww(:,ones(1,m)).*yy(1:n,:);
            //
            mclMline(66);
            b
              = times(
                  mclArrayRef(
                    mwVv(ww, "ww"),
                    colon(),
                    ones(mwVarargin(_mxarray0_, mwVv(m, "m")))),
                  mclArrayRef(
                    mwVv(yy, "yy"), colon(_mxarray0_, mwVa(n, "n")), colon()));
        //
        // 
        // else % even case
        //
        mclMline(68);
        } else {
            //
            // 
            // % Re-order the elements of the columns of x
            // y = [ aa(1:2:n,:); aa(n:-2:2,:) ];
            //
            mclMline(71);
            y
              = vertcat(
                  mwVarargin(
                    mclArrayRef(
                      mwVv(aa, "aa"),
                      colon(_mxarray0_, _mxarray1_, mwVa(n, "n")),
                      colon()),
                    mclArrayRef(
                      mwVv(aa, "aa"),
                      colon(mwVa(n, "n"), _mxarray7_, _mxarray1_),
                      colon())));
            //
            // 
            // % Compute weights to multiply DFT coefficients
            // ww = 2*exp(-i*(0:n-1)'*pi/(2*n))/sqrt(2*n);
            //
            mclMline(74);
            ww
              = _mxarray1_
                * exp(
                    _times_transpose(
                      _mxarray5_,
                      colon(_mxarray4_, mwVa(n, "n") - _mxarray0_),
                      _mxarray1_)
                    * _mxarray6_
                    / (_mxarray1_ * mwVa(n, "n")))
                / sqrt(_mxarray1_ * mwVa(n, "n"));
            //
            // ww(1) = ww(1) / sqrt(2);
            //
            mclMline(75);
            mclIntArrayAssign(
              &ww, mclIntArrayRef(mwVv(ww, "ww"), 1) / sqrt(_mxarray1_), 1);
            //
            // W = ww(:,ones(1,m));
            //
            mclMline(76);
            W
              = mclArrayRef(
                  mwVv(ww, "ww"),
                  colon(),
                  ones(mwVarargin(_mxarray0_, mwVv(m, "m"))));
            //
            // 
            // % Compute DCT using equation (5.92) in Jain
            // b = W .* fft(y);
            //
            mclMline(79);
            b = times(mwVv(W, "W"), fft(mwVv(y, "y")));
        }
    //
    // end
    //
    mclMline(80);
    }
    //
    // 
    // if isreal(a), b = real(b); end
    //
    mclMline(82);
    if (tobool(isreal(mwVa(a, "a")))) {
        b = real(mwVv(b, "b"));
    }
    //
    // if do_trans, b = b.'; end
    //
    mclMline(83);
    if (tobool(mwVv(do_trans, "do_trans"))) {
        b = transpose(mwVv(b, "b"));
    }
    mwValidateOutput(b, 1, nargout_, "b", "images/private/dct");
    mclMlineFunctionReturn()
    return b;
    //
    // 
    //
    mclMlineExitFunctionReturn();
}