const WHITE=0xffFFff, X=0x010101;
fcn houghTransform(image,hx=460,hy=360){
   if(hy.isOdd) hy-=1; // hy argument must be even
   out:=PPM(hx,hy,WHITE);
   rMax:=image.w.toFloat().hypot(image.h);
   dr,dTh:=rMax/(hy/2), (0.0).pi/hx;

   foreach y,x in (image.h,image.w){
      if(image[x,y]==WHITE) continue;
      foreach iTh in (hx){
         th,r:=dTh*iTh, th.cos()*x + th.sin()*y;
	 iry:=hy/2 + (r/dr + 0.5).floor();  // y==0 is top
	 if(out[iTh,iry]>0) out[iTh,iry]=out[iTh,iry] - X;
      }
   }
   out
}