function [xend,yend,xmin,ymin,shapeLength,points,trace] = mtrace(xstart,ystart)
global tetamap;
global edgemap;
trace = zeros(size(edgemap));
xend=0;
yend=0;
xmin=xstart;
ymin=ystart;
shapeLength = 0;
prevLength=0;
isTraceable = true;
currentDegree = 0;
%number of pixels in current degree
pixelCount = 0;
%number of links in list representing the shape
linksCount = 1;
%list describing shape
shape=[];
arcLength=[];
points=[];
%current traced pixel;
curX = xstart;
curY = ystart;
prevX=0;
prevY=0;
%trace direction
prevBlockX=1;
prevBlockY=1;
%trace ahead amount
scanX=0;
scanY=0;
failCount=0;
avgDeg = 0;
%first 2 int point to the next direction of supposed edge
%pixel(clockwise search)
%second 2 ints point to the direction we start looking for an edge at in case we
%did find an edge pixel in that direction
directions = {[1,2,2,1],[1,3,2,1],[2,3,2,1];[1,1,3,2],[2,2],[3,3,2,3];[2,1,3,2],[3,1,2,3],[3,2,2,3]};
clockWise = {[1,2,2,1],[1,3,2,1],[2,3,2,1];[1,1,3,2],[2,2],[3,3,2,3];[2,1,3,2],[3,1,2,3],[3,2,2,3]};
cocountunterClock = {[2,1,2,1],[1,1,2,1],[1,2,2,1];[3,1,3,2],[2,2],[1,3,2,3];[3,2,3,2],[3,3,2,3],[2,3,2,3]};
while(isTraceable)
    trace(curY,curX) = 1;
    %prevent backtracing
    %edgemap(curY,curX)=0;
    
    
    %create adjmatrix for current pixel     
    [curY,curX];
     adjBlock = edgemap(curY-1:curY+1,curX-1:curX+1);
     blockCurY = prevBlockY;
     blockCurX = prevBlockX;
     foundPixel = false;
     adjCount = 1;
      
        
        while( ~foundPixel && adjCount<=8 ) 
            
            if(adjBlock(blockCurY,blockCurX)==1)  
                foundPixel = true;
                
                %set the new current cordinate according to the place of
                %the edge that was found
                curX = curX+(blockCurX-2);
                curY = curY+(blockCurY-2);
                xend = max([xend,curX]);
                yend = max([yend,curY]);
                xmin = min([curX,xmin]);
                ymin = min([curY,ymin]);
                
                pixelDeg = abs(tetamap(curY,curX)*180/pi);
                
                %pixelCount=pixelCount+1;
                shapeLength = shapeLength+1;
                points=[points;shapeLength,currentDegree];
                if(blockCurY-prevBlockY<0)
                    %entry from below
                    prevBlockY = 3;
                    prevBlockX=2;
                elseif(blockCurY-prevBlockY>0)
                    %entry from above
                    prevBlockY=1;
                    prevBlockX=2;
                elseif(blockCurX-prevBlockX>0)
                    %entry from right
                    prevBlockY=2;
                    prevBlockX=1;
                else%entry from left
                    prevBlockY=2;
                    prevBlockX=3;
                end
                    
                    
            else
                %save last entry position                    
                prevBlockY = blockCurY;
                prevBlockX = blockCurX;
                %set pointer to  starting pixel in adjacency block
                %adjBlock(blockCurY,blockCurX)=-1
                newdir = directions{blockCurY,blockCurX};
                 blockCurY=newdir(1);
                blockCurX=newdir(2);
                adjCount = adjCount+1;

            end
        end
        
        %check new degree if a change along 10pixels average to more than X deg , add to our
        %shape contour ignore pixels that suddenly opposite signed degree        
        if(foundPixel)            
                 avgDeg = avgDeg+pixelDeg;
                pixelCount=pixelCount+1;
                
                %fprintf('test\n');
                %fprintf('length: %d  x: %d y:%d  d: %f curDeg: %f, avgDeg: %f\n',shapeLength,curX,curY,pixelDeg,currentDegree,avgDeg/pixelCount);
            %shape{1,linksCount}=[currentDegree,pixelDeg];            
            if(pixelCount>=10 && abs(currentDegree-(avgDeg/pixelCount))>=15)
                d = sqrt(abs(prevX-curX)^2+abs(prevY-curY)^2);
                %if features are close merge them into one.
                if(d<5)
                       %trace(prevY-2:prevY+2,prevX-2:prevX+2)=0;
                        trace(prevY,prevX)=1;
                    
                    avgDeg = avgDeg/pixelCount;
                    diff = abs(currentDegree-avgDeg);
                    median = (shape(size(shape,1))+diff)/2;
                    shape(size(shape,1)) = median;
                    
                    %fprintf('avarging:  length: %d prevLength: %d arclength: %d\n',(shapeLength-prevLength),prevLength,arcLength(size(arcLength)));
                    %arcLength(size(arcLength,1))=arcLength(size(arcLength,1))+shapeLength-prevLength;
                    arcLength(size(arcLength,1))=shapeLength;
                    %shape(size(shape)) = shape(size(shape),1)+shapeLength-prevLength;
                    prevLength=shapeLength;
                    currentDegree= avgDeg;
                else   
                    
                    avgDeg = avgDeg/pixelCount;
                    %fprintf('adding: diff:%f absdiff:%f length: %d prevLength: %d\n',currentDegree-avgDeg,abs(currentDegree-avgDeg),shapeLength-prevLength,prevLength);
                    shape=[shape;abs(currentDegree-avgDeg)];
                    %arcLength=[arcLength;shapeLength-prevLength];
                    arcLength=[arcLength;shapeLength];
                    prevLength=shapeLength;
                    linksCount = linksCount+1;
                    currentDegree=avgDeg;
                end
                trace(curY-2:curY+2,curX-2:curX+2)=1;
                prevY=curY;
                prevX=curX;
                pixelCount=0;                
                avgDeg=0;
            end
            if(pixelCount>25)
                avgDeg=avgDeg/pixelCount;
                pixelCount=1;
                
                %avgDeg=pixelDeg;
            end
            %else
              %  avgDeg = avgDeg+pixelDeg;
               % pixelCount=pixelCount+1;
            %currentDegree = abs(currentDegree-pixelDeg);            
            %end
            
        end
        
        %if no neighbour found
        if(~foundPixel)
            isTraceable = false;
        end
        
          %if reached or passed starting point
         if(curX<=xstart && curY<=ystart)
             isTraceable = false;
         end
end

if(shapeLength>80)
   %imview(trace);
end

%update edge map to prevent retracing
edgemap = edgemap-trace;
%imview(edgemap);
[elem,index] = min(shape);
shape = circshift(shape,[size(shape,1)-index+1,0]);

