SCE Carleton Logo
  Carleton > Engineering > SCE > Faculty > A. Adler > Courses > ELG7173 > Ultrasound Image Processing


Ultrasound Image Processing

Analysis of Ultrasound images

Automatic analysis of medical images is somewhat controversial. The advantages are clear: computer analysis is cheaper and results will not vary for a given image. However, medical images are highly variable. Algorithms are inherently tuned to look for specific features, and are prone to fail spectacularly when given unexpected input data. My opinion is that the best current use of technology is to assist trained technologists. It would attempt to detect candidate features, but would require human validation of all questionable decisions.

Example: Detection of Parameters for glaucoma measurement

Credit: Much of the analysis presented here is from work by Pino Dicorato

Data analysed is from the paper: Elucidating the mechanism of Occludable Angles in patients with Pseudoexfoliation Syndrome by E.A. Ross K.F. Damji W.G. Hodge A. Al-Ani A.A. Al-Ghamdi K. Shah D. Chialant
This work performed by team at U. Ottawa Eye Institute.

Glaucoma is one of the leading causes of blindness in North-America and around the world. In closed angled Glaucoma, high intraocular pressure (fluid pressure in the eye) resulting from inadequate fluid flow between the iris and the cornea causes damage and eventually death of nerve fibres responsible for vision. One important technique to assess patients at risk of glaucoma is to analyze ultrasound images of the eye, for structural changes that reduce the outflow of fluids out of the eye. Sequences of ultrasound images of the eye are analyzed manually; anatomical feature locations are determined by a trained technician, and relevant clinical parameters subsequently measured. We discuss image processing techniques to automate the location of clinically relevant features in ultrasound images of the eye. We describe a feature detection algorithm to identify: 1) the location of the canal of Schlemm (drainage), 2) the location of ciliary body, and 3) the angle formed between the iris and the cornea.

In this study, all image analysis was performed by a trained technologist. We are interested in whether it could be possible to automatically measure these parameters. The expert measurements constitute the gold standard against which computer estimates are validated.

Technical Goal

Measure the following parameters from images

Example Images (Scanned from print-out):

Open Closed

Image Processing

We illustrate steps involved to find Θ1 (the apex of the angle is feature B) and the Scleral Spur (feature A). This illustrates some of the challenges involved in such image processing. Examples are based on this image

Ultra-sound image of Iris: Output from scanner to 3.5" disk.

Locate feature B

This feature is surprisingly hard to accurately obtain. The following code samples are implemented in octave, and require the edge.m code from
Code Image
load and crop image region
im= im(90:180,:);

cmap= [0,0,0;1,1,1;0,0,.8;.8,0,0];

Image is cropped to show core region of interest
Sobel edge detecion
[im1, thresh] = edge( im, 'sobel');
% filt1= [1 2 1;0 0 0;-1 -2 -1];
% filt2= rot90( filt1 )
% combine= sqrt( filt1^2 + filt2^2 )
% Thresh= 23.1 

Problem is that there are holes in edge. A fill will "pour" through holes
Lower Sobel edge threshol
im2 = edge( im, 'sobel', 10);
% Set threshold lower to get more edge
imwrite('img-sobel2.png',im2+1, cmap);

The edge is still unstable and we have artefacts in the black region
Flood fill region of interest
bw= bwfill(im2, 45, 256,4);
% use 4 connected bw fill
fill= 2*bw - im2 + 1;
imwrite('img-s2fill.png',fill, cmap);

Note that the fill region (blue) penetrates into the tissue region
Use Advanced Edge Detection (Canny type)
% Do Sobel edge detection to generate an image at
%  a high (HT) and low threshold (LT)
% Edge extend all edges in the LT image by several
%  pixels, in vertical, horizontal, and 45° dirns.
%  Combine these into edge extended (EE) image
% Dilate the EE image by 1 step
% Select all EE features that are connected to
%  features in the HT image
[im3, thresh] = edge( im, 'andy');

This works fairly well, but there is some edge blurring.
Flood fill region of interest
bw= bwfill(im3, 45, 256,4);
fill= 2*bw - im3 + 1;
imwrite('img-a-fill.png',fill, cmap);

