add all projects

Projects/Dental CT Enhancement/README.md

@@ -0,0 +1,68 @@
+# CT图像处理
+
+> 作者：赵文亮
+>
+> 班级：自64
+>
+> 学号：2016011452
+>
+
+## 运行环境
+
+- Windows 10 x64
+- MATLAB R2018b
+- GUI由MATLAB App designer 编写
+
+## 运行方式
+
+### 打开程序
+
+- 在`app`文件夹下，双击`CT_Process.mlapp`
+- 或者在MATLAB中，切换到`app`目录下，命令行中输入`CT_Process`
+- 程序初始界面为`待处理图像`标签页
+
+### 图像处理流程
+
+#### 脊柱去除
+
+- 切换到`脊柱去除`标签页
+
+
+- 可以修改**滤波半径、亮度调整(Gamma校正)参数**
+- 实时查看脊柱去除效果
+
+#### 固定装置去除
+
+- 切换到`固定装置去除`标签页
+
+
+- 查看固定装置检测于去除结果
+
+#### 处理结果查看
+
+- 切换到`处理后图像`标签页
+- 对比处理前后的图像
+
+## 目录结构
+
+- src/：源代码
+
+- dataset/：输入图像
+
+- results/：处理结果
+
+  - remove_spine.jpg: 去除脊柱后的结果
+  - remove_fixers.jpg：去除固定装置后的结果（即最终结果）
+
+  > 该结果对应的参数为：
+  >
+  > - 滤波半径$D_0=15$
+  > - Gamma校正$a = 1.1, \gamma = 0.7$，校正公式为$I'=aI^\gamma$
+  >
+  > 用户可以在图形界面中修改以上参数以获得更加满足自己喜好的结果。
+
+- app/：图形界面文件`CT_Process.mlapp`
+
+- report.pdf：报告
+
+- README.md：说明文档

Projects/Dental CT Enhancement/src/crop.m

@@ -0,0 +1,5 @@
+% crop image I according the range defined by roi
+function c = crop(I, roi)
+    roi = floor(roi);
+    c = I(roi(3) : roi(4), roi(1) : roi(2));
+end

Projects/Dental CT Enhancement/src/detect.m

@@ -0,0 +1,6 @@
+% detect all the fixers
+function [fixers, masks, lineMasks] = detect(G)
+    [fixers{1}, masks{1}, lineMasks{1}] = detectFixer(G, [1100, 1130, 701, 1000], [0.3, 0.6], 5, false);
+    [fixers{2}, masks{2}, lineMasks{2}] = detectFixer(G, [841, 947, 1427, 1675], [0.4, 0.8], 3, false);
+    [fixers{3}, masks{3}, lineMasks{3}] = detectFixer(G, [1079, 1124, 2061, 2309], [0.4, 0.8], 3, false);
+end

Projects/Dental CT Enhancement/src/detectFixer.m

