Commit 174904aa authored by Nathaniel Callens's avatar Nathaniel Callens
Browse files

deletions

parent 3ee02b10
Loading
Loading
Loading
Loading
+0 −431

File deleted.

Preview size limit exceeded, changes collapsed.

+0 −877

File deleted.

Preview size limit exceeded, changes collapsed.

+76 −13
Original line number Diff line number Diff line
%% Cell type:code id:5bb42c2c tags:

``` python
import numpy as np
from prediction_MSE_Scout import file_extractor, image_extractor, im_distribution
from matplotlib import pyplot as plt
from itertools import product
from sklearn.metrics import mean_squared_error
import os
import sys
from PIL import Image
from scipy.optimize import minimize
from time import time
from numpy import linalg as la
from scipy.stats import gaussian_kde
import seaborn as sns
from collections import Counter
import pandas as pd
import scipy as sp
```

%% Output

    Average Error: 16.119272850205995
    Standard Deviaiton of Mean Errors: 3.4250443465496736
    Average Difference: 38.50191468436905
    Average Time per Image for First: 0.06527793407440186


    Std Deviation of E:  16.59385258507158
    Normal bits:  15
    Encoded Bits:  5.995576969735472
    (258, 322)

%% Cell type:code id:ec24fcba tags:

``` python
def plot_hist(tiff_list):
    """
    This function is the leftovers from the first attempt to plot histograms.
    As it stands it needs some work in order to function again. We will
    fix this later. 1/25/22
    """

    image = tiff_list
    image = Image.open(image)    #Open the image and read it as an Image object
    image = np.array(image)[1:,:]    #Convert to an array, leaving out the first row because the first row is just housekeeping data
    image = image.astype(int)
    A = np.array([[3,0,-1],[0,3,3],[1,-3,-4]]) # the matrix for system of equation
    z0 = image[0:-2,0:-2]   # get all the first pixel for the entire image
    z1 = image[0:-2,1:-1]   # get all the second pixel for the entire image
    z2 = image[0:-2,2::]    # get all the third pixel for the entire image
    z3 = image[1:-1,0:-2]   # get all the forth pixel for the entire image
    # calculate the out put of the system of equation
    y0 = np.ravel(-z0+z2-z3)
    y1 = np.ravel(z0+z1+z2)
    y2 = np.ravel(-z0-z1-z2-z3)
    y = np.vstack((y0,y1,y2))
    # use numpy solver to solve the system of equations all at once
    #predict = np.floor(np.linalg.solve(A,y)[-1])
    predict = np.round(np.round((np.linalg.solve(A,y)[-1]),1))
    # flatten the neighbor pixlels and stack them together
    z0 = np.ravel(z0)
    z1 = np.ravel(z1)
    z2 = np.ravel(z2)
    z3 = np.ravel(z3)
    neighbor = np.vstack((z0,z1,z2,z3)).T
    # calculate the difference
    diff = np.max(neighbor,axis = 1) - np.min(neighbor, axis=1)

    # flatten the image to a vector
    image = np.ravel(image[1:-1,1:-1])
    error = image-predict

    return image, predict, diff, error, A
```

%% Cell type:code id:c2430512 tags:

``` python
class NodeTree(object):
    def __init__(self, left=None, right=None):
        self.left = left
        self.right = right

    def children(self):
        return self.left, self.right

    def __str__(self):
        return self.left, self.right


def huffman_code_tree(node, binString=''):
    '''
    Function to find Huffman Code
    '''
    if type(node) is str:
        return {node: binString}
    (l, r) = node.children()
    d = dict()
    d.update(huffman_code_tree(l, binString + '0'))
    d.update(huffman_code_tree(r, binString + '1'))
    return d


def make_tree(nodes):
    '''
    Function to make tree
    :param nodes: Nodes
    :return: Root of the tree
    '''
    while len(nodes) > 1:
        (key1, c1) = nodes[-1]
        (key2, c2) = nodes[-2]
        nodes = nodes[:-2]
        node = NodeTree(key1, key2)
        nodes.append((node, c1 + c2))
        nodes = sorted(nodes, key=lambda x: x[1], reverse=True)
    return nodes[0][0]
```

