def lpf(image, sigma, mode=2):
    (Mx, My) = image.shape
    if mode == 1:
        kernel = matlab_style_gauss2D(image.shape, sigma)
        kernel /= numpy.max(kernel)

        if Mx == 1 or My == 1:
            fft = numpy.fft.fft(image)
            fft = numpy.fft.fftshift(fft)
            fft *= kernel
            result = numpy.real(numpy.fft.ifft(numpy.fft.fftshift(fft)))
        else:
            fft = numpy.fft.fftshift(numpy.fft.fft2(image))
            fft *= kernel
            result = numpy.real(numpy.fft.ifft2(numpy.fft.fftshift(fft)))
    elif mode == 2:
        new_dim = 2 * array(image.shape)
        kernel = matlab_style_gauss2D((new_dim[0], new_dim[1]), sigma * 2)
        kernel /= numpy.max(kernel)
        kernel = kernel[Mx:, My:]

        image = image.astype(numpy.double)
        if Mx == 1 or My == 1:
            dct = fftpack.dct(image, type=1)
            dct *= kernel
            result = numpy.real(fftpack.idct(dct, type=1))
        else:
            dct = fftpack.dct(fftpack.dct(image.T, type=2, norm='ortho').T, type=2, norm='ortho')
            dct *= kernel
            result = numpy.real(fftpack.idct(fftpack.idct(dct.T, type=2, norm='ortho').T, type=2, norm='ortho'))
    return result

def process_cube(img, weight, quality):
    # TODO: check to make sure that size of img, and q_tables are consistent

    img = img.copy()

    # print('process_cube input: {}'.format(img))

    this_quality = np.round(np.max(weight)*quality)

    if this_quality < 0:
        this_quality = 0
    if this_quality > quality - 1:
        this_quality = quality - 1

    for i in range(img.shape[3]):
        img[:, :, :, i] = cv2.cvtColor(img[:, :, :, i], cv2.COLOR_BGR2LAB)
    img = np.float32(img)

    # print('process_cube pre DCT: {}'.format(img))

    # img_dct = dct(dct(dct(img, axis=0)/4, axis=1)/4, axis=3)/4
    img_dct = dct(dct(img, axis=0)/4, axis=1)/4

    Q_luma = luminance_tables[:, :, :, this_quality].astype(np.float32)
    Q_chroma = chrominance_tables[:, :, :, this_quality].astype(np.float32)

    # Q_luma[:, :, :] = .01
    # Q_chroma[:, :, :] = .01

    # print('Q_luma: {}'.format(Q_luma))
    # print('Q_chroma: {}'.format(Q_chroma))

    # print('dct, pre rounding: {}'.format(img_dct))
    img_dct[:, :, 0, :] /= Q_luma
    img_dct[:, :, 1, :] /= Q_chroma
    img_dct[:, :, 2, :] /= Q_chroma

    img_dct = np.round(img_dct)

    img_dct[:, :, 0, :] *= Q_luma
    img_dct[:, :, 1, :] *= Q_chroma
    img_dct[:, :, 2, :] *= Q_chroma
    # print('dct, post rounding: {}'.format(img_dct))

    # img_processed = idct(idct(idct(img_dct, axis=0)/4, axis=1)/4, axis=3)/4
    img_processed = idct(idct(img_dct, axis=0)/4, axis=1)/4

    # print('process_cube post DCT: {}'.format(img_processed))

    img_processed = np.clip(img_processed, 0, 255)
    img_processed = np.uint8(img_processed)

    for i in range(img.shape[3]):
        img_processed[:,:,:,i] = cv2.cvtColor(img_processed[:,:,:,i], cv2.COLOR_LAB2BGR)

    # print('process_cube output: {}'.format(img))

    # print('pre dct / post_dct: {}'.format(pre_dct / post_dct))
    return img_processed

def idct2(arr):
    """
    @params { np.ndarray } arr
    @return { np.ndarray }
    """
    array = np.float64(arr)
    result = idct(idct(array, axis=0), axis=1)
    return result

