matlab determine a curve

Imagine i have an original plot with superimposed '8' shape curves (the curves are all different but similar to the '8' shape). Does anyone know how to obtain a mean curve having a matrix with the correspondent x,y points from the original plot? I mean, I pretend a medium single curve.

Any code or just ideas would be very very helpful for me. Thank you very very much!

[Image Analyst]

You'd need to upload a picture. Also, do any of the curves overlap along any part of themselves? If so, how can we take that into account unless you have all the individual curves? For example 10 curves perfectly overlapped plus one slightly shifted will appear as just two curves, but obviously you don't want to average just two sets of coordinates.

[Image Analyst]

Do you have the individual curves? You should, if you plotted it.

[joo]

Yes, I have. But it is like this: each curve begins and ends in different points. imagine: the first one begins at (132,152) and ends at (128,162); and the second begins at (153,151) and ends at (132,166). so each curve is different from the other. and this information is took manually by me because what i have is two columns: one with x and one with y. (if you want me to give you the excel file, it would be my pleasure!) thank you very much again! –

[Image Analyst]

Interesting question. I can't think of anyway off the top of my head that's simple. All I can think of are imaging-based approaches like assigning them to pixels, adding together, blurring the image, and trying to find medial axis. Or maybe you can align them all by searching for the left-most point and using circshift to wrap around all the rows of your arrays so that the left-most point is always the top row. That might work for certain situations, like the shapes in your demo.

Maybe you should ask in the newsgroup. Roger Stafford monitors that, and he is an excellent mathematician and he'd probably know of a method right away.

Go ahead and upload your workbook to your favorite filesharing website - one that does not require me to register for anything to download it.

[joo]

i tought about the problem and i concluded i really cannot divide each cycle accurately and without mistake. so, i cannot say to which of the curves a certain point (x, y) belongs to. what i have is only a column with the x points and another with the y points. thank you very very much!

[Walter Roberson]

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

[Image Analyst]

Can you post your code where you use xlsread to read it in and loop over each curve to find the left-most part and record its row? And your code to use circshift to align them? Can you do any of that to help us help you?

Answer 1

[Image Analyst]

You have a little noise on your y-data, as you can see from this example code. You might need to smooth that out first.

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;

folder = 'C:\Users\joo\Documents\Temporary';
fullFileName = fullfile(folder, 'livro1.xlsx');
% 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

[num txt raw] = xlsread(fullFileName);
x = num(:, 1);
y = num(:, 2);

% Plot x
subplot(2,2,1);
plot(x, 'b-');
xlabel('Element Number', 'FontSize', fontSize);
ylabel('X', 'FontSize', fontSize);
title('X', 'FontSize', fontSize);
grid on;

