diff --git a/inteligent_scissors/fLiveWireCalcP.m b/inteligent_scissors/fLiveWireCalcP.m
new file mode 100644
index 0000000..1297b65
--- /dev/null
+++ b/inteligent_scissors/fLiveWireCalcP.m
@@ -0,0 +1,151 @@
+function [iPX, iPY] = fLiveWireCalcP(dImg, iXS, iYS, dRadius)
+%FLIVEWIRECALCP Calculates the path maps in a live-wire implementation .
+
+iLISTMAXLENGTH = 10000;
+
+if nargin < 4, dRadius = 10000; end
+
+[iNRows, iNCols] = size(dImg);
+iNPixelsToProcess = min([int32(pi.*dRadius.^2), iNRows.*iNCols]);
+iNPixelsProcessed = int32(0);
+
+% Initialize live-wire data structures
+iLI = zeros(iLISTMAXLENGTH, 3, 'int32');   % The active list coordinates (x, y, lin)
+dLG = zeros(iLISTMAXLENGTH, 1);            % The active list weights
+iLInd = uint16(0);               % The active list index: Points to the last element of iLI and dLG
+
+lE = false(size(dImg));          % Binary image indicating whether a pixel has peen processed or not
+iPX = zeros(size(dImg), 'int8');% Vector field following the minimum cost path from each processed pixel to the seed
+iPY = zeros(size(dImg), 'int8');% Vector field following the minimum cost path from each processed pixel to the seed
+
+% Build neighbourhoor index lookup tables for faster loop execution
+[iXX, iYY] = meshgrid(1:iNCols, 1:iNRows);
+iXX = int16(iXX); iYY = int16(iYY);
+iNX = cat(3, iXX - 1, iXX    , iXX + 1, iXX - 1, iXX + 1, iXX - 1, iXX    , iXX + 1);
+iNY = cat(3, iYY - 1, iYY - 1, iYY - 1, iYY    , iYY    , iYY + 1, iYY + 1, iYY + 1);
+iNX = permute(iNX, [3 2 1]);
+iNY = permute(iNY, [3 2 1]);
+
+iDirX = int16([1; 0; -1; 1; -1;  1;  0; -1]);
+iDirY = int16([1; 1;  1; 0;  0; -1; -1; -1]);
+%dWeight = [1, 0.7, 1, 0.7, 0.7, 1, 0.7, 1];
+
+%to calculate the dImgP and dImgQ(8 neighbors)
+dImgP=dImg;
+dImgNQ = cat(3,dImg,dImg,dImg,dImg,dImg,dImg,dImg,dImg);
+%dImgNQ = permute(dImgNQ, [3 1 2]);
+id = repmat({':'}, ndims(dImg), 1);
+n = size(dImg, 1);
+m = size(dImg, 2);
+for i = 1:8
+        if iDirX(i) > 0
+            id{1} = [2:n 1];
+        elseif iDirX(i) < 0
+            id{1} = [n 1:n - 1];
+        else 
+            id{1} =  [1:n];  
+        end
+
+        if iDirY(i) > 0
+            id{2} = [2:n 1];
+        elseif iDirY(i) < 0
+            id{2} = [n 1:n - 1];
+        else 
+            id{2} = [1:n];  
+        end    
+        dImgNQ(:,:,i)= dImg(id{:});
+end
+dF = cat(3,dImg,dImg,dImg,dImg,dImg,dImg,dImg,dImg);
+%dF = permute(dF, [3 1 2]);
+for i = 1:8
+    dF(:,:,i) = fLiveWireGetCostFcn(dImgP,dImgNQ(:,:,i));
+end
+% Initialize active list with zero cost seed pixel.
+iLInd = iLInd + 1;
+iLI(iLInd, :) = [iXS, iYS, (iXS - 1).*iNRows + iYS];
+dLG(iLInd) = 0;
+
+    
+% While there are still objects in the active list and pixel limit not
+% reached
+while (iLInd > 0) && (iNPixelsProcessed < iNPixelsToProcess)
+    % Determine pixel q in list with minimal cost and remove from active
+    % list. Mark q as processed.
+    [temp, iInd] = min(dLG(1:iLInd));
+    iQI = iLI(iInd, :);
+    dQG = dLG(iInd);
+    iLI(iInd, :) = iLI(iLInd, :);
+    dLG(iInd, :) = dLG(iLInd, :);
+    iLInd = iLInd - 1;
+    
+    lE(iQI(2), iQI(1)) = true;
+
+    % Generate neighbourhood of q
+    iN = zeros(8, 4, 'int16');
+    iN(:, 1) = iNX(:, iQI(1), iQI(2));
+    iN(:, 2) = iNY(:, iQI(1), iQI(2));
+    iN(:, 3) = iDirX;
+    iN(:, 4) = iDirY;
+    
+    % Mask out pixels that are out of bounds (q is a border pixel)
+    lValid = (iN(:, 1) > 0) & (iN(:, 1) <= iNCols) & ...
+             (iN(:, 2) > 0) & (iN(:, 2) <= iNRows);
+
+    iN = iN(lValid, :);
+    %dThisWeight = dWeight(lValid);
+
+    % Loop over the neighbourhood of q
+    for iI = 1:size(iN, 1);
+        if lE(iN(iI, 2), iN(iI, 1)), continue, end % Skip if neighbourhood pixel already processed
+        iR = iN(iI,:);
+    
+        iLinInd = (int32(iR(2)) - 1).*iNRows + int32(iR(1));
+        % Compute the new accumulated cost to the neighbourhood pixel
+ 
+        %dImgP=dImg;
+        %id = repmat({':'}, ndims(dImg), 1);
+        %n = size(dImg, 1);
+       % if iR(3) > 0
+       %     id{1} = [2:n 1];
+       % elseif iR(3) < 0
+       %     id{1} = [n 1:n - 1];
+       % else 
+       %     id{1} =  [1:n];  
+       % end
+       % n = size(dImg, 2);
+       % if iR(4) > 0
+       %     id{2} = [2:n 1];
+       % elseif iR(4) < 0
+        %    id{2} = [n 1:n - 1];
+       % else 
+       %     id{2} = [1:n];  
+      %  end
+        
+       % dImgQ= dImg(id{:});
+       % dF = fLiveWireGetCostFcn(dImgP,dImgQ);%此处循环次数太多，时间复杂度太高，换成只计算一遍
+        dThisG = dQG + dF(iR(2), iR(1),iI);%.*dThisWeight(iI);
+        
+        % Check whether r is already in active list and if the
+        % current cost is lower than the previous
+        iInd = find(iLI(1:iLInd, 3) == iLinInd);
+        if ~isempty(iInd)
+            if dThisG < dLG(iInd)
+                dLG(iInd) = dThisG;
+                iPX(iR(2), iR(1)) = iR(3);
+                iPY(iR(2), iR(1)) = iR(4);
+            end
+        else
+        % If r is not in the active list, add it!
+            iLInd = iLInd + 1;
+            iLI(iLInd, :) = [int32(iR(1)), int32(iR(2)), iLinInd];
+            dLG(iLInd) = dThisG;
+
+            iPX(iR(2), iR(1)) = iR(3);
+            iPY(iR(2), iR(1)) = iR(4);
+        end 
+        
+    end % FOR loop over the neighborhood of q
+    
+    iNPixelsProcessed = iNPixelsProcessed + 1;
+  
+end % WHILE
diff --git a/inteligent_scissors/fLiveWireGetCostFcn.m b/inteligent_scissors/fLiveWireGetCostFcn.m
new file mode 100644
index 0000000..d7469c5
--- /dev/null
+++ b/inteligent_scissors/fLiveWireGetCostFcn.m
@@ -0,0 +1,31 @@
+function dF = fLiveWireGetCostFcn(dImgP,dImgQ)
+%calculate the local cost of p and q
+
+%if nargin < 2,
+    dWz = 0.43;
+    dWg = 0.43;
+    dWd = 0.14;
+%end
+
+% Calculat the cost function
+% The gradient strength cost Fg
+%dImg = double(dImg);
+%[dY, dX] = gradient(dImg);
+%dFg = sqrt(dX.^2 + dY.^2);
+%dFg = 1 - dFg./max(dFg(:));
+dImgQ = double(dImgQ);
+[GmagQ, GdirQ] = imgradient(dImgQ);%it's to calculate the gradient and direction of image
+dFg = 1 - GmagQ./max(GmagQ(:));
+
+% The zero-crossing cost Fz
+lFz = ~edge(dImgQ, 'zerocross');%zero crossing edge detector, 0 stands for "moving fast"
+
+% The gradient direction Fd ??
+dImgP = double(dImgP);
+[GmagP, GdirP] = imgradient(dImgP);
+GdirQ=abs(GdirQ);
+GdirP=abs(GdirP);
+dFd=(abs(GdirQ-pi)+abs(GdirP-pi))./pi;
+
+% The Sum:
+dF = dWz.*double(lFz)+ dWg.*dFg +dWd.*dFd;
diff --git a/inteligent_scissors/fLiveWireGetPath.m b/inteligent_scissors/fLiveWireGetPath.m
new file mode 100644
index 0000000..1527035
--- /dev/null
+++ b/inteligent_scissors/fLiveWireGetPath.m
@@ -0,0 +1,32 @@
+function [iX, iY] = fLiveWireGetPath(iPX, iPY, iXS, iYS)
+%FLIVEWIREGETPATH Traces the cheapest path from (IXS, IYS)^T through the
+%pathmaps IPX, IPY back to the seed .
+
+iMAXPATH = 1000;
+
+% Initialize the variables
+iPX  = int16(iPX);
+iPY  = int16(iPY);
+iXS = int16(iXS);
+iYS = int16(iYS);
+
+iX = zeros(iMAXPATH, 1, 'int16');
+iY = zeros(iMAXPATH, 1, 'int16');
+
+iLength = 1;
+iX(iLength) = iXS;
+iY(iLength) = iYS;
+
+% While not at the seed point: march back in the direction indicated by the
+% path maps iPX (x-direction) and iPY (y-direction).
+while (iPX(iYS, iXS) ~= 0) || (iPY(iYS, iXS) ~= 0) % We're not at the seed
+    iXS = iXS + iPX(iYS, iXS);
+    iYS = iYS + iPY(iYS, iXS);
+    iLength = iLength + 1;
+    iX(iLength) = iXS;
+    iY(iLength) = iYS;
+end
+
+% revert vectors (to make it a forward path) and don't return the seed point.
+iX = iX(iLength - 1:-1:1);
+iY = iY(iLength - 1:-1:1);
diff --git a/inteligent_scissors/livewire.m b/inteligent_scissors/livewire.m
new file mode 100644
index 0000000..15043db
--- /dev/null
+++ b/inteligent_scissors/livewire.m
@@ -0,0 +1,263 @@
+function varargout = livewire(I, dRadius, dCaptureRadius)
+
+dPointerCross = NaN(16);
+dPointerCross(8   ,1:15) = 2;
+dPointerCross(1:15,   8) = 2;
+
+dPointerCrossMag = ...
+   [NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN; ... 
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN; ... 
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN; ...
+      2,   2,   2,   2,   2,   2,   2,   2,   2,   2,   2,   2,   2,   2,   2,   2; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN,   2, NaN, NaN, NaN,   2, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN,   2, NaN, NaN, NaN,   2, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN,   2, NaN, NaN, NaN,   2, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN,   2, NaN, NaN, NaN,   2, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN,   2, NaN, NaN, NaN,   2,   2,   2, NaN, NaN; ...
+    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN];
+
+% Check if number of output arguments is acceptable.
+if nargout == 2 || nargout > 3
+    error('LIVEWIRE expects either 0, 2 or 3 output arguments.'); 
+end
+
+% Check input argument and try to get a figure/axes/image handle if no
+% input image is supplied (operate on current axes).
+if ~nargin
+    hF = get(0, 'CurrentFigure');
+    if isempty(hF), error('Found no figure. Call LIVEWIRE with an image argument!'); end
+    
+    hA = get(hF, 'CurrentAxes');
+    if isempty(hA), error('Found no axes. Call LIVEWIRE with an image argument!'); end
+    
+    hI = findobj(hA, 'Type', 'image');
+    if isempty(hF), error('Found no image object. Call LIVEWIRE with an image argument!'); end
+    
+    dImg = double(get(hI, 'CData'));
+else
+    if size(I, 3) == 3, I = rgb2gray(I); end
+    dImg = double(I);
+    hI = imshow(dImg, []);
+    hA = gca;
+    hF = gcf;
+end
+
+% Define optional input parameters if omitted
+if nargin < 2
+    SOptions.dRadius = 200;
+    SOptions.dCaptureRadius = 4;
+end
+
+% Bring figure to the front and set callbacks
+figure(hF);
+try
+    set(hF, ...
+        'WindowButtonDownFcn'   , @fButtonDownFcn,...
+        'KeyPressFcn'           , @fKeyPressFcn, ...
+        'KeyReleaseFcn'         , @fKeyReleaseFcn, ...
+        'Pointer'               , 'custom', ...
+        'PointerShapeCData'     , dPointerCrossMag, ...
+        'PointerShapeHotSpot'   , [8, 8], ...
+        'DoubleBuffer'          , 'on'); %'WindowButtonMotionFcn' , @fMotionFcn, ...
+catch
+    set(hF, ...
+        'WindowButtonDownFcn'   , @fButtonDownFcn,...
+        'KeyPressFcn'           , @fKeyPressFcn, ...
+        'Pointer'               , 'custom', ...
+        'PointerShapeCData'     , dPointerCrossMag, ...
+        'PointerShapeHotSpot'   , [8, 8], ...
+        'DoubleBuffer'          , 'on'); %'WindowButtonMotionFcn' , @fMotionFcn, ...
+end
+
+% Create the line objects (use two lines for better visibility)
+hLine1 = line(...
+    'Parent'    , hA, ...
+    'XData'     , [], ...
+    'YData'     , [], ...
+    'Clipping'  , 'off', ...
+    'Color'     , 'g', ...
+    'LineStyle' , '-', ...
+    'LineWidth' , 1.5);
+
+hLine2 = line(...
+    'Parent'    , hA, ...
+    'XData'     , [], ...
+    'YData'     , [], ...
+    'Clipping'  , 'off', ...
+    'Color'     , 'r', ...
+    'LineStyle' , ':', ...
+    'LineWidth' , 1.5);
+
+% Initialize the global variables
+%dF          = fLiveWireGetCostFcn(dImg); % The cost function of the live-wire algorithm, see Ref.
+iPX         = zeros(size(dImg), 'int8'); % The path map that shows the cheapest path to the sast anchor point.
+iPY         = zeros(size(dImg), 'int8'); % The path map that shows the cheapest path to the sast anchor point.
+dXData      = [];                        % The x-coordinates of the path
+dYData      = [];                        % The y-coordinates of the path
+iAnchorList = zeros(200, 1);             % A list of the anchor point indices in dXData and dYData for undo operations
+iAnchorInd  = 0;                         % The index of the list
+lRegularEnd = false;                     % Indicates whether path was successfully drawn or the UI was aborted.
+lControl    = false;                     % Indicates whether the control key is pressed
+
+uiwait(hF);
+
+if ~lRegularEnd
+    switch nargout
+
+        case 1
+            varargout{1} = false(size(dImg));
+
+        case 3
+            varargout{1} = false(size(dImg));
+            varargout{2} = [];
+            varargout{3} = [];
+
+        otherwise
+
+    end
+    return
+end
+
+% Assign the return values
+lMask = poly2mask(dXData, dYData, size(dImg, 1), size(dImg, 2));
+switch nargout
+    
+    case 0
+        figure, imshow(lMask)
+
+    case 1
+        varargout{1} = lMask;
+
+    case 3
+        varargout{1} = lMask;
+        varargout{2} = dXData;
+        varargout{3} = dYData;
+
+    otherwise
+end
+
+
+    % Executed when the mouse button is pressed. Used to determine a
+    % new anchor point or to end the user interaction.
+    function fButtonDownFcn(hObject, eventdata) %#ok<*INUSD> eventdata is repeatedly unused
+        [dX, dY] = fGetAxesPos;
+        if ~dX, return, end
+        
+        if isempty(dXData)
+            % The starting point of the path is selected
+            dXData = dX;
+            dYData = dY;
+        else
+            % A new anchor point and the cheapest path to the last anchor
+            % point is appended to the path
+            [iXPath, iYPath] = fLiveWireGetPath(iPX, iPY, dX, dY);
+            if isempty(iXPath)
+                iXPath = dX;
+                iYPath = dY;
+            end
+            dXData = [dXData, double(iXPath(:)')];
+            dYData = [dYData, double(iYPath(:)')];
+        end
+
+        iAnchorInd = iAnchorInd + 1;
+        iAnchorList(iAnchorInd) = length(dXData); % Save the previous path length for the undo operation
+        
+        % Update the UI and calculate the new path map iP
+        set([hLine1, hLine2], 'XData', dXData, 'YData', dYData);
+        drawnow expose
+        [iPX, iPY] = fLiveWireCalcP(dImg, dX, dY, SOptions.dRadius);
+
+        % If right-click, double-click or shift-click occurred, close path
+        % and return.
+        if ~(strcmp(get(hF, 'SelectionType'), 'normal')) && ~(lControl)
+            [iXPath, iYPath] = fLiveWireGetPath(iPX, iPY, dXData(1), dYData(1));
+            if isempty(iXPath)
+                iXPath = dXData(1);
+                iYPath = dYData(1);
+            end
+            dXData = [dXData, double(iXPath(:)')];
+            dYData = [dYData, double(iYPath(:)')];
+            set([hLine1, hLine2], 'XData', dXData, 'YData', dYData);
+            drawnow expose
+            set(hF, 'WindowButtonMotionFcn', '', 'WindowButtonDownFcn', '', 'KeyPressFcn', '');
+            lRegularEnd = true;
+            uiresume(hF);
+        end
+
+    end
+
+    function [dX, dY] = fGetAxesPos()
+        dPos  = get(hA, 'CurrentPoint');
+        dXLim = get(hA, 'XLim');
+        dYLim = get(hA, 'YLim');
+        dX = dPos(1, 1);
+        dY = dPos(1, 2);
+        if (dX < dXLim(1)) || (dX > dXLim(2)) || ...
+           (dY < dYLim(1)) || (dY > dYLim(2))
+            dX = 0;
+            dY = 0;
+        end
+    end
+
+    function fKeyPressFcn(hObject, eventdata)
+        
+        if strcmp(eventdata.Key, 'control')
+            set(hF, 'PointerShapeCData', dPointerCross);
+            lControl = true;
+        end
+       
+        switch eventdata.Key
+            % ENTER: close the path and return
+            case 'return'
+                [iXPath, iYPath] = fLiveWireGetPath(iPX, iPY, dXData(1), dYData(1));
+                if isempty(iXPath)
+                    iXPath = dXData(1);
+                    iYPath = dYData(1);
+                end
+                dXData = [dXData, double(iXPath(:)')];
+                dYData = [dYData, double(iYPath(:)')];
+                set([hLine1, hLine2], 'XData', dXData, 'YData', dYData);
+                drawnow expose
+                set(hF, 'WindowButtonMotionFcn', '', 'WindowButtonDownFcn', '', 'KeyPressFcn', '');
+                lRegularEnd = true;
+                uiresume(hF);
+
+            % DEL or BACKSPACE: delete path to last anchor
+            case {'delete', 'backspace'}
+                iAnchorInd = iAnchorInd - 1;
+                if ~iAnchorInd
+                    return
+                end
+                
+                dXData = dXData(1:iAnchorList(iAnchorInd));
+                dYData = dYData(1:iAnchorList(iAnchorInd));
+                
+                set([hLine1, hLine2], 'XData', dXData, 'YData', dYData);
+                drawnow expose
+                [iPX, iPY] = fLiveWireCalcP(dImg, dXData(end), dYData(end), SOptions.dRadius);
+                fMotionFcn(hObject, []);
+
+            otherwise
+
+        end
+    end
+
+    function fKeyReleaseFcn(hObject, eventdata)
+        if strcmp(eventdata.Key, 'control')
+            lControl = false;
+            set(hF, 'PointerShapeCData', dPointerCrossMag);
+        end
+    end
+
+    % Closes the figure if requested.
+    function fCloseGUI(hObject, eventdata) %#ok<DEFNU> <-stupid!
+        delete(hObject);
+    end
+
+end
diff --git a/inteligent_scissors/main.m b/inteligent_scissors/main.m
new file mode 100644
index 0000000..7ac6283
--- /dev/null
+++ b/inteligent_scissors/main.m
@@ -0,0 +1,14 @@
+clear;
+close all;
+I = imread('a.png');
+if size(I, 3) == 3, I = rgb2gray(I); end
+dImg = double(I);
+%dF  = fLiveWireGetCostFcn(dImg,dImg-1); % The cost function of the live-wire algorithm, see Ref
+iPX = zeros(size(dImg), 'int8'); % The path map that shows the cheapest path to the sast anchor point.
+iPY = zeros(size(dImg), 'int8'); % The path map that shows the cheapest path to the sast anchor point.
+dX = 50; dY = 50;
+iXPath = dX; iYPath = dY;
+[iPX, iPY] = fLiveWireCalcP(dImg, dX, dY);
+dX = 51; dY = 55;
+[iXPath, iYPath] = fLiveWireGetPath(iPX, iPY, dX, dY);
+[iPX, iPY] = fLiveWireCalcP(dImg, dX, dY);