How can I develop a matlab code which analyses and checks how thick the cracks of a wheel surface are?
Show older comments
Recently I have done a measurement for a wheel to find out how overused it is by using the confocal sensor and the data recorder to check how deep the cracks on the surface are by the light distance.
The measured data looks like this e. g.:
With the green diagram you can see how long the distance of the measured light = how deep the cracks on the surface are. The red and the yellow ones are the high and low voltage level.
Now the challenge is that I need to find out how thin/thick the cracks are using this data. Is there anyway I can use matlab to solve this problem?
I have already written a matlab code that checks how many cracks are there by adding a threshold or a mm light distance. But still have not found out a way yet how to find out how long/thick the cracks on the surface are...
Thanks a lot in advance for the help!
% Create PNRF reader OLE Automation server
FromDisk = actxserver('Perception.Loaders.pNRF');
% Use the open file dialog to select a recording file
[FileName,PathName] = uigetfile('*.pNRF','U:\SP90000_12_12_23');
if isequal(FileName,0)
error('Benutzer hat die Auswahl abgebrochen.');
end
MyRecordingName = fullfile(PathName, FileName);
% Laden der Aufzeichnung
MyData = FromDisk.LoadRecording(MyRecordingName); % Zugriff auf den ersten Recorder
% Initialize channel arrays
channelNames = {}; % Initialisiert als leere Zelle
channel = {}; % Initialisiert als leere Zelle
% Loop through all recorders and channels to get the names of all channels
indexer = 0;
for r = 1: MyData.Recorders.Count
% Get an interface to each recorder
MyRecorder = MyData.Recorders.Item(r);
% Get the channel collection from the recorder
myChannels = MyRecorder.Channels;
% Loop through each channel in the recorder
for i = 1:myChannels.Count
myChannel = myChannels.Item(i);
channelNames{end+1} = string(myChannel.Name);
channel{end+1} = myChannel;
end
end
% Show a modal selection list with all the channels so the operator can select a channel
[indx, tf] = listdlg('ListString',channelNames);
if ~tf
error('Kein Kanal ausgewählt.');
end
myChannel = channel{indx};
SelectedChannelName = myChannel.Name;
% Now we get the DataSrc (data source) interface from the selected channel
ItfData = myChannel.DataSource(3);
% Check if we have synchronous or asynchronous data by using the TimeInfo property
myTimeInfo = ItfData.TimeInfo;
isAsyncData = ~strcmp(myTimeInfo, 'DataSourceTimeInfo_Implicit');
% Get the sweeps collection from the datasource interface
MySweeps = ItfData.Sweeps;
% Get the start and end time of the recorded data of the first channel
dStartTime = MySweeps.StartTime;
dEndTime = MySweeps.EndTime;
% Get the data segments collection between the start and end time
SegmentsOfData = ItfData.Data(dStartTime, dEndTime); % Now the actual data is read.
% Take the first data segment
myFirstSegment = SegmentsOfData.Item(1);
% Check the number of samples
NumberOfSamples = myFirstSegment.NumberOfSamples;
if NumberOfSamples > 500000
NumberOfSamples = 500000;
end
% Show the data in a diagram
figure1 = figure('NumberTitle','off','Name', SelectedChannelName);
% axes1 = axes('Parent',figure1,'Position',[0.2 0.05 0.6 0.9]);
%set(figure1, 'Position', [100 100 1200 800]);
if isAsyncData
[WaveformDateWithTimes, Times] = myFirstSegment.WaveformWithTimes(5, 1, NumberOfSamples, 1);
plot(Times, WaveformDateWithTimes);
title('Not-Equidistant data points');
else
WaveformData = myFirstSegment.Waveform(5, 1, NumberOfSamples, 1);
tEnd = myFirstSegment.StartTime + (NumberOfSamples - 1) * myFirstSegment.SampleInterval;
t = myFirstSegment.StartTime: myFirstSegment.SampleInterval : tEnd;
plot(t, WaveformData);
title('Equidistant data points');
end
legend(SelectedChannelName);
xlabel(['Time (',ItfData.XUnit ,')']);
ylabel([SelectedChannelName, ' (', ItfData.YUnit ,')']);
% Create a Textbox to display data values from the information sheet
% myDataValues = MyData.DataValues;
% cInfo = '';
% for i = 1:myDataValues.Count
% myDataValue = myDataValues.Item(i);
% cInfo = sprintf('%s%s: %s %s\n', cInfo, myDataValue.Name, myDataValue.Value, myDataValue.Units);
% end
% annotation(figure1, 'textbox', [0.66 0.05 0.32 0.9], 'String', cInfo, 'FitBoxToText', 'on');
% Schätze die Bildschirmgröße ab und positioniere das Fenster
screenSize = get(0, 'ScreenSize'); % Dies gibt [left bottom width height]
figureWidth = 1200;
figureHeight = 800;
figureLeft = (screenSize(3) - figureWidth) / 2; % Zentriert horizontal
figureBottom = (screenSize(4) - figureHeight) / 2; % Zentriert vertikal
set(figure1, 'Position', [figureLeft figureBottom figureWidth figureHeight]);
% Ersetzen oder fügen Sie nach der Plot-Erstellung hinzu
if isAsyncData
[WaveformDateWithTimes, Times] = myFirstSegment.WaveformWithTimes(5, 1, NumberOfSamples, 1);
plot(Times, WaveformDateWithTimes);
WaveformData = WaveformDateWithTimes; % Verwendung der asynchronen Daten für die Statistik
else
WaveformData = myFirstSegment.Waveform(5, 1, NumberOfSamples, 1);
tEnd = myFirstSegment.StartTime + (NumberOfSamples - 1) * myFirstSegment.SampleInterval;
t = myFirstSegment.StartTime: myFirstSegment.SampleInterval : tEnd;
plot(t, WaveformData);
end
% Statistikberechnungen
averageValue = mean(WaveformData);
minValue = min(WaveformData);
maxValue = max(WaveformData);
numberOfPositiveValues = sum(WaveformData > 0);
% Anzahl der Werte größer als 0.25 mm
numberOfValuesAboveThreshold = sum(WaveformData > 0.25);
% Anzeige der berechneten Werte
disp(['Durchschnittswert: ', num2str(averageValue)]);
disp(['Minimaler Wert: ', num2str(minValue)]);
disp(['Maximaler Wert: ', num2str(maxValue)]);
disp(['Anzahl positiver Werte: ', num2str(numberOfPositiveValues)]);
disp(['Anzahl der Werte über 0.1 mm: ', num2str(numberOfValuesAboveThreshold)]);
Accepted Answer
Joe Vinciguerra
on 20 May 2024
Edited: Joe Vinciguerra
on 20 May 2024
1 vote
"Is there anyway I can use matlab to solve this problem?"
Yes.
However, you don't seem to have all the information you need to do so. You measured the depth with a confocal sensor, but there is no indication that you have any data that can provide the width (the dimension normal to the one you measured). You need something that measures the relative motion between the sensor and the wheel surface. This can be an encoder or a tachometer for example, depending on your measurement process and your required resolution.
Once you have a sensor that provides distance data [edit to add: or some way to convert time to distance] normal to your confocal sensor you can use matlab to find the rising and falling edges of the cracks, then calcuate the distance between them.
7 Comments
Hello Dear @Joe Vinciguerra. First and foremost, thank you so much for you fast response to my issue. And I know what you mean by the speed/relative motion that is missing. Now what I can tell you is that there are some settings that we did on:
- the side of the sensor (like the measurement range) so that it can provide precise data.
- On the side of the data recorder, you can see the picture how we set the device for the measurement:
3. and the most important thing that I forgot to mention, is that I used different measurement motion values using different speed for the wheel (I placed the wheel on a three-jaw chuck motor like this one: https://www.shop.santool.de/en/clamping-tools1/chucks-zentra/three-jaw-lathe-chucks/din-6350-zentra/one-piece-jaws/steel/zentra-three-jaw-lathe-chuck-o500-mm-steel.html). And I set different speed for the motor using the XENAX Controller, using the SP (Speed) and AC (Acceleration Command):
So as far as I understand, i need to know the used speed of the Wheel so that I can use this formula:
v = d / t -> d = v*t
To find out the distance/thickness of the cracks, right?
If yes there are 2 challenges now: How to use this speed (Inc/s ecoder counter) of the SP command and how to implement it in matlab..
Joe Vinciguerra
on 20 May 2024
I'm not familiar with XENAX, but the protocol looks very similar to Galil's DMC code. If it is, then you can use the "TV" (Tell Velocity) command in your controller software to get the actual velocity. I see you are already use the "TP" command to get position, the "TV" command might work the same way, but the units it outputs depends are your particular setup (typically default to encoder counts per second).
Alternately, if you can read TP fast enough and syncronize it with the confocal data, just use that as your positional data.
Either way, it might go without saying, but don't forget to convert rotational units to linear units by using the wheel radius.
@Joe Vinciguerra Thank you for your comment and the help again. I am gonna try and connect the Xenax to my pc again and see if I can found out through TV what speed equals the used SP Value.
Otherwise I checked (with the help of chatgpt, if I can convert the SP value to the actual velocity /speed of the wheel, but I need to find the diameter of the wheel out so it might not be accurate). Here is what I got:
But either way, I will check this with the "TV" command tomorrow if it this calculation of the speed is accurate.
Now one more question, once I find out the real velocity / speed of the wheel, how do I find out how thick the cracks are that are seen in the green diagram? Does the code do the thing:
% Example data (replace with actual data acquisition code)
encoder_pulses = [100, 105, 110, 95, 90, 100]; % Encoder pulses at each time interval
surface_heights = [0.1, 0.15, 0.2, 0.5, 0.3, 0.1]; % Heights measured by surface sensor
time_interval = 0.001; % Time interval in seconds
% Encoder parameters
pulses_per_revolution = 4000;
wheel_diameter = 0.012;
circumference = pi * wheel_diameter;
distance_per_pulse = circumference / pulses_per_revolution;
% Calculate displacement at each time step
displacement = cumsum(encoder_pulses) * distance_per_pulse;
% Simple threshold-based crack detection
threshold = 0.25; % Threshold height to identify a crack
crack_indices = find(surface_heights > threshold);
% Calculate crack thickness
crack_thickness = diff(displacement(crack_indices)); % Distance between detected cracks
% Display detected cracks and their thickness
for i = 1:length(crack_indices)-1
fprintf('Crack detected at position %.3f m with thickness %.6f m\n', ...
displacement(crack_indices(i)), crack_thickness(i));
end
Joe Vinciguerra
on 20 May 2024
Did you write this? Or did ChatGPT write this? Or did you get it from somewhere else?
Ramo
on 21 May 2024
Your latest codes suggests you have encoder feedback. If that's true then you don't need velocity.
Between using AI to generate your code and not understanding what data-processing steps your project requires, I can't offer any more help. This forum is for matlab specific questions, not to explain line-by-line what the code you had generated does and if it will work for you application.
I'll just leave you with this for once you have the data you need. Here is one way to do it:
% this is just an example, you need to change this as appropriate for your
% application and requirements. I'm make a lot of assumptions for
% demonstration purposes.
% generate random sample data to represent confocal results
displacement_data = 0:0.001:0.2-0.001; % millimeters
confocal_data = zeros(size(displacement_data));
confocal_data(round(0.25*numel(confocal_data)):round(0.75*numel(confocal_data))) = 0.2;
confocal_data = confocal_data + randn(size(confocal_data))/500; % millimeters
% find the width of the crack
dConfocal = abs(diff(confocal_data));
indx = find(islocalmax(dConfocal, MinProminence=(max(dConfocal)-min(dConfocal))/2));
gap = displacement_data(indx);
width = diff(gap); % millimeters
% plot the data
figure(); nexttile(); hold on;
plot(displacement_data, confocal_data)
ylabel('crack depth (m)')
xlabel('crack width (mm)')
xline(gap)
text(0.4, 0.5, "← " + string(width) + " mm →", "units", "normalized")
Best wishes.
More Answers (1)
Image Analyst
on 22 May 2024
You can use regionprops, if you have the Image Processing Toolbox, to measure the number of cracks, their width, and the signal value in the cracks. Just pass it your segmented (thresholded) signal. regionprops works with 1-D signals too, not just 2-D images.
props = regionprops(binarySignal, originalSignal, 'Area', 'MeanIntensity')
numCracks = numel(props);
allCrackWidths = [props.Area];
allSignalsInCracks = [props.MeanIntensity];
You might have to transpose your signal since for 1-D signals I think regionprops requires either column vectors or row vectors (I forget which). So if you don't get a bunch of results, then try transposing your signals.
4 Comments
@Ramo No response so maybe you're waiting for a complete demo. See below:
% Generate test signal as a series of pulses of various widths and spacings.
% Pulses could be between 50 and 200 elements long.
meanSignals = [.3, .9];
numSamples = 2000;
signal = 0.05 * rand(1, numSamples) + meanSignals(1);
pointer = 1;
while pointer < numSamples
index1 = min(numSamples, pointer + 50 + randi(150));
index2 = min(numSamples, index1 + 50 + randi(150));
signal(index1 : index2) = signal(index1 : index2) + (meanSignals(2) - meanSignals(1));
pointer = index2;
end
plot(signal, 'b-')
grid on;
% Now we have our signal of noisy square wave with variable pulse widths
% and we can threshold and find the widths and means.
thresholdValue = 0.625; % Whatever you want.
yline(thresholdValue, 'Color', 'r', 'LineWidth', 3)
binarySignal = signal > thresholdValue;
% Get widths and mean signal in the high parts of the pulses.
props = regionprops(binarySignal, signal, 'Area', 'MeanIntensity');
numCracks = numel(props)
allCrackWidths = [props.Area]
allSignalsInCracks = [props.MeanIntensity]
% Get widths and mean signal in the low parts of the pulses.
% We just need to invert the binary signal.
propsValley = regionprops(~binarySignal, signal, 'Area', 'MeanIntensity');
numValleys = numel(propsValley)
allValleyWidths = [propsValley.Area]
allSignalsInValleys = [propsValley.MeanIntensity]
regionprops gives a one-line solution. Can't get much easier than that. You can get all the data back in a table variable if you want, rather than a structure array.
First and foremost, thanks a lot for the suggestion. My code works without any errors, I am just gonna write it in case someone in the future needs to do something similar:
% Create PNRF reader OLE Automation server
FromDisk = actxserver('Perception.Loaders.pNRF');
% Use the open file dialog to select a recording file
[FileName, PathName] = uigetfile('*.pNRF', 'Select the .pNRF File');
if isequal(FileName, 0)
error('User canceled the selection.');
end
MyRecordingName = fullfile(PathName, FileName);
% Load the recording
MyData = FromDisk.LoadRecording(MyRecordingName);
% Initialize channel arrays
channelNames = {};
channel = {};
% Loop through all recorders and channels to get the names of all channels
for r = 1: MyData.Recorders.Count
MyRecorder = MyData.Recorders.Item(r);
myChannels = MyRecorder.Channels;
for i = 1:myChannels.Count
myChannel = myChannels.Item(i);
channelNames{end+1} = char(myChannel.Name); % Convert to char to handle strings correctly
channel{end+1} = myChannel;
end
end
% Debug: Print all channel names to verify them
%disp('Available Channels:');
%disp(channelNames);
% Find the index of the required channel
distanceIndex = find(strcmp(channelNames, 'Ch B3_s'));
% Ensure that the required channel is found
if isempty(distanceIndex)
error('Distance Channel (Ch B3_s) is missing.');
end
% Extract data from the specified channel
distanceChannel = channel{distanceIndex};
distanceItfData = distanceChannel.DataSource(3);
% Get the sweeps collection from the data source interface
MySweeps = distanceItfData.Sweeps;
dStartTime = MySweeps.StartTime;
dEndTime = MySweeps.EndTime;
% Get the data segments collection between the start and end time
distanceSegments = distanceItfData.Data(dStartTime, dEndTime);
% Calculate the total number of samples in the .pNRF file
totalSamples = 0;
for k = 1:distanceSegments.Count
segment = distanceSegments.Item(k);
totalSamples = totalSamples + segment.NumberOfSamples;
end
% Display the total number of samples
disp(['Total number of samples: ', num2str(totalSamples)]);
% Limit the number of samples for plotting to avoid overloading
%if totalSamples < 500000
% NumberOfSamples = 500000;
%else
NumberOfSamples = totalSamples;
%end
% Take the first data segment and plot the data (if needed)
[distanceData, timeData] = distanceSegments.Item(1).WaveformWithTimes(5, 1, NumberOfSamples, 1);
% Assuming each sample corresponds to 1 mm (or modify accordingly)
mmPerSample = 1; % change this value if each sample corresponds to a different distance in mm
xAxis = (0:NumberOfSamples-1) * mmPerSample;
% Debug: Check if data is loaded correctly
%disp('Loaded Data:');
%disp(['Number of Samples: ', num2str(NumberOfSamples)]);
%disp('First few data points:');
%disp(distanceData(1:10)');
plot(xAxis, distanceData, 'b');
xlabel('Time (s)');
ylabel('Ch B3_s');
grid on;
thresholdValue = 0.6;
yline(thresholdValue, 'Color','r','LineWidth', 2);
binarySignal = distanceData > thresholdValue;
% Get widths and mean signal in the high parts of the pulses.
props = regionprops(binarySignal, distanceData, 'Area', 'MeanIntensity');
numCracks = numel(props);
disp('Number of Cracks:');
disp(numCracks);
allCrackWidths = [props.Area];
disp('Crack Widths:');
disp(allCrackWidths);
allSignalsInCracks = [props.MeanIntensity];
disp(allSignalsInCracks);
% Get widths and mean signal in the low parts of the pulses.
% We just need to invert the binary signal.
propsValley = regionprops(~binarySignal, distanceData, 'Area', 'MeanIntensity');
numValleys = numel(propsValley);
disp(numValleys);
allValleyWidths = [propsValley.Area];
disp(allValleyWidths);
allSignalsInValleys = [propsValley.MeanIntensity];
This is the output that I get:
Total number of samples: 13855280
Number of Cracks:
3
Crack widths:
3961 760 357
0.6240 0.6247 0.6256
4
776407 8640 5713239 7351916
The only thing that is not accurate is that the x-axis which was loaded from the pNRF (Perception Viewer Software) file, has a time unit and it's the simulation time. (See the unit in the simulation above).
Now the only hurdle that I need to check and add in to the code is the encoder data so that the cracks widths are accurate. Or does the command "allCrackWidths = [props.Area];" already does that correctly? I am sorry, I have never evaluate a measurement with matlab, so I am just not sure.
Image Analyst
on 23 May 2024
Be sure to adjust the threshold for your data. I didn't have or use your data so I just made up some and the threshold may not be best for your signal.
Thresholding tells you whether you are in a crack or not. I'm not sure if the cracks are the high signal above the threshold or the low signal below the threshold so make sure you check on that.
regionprops gives you a list of the measurements you specify for all the "runs" of where the binary signal is true/1 (in other words all runs of where you're in a crack). You can get the measurements either in a structure array, where each structure in the array has fields with the measurements you requested, OR you can get the measurements back in a table if you give 'table' as the first argument to regionprops. Assuming you get the results in a structure array, you can either use that by itself or you can extract all the fields out of the structures into a simple double vector doing
allCrackWidths = [props.Area];
If you had the result in a table you'd do
allCrackWidths = props.Area;
instead (if you wanted to extract out the widths into their own vector instead of leaving them in the table). Each element in the vector is the width for one "run" of crack. It's just called area because it's mainly used for images but you have a 1-D signal so the "Area" is really the "length" of that region.
I don't know anything about encoder data but if you want to concatenate some other data onto the width vector, you can certainly do that. You can build a matrix or table with columns representing anything you want.
Categories
Find more on Shifting and Sorting Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