@@ -0,0 +1,87 @@
+% detect a fixer
+function [fixer, mask, lineMasks] = detectFixer(Img, roi, thresh, sigma, draw, show)
+    if nargin < 6
+        show = false;
+    end
+    % crop
+    I = Img(roi(1) : roi(2), roi(3) : roi(4));
+    BW = edge(I, 'canny', thresh, sigma);
+
+    [H, theta, rho]= hough(BW, 'RhoResolution', 0.5, 'ThetaResolution', 0.5);  
+    peak = houghpeaks(H, 2); 
+    lines = houghlines(BW, theta, rho, peak, 'MinLength', 10); 
+
+    % show hough result
+    if show
+        figure();
+        colormap(gca, hot);
+        imshow(imadjust(rescale(H)),'XData',theta,'YData',rho,...
+      'InitialMagnification','fit');
+        xlabel('\theta'), ylabel('\rho');
+        axis on, axis normal, hold on;
+    end
+
+    % adjust the slope of lines
+    fixer = lines;
+    if lines(1).theta ~= lines(2).theta
+        tmp_p = (lines(1).point1 + lines(1).point2) / 2;
+        delta_p = (lines(1).point2 - lines(2).point1) / 2;
+        lines(1).point1 = tmp_p;
+        lines(1).point2 = tmp_p + (lines(2).point2 - lines(2).point1) / (lines(2).point2(2) - lines(2).point1(2)) * delta_p(2);
+    end
+
+    for i = 1 : length(lines)
+        p1 = lines(i).point1 + [roi(3) -  1, roi(1) - 1];
+        p2 = lines(i).point2 + [roi(3) -  1, roi(1) - 1];
+        fixer(i).point1 = (p1 - p2) / (p1(2) - p2(2)) * (1 - p2(2)) + p2;
+        fixer(i).point2 = (p1 - p2) / (p1(2) - p2(2)) * (size(Img, 1) - p2(2)) + p2;    
+    end
+
+    % ensure fixer(1) is the left bound
+    mid1 = (fixer(1).point1 + fixer(1).point2) / 2;
+    mid2 = (fixer(2).point1 + fixer(2).point2) / 2;
+    if mid1(1) > mid2(1)
+        tmp = fixer(1);
+        fixer(1) = fixer(2);
+        fixer(2) = tmp;
+    end
+
+    % mask denotes the fixer area
+    canvas = zeros(size(Img));
+    mask = insertShape(canvas, 'FilledPolygon', ...
+    [fixer(1).point1 fixer(1).point2 fixer(2).point2 fixer(2).point1],...
+    'Color', 'White', 'Opacity', 1, 'LineWidth', 10);%, 'SmoothEdges', false); 
+    mask = mask(:, :, 1);
+
+    % lineMasks denote the fixer bounds
+    lineMasks{1} = zeros(size(Img));
+    lineMasks{2} = zeros(size(Img));
+    lineMasks{1} = insertShape(canvas, 'Line', [fixer(1).point1 fixer(1).point2],...
+        'Color', 'White', 'Opacity', 1, 'LineWidth', 1, 'SmoothEdges', false);
+    lineMasks{2} = insertShape(canvas, 'Line', [fixer(2).point1 fixer(2).point2],...
+        'Color', 'White', 'Opacity', 1, 'LineWidth', 1, 'SmoothEdges', false);
+    lineMasks{1} = lineMasks{1}(:, :, 1);
+    lineMasks{2} = lineMasks{2}(:, :, 1);
+
+    if (show)
+        figure();
+        imshow(Img);
+        hold on;
+        for i = 1 : length(fixer)    
+            xy=[fixer(i).point1; fixer(i).point2];    
+            plot(xy(:,1), xy(:,2), 'LineWidth', 1);    
+        end
+        hold off;
+    end
+
+    if (draw)
+        imwrite(I, 'roi.jpg');
+        imwrite(BW, 'bw.jpg');
+        f = BW;
+        colors = {'Yellow', 'cyan'};
+        for k=1:length(lines)
+            f = insertShape(double(f), 'Line', [lines(k).point1, lines(k).point2], 'Color', colors{k}, 'Opacity', 1, 'LineWidth', 3);
+        end
+        imwrite(f, 'draw.jpg');
+    end
+end

Projects/Dental CT Enhancement/src/drawAllFixers.m

@@ -0,0 +1,21 @@
+function G = drawAllFixers(G, fixers, show)
+    if show
+        imshow(G);
+        hold on;
+        for k = 1 : numel(fixers)
+            for i = 1 : length(fixers{k})
+                xy=[fixers{k}(i).point1; fixers{k}(i).point2];
+                plot(xy(:,1), xy(:,2), 'LineWidth', 1);
+                hold on;
+            end
+        end
+        hold off;
+    else
+        color = {'red', 'green', 'yellow'};
+        for k = 1 : numel(fixers)
+            for i = 1 : length(fixers{k})
+                G = insertShape(G, 'Line', [fixers{k}(i).point1 fixers{k}(i).point2], 'Color', color{k}, 'Opacity', 1, 'LineWidth', 5);
+            end
+        end
+    end
+end