def LSDecompFW(wav, width= 16384, max_nnz_rate=8000.0/262144.0, sparsify = 0.01, taps = 10, 
               level = 3, wl_weight = 1, verbose = False,fc=120):
    
    MaxiterA = 60
   
    length = len(wav)
    
    
    n = sft.next_fast_len(length)
    
    signal = np.zeros((n))
    signal[0:length] = wav[0:length]
     
    h0,h1 = daubcqf(taps,'min')
    L = level
    
    
    #print(n)
    original_signal = lambda s: sft.idct(s[0:n]) + (1.0)*(wl_weight)*idwt(s[n+1:],h0,h1,L)[0]
    LSseparate = lambda x: np.concatenate([sft.dct(x),(1.0)*(wl_weight)*dwt(x,h0,h1,L)[0]],axis=0)
    
    #measurment
    y = signal 
    #FISTA
    ###############################
    cnnz = float("Inf")

    
    c = signal 
    temp = LSseparate(y)
    temp2 = original_signal(temp)
    print('aaa'+str(temp2.shape))
    
    maxabsThetaTy = max(abs(temp))
    
    while cnnz > max_nnz_rate * n:
            
        #FISTA
            tau = sparsify * maxabsThetaTy
            tolA = 1.0e-7
            
            fh = (original_signal,LSseparate)
            
            c = relax.fista(A=fh, b=y,x=LSseparate(c),tol=tolA,l=tau,maxiter=MaxiterA )[0]
            
            cnnz = np.size(np.nonzero(original_signal(c)))
            
            print('nnz = '+ str(cnnz)+ ' / ' + str(n) +' at tau = '+str(tau))
            sparsify = sparsify * 2
            if sparsify == 0.166:
                sparsify = 0.1
    signal_dct = sft.idct(c[0:n])
    signal_dct = signal_dct[0:length]
    signal_wl = (1.0) * float(wl_weight) * idwt(c[n+1:],h0,h1,level)[0] 
    signal_wl = signal_wl[0:length]
   
    return  signal_dct,signal_wl

def laplacian_pca_TV(res, x, f0, lam, gam, iter = 10):
    '''
    TV version of Laplacian embedding
    :param res: resolution of the grid
    :param x: numpy array of data in rows
    :param f0: initial embedding matrix
    :param lam: sparsity parameter
    :param gam: fidelity parameter
    :param iter: number of iterations to carry out
    :return: returns embedding matrix
    '''
    # f0 is an initial projection
    n = res ** 2
    num_data = x.shape[0]

    D = sparse_discrete_diff(res)
    M = 1/(lam*laplacian_eigenvalues(res).reshape(n)+gam)

    f = f0
    y = x .dot(f)
    z = shrink(y .dot(D.T), lam)

    for i in range(iter):

        # Update z
        z_old = z
        z = shrink(y .dot (D.T), lam)

        # Update f
        f_old = f
        u, s, v = la.svd(x.T .dot (y), full_matrices=False)
        f = u .dot(v)

        # Update y
        y_old = y
        q = lam * z .dot (D) + gam * x .dot(f)
        # print('norm of y before is %f' % np.sum(q ** 2))
        y = fftpack.dct(q.reshape((num_data, res, res)), norm='ortho') # Images unraveled as rows
        y = fftpack.dct(np.swapaxes(y,1,2), norm='ortho') # Swap and apply dct on the other side
        # print('norm of y after is %f' % np.sum(y ** 2))
        y = np.apply_along_axis(lambda v: M * v, 1, y.reshape((num_data, n)))
        y = fftpack.idct(y.reshape((num_data, res, res)), norm='ortho')
        y = fftpack.idct(np.swapaxes(y,1,2), norm='ortho')
        y = y.reshape((num_data, n))

        zres = np.sqrt(np.sum((z - z_old) ** 2))
        znorm = np.sqrt(np.sum((z)**2))
        yres = np.sqrt(np.sum((y - y_old) ** 2))
        ynorm = np.sqrt(np.sum((y)**2))
        fres = np.sqrt(np.sum((f - f_old) ** 2))

        value = np.sum(abs(z)) + 0.5*lam*np.sum((z-y .dot(D.T))**2) + 0.5*gam*np.sum((y- x .dot(f) )**2)
        print('Iter %d Val %f Z norm %f Z res %f Ynorm %f Y res %f F res %f' % (i, value, znorm, zres, ynorm, yres, fres))

    return f

 def idct(self, arr,coef_size=0):
     if(coef_size==0):
         idcta2 = fftpack.idct(arr,norm='ortho')
         return idcta2
 
     arrlen = len(arr)
     if coef_size > arrlen:
         new_arr = self.interpolation_with_zeros(arr,coef_size)
     else:
         new_arr = arr[0:int(coef_size)]
     
     idcta2 = fftpack.idct(new_arr,norm='ortho')
     return idcta2

