How can I fit a concave baseline curve to a spectrum?

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.

[John BG]

Mr Longridge

I am half way to what you are asking for.

If I give you working code, not just guidelines, will you switch your accepted answer to mine?

on this question or posting a new question, with same problem I mean.

Appreciating time and attention

awaiting answer

John Bofarull Guix

[e-mail]

Answer 1

[John BG]

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

[e-mail]

Answer 2

[John D'Errico]

Um, the red curves as you have dawn them are not in fact convex. :)

Anyway, the simple answer is...

1. Using the curve as pairs of points in the (x,y) plane, generate the convex hull. Essentially, this will be a set of line segments in the (x,y) plane, a polygon.

2. Then compute the normal vectors for each line segment in that polygon.

3. Ensure the normal vectors are pointing outwards from the convex hull.

4. Drop all line segments with a normal vector that points downwards, so has a negative y-component.

What remains will be a convex polygon that touches the curve at a limited set of points. So it will be a piecewise linear upper envelope on the curve, or a least upper bound function, that has the property of being convex.

That is indeed what you have requested, although not what you have drawn. You can form a smooth approximation to that curve using pchip. Do NOT use a traditional spline here though. In fact, even pchip might be problematic, in the sense that there MAY be some points where it does not exactly satisfy the envelope condition. (I'd need to think about that to be sure.)

Alternatively, you could use my SLM toolbox, which has the capability to fit a least upper bound function, that is also convex. Since it is indeed a spline, it will also be as smooth as possible.

[Jacob Longridge]

Hi John, Thank you for the quick response. Ahhh, convex or concave...it is a matter of perspective :).

I have utilised your SLM engine-to plot up a number of options, but can't seem to get to the "least upper bound" function. I changed the polynomial fit from 1:

to 3:

as expected, there is a better fit of the absorption, but a poor overall smooth "upper bound fit".

I then moved onto using "knots", which has a similar result (increasing knots, have a better overall fit, but not a good overarching curve:

So after little luck with this, I utilized more variables, such that the concave was downward facing using the script:

>> slm = slmengine(B(:,1),B(:,7),'plot','on','concavedown','on','decreasing','on');
Warning: Options LargeScale = 'off' and Algorithm = 'trust-region-reflective' conflict. Ignoring Algorithm and running
active-set algorithm. These settings will error in a future release. 
> In lsqlin (line 294)
  In slmengine>solve_slm_system (line 2490)
  In slmengine>slmengine_cubic (line 2332)
  In slmengine (line 255) 
>> title('concave down and decreasing')

which produced:

A good fit between ~1000-1150 and 2400-2500, but plots well below the "hull" for the rest of the spectrum Is there some way of forcing the curve to touch the upper bounds and somehow 'ignore' the absorption features?

It seems like a simple thing to solve and I am sure there must be a simple way to solve this... I have attached an example of the data-if anyone is up for the challenge...

Thanks for the help

ps. on line 205 of slm_tutorial, 'time' is spelt 'tiem' :)

[Image Analyst]

Are you going to try the convhull() function John suggested to get the convex hull? This will give the "tips" of the peaks.

[John D'Errico]

No. Look carefully at what you have drawn. Pick the second IR curve for example, or the first one.

Do you see regions in that curve where despite having an overall negative second derivative, there are regions where the curvature is a positive number as you have drawn it?

I honestly don't care if you call it convex or concave. It cannot be both. Period.

As for using SLM, Did you read the help? Apparently you did enough to find a typo on my part. But how do you think it will know that you want an upper or lower bound function unless you tell it so????? The mind reading option is not working.

From the help for slmset:

   'envelope' - allows the user to solve for an envelope of the data
            = 'off' --> the curve will be a simple least squares spline
            = 'supremum' --> comute a model such that all residuals
              (yhat - y) are positive. In effect, the curve will be a
              "supremum" (least upper bound) function.
            = 'infimum' --> comute a model such that all residuals
              (yhat - y) are negative. In effect, the curve will be a
              "infimum" (greatest lower bound) function.

            DEFAULT VALUE:  'off'

READ THE HELP, instead of running a spell checker on it.

And finally, the simple solution is as I outlined it, using a convex hull.