Projects/Dental CT Enhancement/src/getMask.m

@@ -0,0 +1,19 @@
+% universial method to get mask
+function mask = getMask(G, x, show)
+    if nargin == 2
+        show = false;
+    end
+    if length(x) == 1
+        mask = insertShape(zeros(size(G)), 'Line', [x, 1, x, size(G, 1)], 'Color', 'White', 'Opacity', 1, 'LineWidth', 1);
+    elseif length(x) == 2
+        mask = insertShape(zeros(size(G)), 'FilledPolygon', [x(1, 1), x(1, 2), x(1, 1), x(2, 2), x(2, 1), x(2, 2), x(2, 1), x(1, 2)], 'Color', 'White', 'Opacity', 1);
+    elseif length(x) == 4
+        mask = insertShape(zeros(size(G)), 'Line', x, 'Color', 'White', 'Opacity', 1, 'LineWidth', 1);
+    elseif length(x) == 8
+        mask = insertShape(zeros(size(G)), 'FilledPolygon', x, 'Color', 'White', 'Opacity', 1);
+    end
+    mask = mask(:, :, 1);
+    if nargin == 3 && show == true
+        imshow(mask);
+    end
+end

Projects/Dental CT Enhancement/src/intensityAdjust.m

@@ -0,0 +1,47 @@
+function adjusted = intensityAdjust(G, mask, lineMasks, offset, roi, gamma, draw)
+    if offset == 0
+        adjusted = G;
+        return;
+    end
+    if nargin == 4 || (numel(roi) == 1 && roi == -1)
+        roi = [1, size(G, 1), 1, size(G, 2)];
+    end
+    if nargin < 6
+        gamma = false;
+    end
+    if nargin < 7
+        draw = false;
+    end
+    mask = cropMask(mask, roi);
+
+    % crop lineMasks to roi
+    lineMasks{1} = cropMask(lineMasks{1}, roi);
+    lineMasks{2} = cropMask(lineMasks{2}, roi);
+
+    % get left and right side masks
+    sm_left = sideMask(lineMasks{1}, 'left', offset);
+    sm_right = sideMask(lineMasks{2}, 'right', offset);
+    fixer_area = G .* mask;
+    neighbor_area = G .* sm_left + G .* sm_right;
+
+    % adjust intensity according to the intensity of side masks
+    I_fixer = sum(fixer_area, [1, 2]) / sum(fixer_area > 0, [1, 2]);
+    I_neighbor = sum(neighbor_area, [1, 2]) / sum(neighbor_area > 0, [1, 2]);
+
+    % use gamma adjust or linear adjust
+    if gamma
+        adjusted = G .* (1 - mask) + G .^ (log(I_neighbor) / log(I_fixer)) .* mask;
+    else
+        adjusted = G .* (1 - mask) + G .* mask * I_neighbor / I_fixer;
+    end
+
+    if draw
+        imwrite(mask, 'intensity_mask.jpg');
+        imwrite(sm_left + sm_right, 'neighbor_mask.jpg');
+    end
+end
+
+function outputMask = cropMask(mask, roi)
+    outputMask = zeros(size(mask));
+    outputMask(roi(1) : roi(2),  roi(3) : roi(4)) = mask(roi(1) : roi(2),  roi(3) : roi(4));
+end

Projects/Dental CT Enhancement/src/removeFixer.m

