%include "irg:imgraph.inc"
%include "utils:region"
%begin

%bytearray image,copy(0:128*512)

%integerarray  lines(0:500,1:4)
%integer numlines,ref,size,w,th,cth,i

%byte col

%string (255) file

!*****************************************************!
!         DISPLAY FRAME                               !
!*****************************************************!
%routine display frame(%integer x,y,%bytename line)
%label ok

  *move.l d0,d2
  *lsl.l #7,d1
  *andi.l #16_FFFFFFF0,d0
  *lsr.l #3,d0
  *add.l d0,d1
  *lea   frame,a1
  *lea   0(a1,d1.l),a1
  *move.l a1,a3
  *adda.l #128,a3
  *lea   colourreg,a2
  *andi.l #15,d2
  *eor.l  #15,d2
  *clr.l  d1
  *bset   d2,d1
  *moveq  #127,d3
z:*move.l  #511,d0
a:*move.b (a0)+,(a2)
  *move.w d1,(a1)
  *ror.w #1,d1
  *btst  #15,d1
  *beq   ok
  *lea    2(a1),a1
ok:
  *dbra   d0,a
  *move.l a3,a1
  *adda.l #128,a3
  *bset   d2,d1
  *dbra   d3,z
%end

!*****************************************************!
!         FILLIMAGE                                   !
!*****************************************************!
%routine fillimage(%bytearrayname image(0:128*512))
%integer i,j,col,row,k,np
i=0;j=0
%for col=1,1,511 %cycle
  %for row=2+i,2,124+i %cycle
    k=row<<9+col
    np=image(k-512)+image(k+512)
    %if j=0 %then np=np+image(k-1) %else np=np+image(k+1)
    %if np=2 %or np=3 %then image(k)=1 %else image(k)=0
  %repeat
    %if j=0 %then j=1 %elsestart
      %if i=1 %then i=0 %else i=1
      j=0
    %finish
%repeat
%end

%routine thinB(%bytearrayname image(0:128*512))
%integer row,col,s,i,j,k,n,past

%constantintegerarray mask(0:15,0:1)=
2_00001011,2_11011011,
2_11010000,2_11011011,
2_01101000,2_01111110,
2_00010110,2_01111110,
2_01001001,2_11111111,
2_00010011,2_11111111,
2_01110000,2_11111111,
2_01010100,2_11111111,
2_00101010,2_11111111,
2_10010010,2_11111111,
2_00001110,2_11111111,
2_11001000,2_11111111,
2_11010110,2_11111111,
2_11111000,2_11111111,
2_01101011,2_11111111,
2_00011111,2_11111111

%routine thinpointB(%bytearrayname image(0:128*512),%integer i,j,   
%integername n)
%integer test,res,a,b,c,d,e,f,g,h
a=i-1
b=i+1
c=j<<9
d=(j+1)<<9
e=(j-1)<<9
  test=(((((((image(a+d))<<1+image(i+d))<<1+image(b+d))<<1+image(a+c))<<1+ %c
               image(b+c))<<1+image(a+e))<<1+image(i+e))<<1+image(b+e)
n=0
%for a=0,1,15%cycle
  res=(test!!mask(a,0))&mask(a,1)
  %if res=0 %start
    n=1
    image(i+512*j)=0
  %finish
%repeat
%end

fillimage(image)

%cycle
  past=0;n=0
    %for row=1,1,126 %cycle
      %for col=1,1,510 %cycle
        k=row<<9+col
        s=image(k)
        %if s=1 %start
         thinpointB(image,col,row,n)
        %finish
        %if n=1 %then past=1
      %repeat
  %repeat
%repeatuntil past=0
%end

!*****************************************************!
!         WORM                                        !
!*****************************************************!
%routine worm(%bytearrayname image(0:128*512),
              %integerarrayname lines(0:500,1:4),
              %integername totnl,%integer w,th,cth)

%integer m,i,j,k,l,x,y,tot,p,
         numpoints,numlines

%integerarray list(0:100,1:2)

%recordformat pointf(%integer x,y)
%recordformat listofpoints(%integer num,%record(pointf) %array pt(0:1024))

%record(listofpoints) %array ptlist(0:128)