Now, we need to fit a triangle to this shape
Find top and bottom edge
im4= fill==3; % left-most filled point
leftpt= find( any(im4) )(1);
im5= im4( :, leftpt:size(im4,2) ); 
bw= bwfill(im5, 1, 1,4);
fill= im5 + bw + 1;
imwrite('img-fill2.png',fill, cmap);

Three regions are defined
Detect Edges
mask= [0,1,0;1,1,1;0,1,0]/5;
edge1= (conv2(fill==2,mask,'same')>0).*im5;
edge2= (conv2(fill==1,mask,'same')>0).*im5;
fill= 2*edge1 +1; fill(find(edge2))= 4;
imwrite('img-edge2.png',fill, cmap);

Red and blue lines mark the boundaries
Extract fit from lines
function [fit,idx]= imagefit( img, poly_degree )
   pt= find( img > 0 ) - 1; % we want 0 base index
   x= floor( pt / size(img,1) );
   y=   rem( pt , size(img,1) );
   fit = polyfit(x,y, poly_degree);
   yfit= polyval( fit, x);
   idx=  1+ round(x*size(img,1) + yfit);

[fit,idx1]= imagefit( edge1,1);
[fit,idx2]= imagefit( edge2,1);

fill= im5+1; fill(idx1)= 3; fill(idx2)= 4;
imwrite('img-edge3.png',fill, cmap);

This doesn't work. We must force lines to meet at point B.
Extract fit. Lines must meet at centre, and only use initial portion of section
function [fit,idx]= imagefit( img, uselim)
   pt= find( img(:,1:uselim) > 0 ) - 1;
   x= floor( pt / size(img,1) ); 
   y=   rem( pt , size(img,1) );
   ymean=  mean( find (img(:,1)) - 1); 
   fit = [x\(y-ymean); 0]; % slope is x\y, intercept is 0 
   yfit= polyval( fit, x) + ymean;
   idx=  1+ round(x*size(img,1) + yfit);

[fit,idx1]= imagefit( edge1, 60);
[fit,idx2]= imagefit( edge2, 60);

fill= im5+1; fill(idx1)= 3; fill(idx2)= 4;
imwrite('img-edge4.png',fill, cmap);

The point position and angle can now be determined
Global Picture
cmap2= [gray(256); [.8,0,0]];
imwrite('img-featb.png',im+1, cmap2);

Note that feature B is still not accurately estimated.
Improvements are left as an excercise for the reader 8-)

Locate feature A

This feature is difficult for many reasons. Most importantly, ultrasound has intrinsically very low SNR, due to Rayleigh scattering. This means we cannot easily use contrast for discrimination of regions in such tissue.

Some code to illustrate the effect of Rayleigh scattering is here.

Code Image
load and crop image region
im= im(110:200,:);

cmap2= [gray(256); [.8,0,0];[0,0,.8]];

Image is cropped to show core region of interest
Select tissue regions
fill= im;
fill( 24:53, 29:91 ) = 256;
fill( 26:51, 31:89 ) = im( 26:51, 31:89 );
fill( 70:87, 49:94 ) = 257;
fill( 72:85, 51:92 ) = im( 72:85, 51:92 );

View Histograms
r1= im( 24:53, 29:91 );
hist( r1(:), 0:5:255, 1);
hold on
r2= im( 70:87, 49:94 );
hist( r2(:), 0:5:255, 1);

Horiz axis is 0 to 255. Regions are faily distinct.
Threshold at crossover point
imt= im>170;

Threshold works quite well in centre of region, but not well in bottom left
Crop and fill from bottom
imt2= imt(:,1:180);
bw= bwfill(imt2, size(imt2,1), 1, 4);
fill= 2*bw - imt2 + 1;

blue region has "holes"
Fill holes in blue region
bw = fill==3;
bw2= bwfill(bw, 'holes', 4);
imwrite('img-t3.png',bw2 + 1,cmap);

Find edge
mask= [0,1,0;1,1,1;0,1,0]/5;
edge= (conv2(bw2,mask,'same')>0).*(1-bw2);
fill= imt2 +1; fill(find(edge))= 4;
imwrite('img-edget.png',fill(30:size(fill,1),:), cmap);

