with(ImageTools):
img := Read("pentagon.png")[..,..,1]:
img_x := Convolution (img, Matrix ([[1,2,1], [0,0,0],[-1,-2,-1]])):
img_y := Convolution (img, Matrix ([[-1,0,1],[-2,0,2],[-1,0,1]])):
img := Array (abs (img_x) + abs (img_y), datatype=float[8]):
countPixels := proc(M)
	local r,c,i,j,row,col:
	row := Array([]);
	col := Array([]);
	r,c := LinearAlgebra:-Dimensions(M);
	for i from 1 to r do
		for j from 1 to c do
			if M[i,j] <> 0 then
				ArrayTools:-Append(row, i, inplace=true):
				ArrayTools:-Append(col, j, inplace=true):
			end if:
		end do:
	end do:
	return row,col:
end proc:
row,col := countPixels(img);
pTheta := proc(acc,r,c,x,y)
	local j, pos:
	for j from 1 to c do
		pos := ceil(x*cos((j-1)*Pi/180)+y*sin((j-1)*Pi/180)+r/2):
		acc[pos,j] := acc[pos,j]+1;
	end do:
end proc:
HoughTransform := proc(img,row,col)
   local r,c,pMax,theta,numThetas,numPs,acc,i:
   r,c := LinearAlgebra:-Dimensions(img);
   pMax := ceil(sqrt(r^2+c^2)):
   theta := [seq(evalf(i), i = 1..181, 1)]:
   numThetas := numelems(theta):
   numPs := 2*pMax+1:
   acc := Matrix(numPs, numThetas, fill=0,datatype=integer[4]):
   for i from 1 to numelems(row) do
   	pTheta(acc,numPs,numThetas,col[i],row[i]):
   end do:
   return acc;
end proc:
result :=HoughTransform(img,row,col);
Embed(Scale(FitIntensity(Create(result)), 1..500,1..500));