%% Cell type:code id:f62c1af6 tags:
%% Cell type:code id:48abcf1e tags:

``` python
def plot_hist_lstsq(tiff_list, lam=.75):

    image = tiff_list
    image = Image.open(image)    #Open the image and read it as an Image object
    image = np.array(image)[1:,:]    #Convert to an array, leaving out the first row because the first row is just housekeeping data
    image = image.astype(int)
    A = np.array([[3,0,-1],[0,3,3],[1,-3,-4]]) # the matrix for system of equation
    z0 = image[0:-2,0:-2]   # get all the first pixel for the entire image
    z1 = image[0:-2,1:-1]   # get all the second pixel for the entire image
    z2 = image[0:-2,2::]    # get all the third pixel for the entire image
    z3 = image[1:-1,0:-2]   # get all the forth pixel for the entire image
    # calculate the out put of the system of equation
    y0 = np.ravel(-z0+z2-z3)
    y1 = np.ravel(z0+z1+z2)
    y2 = np.ravel(-z0-z1-z2-z3)
    y = np.vstack((y0,y1,y2))
    # use numpy solver to solve the system of equations all at once
    #predict = np.floor(np.linalg.solve(A,y)[-1])
    predict = np.round(np.round((np.linalg.solve(A,y)[-1]),1))

    points = np.array([[-1,-1,1], [-1,0,1], [-1,1,1], [0,-1,1]])
    #fit = la.solve(A,y)

    #mse_start = (points@fit).T


    # flatten the neighbor pixels and stack them together
    z0 = np.ravel(z0)
    z1 = np.ravel(z1)
    z2 = np.ravel(z2)
    z3 = np.ravel(z3)
    neighbor = np.vstack((z0,z1,z2,z3)).T

    f, res, rank, s = la.lstsq(points, neighbor.T, rcond=None)

    #mse_finish = (neighbor-mse_start)**2
    #lstsqur = np.sum(mse_finish, axis=1) / 4

    # calculate the difference
    diff = np.max(neighbor,axis = 1) - np.min(neighbor, axis=1)

    # flatten the image to a vector
    image = np.ravel(image[1:-1,1:-1])
    error = image-predict

    return image, predict, res, error, A, diff, (lam*res + (1-lam)*diff)
```

%% Cell type:code id:b973ed91 tags:

``` python
def huffman(image):
    origin, predict, res, error, A, diff, combo = plot_hist_lstsq(image)
    il = res.astype(int)
    num_bins = 5
    data_points_per_bin = len(il) // num_bins   #l is the list of data that you want to create bins for

    sorted_l = il.copy()
    sorted_l.sort()

    bins = [sorted_l[_ * data_points_per_bin: (_+1)*data_points_per_bin] for _ in range(num_bins)]

    image = Image.open(image)
    image = np.array(image)[1:,:]    #Convert to an array, leaving out the first row because the first row is just housekeeping data
    image = image.astype(int)

    boundary = np.hstack((image[0,:],image[-1,:],image[1:-1,0],image[1:-1,-1]))
    boundary = boundary - image[0,0]
    boundary[0] = image[0,0]

    string = [str(i) for i in boundary]
    freq = dict(Counter(string))
    freq = sorted(freq.items(), key=lambda x: x[1], reverse=True)
    node = make_tree(freq)
    encode1 = huffman_code_tree(node)


    mask = res <= np.max(bins[0])
    string = [str(i) for i in error[mask].astype(int)]
    freq = dict(Counter(string))
    freq = sorted(freq.items(), key=lambda x: x[1], reverse=True)
    node = make_tree(freq)
    encode2 = huffman_code_tree(node)


    mask = res > np.max(bins[0])
    new_error = error[mask]
    mask2 = res[mask] <= np.max(bins[1])
    string = [str(i) for i in new_error[mask2].astype(int)]
    freq = dict(Counter(string))
    freq = sorted(freq.items(), key=lambda x: x[1], reverse=True)
    node = make_tree(freq)
    encode3 = huffman_code_tree(node)


    mask = res > np.max(bins[1])
    new_error = error[mask]
    mask2 = res[mask] <= np.max(bins[2])
    string = [str(i) for i in new_error[mask2].astype(int)]
    freq = dict(Counter(string))
    freq = sorted(freq.items(), key=lambda x: x[1], reverse=True)
    node = make_tree(freq)
    encode4 = huffman_code_tree(node)


    mask = res > np.max(bins[2])
    new_error = error[mask]
    mask2 = res[mask] <= np.max(bins[3])
    string = [str(i) for i in new_error[mask2].astype(int)]
    freq = dict(Counter(string))
    freq = sorted(freq.items(), key=lambda x: x[1], reverse=True)
    node = make_tree(freq)
    encode5 = huffman_code_tree(node)

    mask = res > np.max(bins[3])
    string = [str(i) for i in error[mask].astype(int)]
    freq = dict(Counter(string))
    freq = sorted(freq.items(), key=lambda x: x[1], reverse=True)
    node = make_tree(freq)
    encode6 = huffman_code_tree(node)


    new_error = np.copy(image)
    new_error[1:-1,1:-1] = np.reshape(error,(510, 638))
    keep = new_error[0,0]
    new_error[0,:] = new_error[0,:] - keep
    new_error[-1,:] = new_error[-1,:] - keep
    new_error[1:-1,0] = new_error[1:-1,0] - keep
    new_error[1:-1,-1] = new_error[1:-1,-1] - keep
    new_error[0,0] = keep


    #new_error = np.ravel(new_error)

    new_bins = [np.max(bins[0]),np.max(bins[1]),np.max(bins[2]), np.max(bins[3]), np.max(bins[4])]

    # return the huffman dictionary
    return [encode1, encode2, encode3, encode4, encode5, encode6], np.ravel(image), error, new_error, diff, boundary, new_bins, predict, A, res

```

