Image Processing for ultrasound

I would like to do image processing to recognize a needle in an ultrasound image. Such as the one below

This image is scaled down to 20%

Currently I am thinking about blurring, created a region of interest, thresholding the image and then identifiying the most continous blob of white pixels.

Does anyone have a way to create a blur that utilizes the brightest pixels to change the pixels around it?

Other feedback and methods and input would also be appreciated.

Also @Image Analyst might be interetsted

[Image Analyst]

Can you post the original image without the axes on it? Also what do we know about the needle? Is it always about the same size? Is it always vertical? Is it always in about the same position so we can use an ROI masking out the large white blobs at the top and bottom of the scan? Try using imdilate to connect the dots in the needle to make it easier to get it all in one blob.

[Ethan Zhong]

The original images will be .mat files that are directly sent from the matlab software the ultrasound is run on but I only have the jpg screenshot of it right now. I believe this is being used for a needle insertion purpose so we can assume the needle will be mostly verticle but I would like to make it as robust as possible. The needle may be of different sizes I don't know if it'll be drastically different. I believe the ROI is to take out noise for now as the device has some refraction due to the plastic material and the non-alligned mirrors we are using but we are also working on fixing that.

Answer 1

[Mathieu NOE]

let's try some stuff

ok I am not an image processing expert , so use that or not ...

filename = 'image_needle.jpg';

inpict = im2double(rgb2gray(imread(filename)));

% remove mean along dimension 2 (that will remove already a big portion of
% the streaks)
inpict2 = inpict - mean(inpict,2);

spec_img = fftshift(fft2(inpict2));
sze = size(spec_img);

% low pass elliptical filering of the input image (to remove further the
% salt pepper noise) - adjust filter cut off and order to your own
% preferences
sze = size(spec_img);
cutoff1 = 0.5;
cutoff2 = 0.05;
n = 6;
f = EFilter(sze, cutoff1, cutoff2, n);

% apply filter
spec_img = spec_img.*f;

% generated backward the output image by inverse fft
outpict = real(ifft2(ifftshift(spec_img)));

figure
subplot(1,2,1),imshow(inpict)
subplot(1,2,2),imshow(outpict)

%% attempt # 1 to extract the needle segment

% Perform edge detection using the Canny method
edges = edge(outpict, 'canny');

% Perform Hough Transform to detect lines
[H, T, R] = hough(edges);


% Find peaks in the Hough Transform
peaks = houghpeaks(H, 10, 'threshold', ceil(0.2 * max(H(:))));

% Extract lines based on the Hough Transform and peaks
lines = houghlines(edges, T, R, peaks, 'FillGap', 5, 'MinLength', 7);

% Display the original image
figure
imshow(outpict);
hold on;

% Plot the detected lines
for k = 1:length(lines)
    xy = [lines(k).point1; lines(k).point2];
    plot(xy(:,1), xy(:,2), 'LineWidth', 2, 'Color', 'green')
end

%% attempt # 2 to extract the needle segment
% Threshold the image to create a binary image
tmp = abs(gradient(outpict));
binaryImage = tmp > 0.25*max(tmp,[],'all'); % Adjust the threshold value as needed

binaryImage = bwareaopen(binaryImage, 500); % remove small aeras 

% Label connected components
labeledImage = bwlabel(binaryImage);

% Display the original image
figure;
imshow(outpict);
title('Filtered Image');


% Highlight the white segments on the original image
hold on;
boundaries = bwboundaries(binaryImage);
for k = 1:length(boundaries)
    boundary = boundaries{k};
    plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2);
end
title('White Segments Highlighted');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = EFilter(sze, cutoffM, cutoffm, n, varargin)

%EFILTER constructs an elliptical lowpass filter Butterworth filter
%   E = EFILTER(Size, cutoffM, cutoffm, n) designs an Nth order elliptical
%   lowpass digital Butterworth filter where Size is a two element 
%   [rows, cols] vector specifying the size of the filter to construct, 
%   cutoffM and cutoffm, the cutoff freqency on the major and minor
%   axes are 0 < cutoff <= 1.
%
%   If E = EFilter(Size, cutoffM, cutoffm, n, alpha), where alpha is an angle
%   in radians, it will return and plot an elliptical filter rotated
%   counter-clockwise through alpha.
%
%   If E = EFilter(Size, cutoffM, cutoffm, n, alpha, xoff, yoff), where xoff
%   and yoff are offsets in the x and y direction, it will return and
%   plot an eliptical filter which is offset by the specified amount.
%   An offset of 0 corresponds to the center and an offset of 1
%   corresponds to the edge of the filter. A positive offset shifts the
%   filter in the positive direction.
%
%   Calling EFilter(...) without assigning the output variable
%   plots the 3D surface described by the function.

