from Edge detection based on the threshold by Yujia
This edge detection contains threshold selection, morphological processing and image normalization.

ti_qu_bian_yuan.m
% 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

Contact us