%% Cell type:code id:9e91c81d tags:
%% Cell type:code id:0afd3bef tags:

``` python
def compress_rate(image, res, error, diff, bound, list_dic, bins):
    #original = original.reshape(-1)
    #error = error.reshape(-1)

    o_len = 0
    c_len = 0
    im = np.reshape(image,(512, 640))
    real_b = np.hstack((im[0,:],im[-1,:],im[1:-1,0],im[1:-1,-1]))
    original = im[1:-1,1:-1].reshape(-1)
    diff = diff.reshape(-1)

    # calculate the bit for boundary
    for i in range(0,len(bound)):
        o_len += len(bin(real_b[i])[2:])
        c_len += len(list_dic[0][str(bound[i])])


    for i in range(0,len(original)):
        o_len += len(bin(original[i])[2:])
        if res[i] <= bins[0]:
            c_len += len(list_dic[1][str(int(error[i]))])

        elif res[i] <= bins[1] and res[i] > bins[0]:
            c_len += len(list_dic[2][str(int(error[i]))])

        elif res[i] <= bins[2] and res[i] > bins[1]:
            c_len += len(list_dic[3][str(int(error[i]))])

        elif res[i] <= bins[3] and res[i] > bins[2]:
            c_len += len(list_dic[4][str(int(error[i]))])

        else:
            c_len += len(list_dic[5][str(int(error[i]))])


    return c_len/o_len
```

%% Cell type:code id:4f6a5a0d tags:

``` python
def encoder(error, list_dic, diff, bound, bins, res):
    encoded = np.copy(error).astype(int).astype(str).astype(object)

    diff = np.reshape(diff,(510,638))
    res = np.reshape(res,(510,638))
    print(bins)

    for i in range(encoded.shape[0]):
        for j in range(encoded.shape[1]):
            if i == 0 or i == encoded.shape[0]-1 or j == 0 or j == encoded.shape[1]-1:
                encoded[i][j] = list_dic[0][encoded[i][j]]

            elif res[i-1][j-1] <= bins[0]:
                encoded[i][j] = list_dic[1][encoded[i][j]]

            elif res[i-1][j-1] <= bins[1] and res[i-1][j-1] > bins[0]:
                encoded[i][j] = list_dic[2][encoded[i][j]]

            elif res[i-1][j-1] <= bins[2] and res[i-1][j-1] > bins[1]:
                encoded[i][j] = list_dic[3][encoded[i][j]]

            elif res[i-1][j-1] <= bins[3] and res[i-1][j-1] > bins[2]:
                encoded[i][j] = list_dic[4][encoded[i][j]]
            else:
                encoded[i][j] = list_dic[5][encoded[i][j]]


    return encoded
```

