Straightening Image from RF data

I have created an image in matlab from a set of data which I obtained experimentally. The issue I have is the tissue which I imaged is curved as it floats in water, this is seen the image I have.

I wrote a code to threshold the image so that it straightnes out but it does not seem to solve the problem. The thresholding just pushes the image on to the top and does not seem to do well.

Could you please give me some tips to solve this problem. The image screen shot and the thresholding technique I am using is shown below:

section   = input('Size of each tissue section (mm): ');
sectionrows = (section*(rsize))/(newdepth);
sectionrows=round(sectionrows);
nsections = (sectionrows);
nsections =round(nsections);
reducedrsize = ceil(rsize/nsections);
remainder = mod(rsize,nsections);
Row = zeros(reducedrsize,colsize);
threshold =0.02;
newimage = zeros(rsize,colsize);

for j= 1:colsize
    k = find(hil1(:,j)>threshold);
    rowindex = min(k);
    newimage(1:rsize-rowindex+1,j) = hil1(rowindex:end,j);

   for i = 1:reducedrsize-1
      Row(i,j)=max(hil1((i-1)*nsections+1:(i*nsections),j));
  end
  end

Sorry I cant seem to load the image here.

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/7924-where-can-i-upload-images-and-files-for-use-on-matlab-answers

[Image Analyst]

That's no good. We'd first have to do color segmentation. Just simply upload the pure gray scale image. I don't want pseudocoloring, tick marks, etc. Nothing except your original monochrome image.

Anyway, WHY does it need to be straightened? I bet you could make the measurements you need even though the tissue is curved. What do you want to measure? Length, area, intensity, or something else?

Answer 1

[Image Analyst]

Do you have a calibration from RF signal to kilograms, and from pixels to µm? And are you sure you really need to straighten the thing? Why isn't just going column by column enough? If you did, I'd probably start by thresholding and then going column by column to get the midpoint, like this code:

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

% Read in a gray scale demo image.
folder = 'C:\Users\Srikanta\Documents';
baseFileName = 'tissue.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);
% Crop it.
grayImage = rgb2gray(imcrop(grayImage, [90 140 400 110]));
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows columns numberOfColorBands] = size(grayImage);
% 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') 

% Let's compute and display the histogram.
[pixelCount grayLevels] = imhist(grayImage);
% Suppress the 0 bin:
pixelCount(1) = 0;
subplot(2, 2, 2); 
bar(pixelCount);
grid on;
title('Histogram of original image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

% Get the binary image
binaryImage = grayImage > 10;
% Display the original gray scale image.
subplot(2, 2, 3);
imshow(binaryImage, []);
title('Binary Image', 'FontSize', fontSize);

% Smooth it out some and clean it up.
binaryImage = imclose(binaryImage, true(5));
binaryImage = imfill(binaryImage, 'holes')
subplot(2, 2, 4);
imshow(binaryImage, []);
title('Cleaned Binary Image', 'FontSize', fontSize);

% Get the midpoints
topRow = zeros(columns, 1);
bottomRow = zeros(columns, 1);
midpoints = zeros(columns, 1);
for column = 1 : columns
  oneColumn = binaryImage(:, column);
  if sum(oneColumn) > 0
    topRow(column) = find(binaryImage(:, column), 1, 'first');
    bottomRow(column) = find(binaryImage(:, column), 1, 'last');
    midpoints(column) = (topRow(column) + bottomRow(column)) / 2;
  end
end
hold on;
plot(1:columns, midpoints, 'r-', 'LineWidth', 2);

After that it gets tricky, though I've done it. You can go along the midpoint curve fitting a short stretch of it to a slowly varying curve and take perpendicular lines to get the profile. I posted code for that here just within the last two months or so.

[Srikanta Sharma]

hi Image Analyst, Thank you for you time and this very clever, amazing code. This works very well with gray scale images, The RF data I have collected are voltage amplitudes from which I calculate the density. I tried Thresholding as you suggested but it seems to produce a image as shown here (it is of no use) is http://imageshack.us/photo/my-images/42/tissue2.png/ tissue2 is threshold image and Tissue 1 is the RF image. Tissue1 http://imageshack.us/photo/my-images/4/tissue1.png/

Much to my dismay I have not been able to crack it very well. I do not know the pixel values as you very smartly pointed out. My control signal is from a quartz flat which I use to calculate the relative density. I have posted the code which I am using as a thresholding technique to straighten the tissue, but drawing a curve that follows the tissue seems to be a better option after seeing your code.

The RF values are very important to me as my algorithm which I have written to compute the density column wise requires RF values. I do not think I can derive density values from the binary Image itself?

Also could you please tell me as to how I can use the same technique for the RF data and then plot like the way you have done it? (I am sorry if i have asked a dumb question).

The code I have for thresholding and reducing the data by crawling through the rows for all columns at a known length is shown below, Row corresponds to the reduced column image which is obtained after going through 200 rows in steps. I can calculate the length of the row(mm) based on sampling time and x-origin derived from the oscilloscope which is my metadata.

new image corresponds to the threshold image. Many thanks.

section = input('Size of each tissue section (mm): '); sectionrows = (section*(rsize))/(newdepth); sectionrows=round(sectionrows); nsections = (sectionrows); nsections =round(nsections); reducedrsize = ceil(rsize/nsections); remainder = mod(rsize,nsections); Row = zeros(reducedrsize,colsize); threshold =0.025; threshold1 = 0.019; newimage = zeros(rsize,colsize);

for j= 1:colsize k = find(hil1(:,j)>threshold); rowindex = min(k); newimage(1:rsize-rowindex+1,j) = hil1(rowindex:end,j);

 for i = 1:reducedrsize-1
    Row(i,j)=max(hil1((i-1)*nsections+1:(i*nsections),j));
end
end

[Image Analyst]

Try this code. It now works with your second image. It's like the first code except that the cropping is different (I'm not sure why you can't upload just the image and not the figure with tick marks and axes lines), and it also measures the average width of the tissue and gives the average gray level (RF Signal) along the spine.

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

% Read in a standard MATLAB gray scale demo image.
folder = 'C:\Users\Mark\Documents\Temporary';
baseFileName = 'tissue1.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);
% Display the original RGB image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize);
% Convert it to gray scale and Crop it.
grayImage = rgb2gray(grayImage);
grayImage = imcrop(grayImage, [74 84 425 240]);
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows columns numberOfColorBands] = size(grayImage);
% Display the cropped gray scale image.
subplot(2, 2, 2);
imshow(grayImage, []);
title('Cropped 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') 