%routinespec longest(%integerarray points(0:7,1:2),%integer num,
                     %integername res,len)

%routine findpoints(%integer x,y,%integername num,
%integerarrayname points(0:7,1:2))
%integer i,j
num=0
%for i=1,-1,-1 %cycle
  %for j=1,-1,-1 %cycle
     %if image(x+i+(y+j)<<9)=1 %start
printstring("Findpoints : ");write(x+i,3);write(y+j,3);newline
        points(num,1)=x+i
        points(num,2)=y+j
        num=num+1
     %finish
  %repeat
%repeat
%end

%integerfunction length(%integer x,y)
%integer num,res=0,len=0,tlen=0
%integerarray points (0:7,1:2)
printstring("Length : position = ");write(x,3);write(y,3);newline
  plot(x,y)
  image(x+y<<9)=0
  findpoints(x,y,num,points)
printstring("Length : found points = ");write(num,3);newline
  colour(green)
  %if num=1 %then len=1+length(points(0,1),points(0,2)) %c
  %else %if num>1 %start
    longest(points,num,res,tlen)
    len=len+tlen
printstring("Length : longest line = ");write(tlen,3);newline
printstring("Length : total line = ");write(len,3);newline
  %finish
  image(x+y<<9)=1
  %result=len
%end

%routine longest(%integerarray points(0:7,1:2),%integer num,
                 %integername res,len)
%integer n,i
  res=0;len=0
  %for i=0,1,num-1 %cycle
    colour(white)
    n=length(points(i,1),points(i,2))
printstring("Longest : length= ");write(n,3);newline
    %if n>len %start
      %if i#0 %then image(points(res,1)+points(res,2)<<9)=0
      len=n 
      res=i
printstring("Longest : new result= ");write(i,3);newline
    %finish
  %repeat
%end

%routine push(%integer x,y, %record(listofpoints) %name ptlist)
%integer np,res=0,len=0
%integerarray list(0:7,1:2)
  colour(red)
  plot(x,y)
  image(x+y<<9)=0
printstring("push: ");write(x,3);write(y,3);newline
  ptlist_pt(ptlist_num)_x=x
  ptlist_pt(ptlist_num)_y=y
  ptlist_num=ptlist_num+1
  findpoints(x,y,np,list)
printstring("push: foundpoints: ");write(np,3);newline
  %if np>1 %then longest(list,np,res,len)
  %if np#0 %then push(list(res,1),list(res,2),ptlist)
%end

%routine find lines(%record(listofpoints) %arrayname ptlist(0:128),
%integername numlines)
%integer i,j,col,row,k,s,np,numpoints,res,len
%integerarray list(0:7,1:2)
  openoutput(1,"pub:worm.ep")
  selectoutput(1)
  numlines=0
  %for col=1,1,511 %cycle
    %for row=1,1,126 %cycle
      k=row<<9+col
      s=image(k)
      %if s=1 %start
printstring("Initial point found at ");write(col,3);write(row,3);newline
        image(k)=0
        findpoints(col,row,np,list)                        
        %if np#0 %start
          res=0
          ptlist(numlines)_num=1
          ptlist(numlines)_pt(0)_x=col
          ptlist(numlines)_pt(0)_y=row
          %if np>1 %then longest(list,np,res,len)
          push(list(res,1),list(res,2),ptlist(numlines))
          numlines=numlines+1
        %finish
      %finish
    %repeat
  %repeat
  closeoutput
  selectoutput(0)
%end

%integerfunction dx(%integer n,l,w)
  %result=ptlist(l)_pt(n)_x-2*ptlist(l)_pt(n-w)_x+ptlist(l)_pt(n-2*w)_x
%end

%integerfunction dy(%integer n,l,w)
  %result=ptlist(l)_pt(n)_y-2*ptlist(l)_pt(n-w)_y+ptlist(l)_pt(n-2*w)_y
%end

%integerfunction a(%integer n,l,th,w)
  %result=|dx(n,l,w)|+|dy(n,l,w)|-th
%end

%integerfunction tx(%integer n,l,w,p)
  %integer res
  %if p>n-1 %then res=dx(n,l,w) %else res=tx(n-1,l,w,p)+dx(n,l,w)
  %result=res
%end

