Slides

Variable
- Name: begin with letter, case sensitive, cannot use reserved words, built-in functions
- who: show defined variables, whos show variable information, clear clear workspace variables, clc clear window clear x: clear variable x, clf clear figure window
- By default, doubles are displayed with four decimal places
- Use format xxx to change format format = format short 4 decimals, format long 15 decimals, format loose/compact
- Use ; at the end of command to suppress output
- Different types: float(single, double), integer(int8, int16, ...), unsigned integer(uint8, ...), complex(complex single/double). Type conversion may cause overflow, and overflow is handled using saturation
- Numbers are stored by default as double
- Characters are represented using single quotes
Random numbers：the seed remain the same for each time Matlab is started, unless you change the seed
Input and output
- Input as an integer: x = input('Message'), input as a character/string: x = input('Message', 's')
- Output on screen: disp(x), disp function does not allow for formatting output
- Format output: fprintf('format_string', val), use format string to format output, %d integer, %f float, %c single character, %s string
File Input/Out
- save xxx.mat var1 var2 --ascii -append把变量var1和var2按照ascii追加保存到xxx.mat中
- load xxx.mat读取xxx.mat文件
Function
- Function name should be the same as the file name
Some built-in functions
- sum(x), cumsum(x), prod(x), cumprod(x)
- tic, toc will display time elapse in milliseconds

Cell Arrays, Structures, Strings

Cell arrays: store values of different types, commonly used to store strings of different lengths
- Create a cell array: x = cell(2,2)
- Indexing cell arrays: x{r, c}
- Compare if two cell arrays are equal: cellfun(@isequal, x, y) returns a logical matrix, then we can use sum(cellfun(@isequal, x, y), 'all') == numel(x)
- celldisp(x) display all elements in cell array x
Structures: group values that are logically related, but are not the same thing and not necessarily the same type
- Consist of fields, each field contains an array of some Matlab type
- x = struct('key1', val1, 'key2', val2)
- x.key1 get val1
- rmfield(x, 'key1') remove filed key1 from structure x
- isstruct(x) check if x is a structure
- isfield(x, 'key1') check if key1 is a field of structure x
- fieldnames(x) find all field names of structure x, return a cell array

Strings

Strings are treated the same as single-character arrays

x = ['a', 'bc', 'd']; 
[r, c] = size(x); % r = 1, c = 4
disp(x); % 'abcd'
x = ['ab'; 'c']; % An error occured, the length should be the same for each row
x = {'ab', 'c'; 'def', 'g'}; % Use cell arrays to hold strings with different length

x =  {'a', 'bc'};
q1 = char(x); % 2 * 2 char array: 'a ', 'bc' (there is an space after a so that row is length 2)
q2 = cellstr(x); % 1 * 2 cell array: {'a'}, {'bc'}

strcat(s1, s2, ..., sn)concatenate n strings
strcmp(s1, s2) compare two strings

Permutation and combination
- factorial(n) \(n!\)
- prod(n:-1:n-m+1) (n > m) \(P_{n}^{m}=\frac{n!}{(n-m)!}\)
- nchoosek(n, k): \(C_{n}^{k}\)
Vector
- Norm Properties
  - Homogeneous: p(au) = |a| * p(u)
  - Subadditive: p(u + v) <= p(u) + p(v)
  - Positive: p(u) >= 0
  - Zero vector: p(u) = 0 <=> u = 0
- \(|X|_{p} = (\sum\limits_{i = 1}^{n}|x_i|^{p})^{1/p}\)
- \(x\cdot y = |x||y|cos\theta\)
- Orthogonal matrix: \(AA^{T} = A^{T}A = I\)
- Symmetrix matrix: \(A = A^{T}\)
- Skew-symmetric matrix: \(A^{T} = -A\)
Dimensionality Reduction
- Assumpation: data lies on or near a low dimensional subspace, axes of this subspace are effective representation of the data
- Rank is dimensionality
- Why reduce dimensions
  - Discover hidden correlations/topics
  - Remove redundant and noisy features
  - Interpretation and visualization
  - Easier storage and processing data
- SVD Decomposition(singular value decomposition)
  - \(A_{m*n} = U_{m*r}\Sigma_{r*r}(V_{n*r})^{T}\)
  - U, \(\Sigma\), V are unique. U, V are column orthonormal. \(\Sigma\) is diagonal, singular values are positive are are sorted in descending order
  - Choose the first several large singular values to reduce dimensionality
  - Time complexity of SVD: \(min\{O(nm^{2}), O(n^{2}m)\}\)