% Let's compute and display the histogram.
[pixelCount grayLevels] = imhist(grayImage);
% Suppress the 0 bin:
pixelCount(1) = 0;
subplot(2, 2, 3); 
bar(pixelCount);
grid on;
title('Histogram of original image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

% Get the binary image
binaryImage = grayImage > 10;
% Display the original gray scale image.
subplot(2, 2, 4);
imshow(binaryImage, []);
title('Binary Image', 'FontSize', fontSize);

% Smooth it out some and clean it up.
binaryImage = imclose(binaryImage, true(5));
binaryImage = imfill(binaryImage, 'holes');
figure;
subplot(2, 2, 1);
imshow(binaryImage, []);
title('Cleaned Binary Image', 'FontSize', fontSize);

% Find and keep just the largest blob.
labeledImage = bwlabel(binaryImage);
measurements = regionprops(labeledImage, 'area');
allAreas = [measurements.Area];
[sortedAreas indexes] = sort(allAreas, 'descend');
% Get justthe largest blob.
largestBlob = ismember(labeledImage, indexes(1));
binaryImage = largestBlob > 0; % Binarize the labeled image.

subplot(2, 2, 2);
imshow(binaryImage, []);
title('Cleaned Binary Image', 'FontSize', fontSize);

% Get the midpoints
topRow = zeros(columns, 1);
bottomRow = zeros(columns, 1);
midpoints = zeros(columns, 1);
for column = 1 : columns
  oneColumn = binaryImage(:, column);
  if sum(oneColumn) > 0
    topRow(column) = find(binaryImage(:, column), 1, 'first');
    bottomRow(column) = find(binaryImage(:, column), 1, 'last');
    midpoints(column) = (topRow(column) + bottomRow(column)) / 2;
  end
end
hold on;
% Get rid of points that are not in the object.
x = 1:columns;
noObjectIndexes = midpoints == 0;
midpoints(noObjectIndexes) = [];
x(noObjectIndexes) = [];
plot(x, midpoints, 'r-', 'LineWidth', 2);
% Now show over original image
% Display the cropped gray scale image.
subplot(2, 2, 3);
imshow(grayImage, []);
title('Cropped Grayscale Image', 'FontSize', fontSize);
hold on;
plot(x, midpoints, 'r-', 'LineWidth', 2);

% Now make measurements.
area = sum(binaryImage); % Blob area in pixels.
% Calculate the Euclidean distance transform.
edt = bwdist(~binaryImage);
subplot(2, 2, 4);
imshow(edt, []);
title('Euclidean distance transform', 'FontSize', fontSize);
% Get the half-widths by extracting the EDT along the spine.
% Compute the average RF gray level of the original image along the spine also.
numberOfPoints = length(midpoints);
halfWidths = zeros(numberOfPoints, 1);
sumOfGrayLevels = 0;
for k = 1 : numberOfPoints
  row = int32(midpoints(k));
  col = int32(x(k));
  halfWidths(k) = edt(row, col);
  sumOfGrayLevels = sumOfGrayLevels + double(grayImage(row, col));
end
averageWidth = 2 * mean(halfWidths);
averageGrayLevel = sumOfGrayLevels / numberOfPoints;
message = sprintf('The average width = %.4f\nThe average RF signal along the spine = %.4f\n', averageWidth, averageGrayLevel)
uiwait(msgbox(message));

[Image Analyst]

You have the midpoint row. To get the column of the image above that, just do

upperPartOfColumn = grayImage(1:midpoints(columnNumber), columnNumber);

Answer 2

[Darshan]

<http://mathworks.com/matlabcentral/images/surf.gif>>>>