How can I isolate a faint gradient in an image using matlab?

I am trying to trace the border of a shockwave in a series of images as it advances in time.

For reference watch the video below at the 2:30 mark.

http://www.military.com/video/ammunition-and-explosives/explosives/shockwave-captured-on-high-speed/763995636001/

There is a thin almost translucent shock that moves faster than the debris. I have been converting the image to black and white so that I can use edge finding algorithms in matlab, the only trouble that I have is that I can't seem to isolate the shockwave. I have been using the graythresh() function, but it doesn't capture the wave. It captures the more obvious dirt, smoke, and debris instead of the more subtle wave.

I'm thinking that some sort of analysis of the gradient might be effective. Assuming a relatively static background would it be possible to examine two frames and detect the changes in them to isolate the shockwave?

I'm relatively new to image processing, so I'm open to any suggestions.

Thanks, Chris N.

Here are two sample images for reference. They are screenshots of the video separated by a few seconds. The motion and location of the wave becomes apparent when comparing two images. Otherwise it is hard to find the shock location. That is why I feel that some sort of comparison between two frames may work.

Answer 1

[Image Analyst]

What if you simply subtract the two frames?

diffImage = abs(double(frame1) - double(frame2));
imshow(diffImage);

[Chris]

That absolutely amplifies the subtle differences! Especially if I then use the im2bw() function to make the grayscaled image that I'm working with binary. Thank you.

Now I just need to find out how to clean up the image to be able to trace the boundary!

[Sean de Wolski]

bwmorph will be your friend

[Image Analyst]

Post your binary image. I assume it's good? If not you may need to pick a threshold better than im2bw does it (I assume you used graythresh() also????). You can use things like imclearborder() and bwareaopen() and imfill() and other things to clean up the image. Once you have the clean binary image, you can use bwboundaries() to get (x,y) coordinates of the boundary. Or if you want the boundary as an image instead of a list of coordinates, you can use bwperim().

[Chris]

Thank you for the additional replies Sean and Image Analyst!

Sorry for the delay. I did not receive an email about additional posts on this thread. The image subtraction of my images created the image below. The image below has not had any filters applied to it yet.

I noticed that if I use graythresh() on the image using the automatically computed value for the level parameter, the shock completely disappears. The image becomes mostly black. To get the image below I used im2bw(Image,0) where I randomly set the value of the threshold second parameters in the function to 0. It seemed to work.

I will try your suggestions to see if I can clean it up. My end goal is to be able to take a portion and compute the average distance from the center of the expansion.

This might be off topic, but one thing that I was wondering was if the wave shown in the image below represents the path traveled by the actual wave during two frames. Because I am subtracting two frames that represent two different times wouldn't that be interpreted as the distance traveled by the wave between the frames? Just a question. It's more physics instead of matlab.

[Image Analyst]

Are you sure you cast to double before subtraction? And graythresh() only provides good thresholds, in my opinion, for bimodal histogram which are for images with high contrast, not something like what you have. Can you post the histogram so I can look at it? Or post two of the frames?

And this image looks like an overhead view rather than a side view? If that's the case, you can calculate the distance and frame time differences and calculate velocity as you said.

It will be challenging but I think with the right algorithm you can get the distance of the various rings.

[Chris]

I did not cast to double before subtraction. I'll give that a try and see what I get.

The image is looking at it from the side. Like watching a hotdog explode in a microwave while looking at it from the end. I've included the two frames that I used to create the image.

Here are the images that I subtracted and used graythresh() on the resulting subtraction to create the image in my previous post.

Image 1:

Image 2:

Subtracted Image:

The matlab code that I used to amplify the subtraction was:

% Gray out the frames temp1 = rgb2gray(temp1); temp2 = rgb2gray(temp2);

    temp3 = imsubtract(temp2,temp1);
    subplot(1,2,1)

    imshow(temp3)

    level = graythresh(temp3);

    %Convert image to a black and white image using the threshold set in the
    %variable -->level  . The variable level is in the range [0,1]
    BW = im2bw(temp3,0);

    subplot(1,2,2)
    imshow(BW)

