Issiue with implementation of 2D Discrete Cosine Transform in Python
<p>I'm trying to rewrite Zhao Koch steganography method from matlab into python and I am stuck right at the start.</p>

<p>The first two procedures as they are in matlab:</p>

<p>Step 1:</p>

<pre><code>A = imread(casepath); # Reading stegonography case image and aquiring it's RGB values. In my case it's a 400x400 PNG image, so it gives a 400x400x3 array.
</code></pre>

<p>Step 2:</p>

<pre><code>D = dct2(A(:,:,3)); # Applying 2D DCT to blue values of the image
</code></pre>

<p>Python code analog:</p>

<pre><code>from scipy import misc
from numpy import empty,arange,exp,real,imag,pi
from numpy.fft import rfft,irfft

arr = misc.imread('casepath')# 400x480x3 array (Step 1)
arr[20, 30, 2] # Getting blue pixel value

def dct(y): #Basic DCT build from numpy
    N = len(y)
    y2 = empty(2*N,float)
    y2[:N] = y[:]
    y2[N:] = y[::-1]

    c = rfft(y2)
    phi = exp(-1j*pi*arange(N)/(2*N))
    return real(phi*c[:N])


def dct2(y): #2D DCT bulid from numpy and using prvious DCT function
    M = y.shape[0]
    N = y.shape[1]
    a = empty([M,N],float)
    b = empty([M,N],float)

    for i in range(M):
        a[i,:] = dct(y[i,:])
    for j in range(N):
        b[:,j] = dct(a[:,j])

    return b

 D = dct2(arr) # step 2 anlogue
</code></pre>

<p>However, when I try to execute the code I get the following error:</p>

<pre><code>Traceback (most recent call last):
File "path to .py file", line 31, in &lt;module&gt;
D = dct2(arr)
File "path to .py file", line 25, in dct2
a[i,:] = dct(y[i,:])
File "path to .py file", line 10, in dct
y2[:N] = y[:]
ValueError: could not broadcast input array from shape (400,3) into shape (400)
</code></pre>

<p>Perhaps someone could kindly explain to me what am I doing wrong?</p>

<p>Additional Info:
OS: Windows 10 Pro 64 bit
Python: 2.7.12
scipy:0.18.1
numpy:1.11.2
pillow: 3.4.1</p>

<p>Your code works fine, but it is designed to only accept a 2D array, just like <a href="https://uk.mathworks.com/help/images/ref/dct2.html" rel="nofollow"><code>dct2()</code></a> in Matlab. Since your <code>arr</code> is a 3D array, you want to do</p>

<pre><code>D = dct2(arr[...,2])
</code></pre>

<p>As mentioned in my comment, instead or reinventing the wheel, use the (fast) built-in <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.dct.html" rel="nofollow"><code>dct()</code></a> from the scipy package.</p>

<p>The code from the link in my comment effectively provides you this:</p>

<pre><code>import numpy as np
from scipy.fftpack import dct, idct

def dct2(block):
    return dct(dct(block.T, norm='ortho').T, norm='ortho')

def idct2(block):
    return idct(idct(block.T, norm='ortho').T, norm='ortho')
</code></pre>

<p>But again, I must stress that <strong>you have to call this function for each colour plane individually</strong>. Scipy's <code>dct()</code> will happily accept any N-dimensional array and will apply the transform on the last axis. Since that's your colour planes and not your rows and columns of your pixels, you'll get the wrong result. Yes, there is a way to address this with the <code>axis</code> input parameter, but I won't unnecessarily overcomplicate this answer.</p>

<hr>

<p>Regarding the various DCT implementations involved here, your version and scipy's implementation give the same result if you omit the <code>norm='ortho'</code> parameter from the snippet above. But with that parameter included, scipy's transform will agree with Matlab's.</p>

python matlab dct