How do I calculate the radius of this culvert from the image hyberbola?

Hi everyone,

I took a GPR (Ground-Penetrating Radar) B-scan of a roadway section with a round-shaped steel culvert in it.

I am trying to use MATLAB to detect the curve and/or calculate the approximate radius from the image, the actual diameter is 15 inches.

I tried using " imfindcircles" to do it, tried different ways such as increasing the intensity, turning to grayscale, binarizing it (in order to use objectpolarity option), as well as changing the [rmin , rmax] values, but that doesnt seem to work. The figure file, and the raw image texi file are also attached.

I_scan11 = imread("Scan11_color.png");
%imshow(I_scan11);
I_scan11_gray = im2gray(I_scan11);
imshow(I_scan11_gray)
I_scan11_gray1 = imadjust(I_scan11_gray);
imshow(I_scan11_gray1);
%imhist(I_scan11_gray)

 %I_scan11_BW = I_scan11_gray > 210;
 I_scan11_BW = imbinarize(I_scan11_gray);
 imshow(I_scan11_BW);
% 
% SE = strel("disk",14);  % 20
% I_scan11_BW1 = imdilate(I_scan11_BW,SE);
% %I_scan11_BW1 = imclose(I_scan11_BW,SE);
% %I_scan11_BW1 = bwareaopen(I_scan11_BW1,5);
% imshow(I_scan11_BW1);

[center, radius] = imfindcircles(I_scan11_BW,[400,1120],"Sensitivity",0.99,"ObjectPolarity","dark")
viscircles(center,radius);
% hold on;
% plot(center(:,1),center(:,2),'yx','LineWidth',2);
% hold off;

The following code generates the image as seen above: (exclusive for GPR)

load MECHANICAL_LAB_EXP_05_11_21_022_P_1A.txt;  % Image text file derived from the device
Bscan = MECHANICAL_LAB_EXP_05_11_21_022_P_1A; % This is scan 11-1. 
clear MECHANICAL_LAB_EXP_05_11_21_022_P_1A;
% Apparent/input dielectric constant = 6.0, subsurface material is
% non-magnetic ( magnetic permeability = 1)
[a, b] = size(Bscan);   % a = cross-range (rx), b = range(r)   (rx,r) = (a,b)
dr = (9 * 12) / (b-1);    % Range resolution, in inches % Penetration depth is 9 ft
d = 4.11 * 12;            % GPR cart rolling distance, in inches  % the rolling distance is 4.11 ft
drx = d / (a - 1); % Cross-range resolution, in inches

r = (0:b-1) * dr;     % Range, in inches
rx = (0:a-1) * drx; % Cross-range, in inches


