function [data,fitCurves,pixClusters] = preprocessDebleach(dataOrg,datMean)
    [H,W,T] = size(dataOrg);
    dataVec = reshape(dataOrg,[],T);

    %% Supper Pixel
    [L,N] = superpixels(datMean,1000,'Compactness',1,'NumIterations', 20);
    BW = boundarymask(L);
    pixClusters = label2idx(L);
%     figure;imshow(imoverlay(datMean,BW,'cyan'),'InitialMagnification',67);

    %% Model
    f = @(b,x) (b(1)*exp(-b(2).*x.^(b(4)))+b(3));
    options = optimoptions('lsqcurvefit');
    options.Display = 'off';

    %% Fit
    datFit = zeros(size(dataVec));  % the bleach trend

%     tic;
    fitCurves = cell(N,1);
    curCurves = cell(N,1);
    x = [1:T]';

    for i = 1:N
        pix = pixClusters{i};
        curCurves{i} = mean(dataVec(pixClusters{i},:))';
    end

    for i = 1:N
        z = double(curCurves{i});
        x0 = [max(z)-min(z),1e-3,min(z)/(max(z)-min(z)),1];
        lb = [0,1e-5,0,1];
        ub = [max(z),1,1000,5];
        fitted = lsqcurvefit(f,x0,x,z,lb,ub,options);
        fitcurve = f(fitted,x);
        fitcurve = fitcurve/fitcurve(1);
        fitCurves{i} = fitcurve;
    end

    for i = 1:N
        pix = pixClusters{i};
        datFit(pix,:) = repmat(fitCurves{i}',numel(pix),1);
    end
    data = dataVec./datFit;
%     toc;

end