% Plot y
subplot(2,2,2);
plot(y, 'b-');
xlabel('Element Number', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
title('Y', 'FontSize', fontSize);
grid on;

% Plot curves
subplot(2,2,3);
plot(num(:, 1), num(:, 2), 'b-');
title('X', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
grid on;

% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);

% Calculate distance from (50,360) of every point.
distances = sqrt((x - 50).^2 - (y - 360).^2)
subplot(2,2,4);
plot(distances, 'b-');
title('Distance from (50, 360)', 'FontSize', fontSize);
xlabel('index', 'FontSize', fontSize);
ylabel('Distances', 'FontSize', fontSize);
grid on;

% Find valleys in distances to find out where it turned around.
[valleys, indexesAtMins] = findpeaks(-distances)

% Plot them.
hold on;
plot(indexesAtMins, -valleys, 'r*');

[joo]

hello! thank you so much. i am sorry but i don't quite understand (1)why you Calculate distance from (50,360) of every point. and then (2)what is the final objective of [valleys, indexesAtMins] = findpeaks(-distances) ? can you also tell me (3)what should be the reasoning for the next steps? //sorry to bother you so much but as i am new with matlab i have many difficulties, besides the code also in thinking about the steps to find the solution ... thank you so so much. this is very important for me.

[Image Analyst]

I did that because hte curves were tilted so taking the right most or left most x value didn't quite get the end of the curve. That may have been because of the noise though. Look at the plot of y alone - see that little glitch at the bottom of the y curves? That needs fixing, or at least you need a peak finder that is more robust about ignoring small noise spikes like that.

[joo]

and why do you Calculate distance from (50,360) of every point? thank you very much!

[Image Analyst]

I just answered that. Run my code and look at the graph and I think you'll see why.

Answer 2

[joo]

image analyst i am sorry to insist but i have been working all these days in this and still i can't understand. i see that in the plot you can vizualize (50,360). but why are these numbers and not for example (55,365)-with this pair the results would be different... i really don't understand where do they come from since i investigated all the results and code and nothing related with (50,360)... i am stucked in this for days. i am so grateful for your help. this is very important for me. thank you so much for your attention.

[Image Analyst]

Sure, you can use that point. It doesn't really matter. Of course the distances will be different but that doesn't matter when all you need to do is to determine what element the loop turns around at. You can pick any point that's sort of off the end of the loop end. If you want, simply fit all the points to a line with polyfit() and then use polyval() to find a point exactly along the mean axis of all the curves. You can do that but you'd probably get pretty much the same elements that define where the curve turns around at. All right, here -- I did the whole thing for you.

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;

folder = 'C:\Users\Mark\Documents\Temporary';
baseFileName = 'livro1.xlsx';
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

[num txt raw] = xlsread(fullFileName);
x = num(:, 1);
y = num(:, 2);

% Plot x
subplot(2,2,1);
plot(x, 'b-');
xlabel('Element Number', 'FontSize', fontSize);
ylabel('X', 'FontSize', fontSize);
title('X', 'FontSize', fontSize);
grid on;

% Plot y
subplot(2,2,2);
plot(y, 'b-');
xlabel('Element Number', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
title('Y', 'FontSize', fontSize);
grid on;

% Plot curves
subplot(2,2,3);
plot(num(:, 1), num(:, 2), 'b-');
title('X', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
grid on;

% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);

% Calculate distance from (50,360) of every point.
distances = sqrt((x - 450).^2 + (y - 550).^2)
subplot(2,2,4);
plot(distances, 'b-');
title('Distance from (450, 550)', 'FontSize', fontSize);
xlabel('index', 'FontSize', fontSize);
ylabel('Distances', 'FontSize', fontSize);
grid on;

% Find valleys in distances to find out where it turned around.
[peakValues, indexesAtPeaks] = findpeaks(distances)
% Take just the "good" peaks = those that are higher than
goodPeaks = peakValues > 300;
peakValues = peakValues(goodPeaks);
indexesAtPeaks = indexesAtPeaks(goodPeaks);
% Plot them.
hold on;
plot(indexesAtPeaks, peakValues, 'r*');

% We only want to go between the peaks because the
% data before the first peak and after the last peak
% are only partial peaks.
% Average the x and y between the peaks
figure;
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);

sumX = zeros(1,100);
sumY = zeros(1,100);
count = 0;
for peak = 1 : length(indexesAtPeaks)-1
    index1 = indexesAtPeaks(peak);
    index2 = indexesAtPeaks(peak+1)-1;
    % Extract the curve between those indexes.
    thisX = x(index1:index2);
    thisY = y(index1:index2);
    % Skip bad (short) ones;
    if length(thisX) < 5
        break;
    end
    % Plot them.
    hold on;
    plot(thisX, thisY, 'r*');
    grid on;
    % These may not all be the same number of elements for every peak.
    % So we need to interpolate before we can average.
  originalIndexes = 1 : length(thisX);
  interpolatedIndexes = linspace(originalIndexes(1), originalIndexes(end), 100);  % Interpolate 100 points in our curves.
    interpX = interp1(originalIndexes, thisX, interpolatedIndexes, 'spline');
    interpY = interp1(originalIndexes, thisY, interpolatedIndexes, 'spline');
  % Plot the interpolated curve.
    plot(interpX, interpY, 'b-');
  % Sum them up
  sumX = sumX + interpX;
  sumY = sumY + interpY;
  count = count + 1; % # times through the loop.
end
meanX = sumX / count;
meanY = sumY / count;
% Plot the mean curve.
plot(meanX, meanY, 'g-', 'LineWidth', 3);

[joo]

could you please tell me why do you break the cycle at 'length(thisX) < 5' instead of only disconsider this correspondent iteration, and continue to run the remaining peaks (peaks=6,7,8,9,10) from the beginning? i understood everything besides this... i would be very grateful if you could answer me this last question.

[Image Analyst]

If you look at the curves I plotted, you'll see there is a little tiny peak at the bottom of the curve. This small peak is at the wrong end plus it's so small that it's basically noise. So I put in that code to ignore those small spurious peaks.

[joo]

yes, i understood that. but i don't get why after finding a 'length(thisX) < 5' true, you break (and so only plot the cycles found before); instead of not considering this 'thisX' and continue to run from the beginning the remaining peaks whose 'length(thisX) < 5' is false.

sorry for the bothering but i find the code so perfect, but i don't really understand this part... thank you so much

[Image Analyst]

Sorry - you're right. It should have been continue instead of break so that it would get the rest of the peaks instead of quitting after peak 5.

[joo]

ok. thank you very very much!

[joo]

in your code you filter x,y data responsible for the noise/small peaks. you discover this peaks by x^2 +y^2. how can you justify this? shouldn't the noise data be seen only in x and y data alone? how can i justify this filtering?

i am not at all declining your code... i just need to justify the steps, and this is a doubt of mine... can you clear me this last doubt? thank you very much.

[Image Analyst]

I'm calculating the Euclidean distance from some point way off the end of the curve's travel to the points on the curve using the Pythagorean Theorem, which is the sqrt of the sum of the squares. That distance oscillates as the point on the curve gets closer and farther away from the fixed, distant point as the curve point travels.