% Script file: Edge segmentation from DXA image
% Purposes: This program detects for the edge from femur using Sobel
% combined with morphorgical method. The position is defined by the
% operator, which are narror neck, intertrochanteric, and shaft.
% Define variables:
% I_grey: gray-scale image of Normal_NoFrac_Case
% I_double: double format of I_grey
% M: mean value of grey scale image
% VAR: standard deviation value of grey sacle image
clc; clear all; close all;
I = imread('Precision_1_before.jpg'); % input color image;read the image
I_grey = rgb2gray(I); % convert the image into gray form
I_double = double(I_grey); % convert the image into double format
[height width] = size(I_double);
for i = 2:height-1
for j = 2:width-1
sortseq = [I(i-1,j-1),I(i-1,j),I(i-1,j+1),I(i,j-1),I(i,j),I(i,j+1),I(i+1,j-1),I(i+1,j),I(i+1,j+1)];
s = sort(sortseq);%sort the element of the array
I_double(i,j) = s(5);%select the median
end
end
% I_double = wiener2(I_double, [5,5]); %Wiener low-pass noise-removal
axis on;
title('\fontsize{15}De-noise DXA image','Color','r');
M = mean2(I_double);
VAR = std2(I_double);
% Mo and Vo are the specified standard values
M0 = 195;
V0 = 12;
[row,col] = size(I_double);
for i = 1:row;
for j = 1:col;
if(I_double(i,j)>M)
Nor(i,j) = M0+(V0*(I_double(i,j)-M)^2/VAR)^(1/2);
else
Nor(i,j) = M0-(V0*(I_double(i,j)-M)^2/VAR)^(1/2);
end
end
end
% Sobel operator;
T = 58; % threshold
X0 = I_double;
[row,col] = size(X0);
g = zeros(row,col);
X = Nor;
for i = 2:row-1
for j = 2:col-1
Dx = (X(i+1,j-1)-X(i-1,j-1)+2*(X(i+1,j)-X(i-1,j))+X(i+1,j+1)-X(i-1,j+1))^2;
Dy = (X(i-1,j+1)-X(i-1,j-1)+2*(X(i,j+1)-X(i,j-1))+X(i+1,j+1)-X(i+1,j-1))^2;
g(i,j) = round(sqrt(Dx+Dy));
end
end
for i=1:row
for j=1:col
if g(i,j)>T
X(i,j)=255;
else
X(i,j)=0;
end
end
end
A = strel('disk',5);
X = imdilate(X,A);
X1 = bwmorph(X,'majority',500);
X2 = bwmorph(X,'thin',Inf);
imshow(X2,[]);
axis on;
grid on
title('\fontsize{15}Sobel operator combined with morphological image','Color','r');
% I_plus = uint8(X2) + uint8(I_grey);
% figure, imshow(I_plus)
% figure, imshow(I_grey)
% [c_col c_row c] = improfile;
% improfile();
% index = find(c==max(c));
% x_1 = c_col(index(1));
% y_1 = c_row(index(1));
% x_2 = c_col(index(2));
% y_2 = c_row(index(2));
% improfile;
% [x_1,y_1] = ginput(1);
% [x_2,y_2] = ginput(1);
% delta_pixel_x = x_2-x_1;
% real_distance_x_mm = 0.264583*delta_pixel_x