function [edges, theta, rho] = sobel(im)
    gx = [1 0 -1; 2 0 -2; 1 0 -1];
    gy = [1 2 1; 0 0 0; -1 -2 -1];
    [height, width] = size(im);
    edges = zeros(height, width);
    theta = edges;
    theta(1:height, 1:width) = inf;
    rho = theta;
    for i = 2:height - 1
      for j = 2:width - 1
           new_x = 0;
           new_y = 0;
           for hw = -1 :1
               for wi = -1:1
                   c = im(j+wi,i+hw);
                   new_x = new_x + gx(hw + 2, wi + 2) * c;
                   new_y = new_y + gy(hw + 2, wi + 2) * c;
               end
           end
           grad = sqrt(new_x^2 + new_y^2);
           rho(j, i) = grad;
           if new_x == 0
               theta(j, i) = 0;
           else
               theta(j, i) = posori(atan(new_y / new_x));
           end
           if grad > 0
               edges(j, i) = 1;
           end
      end
    end
end