def razafindradina_embed(grayscale_container_path, grayscale_watermark_path, watermarked_image_path, alpha):
    """    
    Razafindradina embedding method implementation. 
    
    Outputs the resulting watermarked image
    
    23-July-2015
    """

    grayscale_container_2darray = numpy.asarray(Image.open(grayscale_container_path).convert("L"))
    grayscale_watermark_2darray = numpy.asarray(Image.open(grayscale_watermark_path).convert("L"))

    assert (
        (grayscale_container_2darray.shape[0] == grayscale_container_2darray.shape[1])
        and (grayscale_container_2darray.shape[0] == grayscale_watermark_2darray.shape[0])
        and (grayscale_container_2darray.shape[1] == grayscale_watermark_2darray.shape[1])
    ), "GrayscaleContainer and GrayscaleWatermark sizes do not match or not square"

    # Perform DCT on GrayscaleContainer

    # print grayscale_container_2darray

    gcdct = dct(dct(grayscale_container_2darray, axis=0, norm="ortho"), axis=1, norm="ortho")

    # print grayscale_container_2darray

    # Perform SchurDecomposition on GrayscaleWatermark

    gwsdt, gwsdu = schur_decomposition(grayscale_watermark_2darray)

    # alpha-blend GrayscaleWatermark TriangularMatrix into GrayscaleContainer DCT coeffs with alpha

    gcdct += gwsdt * alpha

    # Perform IDCT on GrayscaleContainer DCT coeffs to get WatermarkedImage

    watermarked_image_2darray = idct(idct(gcdct, axis=0, norm="ortho"), axis=1, norm="ortho")

    watermarked_image_2darray[watermarked_image_2darray > 255] = 255
    watermarked_image_2darray[watermarked_image_2darray < 0] = 0

    watermarked_image = Image.fromarray(numpy.uint8(watermarked_image_2darray))

    # watermarked_image.show()

    # Write image to file

    watermarked_image.save(watermarked_image_path)

    return

def inverse_transform(data):
    result = []
    for i in range(len(data[0])):
        result.append([])

    for i in range(len(data)):
        partial_result = idct(data[i])
        for j in range(len(partial_result)):
            result[j].append(partial_result[j])

    final_result = []
    for i in range(len(result)):
        partial_result = idct(result[i])
        final_result.append(partial_result)
    return final_result

def outlier_removal_smoothing_dct(tsData, n):
    """
    >>> data = pd.read_csv('input.csv')
    >>> outlier_removal_smoothing_dct(data, 15)
    {'buy': 28833387.903209731, 'sale': 40508532.108399086, 'method': 'outlier_removal_smoothing_dct'}
    >>> outlier_removal_smoothing_dct(data, 20)
    {'buy': 30315166.377325296, 'sale': 42088164.543017924, 'method': 'outlier_removal_smoothing_dct'}
    """

    tsData['ingestdatetime'] = tsData.apply(lambda row: datetime.strptime(str(row['ingestdate']), "%Y%m%d"), axis=1)
    tsData = tsData.sort(['ingestdate'], ascending=[1])

    N = len(tsData)
    t = range(len(tsData))
    x = [float(elem) for elem in tsData['qtyavail']]
    y = dct(x, norm='ortho')
    window = np.zeros(N)
    window[:n] = 1
    yr = idct(y*window, norm='ortho')
    
    diffs = np.diff(yr)
    result_sale = abs(sum(filter(lambda x: x < 0, diffs)))
    result_buy = abs(sum(filter(lambda x: x > 0, diffs)))

    result = {}
    result['sale'] = result_sale
    result['buy'] = result_buy
    result['method'] = inspect.stack()[0][3]
    return result

    def decompress(self, qmat, block_size,
                   point_count, interval, min_ts):
        """Decompresses the bitarray and returns a Dataset object."""
        coeffs = zeros(block_size)
        values = zeros(point_count)
        ind = 0
        i = 0
        block_num = 0
        while ind < self.bits.length():
            if qmat[i] == 52:
                coeffs[i] = 0.
            elif ind > len(self.bits) - (64 - qmat[i]):
                # File is over
                break
            else:
                v = self.bits[ind:ind + 64 - qmat[i]]
                v.extend(qmat[i] * (False,))
                coeffs[i] = struct.unpack(">d", v.tobytes())[0]
                ind += 64 - qmat[i]

            i += 1

            if i >= block_size:
                values[block_num * block_size:block_num * block_size + block_size] = idct(coeffs, norm='ortho')
                i = 0
                block_num += 1
                # We pad out to a full byte at the end of the block
                if ind % 8 != 0:
                    ind += 8 - (ind % 8)
                if i > 1:
                    raise SystemExit

        return Dataset(point_count, interval, min_ts, values)