% Katie Streit   [e-mail]
% ELEC 301
% Rice University
%
% December 2001

% Much of this code was based on Peter Kovesi's  ([e-mail])
% Matlab function for a lowpass Butterworth filter.

if nargin == 4
alpha = 0;
offx = 0;
offy = 0;
elseif nargin == 5
offx = 0;
offy = 0;
alpha = varargin{1};
elseif nargin == 7
alpha = varargin{1};
offx = varargin{2};
offy = varargin{3};
else
	error('Invalid number of input arguments');
end

    if nargout > 1
        error('Invalid number of output arguments');
    end

    if cutoffM < 0 | cutoffM > 1
	error('cutoffM frequency must be between 0 and 1');
    end
    
    if cutoffm < 0 | cutoffm > 1
        error('cutoffm frequency must be between 0 and 1');
    end  
      
    if rem(n,1) ~= 0 | n < 1
	error('n must be an integer >= 1');
    end
    

%%extracts the sizes from sze 
rows = sze(1);
cols = sze(2);

%x and y matrices normalized to +/-.5 and an offset of offx or offy
x =  ((((ones(rows,1) * [1:cols])-offx*rows/2)  - (fix(cols/2)+1))/cols);
y =  ((([1:rows]' * ones(1,cols))-offy*rows/2) - (fix(rows/2)+1))/rows;

%applies a linear transformation to rotate through alpha. Note that it takes
% uses negative alpha, which is caused by x and y being independent matrices.
x2 = (x*cos(alpha) - y*sin(-alpha));
y2 = (x*sin(-alpha) + y*cos(alpha));

%constructs an elliptical cone (defined by a and b) of height r on each at 
%each elliptical ring. (r is effectively the "radius")
%r = sqrt(((x2/a).^2 + (y2/b).^2));

%Designs the filter
%f = 1./(1.0 + (r./cutoff).^(2*n));

a = cutoffM/2;
b = cutoffm/2;

f = 1./(1+((x2/(a)).^2 + (y2/(b)).^2).^n);

if nargout > 0
  varargout{1} = f;
else
  %Plots a normalized (+/- 1), interpolated 3D image of the filter
  surf([-1:2/(cols-1):1],[-1:2/(rows-1):1], f);
  shading interp;
  title('Elliptical Butterworth filter');
  xlabel('x');
  ylabel('y');
  zlabel('intensity');
  grid on;
end

end

[Mathieu NOE]

yet another try, without any image processing function !!!

result :

code (using the attached functions)

filename = 'image_needle.jpg';
inpict = im2double(rgb2gray(imread(filename)));

% remove mean along dimension 2 
inpict = inpict - mean(inpict,2);

% low pass elliptical filering of the input image (to remove further the
% salt pepper noise) - adjust filter cut off and order to your own
% preferences
spec_img = fftshift(fft2(inpict));

sze = size(spec_img);
cutoff1 = 0.5;
cutoff2 = 0.05;
n = 6;
f = EFilter(sze, cutoff1, cutoff2, n);

% apply filter
spec_img = spec_img.*f;

% generated backward the output image by inverse fft
outpict = real(ifft2(ifftshift(spec_img)));

figure
subplot(2,1,1),imshow(inpict)
subplot(2,1,2),imshow(outpict)

%% attempt to extract the white segment
% Threshold the image to create a binary image
binaryImage = outpict > 0.3*max(outpict(:)); % Adjust the threshold value as needed

% Display the binary image
figure;
imshow(binaryImage);
title('Binary Image (White Segments)');
hold on;

% find the boundary points
[y,x] = ind2sub(size(binaryImage),find(binaryImage>0.5));
[y_selec,x_selec] = myboundary(y,x);
plot(x,y, '.', x_selec, y_selec, '.r')
title('Boundary points Highlighted');

% last round !!!
% let say we don't want to keep line objects with width > tol (in pixels)
tol = 0.03*sze(2); % here the tol is 3% of the picture width
[y_selec_unic,ia,ic] = unique(y_selec);
                   
% "scroll" the image along the y direction and look for narrow  profiles 
m = 0;
for k = 1:numel(y_selec_unic)
    ind = ic == k;
    x_selected = x_selec(ind);
    dx = max(x_selected) - min(x_selected);
    if dx < tol % we keep it
        m = m+1;
        yfinal(m) = y_selec_unic(k);
        xfinal(m) = mean(x_selected);
    end
end

tmp = abs(diff(xfinal));
ind = find(tmp>0.1*max(tmp));
[v,ii] = max(diff(ind));
iii = ind(ii):ind(ii+1);
xxx = xfinal(iii);
yyy = yfinal(iii);

plot(xxx, yyy, 'dg')
hold off;

%% FINAL PLOT !!!!!!!!!

figure
imshow(inpict)
hold on 
plot(xxx, yyy, '.r')
hold off;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = EFilter(sze, cutoffM, cutoffm, n, varargin)

%EFILTER constructs an elliptical lowpass filter Butterworth filter
%   E = EFILTER(Size, cutoffM, cutoffm, n) designs an Nth order elliptical
%   lowpass digital Butterworth filter where Size is a two element 
%   [rows, cols] vector specifying the size of the filter to construct, 
%   cutoffM and cutoffm, the cutoff freqency on the major and minor
%   axes are 0 < cutoff <= 1.
%
%   If E = EFilter(Size, cutoffM, cutoffm, n, alpha), where alpha is an angle
%   in radians, it will return and plot an elliptical filter rotated
%   counter-clockwise through alpha.
%
%   If E = EFilter(Size, cutoffM, cutoffm, n, alpha, xoff, yoff), where xoff
%   and yoff are offsets in the x and y direction, it will return and
%   plot an eliptical filter which is offset by the specified amount.
%   An offset of 0 corresponds to the center and an offset of 1
%   corresponds to the edge of the filter. A positive offset shifts the
%   filter in the positive direction.
%
%   Calling EFilter(...) without assigning the output variable
%   plots the 3D surface described by the function.

% Katie Streit   [e-mail]
% ELEC 301
% Rice University
%
% December 2001

% Much of this code was based on Peter Kovesi's  ([e-mail])
% Matlab function for a lowpass Butterworth filter.

if nargin == 4
alpha = 0;
offx = 0;
offy = 0;
elseif nargin == 5
offx = 0;
offy = 0;
alpha = varargin{1};
elseif nargin == 7
alpha = varargin{1};
offx = varargin{2};
offy = varargin{3};
else
	error('Invalid number of input arguments');
end

    if nargout > 1
        error('Invalid number of output arguments');
    end

    if cutoffM < 0 | cutoffM > 1
	error('cutoffM frequency must be between 0 and 1');
    end
    
    if cutoffm < 0 | cutoffm > 1
        error('cutoffm frequency must be between 0 and 1');
    end  
      
    if rem(n,1) ~= 0 | n < 1
	error('n must be an integer >= 1');
    end
    

%%extracts the sizes from sze 
rows = sze(1);
cols = sze(2);

%x and y matrices normalized to +/-.5 and an offset of offx or offy
x =  ((((ones(rows,1) * [1:cols])-offx*rows/2)  - (fix(cols/2)+1))/cols);
y =  ((([1:rows]' * ones(1,cols))-offy*rows/2) - (fix(rows/2)+1))/rows;

%applies a linear transformation to rotate through alpha. Note that it takes
% uses negative alpha, which is caused by x and y being independent matrices.
x2 = (x*cos(alpha) - y*sin(-alpha));
y2 = (x*sin(-alpha) + y*cos(alpha));

%constructs an elliptical cone (defined by a and b) of height r on each at 
%each elliptical ring. (r is effectively the "radius")
%r = sqrt(((x2/a).^2 + (y2/b).^2));

%Designs the filter
%f = 1./(1.0 + (r./cutoff).^(2*n));

a = cutoffM/2;
b = cutoffm/2;

f = 1./(1+((x2/(a)).^2 + (y2/(b)).^2).^n);

if nargout > 0
  varargout{1} = f;
else
  %Plots a normalized (+/- 1), interpolated 3D image of the filter
  surf([-1:2/(cols-1):1],[-1:2/(rows-1):1], f);
  shading interp;
  title('Elliptical Butterworth filter');
  xlabel('x');
  ylabel('y');
  zlabel('intensity');
  grid on;
end

end

[Ethan Zhong]

Thank you so much!! I'll give this a try and see how it holds up with some other images too. It suppose to run real time so I'll also post updates on lag and everything!

[Mathieu NOE]

as always, my pleasure !

See more (8)