[Image Analyst]

That looks much better. Will that work for you?

[Chris]

I'm having trouble utilizing the bwmorph() functions to clean up and trace the boundary of the shock. Any pointers on where to start? I tried using the 'clean' feature of bwmorph() to try and remove some of the grainy-ness.

[Image Analyst]

Can you indicate what you would like for the ideal output? Would you like the outer boundary/envelope of the flower-like shape? Maybe try activecontour. I've attached a demo for you.

[Chris]

The ideal output would be the pixel locations of the location of the shock front. The shock front is the smooth circular bubble that surrounds the flower-like shape.

[Image Analyst]

Why don't you get the average radial profile? Just scan the image and figure out its distance from the location of the center and add in the intensity. Then plot the intensity and look for something unusual in the profile.

[Chris]

So something like setting the center by finding the pixel location of the center of the expansion, then moving radially out and making an intensity plot of the pixels along a line moving radially out from the center. Then finding where the variation occurs and calling it the 'location' of the shock?

I will definitely try that. It seems like it may be much easier to implement than tracing the boundary. My motivation for doing this is to eventually be able to have the code automatically track the location with me having to tell it where the shock is in every frame.

[Image Analyst]

Chris: Here is some code to do what I said. I also added a Savitzy-Golay filter (in case you have the Signal Processing Toolbox and want to smooth it some before you try to find the anomoloies).

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
  % User does not have the toolbox installed.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end

% Read in a standard MATLAB gray scale demo image.
folder = 'C:\Users\Chris\Documents\Temporary';
baseFileName = 'capture.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  % File doesn't exist -- didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(grayImage);
if numberOfColorBands > 1
  % It's not really gray scale like we expected - it's color.
  % Convert it to gray scale by taking only the green channel.
  grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the original gray scale image.
subplot(2, 1, 1);
imshow(grayImage, []);
axis on;
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') 

promptMessage = sprintf('Locate the center of the shock wave,\nor Cancel to abort processing?');
titleBarCaption = 'Continue?';
button = questdlg(promptMessage, titleBarCaption, 'Continue', 'Cancel', 'Continue');
if strcmpi(button, 'Cancel')
  return;
end
[x, y] = ginput(1);
hold on;
plot(x, y, 'r+', 'LineWidth', 3, 'MarkerSize', 15);

% Scan every pixel finding the distance
maxDistance = ceil(sqrt((columns)^2 + (rows^2)))
profile = zeros(1, maxDistance);
for col = 1 : columns
  for row = 1 : rows
    distance = sqrt((col - x)^2 + (row - y)^2) + 1; % Add 1 so you don't get zero.
    index = int32(distance);
    profile(index) = profile(index) + double(grayImage(row, col));
  end
end
% Crop it to actual data
lastIndex = find(profile > 0, 1, 'last');
profile = profile(1:lastIndex);
subplot(2,1,2);
plot(profile, 'LineWidth', 2);
grid on;
title('Average Radial Profile', 'FontSize', fontSize);

% Now smooth with Savitzky-golay filter
k = 1; % Order of the polynomial
windowSize = 7;
smoothedProfile = sgolayfilt(profile, k, windowSize, [], 2);
hold on;
plot(smoothedProfile, 'r-', 'LineWidth', 5);
legend('Original Profile', 'Smoothed Profile');

[Chris]

Thank for your help with this issue Image Analyst. I was eventually able to achieve everything that I wanted for my shock tracking code.

[Sampada J]

@Chris, @Image Analyst I'm Having the same problem. Even after amplifying the subtraction results, the exact boundary is not visible and visible only as distant dots (would seem noise to a layman). Would be great if you guys can share the code.

[Image Analyst]

Yes, because not much changed between those two frames. You might try a larger time difference between frames if you want to see a bigger difference.

Answer 2

[Rui Yang]

Hello Chris, I have encountered a similar problem. I have read your discussion, but can not understand. Could I have a copy of your code and pictures? My email is [e-mail].

Thank you very much for your sharing.

[Sampada J]

@Rui Yang Please share the code or your own. Would be a great help.