def idctii(x, axes=None):
    """
    Compute a multi-dimensional inverse DCT-II over specified array axes.
    This function is implemented by calling the one-dimensional inverse
    DCT-II :func:`scipy.fftpack.idct` with normalization mode 'ortho'
    for each of the specified axes.

    Parameters
    ----------
    a : array_like
      Input array
    axes : sequence of ints, optional (default None)
      Axes over which to compute the inverse DCT-II.

    Returns
    -------
    y : ndarray
      Inverse DCT-II of input array
    """

    if axes is None:
        axes = list(range(x.ndim))
    for ax in axes[::-1]:
        x = fftpack.idct(x, type=2, axis=ax, norm='ortho')
    return x

def write_to_image(path, text):
    x1, x2, x3, y1, y2, y3 = [4, 5, 3, 3, 5, 4]
    index = 0
    D = 5
    img = Image.open(path)
    img.getdata()
    bitext = text_to_binary(text)
    bitext = bitext + '0000000000000000'
    r, g, b = [np.array(x) for x in img.split()]
    lx, ly = r.shape()
    for x in xrange(0, lx - 2, 8):
        for y in xrange(0, ly - 2, 8):
            if index == len(bitext) - 1:
                break
            metric = r[x:x + 8, y:y + 8].astype('float')
            metric = dct(metric, norm='ortho')
            if bitext[index] == 1:
                metric[x1, y1] = max(metric[x1, y1], metric[x3, y3] + D + 1)
                metric[x2, y2] = max(metric[x2, y2], metric[x3, y3] + D + 1)
            else:
                metric[x1, y1] = min(metric[x1, y1], metric[x3, y3] - D - 1)
                metric[x2, y2] = min(metric[x2, y2], metric[x3, y3] - D - 1)
            index = index + 1
            metric = idct(metric, norm='ortho')
            r[x:x + 8, y:y + 8] = metric.astype('uint8')
    im = Image.merge("RGB", [Image.fromarray(x) for x in [r, g, b]])
    im.save('%s_writed' % path)

def spectrogramToSnippet(spectrogram):
    snippet = np.zeros(shape = (spectrogram.shape[0]* spectrogram.shape[1]))
    num_windows = spectrogram.shape[0]
    for i in range(num_windows):
        idx = [i*window_len, (i+1)*window_len]
        snippet[idx[0], idx[1]] = idct(windows[i, :])
    return snippet

def partialChebyshev(field, order=1, halfLength=None):

    if halfLength is None:
        halfLength = 1.0;

    if(len(shape(field)) == 2):
        length,breadth = shape(field);
    else:
        length = len(field);
        breadth = 1;
    print length;


    trans = dct(field,type=1,axis=0);
    coeff = arange(length,dtype=float);
    temp = zeros(shape(field), dtype=float);

    trans = 2.0*trans*coeff;

    temp[:,length-2] = trans[:,length-1];

    for i in arange(length-2,0,-1):
        temp[:,i-1] = temp[:,i+1] + trans[:,i]


    return idct(temp, type=1,axis=0)/(2.0*(length-1)*halfLength);

