MATLAB: How to fit a concave baseline curve to a spectrum

continuum removalconvex curve fittinghull fittinginfrared spectroscopy

I have a number of reflectance spectra with a consistent wavelength range (350-2500nm) (x). The reflectance (y) shows an overall convex shape that I want to fit for each spectrum, similar to a baseline. Here are two examples
Spectrum 1:
Spectrum 2:
Ideally the fitted curves would look something like this:
Spectrum 1 ideal:
Spectrum 2 ideal:
I have tried the "continuum removed" function: "CR = continuum_removed(v,'Sampling',x,'Plots','yes')
The result of which returned the error: "Error using horzcat Dimensions of matrices being concatenated are not consistent.
Error in continuum_removed (line 56) x_ext = [x(1)-small x x(end)+small];
Error in Untitled3 (line 12) CR = continuum_removed(v,'Sampling',x,'Plots','yes');"
Any help in fitting a curve would be greatly appreciated.

Best Answer

I politely advocate less of Byzantine discussions and more coding of working scripts answering questions:
What Mr Longridge asked for:
to remove the 'sticks' just ctrl+R line 74 of the the working script (also attached):
clearvars
A=imread('IR1.jpg');
% figure(1);imshow(A);
A1=imadjust(A(:,:,1));A2=imadjust(A(:,:,2));A3=imadjust(A(:,:,3));     % improve image contrast
B=zeros(size(A));
B(:,:,1)=A1;B(:,:,2)=A2;B(:,:,3)=A3;
% figure(2);imshow(B);
threshold1=125;threshold2=125;                            % turn image into binary map, it misses some pixels
B(find(B>threshold1))=255; 
B(find(B<=threshold2))=0;                                    % basic image clean-up, converting to binary map
B(:,:,2)=and(B(:,:,2),B(:,:,3));
B(:,:,3)=[];
B=and(B(:,:,1),B(:,:,2));
B=logical(~B);                 
% figure(3);imshow(B);
[sz1B sz2B]=size(B);                                              % thin line to 1 pixel only
L0=zeros(sz1B,1);
y=[]
for k=1:1:sz2B 
    L=B(:,k)';
    px=find(L==1); 
    if length(px>1) 
        px=floor(median(px));
        y=[y px];                                                    % build signal
        L1=L0;
        L1(px)=1;
        C(:,k)=L1;
    end;  
end
% figure(4);imshow(C);
[y x]=ind2sub(size(C),find(C==1));
y=abs(y-sz1B);
plot(y)
%%your start point
e=[x(1) y(1)];
env=e;                                                                                   % envelope points
r1=[x (zeros(1,sz2B)+sz1B)'];                                                     % range for aim points
r2=[(ones(sz1B,1)*sz2B) [sz1B:-1:1]'];
r3=[x([end:-1:1]) zeros(1,sz2B)'];                                               % x axis points
r1(1,:)=[];r3(1,:)=[];                                                                % avoid vertical angle
range1=[r1;r2;r3];                                                                    % all the points that e_aim may take while seeking curve points
s=1;                                                                                        % s: index pointer within aim range
e_aim=range1(s,:);
k=1
while e(1)<length(x)
    alpha=atand((e_aim(2)-e(2))/(e_aim(1)-e(1)));
    seg_x=[e(1):1:e_aim(1)];                                                    % no aiming backwards, or +90º, or -90º
      if alpha>=0
          seg_y=linspace(e(2),e_aim(2),numel(seg_x)); seg_y=floor(seg_y);    end
      if alpha<0
          seg_y=linspace(e_aim(2),e(2),numel(seg_x)); seg_y=floor(seg_y);    
          seg_y=seg_y([end:-1:1]);   end;
       Lx=unique(seg_x); Lx=Lx([2:end]); Ly=seg_y([2:end]);        % reduce x and remove first point 
      x_intersect=[];y_intersect=[];
      hold on;plot([e(1) e_aim(1)],[e(2) e_aim(2)]);
      for k=2:1:length(Lx)
           y2=y(Lx);
          if (Ly(k)==y2(k))
               x_intersect=[x_intersect Lx(k)];y_intersect=[y_intersect y2(k)];
              y_intersect=unique(y_intersect);                                % take intersect with higher y
              xn=unique(find(y_intersect==max(y_intersect)));                               
              y_intersect=y_intersect(xn);
              x_intersect=x_intersect(xn);
              e(1)=x_intersect;                                                   % shift pivot point to intersect point
              e(2)=y_intersect;   
              env=[env;[x_intersect y_intersect]] ;                        % update envelope
          end
      end
      s=s+1;
      e_aim=range1(s,:);
end
hold all;plot(env(:,1)',env(:,2)','*r')                                       % show the envelope points
If you find this answer of any help solving your question,
please mark my answer as the accepted one
thanks in advance
John
jgb2012@sky.com
Related Question