%% Cell type:code id:4b65c7e9 tags:

``` python
def decoder(A, encoded_matrix, list_dic, bins):
    """
    Function that accecpts the prediction matrix A for the linear system,
    the encoded matrix of error values, and the encoding dicitonary.
    """

    the_keys0 = list(list_dic[0].keys())
    the_values0 = list(list_dic[0].values())

    the_keys1 = list(list_dic[1].keys())
    the_values1 = list(list_dic[1].values())

    the_keys2 = list(list_dic[2].keys())
    the_values2 = list(list_dic[2].values())

    the_keys3 = list(list_dic[3].keys())
    the_values3 = list(list_dic[3].values())

    the_keys4 = list(list_dic[4].keys())
    the_values4 = list(list_dic[4].values())

    the_keys5 = list(list_dic[5].keys())
    the_values5 = list(list_dic[5].values())

    error_matrix = np.zeros((512,640))

    for i in range(error_matrix.shape[0]):
        for j in range(error_matrix.shape[1]):
            if i == 0 and j == 0:
                error_matrix[i][j] = int(the_keys0[the_values0.index(encoded_matrix[i,j])])

            elif i == 0 or i == error_matrix.shape[0]-1 or j == 0 or j == error_matrix.shape[1]-1:
                error_matrix[i][j] = int(the_keys0[the_values0.index(encoded_matrix[i,j])]) + error_matrix[0][0]
            else:
                z0 = error_matrix[i-1][j-1]
                z1 = error_matrix[i-1][j]
                z2 = error_matrix[i-1][j+1]
                z3 = error_matrix[i][j-1]
                y0 = int(-z0+z2-z3)
                y1 = int(z0+z1+z2)
                y2 = int(-z0-z1-z2-z3)
                y = np.vstack((y0,y1,y2))
                difference = max(z0,z1,z2,z3) - min(z0,z1,z2,z3)
                predict = np.round(np.round(np.linalg.solve(A,y)[-1][0],1))

                if difference <= bins[0]:
                    error_matrix[i][j] = int(the_keys1[the_values1.index(encoded_matrix[i,j])]) + int(predict)
                elif difference <= bins[1] and difference > bins[0]:
                    error_matrix[i][j] = int(the_keys2[the_values2.index(encoded_matrix[i,j])]) + int(predict)
                elif difference <= bins[2] and difference > bins[1]:
                    error_matrix[i][j] = int(the_keys3[the_values3.index(encoded_matrix[i,j])]) + int(predict)
                elif difference <= bins[3] and difference > bins[2]:
                    error_matrix[i][j] = int(the_keys4[the_values4.index(encoded_matrix[i,j])]) + int(predict)
                else:
                    error_matrix[i][j] = int(the_keys5[the_values5.index(encoded_matrix[i,j])]) + int(predict)


    return error_matrix.astype(int)
```

%% Cell type:code id:280aafd3 tags:

``` python
scenes = file_extractor()
images = image_extractor(scenes)
list_dic, image, error, new_error, diff, bound, bins, predict, A, res = huffman(images[0])
encoded_matrix = encoder(np.reshape(new_error,(512,640)), list_dic, diff, bound, bins, res)

rate = compress_rate(image, res, error, diff, bound, list_dic, bins)
```

%% Output

    [6, 28, 73, 170, 8066]

%% Cell type:code id:d342f424 tags:
%% Cell type:code id:329cc11b tags:

``` python
rate = []
rate_nb = []
rate_u = []
for i in range(len(images)):
    list_dic, image, error, new_error, diff, bound, bins, predict, A_, res = huffman(images[i])
    r = compress_rate(image, res, error, diff, bound, list_dic, bins)
    rate.append(r)
    """encoding, error, image = huffman_nb(images[i])

    r = compress_rate_nb(image, error, encoding)
    rate_nb.append(r)"""

print(f"Compression rate of huffman with different bins: {np.mean(rate)}")
#print(f"Compression rate of huffman without bins: {np.mean(rate_nb)}")
```

