How can I develop a matlab code which analyses and checks how thick the cracks of a wheel surface are?

15 次查看（过去 30 天）
Ramo
Ramo 2024-5-20
评论： Image Analyst 2024-5-23
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)]);

请先登录，再回答此问题。

采纳的回答

Joe Vinciguerra
Joe Vinciguerra 2024-5-20
编辑：Joe Vinciguerra 2024-5-20
"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 个评论
显示 5更早的评论
Ramo
Ramo 2024-5-22
编辑：Ramo 2024-5-22
@Joe Vinciguerra I am still not sure whether the code does what it have to do. My actual code still don't show the number of the cracks that I am expecting. Something is still missing.., and the TV command gives me the velocity that I selected or put with putty. which I already know.. I am kind a stuck somehow. The code works, but I don't understand the results that it shows:
% 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 Recording 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 channels
highVoltageIndex = find(strcmp(channelNames, 'Ch B1_A'));
lowVoltageIndex = find(strcmp(channelNames, 'Ch B2_B'));
distanceIndex = find(strcmp(channelNames, 'Ch B3_s'));
% Ensure that all required channels are found
if isempty(highVoltageIndex)
error('High Voltage Channel (Ch B1_A) is missing.');
end
if isempty(lowVoltageIndex)
error('Low Voltage Channel (Ch B2_B) is missing.');
end
if isempty(distanceIndex)
error('Distance Channel (Ch B3_s) is missing.');
end
% Extract data from the specified channels
highVoltageChannel = channel{highVoltageIndex};
lowVoltageChannel = channel{lowVoltageIndex};
distanceChannel = channel{distanceIndex};
highVoltageItfData = highVoltageChannel.DataSource(3);
lowVoltageItfData = lowVoltageChannel.DataSource(3);
distanceItfData = distanceChannel.DataSource(3);
% Get the sweeps collection from the datasource interface
MySweeps = distanceItfData.Sweeps;
dStartTime = MySweeps.StartTime;
dEndTime = MySweeps.EndTime;
% Get the data segments collection between the start and end time
highVoltageSegments = highVoltageItfData.Data(dStartTime, dEndTime);
lowVoltageSegments = lowVoltageItfData.Data(dStartTime, dEndTime);
distanceSegments = distanceItfData.Data(dStartTime, dEndTime);
% Check the number of samples
NumberOfSamples = distanceSegments.Item(1).NumberOfSamples;
if NumberOfSamples > 500000
NumberOfSamples = 500000;
end
% Take the first data segment
highVoltageData = highVoltageSegments.Item(1).Waveform(5, 1, NumberOfSamples, 1);
lowVoltageData = lowVoltageSegments.Item(1).Waveform(5, 1, NumberOfSamples, 1);
distanceData = distanceSegments.Item(1).Waveform(5, 1, NumberOfSamples, 1);
% Plot the data to visualize
figure;
subplot(3,1,1);
plot(highVoltageData);
title('High Voltage (Ch B1_A)');
subplot(3,1,2);
plot(lowVoltageData);
title('Low Voltage (Ch B2_B)');
subplot(3,1,3);
plot(distanceData);
title('Distance (Ch B3_s)');
% Debug: Print range and sample values from distanceData
disp(['Distance Data Range: min = ', num2str(min(distanceData)), ', max = ', num2str(max(distanceData))]);
disp('Sample Distance Data Values:');
disp(distanceData(1:10)); % Display first 10 values for inspection
% Encoder parameters
pulses_per_revolution = 2000;
wheel_diameter = 0.09;
circumference = pi * wheel_diameter;
distance_per_pulse = circumference / pulses_per_revolution;
% Calculate displacement at each time step based on distance data
displacement = cumsum(distanceData) * distance_per_pulse;
% Simple threshold-based crack detection
threshold = 0.6; % Threshold height to identify a crack
crack_indices = find(distanceData > 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
% Highlight detected cracks on the plot
figure;
plot(displacement, distanceData);
hold on;
yline(threshold, 'r--', 'Threshold'); % Plot a horizontal line at the threshold
for i = 1:length(crack_indices)
plot(displacement(crack_indices(i)), distanceData(crack_indices(i)), 'ro'); % Mark cracks with red circles
end
hold off;
title('Distance Data with Detected Cracks');
xlabel('Displacement (m)');
ylabel('Distance (mm)');
% Perform statistics calculations
averageValue = mean(distanceData);
minValue = min(distanceData);
maxValue = max(distanceData);
numberOfPositiveValues = sum(distanceData > 0);
numberOfValuesAboveThreshold = sum(distanceData > threshold);
% Display calculated values
disp(['Average Value: ', num2str(averageValue)]);
disp(['Minimum Value: ', num2str(minValue)]);
disp(['Maximum Value: ', num2str(maxValue)]);
disp(['Number of Positive Values: ', num2str(numberOfPositiveValues)]);
disp(['Number of Values Above Threshold: ', num2str(numberOfValuesAboveThreshold)]);
% Display the crack widths in mm
crackWidths_mm = crack_thickness;
disp('Crack Widths (in mm):');
disp(crackWidths_mm);
Joe Vinciguerra
Joe Vinciguerra 2024-5-22
Ran in:
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.

请先登录，再进行评论。

更多回答（1 个）

Image Analyst
Image Analyst 2024-5-22
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 个评论
显示 2更早的评论
Ramo
Ramo 2024-5-23
编辑：Ramo 2024-5-23
@Image Analyst
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
Image Analyst 2024-5-23
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Sensor Models 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by