- CUR Decomposition
  - \(A = CUR\)
  - When A is huge and sparse. In SVD, \(\Sigma\) is sparse and small, U and V are big and dense. In CUR, U is dense and small, C and R are big and sparse
  - CUR is lesss accurate, but faster than SVD
- SVD is limisted to linear projectsions, for non-linear methods, use isomap

PCA

PCA is a statistical technique for reducing the dimensionality of a dataset while preserving the maximum amount of information
One-dimensional PCA: for \(x_1, ..., x_n\), find line \(x(t) = tv + b\) defined by vector v and b such that 1D projection \(a_i = v^{t}(x_i - b)\) has the largest possible variance.
- The solution is not unique, since there can be parallel lines(different b)
- Let \(b = \frac{1}{n}\sum\limits_{i=1}^{n}{x_i}\)
- Find v, s.t. \(|v| = 1, v = argmax\sum {a_i}^2, a_i = v^{t}(x_i - \bar{x})\)
- \(v = argmax \;v^{t}Cv\), where C is square, symmetric, and positive semidefinite(nonnegative eigenvalues)
- Theorem: Given set of points \(x1,...,x_n\), the optimal direction for projecting the data is the largest v is the largest eigenvector of the covariance matrix \(C = \sum{(x_i-\bar{x})(x_i - \bar{x})^{T}} = X^{T}X, X \in \mathbb{R}^{n*d}\)
- Do SVD: \(X = U\Sigma V^{T}\), then \(C = V(\Sigma^{T}\Sigma V^{T}x\)

One example

% rows of x are observations, columns of x are dimensions to be reduced
% X is m * n: m obs, n dims
x = [80 85 60 55
    90 85 70 45
    95 80 40 50];
% rows of coef are dimensions, columns of coef are principals
% coef is n * k: we want to reduce dimensions from n to k
% data in coef is principal component coefficients: the coefs 
% to transform from the dimension space to the principal component space

% rows of score are observations, columns of score are principal components
% score is: m * k, m obs, k pcs
% data in score is obs values in the principal component space
[coef, score] = pca(x); 

% An equivalent way to get the scores
colmean = mean(x, 1); % calculate column mean
meanmat = repmat(colmean, size(x, 1), 1); % generate mean matrix for x
xtilde = x - meanmat; % center obs values by subtracting column means
myscore = xtilde * coef; % translate obs values in the dim space to the pc space

For non-linear scenarios, use manifold learning, isomap

Linear programming

Optimization
- Criterion: object function and directino(max/min, in Matlab, the direction is always min)
- Variables
- Constraint(solution space)
Linear programming, the degree of freedom is 1 and the solution space is convex: \(x, y \in S\), \(\forall t \in (0,1),\; tx + (1-t)y \in S\)
Use x = linprog(f,A,b) so solve \(x = argmin \{f\}, s.t. Ax \le b\)
Use x = linprog(f,A,b,Aeq,Beq) to solave \(x = argmin \{f\}, s.t. Ax \le b, Aeqx = beq\)

Use x = linprog(f,A,b,Aeq,Beq,lb,ub) to solve \(x = argmin \{f\}, s.t. Ax \le b, Aeqx = beq, lb \le x \le ub\)

%  max x1 + x2/3
%  x1 + x2 <= 2
%  x1 + x2/4 <= 1
%  x1 - x2 <= 2
% -x1/4 - x2 <= 1
% -x1-x2 <= -1
% -x1 + x2 <= 2
% x1 + x2/4 = 1/2
% -1 <= x1 <= 3/2
% -0.5 <= x2 <= 5/4
A = [1, 1;
    1, 1/4;
    1, -1;
    -1/4, -1;
    -1, -1;
    -1, 1];
b = [2, 1, 2, 1, -1, 2];
Aeq = [1, 1/4];
beq = 1/2;
lb = [-1, -1/2];
ub = [3/2, 5/4];
f = [-1, -1/3];
x = linprog(f,A,b,Aeq,beq,lb,ub);

The assignment problem

Description: Assignment problem
Solution: Hungarian algorithm

Matlab: Hungarian Algorithm in Matlab

clear;clc;
costmat = [14,5,8,7;
    2,12,6,5;
    7,8,3,9;
    2,4,6,10];
[assign,cost] = munkres(costmat);
disp(assign);  % 2 4 3 1, m1 assigned to t2, m2 assigned to t4, m3 assigned to t3, m4 assigned to t1 
disp(cost);    % 15

Basic image processing
- Show image information: info = imfinfo(filename)
- Read input image: img = imread(filename)
- Write output image: imwrite(img, filename, fmt, 'Quality', qualityvalue) refer to the documentation for more arguments
- Rgb to grayscale: img_gray = rgb2gray(img)
- Grayscale to rgb without colormap: img_rgb = cat(3, img_gray, img_gray, img_gray);
- Grayscale to binarize: img_bin = imbinarize(img_gray), can add threshold, refer to the documentation for more details
- Save histogram of an image
  1
  2
  3
  4
  5
  figure;
  imhist(img_gray);
  title('Histogram of the original image');
  % gcf(get current figure)
  saveas(gcf, "histogram.png");
- Histogram equalization of an grayscale image img_heq = histeq(img_gray)
- Adjust image intensity: J = imadjust(I), this increase the contrast, by default, out is [0.01, 0.99]
- Image enhancement: J = imadjust(I,[low_in high_in],[low_out high_out],gamma)
  - If gamma < 1: nonlinear mapping, weighted toward higher(brighter) output values
  - If gamma = 1: linear mapping, default value of gamma
  - If gamma > 1: nonlinear mapping, weighted toward lower(darker) output values
- Find limits to contrast stretch image: lowhigh = stretchlim(I), often used with imadjust: J = imadjust(I,stretchlim(I), []);
- Convert img to double precision: I2 = im2double(I)
Point processing
- Negative/inverted image: J = 255 - I, or J = imadjust(I, [0, 1], [1, 0])
- Darken image: J = I - 128
- Linearly lower contrast: J - I / 2
- Nonlinearly lower contrast: J1 = im2double(I); J = J1.^(1/3)
- Lighten image: J = I + 128
- Linearly raise contrast: J = I * 2
- Nonlinearly raise contrast: J1 = im2double(I); J = J1.^3

Image Pyramid

clear;clc;close all;
img = imread('cameraman.tif');
img = im2double(img);
pyramid = cell(1,4);
residual = cell(1,4);
pyramid{1} = img;
for i = 2:numel(pyramid)
    tmp = imgaussfilt(pyramid{i - 1});
    residual{i - 1} = pyramid{i - 1} - tmp;
    pyramid{i} = imresize(tmp, 0.5,"lanczos2");
end

tmp = imgaussfilt(pyramid{end});
residual{end} = pyramid{end} - tmp;

figure;

[m0, n0] = size(pyramid{1});
for i = 1:numel(pyramid)
    [m, n] = size(pyramid{i});
    if i > 1
    tmp = pyramid{i};
    pyramid{i} = padarray(tmp, [m0 - m, n0 - n]);
    end
    subplot(2,2,i)
    imshow(pyramid{i});
    if i == 1
        title_str = strjoin(["Original image of size [", m, 'x', n, ']' ]);
    else
        title_str = strjoin(["Filtered image of size [", m, 'x', n, ']' ]);
    end
    title(title_str);
end

figure;
for i = 1:numel(pyramid)
    [m, n] = size(residual{i});
    if i > 1
    tmp = residual{i};
    residual{i} = padarray(tmp, [m0 - m, n0 - n]);
    end
    subplot(2,2,i)
    imshow(residual{i});
    title_str = strjoin(["Residual of size [", m, 'x', n, ']' ]);
    title(title_str);
end

Spatial 2d image filter
- Box filter(separable)
  - Example \[ \frac{1}{9}\begin{bmatrix} 1 & 1 & 1\\ 1 & 1 & 1\\ 1 & 1 & 1\end{bmatrix}\]
  - Matlab code
    1
    2
    3
    4
    % Create 3 * 3 box average filter
    box = fspecial("average");
    %Filtered image
    img_box = imfilter(img, box);
- Gaussian filter(separable)
  - Example \[ \frac{1}{16}\begin{bmatrix} 1 & 2 & 1\\ 2 & 4 & 2\\ 1 & 2 & 1\end{bmatrix}\]
  - Matlab code
    1
    2
    3
    4
    5
    % Create 3 * 3 gaussian filter with sigma = 0.5
    gau = fspecial("gaussian");
    img_gau = imfilter(img, gau);
    % This way is recommended
    img_gau2 = imgaussfilt(img);
- Laplacian of gaussian filter
  1
  2
  3
  % Create 5 * 5 Laplacian of Gaussian filter with sigma = 0.5
  logau = fspecial("log");
  img_log = imfilter(img, logau);
- Derivative of gaussian filter
  1
  2
  3
  4
  H1 = fspecial('gaussian', 4, 1);
  H2 = fspecial('gaussian', 4, 4);
  dog = H2 - H1;
  img_dog = imfilter(rgb2gray(img), dog);
- Laplacian filter
  1
  2
  3
  % Create 3 * 3 laplacial filter with alpha = 0.2
  lap = fspecial("laplacian");
  img_lap = imfilter(img, lap);
- Other filters
  - Unchanged \[\begin{bmatrix} 0 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 0\end{bmatrix}\]
  - Shift to left by one \[\begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0\end{bmatrix}\]
  - Sharpening: do nothing for flat areas, stress intensity peaks \[\begin{bmatrix} 0 & 0 & 0\\ 0 & 2 & 0\\ 0 & 0 & 0\end{bmatrix} - \frac{1}{9}\begin{bmatrix} 1 & 1 & 1\\ 1 & 1 & 1\\ 1 & 1 & 1\end{bmatrix} \]
- Sober filter(separable)
  - Vertical sober filter
    - Example \[\begin{bmatrix} 1 & 2 & 1\\ 0 & 0 & 0\\ -1 & -2 & -1\end{bmatrix}\]
    - Matlab cod
      1
      2
      3
      % Create 3 * 3 vertical sober filter that emphasizes horizontal edges
      vs = fspecial("sobel");
      img_vs = imfilter(img, vs);
  - Horizontal sober filter
    - Example \[\begin{bmatrix} 1 & 0 & -1\\ 2 & 0 & -2\\ 1 & -2 & -1\end{bmatrix}\]
    - Matlab cod
      1
      2
      3
      % Create 3 * 3 horizontal sober filter that emphasizes vertical edges
      hs = fspecial("sobel")';
      img_hs = imfilter(img, hs);

Frequency basics

Given \(f(x,y)\) on \(x = 0,...,M - 1\) and \(y = 0,...,N-1\), DFT \(F(u,v) = \sum\limits_{x = 0}^{M - 1}\sum\limits_{y=0}^{N - 1}f(x,y)e^{2\pi i (ux/M + vy/N)}\) for \(u=0,...,M-1\) and \(v=0,...,N-1\). IDFT of \(F(u,v)\) is \(f(x,y) = \frac{1}{MN}\sum\limits_{u = 0}^{M - 1}\sum\limits_{v=0}^{N - 1}F(u,v)e^{2\pi i (ux/M + vy/N)}\)
Fourier spectrum \(|F(u,v)| = [R^{2}(u,v) + I^{2}(u,v)]^{\frac{1}{2}}\), phase angle \(\phi(u,v) = arctan[\frac{I(u,v)}{R(u,v)}]\)
If \(f(x,y)\) is real, then \(F(u,v) = F(-u, -v)\), and \(|F(u,v)| = |F(-u,-v)|\)
FFT of f: F = fft2(f), IFFT of F: f = ifft2(F). If f is real,then the result of ifft2 should be real, but for some MATLAB versions this is not the case, so you can use f = real(ifft2(F))
Move the origin of the transform to the center of the frequency rectangle: Fc = fftshift(F), reverse centering: F = ifftshift(Fc)
Get the phase angle: phi = atan2(imag(F), real(F)), or phi = angle(F)
Fourier spectrum: S = abs(F), F = S.*exp(i*phi)

Visualize the Fourier spectrum

img = imread('q1_origin.png');

% Convert image to grayscale
img_gray = rgb2gray(img);

% Compute Fourier transform of the image
F = fft2(img_gray);

% Move the origin of the transform to the center of the frequency rectangle
Fc = fftshift(F);

% Compute Fourier spectrum
F_mag = abs(Fc);

% Increase in visualization details
F_log = log(1 + F_mag);
F_log = F_log/max(F_log(:));

figure;
imshow(F_log, []);
title("Fourier spectrum");

General steps in DFT filtering
1. Convert to grayscale f = rgb2gray(f);
2. Convert the image to floating point [f, revertclass] = tofloat(f);
3. Obtain size [N, N] = size(f)
4. Obtain FFT F = fft2(f)
5. Generate a filter H, of size M * N, if the filter is centered instead, let H = ifftshift(H) before using the filter. But when you need to show the image, you need to show the centered one, use H = fftshift(H) if H is not centered
6. Multiply the transform by the filter G = H .*F
7. Obtain the inverse FFT of G, ie, the filtered image g = ifft2(G)
8. Convert the filtered image to the class of the input image g = revertclass(g)

Frequency gaussian filter

img = imread("filter_input.jpg");
img = rgb2gray(img);
%  tofloat.m, lpfilter.m(requires dftuv.m) are from the textbook
[f, revertclss] = tofloat(img);
[M, N] = size(f);
F = fft2(f);

% Generate gaussian filter
std = 20;
filter_gaussian = lpfilter('gaussian', M, N, std);
% Center the filter for display
Fc_gaussian = fftshift(filter_gaussian);

G_gaussian = filter_gaussian.*F;
% This is the filtered image
gaussian_filtered_img = ifft2(G_gaussian);
gaussian_filtered_img = revertclss(gaussian_filtered_img);

Frequency box filter

% Generate box filter of the same size of image
filter_box = zeros(M, N);
for i = 1:M
    for j = 1:N
        cond1 = i >= M / 2 - 30 && i <= M / 2 + 30;
        cond2 = j >= N / 2 - 30 && j <= N / 2 + 30;
        if cond1 && cond2
            filter_box(i, j) = 1;
        end
    end
end
% The created box filter is already centered
Fc_box = filter_box; 
% The filter is centered so ifftshift is needed
filter_box = ifftshift(filter_box);

G_box = filter_box.*F;
g_box = real(ifft2(G_box));
g_box = revertclss(g_box);

Frequency lowpass/highpass filter

% Generate lowpass/highpass filter
% This uses lpfilter.m
std = 5;
filter_lp = lpfilter('ideal', M, N, std);
Fc_lp = fftshift(filter_lp);

G_lp = filter_lp.*F;
g_lp = ifft2(G_lp);
g_lp = revertclss(g_lp);

filter_hp = 1 - filter_lp;
Fc_hp = fftshift(filter_hp);
G_hp = filter_hp.*F;
g_hp = ifft2(G_hp);
g_hp =revertclss(g_hp);

Frequency bandpass/bandreject filter

% Generate bandpass/bandreject filter
% This uses bandfilter.m
rad = 50;
w = 30;
filter_bp = bandfilter('ideal', 'pass', M, N, rad, w);
Fc_bp = fftshift(filter_bp);

G_bp = filter_bp.*F;
g_bp = ifft2(G_bp);
g_bp = revertclss(g_bp);

filter_br = 1 - filter_bp;
Fc_br = fftshift(filter_bp);
G_br = filter_br.*F;
g_br = ifft2(G_br);
g_br = revertclss(g_br);

Frequency notchpass/notchreject filter

% Generate notchpass/notchreject filter
% This uses cnotch.m, iseven.m, isodd.m
% You should change C1 at your needs
w = 4;
C1 = zeros(M, 2);
for i = 1:M
    C1(i,1) = i;
    C1(i, 2) = 247;
end
C1(260:290,:) = [];
filter_np = cnotch('gaussian','pass', M, N, C1, w);
Fc_np = fftshift(filter_np);

G_np = filter_np.*F;
g_np = ifft2(G_np);
g_np = revertclss(g_np);


filter_nr = 1 - filter_np;
Fc_nr = fftshift(filter_nr);

G_nr = filter_nr.*F;
g_nr = ifft2(G_nr);
g_nr = revertclss(g_nr);

Image transformation

Translation

Transformation matrix \[\begin{bmatrix} 1 & 0 & t_x\\ 0 & 1 & t_y\\ 0 & 0 & 1\end{bmatrix}\]

Matlab code

img = imread("q1_origin.png");
% Get spatial referencing information 
img_ref = imref2d(size(img));

% Set tx and ty, and generate the transformation matrix
tx = 15;
ty = 20;
A = [1, 0, tx;
    0, 1, ty;
    0, 0, 1];

tform = transltform2d(A);
% Adjust the spatial referencing information of the translated image
translated_ref = img_ref;
translated_ref.XWorldLimits(2) = translated_ref.XWorldLimits(2) + tx;
translated_ref.YWorldLimits(2) = translated_ref.YWorldLimits(2) + ty;
[translated_img, translated_ref] = imwarp(img, tform,OutputView=translated_ref);

% Display the original and translated images
figure;
subplot(1,2,1)
imshow(img, img_ref)
subplot(1,2,2)
imshow(translated_img, translated_ref);
sgtitle('\color{red} Original image(left), translated image(right)');

Euclidean(rigit): rotation + translation

Transformation matrix \[\begin{bmatrix} cos(\theta) & -sin(\theta) & t_x\\ sin(\theta) & cos(\theta) & t_y\\ 0 & 0 & 1\end{bmatrix}\]

Matlab code

img = imread("q1_origin.png");

% Set theta, tx and ty to generate the transformation matrix
theta = 30;
tx = 10;
ty = 20;
A = [cosd(theta), -sind(theta), tx;
    sind(theta), cosd(theta), ty;
    0, 0, 1];

tform = rigidtform2d(A);
translated_img = imwarp(img, tform);

figure;
imshowpair(img, translated_img, "montage");
title('\color{red}Original Image(left), translated image(right)');

Similarity: uniform scaling + rotation + translation

Transformation matrix \[\begin{bmatrix} s*cos(\theta) & -s*sin(\theta) & t_x\\ s*sin(\theta) & s*cos(\theta) & t_y\\ 0 & 0 & 1\end{bmatrix}\]

Matlab code

img = imread("q1_origin.png");

% Set scale, theta, tx and ty to generate the transformation matrix
scale = 2;
theta = 30;
tx = 10;
ty = 20;
A = [scale * cosd(theta), -scale * sind(theta), tx;
    scale * sind(theta), scale * cosd(theta), ty;
    0, 0, 1];

tform = simtform2d(A);
translated_img = imwarp(img, tform);

Affine: uniform scaling + shearing + rotation + translation

Transformation matrix \[A = s * \begin{bmatrix} cos(\theta) & -sin(\theta) \\ sin(\theta) & cos(\theta) \end{bmatrix} * \begin{bmatrix} 1 & h_1\\ h_2& 1\end{bmatrix} \] \[\begin{bmatrix} A(1,1) & A(1,2) & t_x\\ A(2,1) & A(2,2) & t_y\\ 0 & 0 & 1\end{bmatrix}\]

Matlab code

img = imread("q1_origin.png");

% Set scale, theta, tx, ty h1, h2 to generate the transformation matrix
tx = 10;
ty = 20;
scale = 3;
theta = 30;
h1 = 3;
h2 = 2;
r = [cosd(theta), -sind(theta);
    sind(theta), cosd(theta)];
h = [1, h1;
    h2, 1];
top = scale * r * h;
A = [0, 0, tx;
    0, 0, ty;
    0, 0, 1];
A(1:2, 1:2) = top;


tform = affinetform2d(A);
translated_img = imwarp(img, tform);

Projection: affine transformation + projective wrap
- Transformation matrix \[\begin{bmatrix} a_{11} & a_{11} & a_{11}\\ a_{21} & a_{22} &a_{23}\\ a_{31} & a_{32} & a_{33}\end{bmatrix}\]
- Matlab code
  1
  2
  3
  4
  5
  6
  7
  8
  9
  img = imread("q1_origin.png");
  
  % Set the transformation matrix
  A = [2 0 0;
  0.33 1 0;
  0 0 1];
  
  tform = affinetform2d(A);
  translated_img = imwarp(img, tform);

Hough

Line fitting: (x,y) => (m,b) / (\(\theta, \rho\))
Hough transformation is a generic framework for detecting a parametric model
A line becomes a point, a point becomes a line

Matlab code

% Reference for hough: https://www.mathworks.com/help/images/ref/hough.html
img = imread("filter_input.jpg");

img = rgb2gray(img);

% If method is not specified, the default one is "Sobel"
% Reference for edge:https://www.mathworks.com/help/images/ref/edge.html#d124e102512
% Default method for edge is sobel,
BW = edge(img);
[H,theta,rho] = hough(BW);
subplot(2,1,1);
imshow(img);
title('Original image');
subplot(2,1,2);
% rescale(H) rescales H to [0,1]
imshow(imadjust(rescale(H)),'XData',theta,'YData',rho,'InitialMagnification','fit');
title('Hough transform of original image');
xlabel('\theta'), ylabel('\rho');
axis on, axis normal, hold on;
% Add color map
colormap(gca,hot);

1-D digital filter

Rational transfer function: \(Y(z) = \frac{b_1 + b_2 z^{-1}+...+b_{n_b + 1}z^{-n_{b}}}{1 + a_2 z^{-1} + ...+a_{n_a + 1}z^{-n_a}}X(z)\)
y = filter(b, a, x) filters the input data x using a rational transfer function defined by the numerator and denominator coefficients b and a
- If x is a vector, then filter returns the filtered data as a vector of the same size as x
- If x is a matrix, then filter acts along the first dimension and returns the filtered data for each column
- If x is a multidimensional array, then filter acts along the first array dimension whose size does not equal 1

Example

% Reference: https://www.mathworks.com/help/matlab/ref/filter.html#buagwwg-2
t = linspace(-pi,pi,100);
rng default  %initialize random number generator
x = sin(t) + 0.25*rand(size(t));

windowSize = 5; 
b = (1/windowSize)*ones(1,windowSize);
a = 1;

y = filter(b,a,x);

plot(t,x);
hold on
plot(t,y);
legend('Input Data','Filtered Data');

Matlab Basics

Ctrl + R注释, Ctrl + T取消注释
Matlab使用三个...表示连接下一行
特殊值：i, j 虚数单位; eps浮点数相对精度; Inf/inf无穷大; intmax/intmin整数范围; realmax/realmin实数范围; pi; NaN/nan(NaN参与的结果都是NaN);
复数: x = a + bi, real(x) = a, imag(x) = b
求n次根号不要用幂，幂只能得到第一象限的根，用nthroot(a, n)
Matlab中标量是1*1数组，空数组是有一个维度为0，空数组用[]表示，可以把[]赋给值表示删除
数组结构: ndims(x)x的维度，size(x)x的大小，size(x, nd)x在第nd维度上的大小，length(A)所有维度的最大值，numel(x)x中总元素个数
自变量数组创建 a:inc:b初始，step，结束，并且abs(end) < abs(b), inc可以是负数; linspace(a,b,n)初始，end，个数
二维数组创建: 行内空格或者逗号，行间分号或者回车；中等二维数组使用Matlab的New Variable功能创建，；大的二维数组使用.m文件创建
一些内置数组函数
- diag(x)生成对角矩阵，或者取矩阵的对角线
- ones(a, b)全1矩阵
- eye(n)n阶单位阵
- zeros(a, b)全0矩阵
- magic(n)n阶数组，行列对角线的和都一样，\(n > 2\)
- rand(a, b)a*b阶元素在(0,1)均匀分布矩阵
- randn(a, b)a*b阶标准正态分布矩阵
- randi(max, n)1:max均匀分布的n * n随机整数矩阵, randi([min, max], m, n)生成在[min, max]内均匀分布的m*n阶随机矩阵
一些操作函数: permute(A, [dims])按照A维度index的顺序重排A, repmat(A, x, y)把A当单位，生成xy阶单位的矩阵; reshape(A, x, y)改变行列数; flipud(A)水平中心为轴，上下翻转; fliplr(A)垂直中心为轴，左右翻转; rot90(A, n)逆时针旋转A角度90 n
数组下标: 二维数组默认按列填充，序号坐标按列填充。转换函数[row, col] = ind2sub(size(A), idx), idx = sub2ind(size(A), row, col)。end可以代表最后一行，最后一列，或者最后一个
数组访问: A(r, c), A(r, :), A(:, c), A(idx), A(:)访问所有，返回列向量。注意序号访问返回结果和序号相同: A([idx])返回行向量, A([idx]')返回列向量
数组运算: 数组之间是对应元素的运算，数组和标量之间是标量和数组每一个元素的运算
逻辑运算符: A&B(A&&B is short-circuiting), A|B(A||B is short-circuiting), ~A, xor(A,B), all(A)是否全部非零, any(A)是否存在非零
初等函数都是对数组每一个元素进行运算的
矩阵乘除法
- a*B标量和每一个矩阵元素相乘
- A*B矩阵乘法
- X = A\B给出方程AX = B的解
- X = B/A给出方程XA = B的解
- x = A\b给出方程Ax = b, row(A) > col(A)的最小二乘解, b是列向量
- x = A\b给出方程Ax = b, row(A) < col(A)的最小二乘基础解, b是列向量
矩阵的标量特征参数: rank(A), trace(A), det(A), inv(A)逆矩阵
特征值和特征向量
- vals = eig(A)返回列向量特征值, vals = eig(A, 'matrix')返回矩阵形式的特征值矩阵
- [V, D] = eig(A)，其中V是特征向量矩阵，每一列是一个特征向量，D是特征值矩阵，D除了对角线是特征值外都是0。A * V = V *D
除法解方程: Ax = b => x = A\b, xC = d => x = d/C
二项分布: pk = binopdf(k, N, p)一共N次发生k次的概率, Fk = binocdf(k, N, p)一共N次发生不大于k次的总概率, E(k)=Np, D(k)=Npq。r = binornd(N, P, m, n)产生m*n服从B(N,p)的随机数组
正态分布: px = normpdf(x, mu, sigma)取值x的概率密度, Fx = normcdf(x, mu, sigma)不大于x的总概率, r = normrnd(mu, sigma, m, n)产生m*n服从N(mu, sigma^2)的随机数组。E(x)=mu, D(x)=sigma^2
统计分析命令:
- max(X)/min(X)每一列的最大/最小值, max(x, [], 'all')求矩阵最大值
- mean(X)每一列的平均值, mean(X, 1)每一列平均值, mean(X, 2)每一行平均值, mean(X, 'all')矩阵平均值
- std(x)每一列标准差, std(x, 0, 1)每一列标准差, std(x, 0, 2)每一行标准差, std(x, 0, 'all')全部的标准差
- var(x)每一列方差, var(x, 0, 1)每一列方差, var(x, 0, 2)每一行方差, var(x, 0, 'all')全部的方差
- cov(x)不同列的协方差, diag(cov(x))' = var(x)
- corrcoef(x)不同列的协相关系数
Matlab绘制连续曲线的时候, 会根据用户指定的离散点，自动进行插值计算，从而绘制出连续的曲线

基本plot

plot(x, y, key, val)y后面是参数，比如Color, LineStyle, LineWidth, Marker, MarkerSize, MarkerEdgeColor, MarkerFaceColor
plot(x)当x是一维数组的时候，下标是x轴，x是y轴。当x是二维数组的时候，行是x轴坐标1:n，列是y轴，画出列数条曲线
先画图，再设置其他参数
axis([x1, x2, y1, y2])二维绘图坐标轴范围
grid on 画出分格线, box on坐标呈现封闭形式, title(xxx), xlabel(xxx), ylabel(xxx), text(x, y, xxx)在(x, y)处写文字, legend(s1, s2)图例
hold on同一张图不刷新，从而在一张图上面画多个曲线
subplot(m, n, k)在mn的图的第k个子图, subplot('position', [left, bottom, width, height]),图总体是[0,1] [0,1]并且原点是左下角,这个可以指定位置画图，比如一个图占一行/一列
plot3(x, y, z)三维画图，画曲线

三维曲面图流程

x = linespace(x1,x2,n);
y = linespace(y1,y2,m);
[X, Y] = meshgrid(x, y);
Z = fun(X., Y.) % 这里注意数组运算需要加.
surf(X, Y, Z); % 三维曲面图
hold one;
stem3(X, Y, Z); % 杆状图，指向曲面上的点

其他画图
- bar(x)x是一维数组的画, x轴是idx, y轴是x值。x是二维数组的话是分组直方图，每一列是一组。bar(x, 'stacked')对二维数组画分组直方图。
- hist(x)x是一维数组, x轴是idx, y轴是x值。x是二维数组的话是分组频率图，每一列是一组。hist(x, nbins)
- pie(x, explode, labels)x是一维数组, explode是扇形间距, labels是cell标签(大括号)
控制流
- if expr stat elseif expr stat else stat end中的表达式可以是逻辑数组，这时候只有全部是逻辑1才是true从而执行对应的statement
- switch expr case val stat otherwise stat end这里不用在case最后写break，但是执行中会break。这里的expr是标量数值或者标量字符串，数值用==进行比较，字符串用strcmp(s1, s2)进行比较
- for idx = arr stat end, while expr stat end

Assignments

Solve Ax = b

Use x = A\b, where b is a column vector

A = [1, 1;
    2,4];
b = [30; 80];
res = A\b;

Use symbolic function

syms x y
eq1 = x + y == 30;
eq2 = 2 * x + 4 * y == 80;
[resx, resy] = solve([eq1, eq2], [x, y]);

Symbolic function value

syms x
f = x^2 + 2 * x + 1;
x = 2;
subs(f); % 9

In triangle ABC, \(\frac{a}{sinA} = \frac{b}{sinB} = \frac{c}{sinC}\), \(c^2 = a^2 + b^2 - 2abcosC\)
- sin(x) for values, sind(x)for degrees, asin(x) for values, asind(x) for degrees

My Blog

COSI 177A Scientific Data Processing in Matlab

Slides

Matlab Basics

Assignments