function [datCorrect,weights,hfc,vesselRegion] = correctHemodynamic(data,vesselRegion,window)
   

    %% select dark region
%     vesselRegion = datMean<0.2;
    % The vessel region used bwareaopen with param 1000, but Daniel set it
    % to less. 
%     vesselRegion = bwareaopen(vesselRegion,300,4);
%     vesselRegion = bwareaopen(vesselRegion,1000,4); 
%     vesselRegion = myGrow2D(1-datMean,vesselRegion,2,150);
%     BW = boundarymask(vesselRegion);
%     figure;imshow(imoverlay(datMean,BW,'blue'),'InitialMagnification',67);
    vessel = find(vesselRegion>0);
    vesselCurve = mean(data(vessel,:),1);

    %% seperate low-frequency part and high-frequency part
    [~,hfc_New] = separate_lfc_hfc_components(vesselCurve,window);
    [~,datHfc] = separate_lfc_hfc_components(data,window);
    IterTimes = 10;
%     errors = zeros(IterTimes,1);

    %% linear regression
    % Y ~ F ` (aX1 + b)
    for i = 1:IterTimes
%         tic;
        hfc = hfc_New;
        X = [log(hfc)',ones(numel(hfc),1)];
        weights = inv(X'*X)*X' * log(datHfc');
        %% update hfc
%         hfc_New = weights(1,:)*datHfc/sum(weights(1,:));
        vesselCurve = weights(1,:)*data/sum(weights(1,:));
        [~,hfc_New] = separate_lfc_hfc_components(vesselCurve,50);
        if(max(abs(hfc_New - hfc)./hfc)<1e-3)
            break;
        end
%         toc;
    end
    datFit = exp(X*weights)';
    datCorrect = data./datFit;

end