%integerfunction ty(%integer n,l,w,p)
  %integer res
  %if p>n-1 %then res=dy(n,l,w) %else res=ty(n-1,l,w,p)+dy(n,l,w)
  %result=res
%end

%integerfunction c(%integer n,w,p,cth)
  %result=|tx(n,l,w,p)|+|ty(n,l,w,p)|-cth
%end

%integerfunction ex(%integer n,l,r)
  %result=(ptlist(l)_pt(n)_x+ptlist(l)_pt(r)_x)>>1 - ptlist(l)_pt((n+r)>>1)_x
%end

%integerfunction ey(%integer n,l,r)
  %result=(ptlist(l)_pt(n)_y+ptlist(l)_pt(r)_y)>>1 - ptlist(l)_pt((n+r)>>1)_y
%end

%integerfunction dev(%integer n,w,r)
  %result=|ex(n,l,r)+ey(n,l,r)|-w
%end

findlines(ptlist,numlines)

p=w*2
k=w*2
totnl=0

%for l=0,1,numlines-1 %cycle
  %if ptlist(l)_num>w*2 %start
    numpoints=1
    list(0,1)=ptlist(l)_pt(0)_x
    list(0,2)=ptlist(l)_pt(0)_y
    %for i=w*2,1,ptlist(l)_num-1 %cycle
      j=a(i,l,th,w)
      %if j>0 %start
        list(numpoints,1)=ptlist(l)_pt(i-w)_x
        list(numpoints,2)=ptlist(l)_pt(i-w)_y
        numpoints=numpoints+1
      %finish
      %if k=i %start
        j=c(k,w,p,cth)
        %if j>0 %start
          list(numpoints,1)=ptlist(l)_pt(i-w)_x
          list(numpoints,2)=ptlist(l)_pt(i-w)_y
          numpoints=numpoints+1
          p=k-w
          k=p+2*w
        %finish
      %finish
      %if k=i %then k=k+1
    %repeat
    list(numpoints,1)=ptlist(l)_pt(ptlist(l)_num-1)_x
    list(numpoints,2)=ptlist(l)_pt(ptlist(l)_num-1)_y
    numpoints=numpoints+1
    %for i=0,1,numpoints-2 %cycle
      lines(totnl,1)=list(i,1)
      lines(totnl,2)=list(i,2)
      lines(totnl,3)=list(i+1,1)
      lines(totnl,4)=list(i+1,2)
      totnl=totnl+1
    %repeat
  %finish
%repeat
%end

%routine dump bitmap(%bytearrayname image(0:128*512))
%integer ref,size
  open output(1,file.".thn")
  select output(1)
  printsymbol('*')
  closeoutput
  selectoutput(0)
  accessfile(file.".thn",1,ref,size)
  writeregion(ref,0,128*512,image(0)) 
  deaccessfile(ref)
  printstring("Done");newline
%end

clear
enable(7)

prompt("Input File:");read(file)

accessfile(file,0,ref,size)
readregion(ref,0,size,image(0))
deaccessfile(ref)

%for i=0,1,128*512 %cycle
  copy(i)=image(i)
%repeat

displayframe(0,0,copy(0))
thinb(copy)
!dumpbitmap(copy)
%for i=0,1,128*512 %cycle
  image(i)=copy(i)
%repeat

%cycle
%for i=0,1,128*512 %cycle
  copy(i)=image(i)
%repeat
prompt("Worm length:");read(w)
prompt("Threshold:");read(th)
prompt("Curve Threshold:");read(cth)

setcolour(black)
fill(0,135,512,270)
displayframe(0,0,image(0))

worm(copy,lines,numlines,w,th,cth)
openoutput(1,"w.lis")
selectoutput(1)
col=blue
write(numlines,3);newline
%for i=0,1,numlines-1 %cycle
  colour(col)
  write(lines(i,1),3);write(lines(i,2),3)
  write(lines(i,3),3);write(lines(i,4),3);newline
!  line(lines(i,1),lines(i,2),lines(i,3),lines(i,4))
  line(lines(i,1),135+lines(i,2),lines(i,3),135+lines(i,4))
  %if col=blue %then col=green %else col=blue
%repeat
closeoutput
selectoutput(0)
%repeat
%endofprogram
