%% Data
% Imports the array data (x and y) to the workspace by loading Excel data
% that has been imported and converted to *.mat format.
load('900day_r') % y data.
load('x_degreesb'); % x axis data
%% Remove non linear trends
opol = 0;% Degree of filtering on original profile
[p,s,mu] = polyfit(x_degreesb,x900day_r,opol);
f_y = polyval(p,x_degreesb,[],mu);
x900day_r = x900day_r - f_y;
max_x = max(x900day_r);% Find maximum in array
x900day_r = x900day_r/max_x;% Normalize height to max
min_x = min(x900day_r);% Find minimum in array
x900day_r = x900day_r - min_x;% Shift profile (lowest value in array = 0)
%% Find Peaks & Valleys
[pks, locs] = findpeaks(x900day_r); % returns peaks & locations
x900day_r_Inv = max(x900day_r)-x900day_r; % invert y data
vlys = max(vlys)-vlys; % invert data for valley markers
%% Plot Profile
% Plot profile and markers for peaks
plot(x_degreesb,x900day_r,'b',x_degreesb(locs),pks+0.04,'v','markersize',5,'markerfacecolor','b','linewidth',1);
hold on
% Plot profile and markers for valleys
plot(x_degreesb(min_locs),vlys-0.04,'^','markersize',5,'markerfacecolor','b','linewidth',1);
% Plot characteristics
axis('tight') % Makes the graph fill the figure.
grid on % Turns the grid on (major not minor).
% Sets up the figure (fig) for display
fig = gca;
set(fig,'fontname','Bodoni 72','fontsize',20);
% Set y limits to auto. Set x limits and x tick marks
ylim('auto'); % Sets the y limits
xlim([0, 360]); % Sets the x limits from 0 to 360
set (fig, 'XTick', [0, 45, 90, 135, 180, 225, 270, 315, 360]) % Sets x axis tick marks
% Set fig for presentation appearance
x = xlabel('$Location On Cylinder Perimiter {[degrees]}$'); % x label.
y = ylabel('$Curve Height {[mm]}$'); % y label.
t = title('$900 Days Cylinder$ r'); % Title (presentation fraction).
% Set vertical lines at quadrant boundaries
hold on
x1 = 90;
x2 = 180;
x3 = 270;
x4 = 360;
y1 = get(gca, 'ylim');
plot ([x1 x1],y1,'k')%Caudal Medial(0:90)
plot ([x2 x2],y1,'k')%Cranial Medial(90:180)
plot ([x3 x3],y1,'k')%Cranial Lateral(180:270)
plot ([x4 x4],y1,'k')%Cadual Lateral(270:360)
hold on
% Interpretation of text characters for presentation
set(t,'Interpreter','Latex');
set(x,'Interpreter','Latex');
set(y,'Interpreter','Latex');
%% Isolate Cranial Medial & Lateral Section of Profile
x = x_degreesb;% Quadrant data
y = x900day_r;% Profile data
indx = (x >= 90) & (x <= 270);% Index section
[pks, locs, widths, proms] = findpeaks(y(indx));% Find peaks in section
Rmax = max(pks);% Find maximum peak
Rmin = min(y(indx));% Find minimum in section
CL = (Rmax + Rmin)/2;% Center line of sectioned profile
%% Plot Center Line
hold on
x1 = 90;
x2 = 270;
y1 = CL;
plot ([x1 x2],[y1 y1],'r','linewidth',1.5);
%% Plot Rmax
hold on
x1 = 90;
x2 = 270;
y1 = Rmax;
plot ([x1 x2],[y1 y1],'k','linewidth',1.5);
%% Plot Rmin
hold on
x1 = 90;
x2 = 270;
y1 = Rmin;
plot ([x1 x2],[y1 y1],'k','linewidth', 1.5);
%% Highlight Region of Interest
%subplot(2,1,2)
x = x_degreesb;% Quadrant data
y = x900day_r;% Profile data
indx = (x >= 90) & (x <= 270);% Index section
plot(x(indx),y(indx),'k','linewidth',2)% Plot 90:270 curve in black
level = CL;% set the centerline as the level for area shading
% Shade the area of the curve in the indexed section grey [.7 .7 .7] above the level CL
area(x(indx), max(y(indx), level),level, 'EdgeColor', 'none', 'FaceColor',[.7 .7 .7],'showbaseline', 'off');
% Shade the area of the curve in the indexed section dark grey below the level CL
area(x(indx), min(y(indx), level),level, 'EdgeColor', 'none', 'FaceColor',[.5 .5 .5],'showbaseline','off');
Does anyone know how I can find the area above (light grey) and below (dark grey) the centerline (red) for the specific range between 90:270 using MATLAB? I have been trying to use trapz, setting the level (red line in Data Plot picture) but can't seem to get trapz to calculate just the highlighted areas. I posted the code, but not the data, as its a rather large set of data that makes up the curve. Any help would be greatly appreciated!
RP
@Some Guy: Thanks for the quick responses. Here is an example script that you can run. You can change the ratio of blue area to grey area by changing the level value. In pic 1 with level set to 2 they should be equal. In pic 2 blue should be more than grey with level set to 1.6. Thats what I was trying to do. Any thoughts?
%% Example
x = 0:.01:4*pi;% x data
y = sin(x)+2;% y data
level = 1.6;% level
plot(x, y)
hold on
x_interest = 0:.01:x(length(y));
y_interest = sin(x_interest)+2;
xlim ([0 x(length(y))])
% Shaded area above level
area(x_interest, max(y_interest, level), level, ...
'EdgeColor', 'none', 'FaceColor', [.6 .7 .8], ...
'ShowBaseLine', 'off');
% Shaded area below level
area(x_interest, min(y_interest, level), level, ...
'EdgeColor', 'none', 'FaceColor', [.5 .5 .5], ...
'ShowBaseLine', 'off');
y_above = max(y_interest - level,0); % Take only the part of curve above y_split
y_below = max(-y_above,0); % Take the part below y_split
A_above = trapz(y_above)
A_below = trapz(y_below)
y_aboveyou should take the max with zero only after calculating y_below. If you do that you will get correct results. - Some Guy