@@ -0,0 +1,33 @@
+% remove fixer decided by lineMask
+function outputImg = removeFixer(Img, lineMask, blurLen, threshold)
+    if nargin == 3
+        threshold = -1;
+    end
+    [yy, xx] = ind2sub(size(Img), find(lineMask > 0));
+    outputImg = Img;
+
+    % blur the bound of fixer
+    for i = 1 : length(yy)
+        y = yy(i);
+        x = xx(i);
+        if ((size(threshold, 2) == 1 && y < threshold) || (size(threshold, 2) == 2 && (y > threshold(2) || y < threshold(1))))
+            continue;
+        end
+        for d = 1 : blurLen
+             outputImg(y, x - d : x + d) = (outputImg(y, x - d : x + d) + interpolate(x - d, x + d, Img(y, x - d), Img(y, x + d))) / 2;
+        end
+    end
+end
+
+function y = interpolate(x1, x2, y1 ,y2)
+    if (x2 == x1)
+        y = y1;
+        return;
+    end
+    b = -(y2 - y1) / (x2 - x1) * 3 / 2;
+    a = 4 * b / (3 * (x2 - x1) ^ 2);
+    x = x1 : x2;
+    x_mean = mean([x1, x2]);
+    f = @(x)(a * (x - x_mean) .^3 - b * (x - x_mean) + mean([y1, y2]));
+    y = f(x);
+end

Projects/Dental CT Enhancement/src/removeFixers.m

@@ -0,0 +1,10 @@
+% remove all the fixers
+function J = removeFixers(I, lineMasks, offset, threshold)
+    J = I;
+    for u = 1 : 3
+        for v = 1 : 2 
+            J = removeFixer(J, lineMasks{u}{v}, offset(u, v), -1);
+            J = removeFixer(J, lineMasks{u}{v}, 30, threshold(u, v));
+        end
+    end
+end

Projects/Dental CT Enhancement/src/removeSpine.m

@@ -0,0 +1,41 @@
+function outputImg = removeSpine(I, D0, gamma, a)
+    if nargin < 3
+        gamma = 0.8;
+        a = 1.1;
+    elseif nargin < 4
+        a = 1.1;
+    end
+    % high pass filter
+    roi =  [500, 1; 2600, size(I, 1)];
+    I_crop = crop(I, roi);
+
+    I_fft = fftshift(fft2(I_crop));
+    H = zeros(size(I_fft));
+    [m, n] = size(H);
+    order = 2;
+    for u = 1 : m
+        for v = 1 : n
+            D = sqrt((u - (m + 1) / 2) ^ 2 + (v - (n + 1) / 2) ^ 2);
+            H(u, v) = 1 - 1 /(1 + (sqrt(2) -  1) * (D / D0) ^ (2 * order));
+        end
+    end
+
+    I_fft = I_fft .* H;
+    I_crop = real(ifft2(ifftshift(I_fft)));
+    I_crop = (I_crop - min(min(I_crop)))./ (max(max(I_crop)) - min(min(I_crop))) * (max(max(I)) - min(min(I)) + min(min(I)));
+    J = I;
+    J(roi(3) : roi(4), roi(1) : roi(2)) = I_crop  .^ gamma * a;    
+    % blend
+    sigma = 300;
+    x = 1 : size(I, 2);
+    d = 500;
+    c = 1547;
+    mask1 = meshgrid(gaussmf(x, [sigma, c - d / 2]), 1 : size(I, 1));
+    mask2 = meshgrid(gaussmf(x, [sigma, c + d / 2]), 1 : size(I, 1));
+    mask1(:, c - d / 2 : end) = 0;
+    mask2(:,1 : c + d / 2) = 0;
+    mask = mask1 + mask2;
+    mask(:, c - d / 2: c + d / 2) = 1;
+    mask(mask > 1) = 1;
+    outputImg = J .* mask + I .* (1 - mask);
+end

Projects/Dental CT Enhancement/src/sideMask.m

@@ -0,0 +1,17 @@
+% generate a sidemask according to the lineMask and direction
+function mask = sideMask(lineMask, direction, offset)
+    mask = zeros(size(lineMask));
+    if strcmp(direction, 'left')
+        factor = -1;
+    else
+        factor = 1;
+    end
+    index = find(lineMask == 1);
+    for i = 1 : offset
+        try
+            mask(i * size(lineMask, 1) * factor + index) = 1;
+        catch
+            break;
+        end
+    end
+end