Fit line to straight segment
function [fit,idx]= imagefit( img, poly_degree )
   pt= find( img > 0 ) - 1; % we want 0 base index
   x= floor( pt / size(img,1) );
   y=   rem( pt , size(img,1) );
   fit = polyfit(x,y, poly_degree);
   yfit= polyval( fit, x);
   idx=  1+ round(x*size(img,1) + yfit);

[fit,idx1]= imagefit( edge(:,38:149),1);
fill2= fill;
fill2(idx1 + size(fill,1)*38)= 3;
imwrite('img-edgef.png',fill2, cmap);

Find edge point maximum depth below line
for i=39:149
   edge_y(i)= mean(find(fill(:,i)==4));
   line_y(i)= mean(find(fill2(:,i)==3));
dif=  edge_y - line_y;
ptx= mean( find( dif == max(dif) ) );
fill2(:,round(ptx))= 3;
imwrite('img-edgep.png',fill2, cmap);

On the original image this gives
pty= mean(find(fill(:,round(ptx))==4));
fill= im;
fill( :, round(ptx) ) = 256;
fill( round(pty), : ) = 256;

While this works, for this particular image, this processing is fairly unreliable, because:

  • The threshold value needs to be accurately chosen
  • The horizontal limits for the line fit need to be accurately chosen
In the following section, we focus on a technique to better separate the regions to find an appropriate threshold. The essense of these techniques is to use texture analysis.

NOTE: This material is for explanation of the concepts of texture analysis. Correct implementation of these ideas uses techniques such as Markov random fields.
Code Image
Analyse region autocorrelations
sz= 8; sz2= 2*sz+1;

xc1= xcorr2(r1, 'unbiased');
xctr= ceil(size(xc1,2)/2);
yctr= ceil(size(xc1,1)/2);
xc1=xc1( yctr+(-sz:sz), xctr+(-sz:sz));

xc2= xcorr2(r2, 'unbiased');
xctr= ceil(size(xc2,2)/2);
yctr= ceil(size(xc2,1)/2);
xc2=xc2( yctr+(-sz:sz), xctr+(-sz:sz));

minfill= min([xc1(:);xc2(:)]);
fill= [xc1-minfill, zeros(sz2,2), xc2-minfill];

Xcorr of R1 on left and R2 on right
R1 is more spatially uniform than R2
Develop a "sample" of regions
sz2= 2*4+1;

f1e= interpft( fft2(r1), 2*sz2, 2*sz2 );
i1e= real(ifft2(f1e)); i1e=i1e(1:sz2, 1:sz2);

f2e= interpft( fft2(r2), 2*sz2, 2*sz2 );
i2e= real(ifft2(f2e)); i2e=i2e(1:sz2, 1:sz2);

Original histogram of regions

Selectively Enhance with region 1
r1h= xcorr2(r1,i1e,'unbiased')(:);
r2h= xcorr2(r2,i1e,'unbiased')(:);
scale= 255/max([r1h;r2h]);
hold off; hist(scale*r1h, 0:5:250, 1); 
hold on; axis ([0, 250, 0, .15]);
hist(scale*r2h, 0:5:250, 1); 

Region 2 now overlaps Region 1
Selectively Enhance with region 2
r1h= xcorr2(r1,i2e,'unbiased')(:);
r2h= xcorr2(r2,i2e,'unbiased')(:);
scale= 255/max([r1h;r2h]);
hold off; hist(scale*r1h, 0:5:250, 1); 
hold on; axis ([0, 250, 0, .15]);
hist(scale*r2h, 0:5:250, 1); 

Greater separation between regions
Enhance image with region 2
ime2= xcorr2(im, i2e, 'unbiased');
ime2= ime2-min(ime2(:));
ime2= 255*ime2/max(ime2(:));
ime2crop= ime2(20:size(ime2,1)-10,5:size(ime2,2)-4);
fill= [ ime2crop; im(25:size(im,1)-6,:)];

Greater separation between regions

Last Updated: $Date: 2004-04-12 09:48:56 -0400 (Mon, 12 Apr 2004) $