Image illumination correction of arterial contour

I have a code which which I am a bit confused on how to optimize in order to make it work for an image http://tinypic.com/view.php?pic=2jdx6g&s=7 since I am not able to modify the function imopen to remove the background illumination in order to obtain clear picture of the coronary vessels. Any suggestion?

Answer 1

[Image Analyst]

That is not what imopen is supposed to do. If anything it would be imclose() which does a dilation (smear bright stuff over the dark vessels) followed by an erosion (to return enlarged bright things to their original size). But I'd look at the literature to see what people are doing successfully with angiograms rather than spend time implementing something that ends up being really primitive and ineffective. http://iris.usc.edu/Vision-Notes/bibliography/contentsmedical.html#Medical%20Applications,%20CAT,%20MRI,%20Ultrasound,%20Heart%20Models,%20Brain%20Models

[mona]

Well you said that, I actually need to implement something that is as primitive as extracting the background and then subtracting the high contrast image from the background or smear the background only thereby obtaining the vessel contour. I am specifically confused about is how I would work with the background which seem to mix with the image at certain points.

[Image Analyst]

That could work. I'm not saying it's bad, just that it's basic, rudimentary, or primitive. But I doubt such a simple concept is the current state of the art - I bet there are better, more sophisticated, but unfortunately more complicated methods out there. But this simple method might work for your needs. All I'm saying is that maybe you should see what others are doing since they're probably using something better than what we're discussing.

[mona]

What I need to know at this very level is, how could this basic thing be implemented since I am trying to do the same with the code and after that I would move on to next level.

[Image Analyst]

Like I said, use imclose() to get the background image. Then, since it's radiology, use background subtraction (rather than division) to get the background corrected image. Cast both images to double before you do the subtraction.

[mona]

You are right. But I have tried the code while taking log and the exponential after which I have with imclose command. The difference between the image is not as desired due to some reason. It seems I am missing something significantly large. Little help here please,

B=log(1+I);
coeff=fir1(64,1/10000);
h=coeff'*coeff;
figure(2),freqz2(h)
g1 = imfilter(B, h);
g2=exp(g1)-1;

[Image Analyst]

I don't see imclose() and don't know the purpose of log() and fir1() or any of that code. I thought you just wanted to background correct the images. If you want to discover (-mu*x) of your contrast agent, you just take the log of both images (with and without contrast) and subtract them. If you don't have the image without contrast agent (for some strange reason) then you can estimate it by doing a morphological closing.

The theory says, for the background image

image1 = initialIntensity * exp(-mu1 * x1); % x = thickness

For the image with contrast

image2 = initialIntensity * exp(-mu1 * x1 - mu2 * x2); % x = thickness

So the steps would be

image2 = imclose(image1, close(15)); % May need to adjust the 15
mu2TimesThickness2 = log(image1+1) - log(image2+1);

or something like that. Don't you agree with my math?

[mona]

Yes, I do agree with you on this context. Although I don't understand what you assumed related to "background correct" but here's a simple case scenario of what I intend.

I just need the vessel to be darkened enough so that it is the most prominent feature of the image and for that purpose, the patches in the background should be suppressed as much as possible. Therefore, in a very simple coding,

I =double(rgb2gray(my image));
se = strel('disk',500);
g2 = imclose(I,se);
g2=mat2gray(g2);
I2 = I - g2;
figure(), imshow(I2,[])

which ofcourse doesn't work but atleast portrays the idea.

[Image Analyst]

You didn't do what I said. Try this:

grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows columns numberOfColorBands] = size(grayImage);
if numberOfColorBands > 1
  grayImage = rgb2gray(grayImage);
end
% Display the original gray scale image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Give a name to the title bar.
set(gcf,'name','Demo by ImageAnalyst','numbertitle','off') 

% Close the image
se = strel('disk', 21);
closedImage = imclose(grayImage, se);
subplot(2, 2, 2);
imshow(closedImage, []);
title('Closed Image', 'FontSize', fontSize);

% Subtract to get the vessels alone.
vesselImage = log(1+double(closedImage)) - log(1+double(grayImage));
subplot(2, 2, 3);
imshow(vesselImage, []);
title('Subtracted Image', 'FontSize', fontSize);

I think the output looks pretty good. If you don't want the vessels to be white (which is normal for "foreground" objects, especially if you want to use regionprops()), then you can just reverse the subtraction:

vesselImage = log(1+double(grayImage)) - log(1+double(closedImage));