figure; 
imagesc(rx, r, Bscan'); colormap jet; colorbar; caxis([-2.5e6,3e6]);
xlabel('Cross-range, {\it r_x} (in)'); ylabel('Range, {\it r} (in)');
title(scan_title);
xlim([d-33 d]); ylim([0,22]);
set(gcf, "Position", [744 757 560 420]);
set(gca, 'fontsize', 12); 

Please help me! Thanks!

[DGM]

I'm not sure I understand what I'm looking at. If the projection of a circular object is a hyperbola, then trying to work off of the curvature might not be as easy as trying to find the asymptotes and minimum range. I don't know if that's appropriate for the available data though. I can't really tell if the image shows a part of a circle or part of a hyperbola.

[Koosha]

Since the GPR sends EM waves and a steel culvert is a conductor, we have 100% reflection and nothing penetrates below, and metallic objects (culvert, rebar..etc) are depicted as hyperbolas. The range axis depicts the B scan slice depth and the cross-range is the range I rolled the device on the ground, perpendicular to the culvert.

To summarize, the Image is showing a hyperbola, probably thats why the imfindcircles did not work.

Now I am trying to empirically estimate the culvert radius based on the projected hyperbola curve.

Answer 1

[William Rose]

@DGM, this is an inteesting and non-trivial problem. Here are suggesitons, based on my experience using the image processing tools in Labview to determine vl=blood vessel radius from ultrasound images. I have not used the Matlab image processing tools, so my suggesitons are not detailed.

Use the values from the b&w image. Let the user set a depth below which to analyze. Set all array values at shallower depths to zero to avoid getting confused later by the high reflactance values near the surface. In the image you provided, a dpeth of 15" or 16" might be a reasonable depth setting, since the top of the pipe seems to be just below 16". THen you find the depth coordinate for each "column" of the array which is maximum for that column. This is the vertical coordinate of the hyperbola, for that x coordinate. One you have done all that, you have one depth for every x-coordinate. The yu use least squares to find the coefficients of hte hyperbola that give the best fit to the x,y set of points. Then you convert the hyperbola coefficients to a radius.

Good luck.

[William Rose]

[edit- corrected typos]

The color figure which you uploaded is a colorized version of an array of intensities. I recommend that you work with the intensity array itself, instead of converting the colorized version back to black and white, which can introduce inconsistencies. If you upload the intensity array, I will show you what I mean in my earlier answer.

Are you sure the GPR does not do a correction for angle, as part of its internal processing of a B-scan? I assume that the image is generated by a phased emission of a signal from a flat transducer placed on the surface. The software adjusts the phasing of an array of individual transmitters to steer the beam back and forth across a range of angles. If the direction of "down" is known (and it is - its perpendicular to the transmitter plate, usually), then it is possible to do an inverse cosine correction for the instantaneous beam angle, so that the image presented is already depth-corrected. Are you confident that this has not already been done in the internal software? Ask the GPR manufcturer, if you don't know.

[Koosha]

William, this is not a GPR problem, it is merely a MATLAB problem.

As for the depth, the actual culvert is located at 25 inches, and the Image is very close to that (18~20), since the apparent dielectric constant input in the device is very close to that of the medium

d = DeltaT.c/2.sqrt(e') , supposed that material is non-magnetic (magnetic permeability = 1). 2d/deltaT = phase velocity << c (speed of light in open space), the 2 factor is for incidence and reflection time.

I am a researcher in NDT/E (nondestructive testing and evaluation), and this problem has not been solved anywhere, neither does the manufacturer know.

The text matrix of the images were transfered to MATLAB and generated into an image there.

I hope that answered your question.

[William Rose]

@Koosha Raisi, thank you. Can you attach the array of intensity values?

Have you worked out the analytical mathematics of the transformation or a circle into its equivalent hyperbola? YOu need to do that to nvert the projection.

[Koosha]

Thank you for your reply, please find the array attached as requested.

[Koosha]

@William Rose Thank you, the code works perfectly.

[William Rose]

@SERG,

Great! I'm glad to hear it. I included a lot of comments in the code, to explain what it is happening.

The part which is hardest to understand, in my opinion, is the section where I use logical indexing to only keep ridge points that are not on the bottom edge of the image: lines 3,4,5 below.

[Bpk,jpk]=max(Bscan1');  %find max point in each column of Bscan1'
[row,col]=size(Bscan1);
Bpk=Bpk(jpk<col);   %discard peak values along the bottom edge
xpk=x1(jpk<col);    %x values that go with the peak values
jpk=jpk(jpk<col);
ypk=(jpk+jlo-2)*dy; %y  values that go with the peak values

I want to eliminate peak points on the bottom edge of the image, because they probably do not represent ridge points but rather points on the side of a hill. Logical indexing is a way to do it. It has taken me a while to understand how to use logical indexing in Matlab. This is one example of how logical indexing is useful.

[Koosha]

Thank you @William Rose ! it was a very pristine and well-explained code. I am actually busy for a few deadlines but I will come back to you if I had any problems understanding the indexing part.

Either way, you have solved my problem and thank you so much!

[William Rose]

You're welome @SERG!

Answer 2

[Image Analyst]

Do you just want to find the peak/ridge line of the curved thing at the bottom of the grayscale image? I could probably find that and then to get the radius of curvature at any point, you just have to fit the data of the ridgeline (x,y) coordinates into a circle fitting function. Of course the radius of curvature changes as you move along the ridgeline so that's why you have to specify both the center location and the window width to get the radius of curvature just within that window.

Like @William Rose asked, please attach your actual data (not a pseudocolored image) in a .mat file with the paperclip icon.

[Koosha]

Thank you for your answer, please find the text file of the Image attached.

[William Rose]

@Koosha Raisi,

I downloaded the text file and loaded into the Matlab workspace. When I plot it, it does not look like the image you provided earlier. There is no hyperbola or circle, and the range of values indicated on the color bar is different by about a factor of 10.

a=load('MECHANICAL_LAB_EXP_05_11_21_022_P_1A.txt','-ascii');
surf(a,'EdgeColor','none'); view(0,90); colorbar

This indicates that there is something incorrect about the format of the text file, or it indicates that the data in this text file is not the data used to generate the original image you posted. The text file has 84 rows and 512 columns. Is that the resolution of your GPR image? Maybe if we reformat the array into a different number of rows and columns, the image will look right. Do you have a suggesiton?

[Koosha]

Dear William, sorry for not saying earlier, the code used for image generation is as follows:

[a, b] = size(Bscan);   % a = cross-range, b = range  
dr = (9 * 12) / (b-1);    % Range resolution, in inches
d = gpr_roll * 12;            % GPR cart rolling distance, in inches
drx = d / (a - 1); % Cross-range resolution, in inches

r = (0:b-1) * dr;     % Range, in inches
rx = (0:a-1) * drx; % Cross-range, in inches


figure; 
imagesc(rx, r, Bscan'); colormap jet; colorbar; caxis([-2.5e6,3e6]);
xlabel('Cross-range, {\it r_x} (in)'); ylabel('Range, {\it r} (in)');
title(scan_title);
xlim([d-33 d]); ylim([0,22]);
set(gcf, "Position", [744 757 560 420]);
set(gca, 'fontsize', 12); 

This is the code used for image generation of GPR. Bscan is basically the image text file.

[William Rose]

@Koosha Raisi,

You said in your comment "this problem has not been solved anywhere". Are you sure? Have you read "A Review on Migration Methods in B-Scan Ground Penetrating Radar Imaging", 2014? Have you read Schneider (1978), “Integral formulation for migration in two and three dimensions”, or some of the 1100 papers which have cited it? I recommend looking at equation 1 in the 2014 paper, as it applies to your problem.

I am not familiar with GPR, but by glancing at the 2014 paper, I see the reciever moves along the ground to generate a B-mode image. The solution for the B-mode image of a circular reflective object is not complicated. It is probably better if you work it out, so you understand the math. Draw a picture and do the geometrical calculations to find the distance to the reflective circle, as a function of the x-position of the transceiver. This will yield an equation for a hyperbola which involves the pipe radius R.

[William Rose]

@Koosha Raisi, please save Bscan by doing

save('Bscan');

This will generate file Bscan.mat. Then attach Bscan.mat to your response. Thank you.

[Koosha]

What is the problem ? load the text file and assign it to a variable (e.g., bscan) and use the code I provided to generate the Image, you can play with the intensity of the colorbar as well, if needed. thank you for pointing out the paper, I will look it up. meanwhile, i will appreciate further assistance or ideas.

[William Rose]

@Koosha Raisi,

I did load the text file. It makes the image I showed you, which looks nothing like the image you presented. The image made by the text fie does not have a hyperbolic ridge so it cannot be used to fit a hyperbola or estimate a radius. That is the problem.

[William Rose]

@Koosha Raisi, when I run your code I get an error:

"Unrecognized function or variable 'gpr_roll'."

[yanqi liu]

may be set gpr_roll can get different view，such as

Bscan=load('MECHANICAL_LAB_EXP_05_11_21_022_P_1A.txt','-ascii');

[a, b] = size(Bscan);   % a = cross-range, b = range  
dr = (9 * 12) / (b-1);    % Range resolution, in inches
d = 2 * 12;            % GPR cart rolling distance, in inches
drx = d / (a - 1); % Cross-range resolution, in inches

r = (0:b-1) * dr;     % Range, in inches
rx = (0:a-1) * drx; % Cross-range, in inches


figure; 
imagesc(rx, r, Bscan'); colormap jet; colorbar; caxis([-2.5e6,3e6]);
xlabel('Cross-range, {\it r_x} (in)'); ylabel('Range, {\it r} (in)');
% title(scan_title);
xlim([0 d]); ylim([0,22]);
% set(gcf, "Position", [744 757 560 420]);
set(gca, 'fontsize', 12); 

what is target？get the curve by fit？may be just draw on the figure to explanation

[William Rose]

@Koosha Raisi,

You wrote "meanwhile, i will appreciate further assistance or ideas". What did you do with the ideas you have already recieved?

1. @Image Analyst had an exclelent suggesiton and gave you some code. What have you done with it?
2. I provided pseudo-code in my initial answer. Have you tried to implement it?
3. Decide which of the reconstruction algorithms in the 2014 paper you want to use. Write some code to implement the algorithm. Then people on this site will be more than happy to assist you.
4. I recommended that you derive the equation for the hyperbola that results from a B-mode scan of a circular object. It is an exericse in advanced high school analytic geometry. Did you do it? It is a key step in solving your original request, which is to esitmate the radius.

[Koosha]

The GPR rolling distance is 4.11 ft, I forgot to put it. i fixed the first post.

[William Rose]

@Koosha Raisi,

Becaue your code gave an error due to gpr_roll being udefined, I assigned a value to gpr_roll. It is evident that gpr_roll is a distance in feet. I tried a value of 4. Then the code has another error:

Unrecognized function or variable 'scan_title'.

which I fixed by defining scan_title. Then I get an image like the one you originall posted.

[William Rose]

@Koosha Raisi,

The data file you provided has much more data that what your plot shows. The dynamic range and the vertical and horizontal extents of the text file data are all significantly much larger than what your plot shows. Here is the full image, on the left, and the image you generate, on the right. For the full image, I commented out commands caxis(), xlim(), and ylim(). I added a red rectangle to the full image to show the region covered by your image. Ths explains the earlier discrepancy. Don;t make the people who are trying to help you do a bunch odetective work, please.

Knowing this, we can extract the subset of data from Bscan that correspond to the region of interest. Then we can take the next steps.

[William Rose]

Forgot to attach the image I mentioned. Sorry about that. Here it is. See comment above for explanation.

[William Rose]

@Koosha Raisi,

The attached edited script extracts the data for the region of interest and finds the points along the ridge within that region. It generates a plot of the data pouints found. See below. It is a 3D plot, so if you run the code you will be able to roatet the plot in 3D and you will see that the red points are along the ridge.

The next step is to find the hyperbola that gives the best fit to these points, and use the fit to esitmate the pipe radius. You are almost there. Have you figured out the equation for the hyperbola? Assume the transceiver moves along the horizontal axis, at y=0. Assume the pipe center is at x0, y0 and has radius R. What is the distance r(x) from the transceiver to the closest point on the pipe, as the transceiver moves along? The is the range that will be recorded by the transceiver when it is at position x. (The transceiver's y cooridnate is always 0, by our choice of coordinate system.) Hint: the equation relating r(x) and x is the equaltion of a hyperbola, and it involves the constants x0, y0, and R.

The final step is to find the values of x0, y0, and R that give the best fit of the hyperbola to the ridge points. Then we will know R, and as a bonus, we will know the pipe center coordinates. We can use fmincon() or lsqcurvefit() to estimate x0, y0, R.

[William Rose]

@Koosha Raisi,

Script finds peak points in the B-mode scan region of interest.

Script finds the hyperbola that gives the best fit to the peak points in the region of interest.

Script estimates the pipe radius and location that would produce the best fit hyperbola.

Script makes three plots, including the plot below.

[Rik]

[removed by request; water under the bridge]