def convolveGaussianDCT(x, sigma, pad_sigma=4, mode="same", cache={}):
    """
    1D convolution of x with Gaussian of width sigma pixels
    If pad_sigma>0, pads ends with zeros by int(pad_sigma*sigma) pixels
    Otherwise does unpadded fast cosine transform, hence reflection from the ends
    """

    fill = int(pad_sigma * sigma)
    actual_size = x.size + fill * 2
    if fill > 0:
        s = nearestFFTnumber(actual_size)
        fill2 = s - x.size - fill
        padded_x = np.pad(x, (fill, fill2), mode="constant")
    else:
        padded_x = x

    s = padded_x.size
    hnorm = sigma / float(s)
    gauss = cache.get((s, hnorm))
    if gauss is None:
        gauss = np.exp(-(np.arange(0, s) * (np.pi * hnorm)) ** 2 / 2.0)
        cache[(s, hnorm)] = gauss
    res = fftpack.idct(fftpack.dct(padded_x, overwrite_x=fill > 0) * gauss, overwrite_x=fill > 0) / (2 * s)
    if fill == 0:
        return res
    if mode == "same":
        return res[fill:-fill2]
    elif mode == "valid":
        return res[fill * 2 : -fill2 - fill]
    else:
        raise ValueError("mode not supported for convolveGaussianDCT")

def icosine_transform(data):
    w, h = data.shape
    image = data.copy()
    for i in range(w):
        image[i][:] = idct(data[i][:], norm='ortho')

    return image

    def on_adjust(self, event=None):
        """
        Executed when slider changes value
        """
        width, height = self.lena.shape
        quant = self.sld.GetValue()

        # reconstructed
        rec = np.zeros(self.lena.shape, dtype=np.int64)

        for y in xrange(0,height,8):
            for x in xrange(0,width,8):
                d = self.lena[y:y+8,x:x+8].astype(np.float)
                D = dct(dct(d.T, norm='ortho').T, norm='ortho').reshape(64)
                Q = np.zeros(64, dtype=np.float)
                Q[self.unzig[:quant]] = D[self.unzig[:quant]]
                Q = Q.reshape([8,8])
                q = np.round(idct(idct(Q.T, norm='ortho').T, norm='ortho'))
                rec[y:y+8,x:x+8] = q.astype(np.int64)

        diff = np.abs(self.lena-rec)

        im = self.axes.flat[0].imshow(self.lena, cmap='gray')
        self.axes.flat[0].set_title('Original Image')

        im = self.axes.flat[1].imshow(rec, cmap='gray')
        self.axes.flat[1].set_title('Reconstructed Image')

        self.hinton(ax=self.axes.flat[2])
        self.axes.flat[2].set_title('DCT Coefficient Mask')

        im = self.axes.flat[3].imshow(diff, cmap='hot', vmin=0, vmax=255)
        self.axes.flat[3].set_title('Error Image')

        for ax in self.axes.flat:
            ax.axis('off')

        self.fig.subplots_adjust(right=0.8)
        cbar_ax = self.fig.add_axes([0.85, 0.15, 0.05, 0.7])
        self.fig.colorbar(im, cax=cbar_ax)

        p = self.psnr(self.lena, rec)
        s = self.compute_ssim(self.lena, rec)

        self.statusbar.SetStatusText('PSNR :%.4f SSIM: %.4f' % (p, s))

        self.canvas.draw()

def idct(c):
    """
    Apply inverce DCT to all elements in b

    :param c:
    :return:
    """
    return fftpack.idct(c, norm='ortho')

def background_reconstruction_dct(I):
    """
    background_reconstruction_dct(I)      where I is a BW image matrix
    Reconstructs the background image by the DCT method proposed in the paper by 
    L.Chen et al (2008) :
    'Automatic TFT-LCD mura defect inspection using discrete cosine transform-based background filtering and ‘just noticeable difference’ quantification strategies'
    
    """
    # Type II DCT is used, the inverse is the Type-III DCT according to SF2 handout
    M = I.shape[0] # number of columns in image matrix
    N = I.shape[1] # number of rows in image matrix
    
    X = fftpack.dct(fftpack.dct(I,  type = 2, norm = 'ortho').T, type = 2, norm = 'ortho').T
    X[1:,:][:,1:]= 0;
    Y =fftpack.idct( fftpack.idct(X, type = 2, norm = 'ortho').T, type =2, norm = 'ortho').T
    Y = np.uint8(Y)
    #cv2.imwrite('stage_3.png', Y)
    return Y

 def foreach_idct(self):
     """
     Выполняет обратное дискретное преобразование косинусов
     :return:
     """
     for i in range(self.nrows):
         for j in range(self.ncols):
             self.image_blocks[i, j, ...] = idct(np.transpose(self.image_blocks_magnitude[i, j, ...]),
                                                 **self.dct_params)