%% Output

    Compression rate of huffman with different bins: 0.42801850527808777

%% Cell type:code id:c4242b52 tags:

``` python
reconstruct_image = decoder(A, encoded_matrix, list_dic, bins)
np.allclose(image.reshape(512,640), reconstruct_image)
```

%% Cell type:code id:a9502e22 tags:

``` python
x = np.abs(error.copy())
y = diff.copy()
plt.hexbin(x,y,cmap="rocket")
plt.colorbar()


def rel_freq(x):
    freqs = [x.count(value) / len(x) for value in set(x)]
    return freqs
print(sp.stats.entropy(rel_freq(list(diff))))

def entropy_check(x, y):
    #freq = rel_freq(list(np.ravel(o)))
    means = []
    for i in range(len(images)):
        p, d, o, e, A = predict(images,0)
        d = d.reshape((510,638))
        x = np.abs(np.ravel(e))
        y = np.ravel(d)

        mask1 = y <= 25
        x_masked1 = x[mask1]

        mask2 = y > 25
        x_masked2 = x[mask2]
        mask2 = y[mask2] <= 40
        x_masked2 = x_masked2[mask2]

        mask3 = y > 40
        x_masked3 = x[mask3]
        mask3 = y[mask3] <= 75
        x_masked3 = x_masked3[mask3]

        mask4 = y > 75
        x_masked4 = x[mask4]


        e_m1 = sp.stats.entropy(rel_freq(list(x_masked1)))
        e_m2 = sp.stats.entropy(rel_freq(list(x_masked2)))
        e_m3 = sp.stats.entropy(rel_freq(list(x_masked3)))
        e_m4 = sp.stats.entropy(rel_freq(list(x_masked4)))
        means.append([e_m1, e_m2, e_m3, e_m4])
    return np.mean(np.array(means).reshape(len(images),4), axis=0)

#print(entropy_check(x, y))
```

%% Output

    4.76234755148326


%% Cell type:code id:c0bb307b tags:
%% Cell type:code id:a2582804 tags:

``` python
new_im = image.copy().reshape(512,640)
z0, z1, z2, z3 = new_im[0,0], new_im[0,1], new_im[0,2], new_im[1,0]
y0 = -z0+z2-z3
y1 = z0+z1+z2
y2 = -z0-z1-z2-z3
b = np.vstack((y0,y1,y2))

xs = [-1,-1,-1,0]
ys = [-1,0,1,-1]
zs = [z0,z1,z2,z3]


fit= la.solve(A,b)


plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')

xlim = ax.get_xlim()
ylim = ax.get_ylim()
X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),
                  np.arange(ylim[0], ylim[1]))
Z = np.zeros(X.shape)
for r in range(X.shape[0]):
    for c in range(X.shape[1]):
        Z[r,c] = fit[0] * X[r,c] + fit[1] * Y[r,c] + fit[2]
ax.plot_wireframe(X,Y,Z, color='k')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()


def MSE(xs, ys, zs, fit):
    total = 0
    for i in range(4):
        z_ = fit[0] * xs[i] + fit[1] * ys[i] + fit[2]
        total += (zs[i] - z_)**2
    return total / 4

def distance_to_plane(xs, ys, zs, fit):
    distances = []
    for i in range(4):
        num = np.abs(fit[0]*xs[i] + fit[1]*ys[i] + fit[2])
        denom = np.sqrt(fit[0]**2 + fit[1]**2 + fit[2]**2)
        distances.append(num/denom)
    return distances

print(distance_to_plane(xs,ys,zs,fit))

points = np.array([[-1,-1,1], [-1,0,1], [-1,1,1], [0,-1,1]])
my_mse = MSE(xs, ys, zs, fit)
my_pred = points@fit.reshape(3,)
print(mean_squared_error(zs, my_pred))
print(my_mse)
```

%% Cell type:code id:671f7847 tags:
%% Cell type:code id:487fc2f2 tags:

``` python
f, res, rank, s = np.linalg.lstsq(points, np.array(zs), rcond=None)
print(f[0]*0 + f[1]*0 + f[2])
print(new_im[1,1])

plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')

xlim = ax.get_xlim()
ylim = ax.get_ylim()
X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),
                  np.arange(ylim[0], ylim[1]))
Z = np.zeros(X.shape)
for r in range(X.shape[0]):
    for c in range(X.shape[1]):
        Z[r,c] = f[0] * X[r,c] + f[1] * Y[r,c] + f[2]
ax.plot_wireframe(X,Y,Z, color='k')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
```

%% Cell type:code id:eec0746a tags:
%% Cell type:code id:b4998aef tags:

``` python
def plot_hist_lstsq(tiff_list):
def predict_pix_lstsq(tiff_list):

    image = tiff_list
    image = Image.open(image)    #Open the image and read it as an Image object
    image = np.array(image)[1:,:]    #Convert to an array, leaving out the first row because the first row is just housekeeping data
    image = image.astype(int)
    A = np.array([[3,0,-1],[0,3,3],[1,-3,-4]]) # the matrix for system of equation
    z0 = image[0:-2,0:-2]   # get all the first pixel for the entire image
    z1 = image[0:-2,1:-1]   # get all the second pixel for the entire image
    z2 = image[0:-2,2::]    # get all the third pixel for the entire image
    z3 = image[1:-1,0:-2]   # get all the forth pixel for the entire image
    # calculate the out put of the system of equation
    y0 = np.ravel(-z0+z2-z3)
    y1 = np.ravel(z0+z1+z2)
    y2 = np.ravel(-z0-z1-z2-z3)
    y = np.vstack((y0,y1,y2))

    # use numpy solver to solve the system of equations all at once
    predict = np.round(np.round((np.linalg.solve(A,y)[-1]),1)) #round the solution to the nearest integer so that encoding/decoding is easier

    points = np.array([[-1,-1,1], [-1,0,1], [-1,1,1], [0,-1,1]]) #Matrix system of points that will be used to solve the least squares fitting hyperplane



    # flatten the neighbor pixels and stack them together
    z0 = np.ravel(z0)
    z1 = np.ravel(z1)
    z2 = np.ravel(z2)
    z3 = np.ravel(z3)
    neighbor = np.vstack((z0,z1,z2,z3)).T

    f, res, rank, s = la.lstsq(points, neighbor.T, rcond=None)


    # calculate the difference
    diff = np.max(neighbor,axis = 1) - np.min(neighbor, axis=1)

    # flatten the image to a vector
    image = np.ravel(image[1:-1,1:-1])
    error = image-predict

    return image, predict, res, error, A, diff

i, p, l, e, A, d = plot_hist_lstsq(images[0])

plt.hexbin(np.abs(e), l, mincnt=1, bins="log")
plt.xlabel("Error")
plt.ylabel("Lst Sqr Residual")
plt.colorbar()
plt.show()
```

%% Output


    'plt.hexbin(np.abs(e), diff, mincnt=1, bins="log")\nplt.colorbar()\nplt.show()\n\nplt.hexbin(diff, l, mincnt=1, bins="log")\nplt.show()\n\nplt.hexbin(np.abs(e), .9*diff + .1*l, mincnt=1, bins="log")\nplt.colorbar()\nplt.show()'

%% Cell type:code id:700f6e7f tags:
%% Cell type:code id:db376cb9 tags:

``` python
il = l.astype(int)
bins = []
num_bins = 5
data_points_per_bin = len(il) // num_bins   #l is the list of data that you want to create bins for

sorted_l = il.copy()
sorted_l.sort()

bins = [sorted_l[_ * data_points_per_bin: (_+1)*data_points_per_bin] for _ in range(num_bins)]

for b in bins:
    print(np.max(b))


```

%% Cell type:code id:0c297da9 tags:
%% Cell type:code id:7575133b tags:

``` python
imm, p, res, e, A, d = plot_hist_lstsq(images[0])
res = res.astype(int)
uni = np.unique(res)
uni_n = len(np.unique(res))
entropy = []
for i in range(uni_n):
    mask = res == uni[i]
    mask_error = e[mask]
    entropy.append(sp.stats.entropy(rel_freq(list(mask_error))))
print(np.mean(entropy))
print(sp.stats.entropy(rel_freq(list(imm))))
print(len(bin(imm[0])[2:]))
```

%% Output

    3.0392286797151855
    5.304096657944529
    15

%% Cell type:code id:d7fc288d tags:
%% Cell type:code id:dcc26973 tags:

``` python
fre = rel_freq(list(res))
print(np.array(fre)@np.array(entropy))
print(3.93012/15)
def predict_pix(tiff_image, difference = True):
    """
    This function predict the pixel values excluding the boundary.
    Using the 4 neighbor pixel values and MSE to predict the next pixel value
    (-1,1) (0,1) (1,1)  => relative position of the 4 other given values
    (-1,0) (0,0)        => (0,0) is the one we want to predict
    take the derivative of mean square error to solve for the system of equation
    A = np.array([[3,0,-1],[0,3,3],[1,-3,-4]])
    A @ [a, b, c] = [-z0+z2-z3, z0+z1+z2, -z0-z1-z2-z3] where z0 = (-1,1), z1 = (0,1), z2 = (1,1), z3 = (-1,0)
    and the predicted pixel value is c.

    Input:
    tiff_image (string): path to the tiff file

    Return:
    image   (512 X 640): original image
    predict (325380,): predicted image excluding the boundary
    diff.   (325380,): IF difference = TRUE, difference between the min and max of four neighbors exclude the boundary
                       ELSE: the residuals of the four nearest pixels to a fitted hyperplane
    error   (325380,): difference between the original image and predicted image
    A       (3 X 3): system of equation
    """
    image = Image.open(tiff_image)    #Open the image and read it as an Image object
    image = np.array(image)[1:,:]    #Convert to an array, leaving out the first row because the first row is just housekeeping data
    image = image.astype(int)
    print(image.shape)
    # use
    A = np.array([[3,0,-1],[0,3,3],[1,-3,-4]]) # the matrix for system of equation
    # where z0 = (-1,1), z1 = (0,1), z2 = (1,1), z3 = (-1,0)
    z0 = image[0:-2,0:-2]   # get all the first pixel for the entire image
    z1 = image[0:-2,1:-1]   # get all the second pixel for the entire image
    z2 = image[0:-2,2::]    # get all the third pixel for the entire image
    z3 = image[1:-1,0:-2]   # get all the forth pixel for the entire image
    # calculate the out put of the system of equation
    y0 = np.ravel(-z0+z2-z3)
    y1 = np.ravel(z0+z1+z2)
    y2 = np.ravel(-z0-z1-z2-z3)
    y = np.vstack((y0,y1,y2))
    # use numpy solver to solve the system of equations all at once
    #predict = np.floor(np.linalg.solve(A,y)[-1])
    predict = np.round(np.round((np.linalg.solve(A,y)[-1]),1))

    #Matrix system of points that will be used to solve the least squares fitting hyperplane
    points = np.array([[-1,-1,1], [-1,0,1], [-1,1,1], [0,-1,1]])

    # flatten the neighbor pixlels and stack them together
    z0 = np.ravel(z0)
    z1 = np.ravel(z1)
    z2 = np.ravel(z2)
    z3 = np.ravel(z3)
    neighbor = np.vstack((z0,z1,z2,z3)).T

    if difference:
        # calculate the difference
        diff = np.max(neighbor,axis = 1) - np.min(neighbor, axis=1)

    else:
        #Compute the best fitting hyperplane using least squares
        #The res is the residuals of the four points used to fit the hyperplane (summed distance of each of the
        #points to the hyperplane), it is a measure of gradient
        f, diff, rank, s = la.lstsq(points, neighbor.T, rcond=None)

    # calculate the error
    error = np.ravel(image[1:-1,1:-1])-predict

    return image, predict, diff, error, A
```

%% Output

    3.9301269429408086
    0.262008
+0 −717

File deleted.

Preview size limit exceeded, changes collapsed.

+0 −280

File deleted.

Preview size limit exceeded, changes collapsed.

Loading