function avg_val = cv_sub_pixel(num_rows, num_cols, x_side, y_side, ...
   x_center, y_center, theta, low_val, high_val)
%
line_eps = 0.01;
xmin = -x_side/2;
xmax = -xmin;
x_step = x_side/num_cols;
v_line_x = xmin:x_step:xmax;

ymin = -y_side/2;
ymax = -ymin;
y_step = y_side/num_rows;
v_line_y = ymin:y_step:ymax;

num_lines_x = length(v_line_x);
num_lines_y = length(v_line_y);
y_vals1 = v_line_y(1)*ones(num_lines_x,1);
y_vals2 = v_line_y(2)*ones(num_lines_x,1);
v_lines = [v_line_x' y_vals1 v_line_x' y_vals2];

x_vals1 = v_line_x(1)*ones(num_lines_y,1);
x_vals2 = v_line_x(2)*ones(num_lines_y,1);
h_lines = [x_vals1 v_line_y' x_vals2 v_line_y'];

[num_v_lines,dummy] = size(v_lines);
[num_h_lines,dummy] = size(h_lines);

rtheta = pi*theta/180;
ideal_line = [x_center, y_center 0; x_center+cos(rtheta) y_center+sin(rtheta) 0];

intercept_points = zeros(4,3);
[intercept_points(1,1:3), left_status] = cv_int_line_line(...
   [v_lines(1,1:2) 0; v_lines(1,3:4) 0], ideal_line, line_eps);
[intercept_points(2,1:3), top_status] = cv_int_line_line(...
   [h_lines(num_h_lines,1:2) 0; h_lines(num_h_lines,3:4) 0], ideal_line, line_eps);
[intercept_points(3,1:3), right_status] = cv_int_line_line(...
   [v_lines(num_v_lines,1:2) 0; v_lines(num_v_lines,3:4) 0], ideal_line, line_eps);
[intercept_points(4,1:3), bot_status] = cv_int_line_line(...
   [h_lines(1,1:2) 0; h_lines(1,3:4) 0], ideal_line, line_eps);
intercept_points

in_bounds = zeros(4,1);
in_bounds(1) = intercept_points(1,1) <= xmax & intercept_points(1,1) >= xmin ...
             & intercept_points(1,2) <= ymax & intercept_points(1,2) >= ymin;
in_bounds(2) = intercept_points(2,1) <= xmax & intercept_points(2,1) >= xmin ...
             & intercept_points(2,2) <= ymax & intercept_points(2,2) >= ymin;
in_bounds(3) = intercept_points(3,1) <= xmax & intercept_points(3,1) >= xmin ...
             & intercept_points(3,2) <= ymax & intercept_points(3,2) >= ymin;
in_bounds(4) = intercept_points(4,1) <= xmax & intercept_points(4,1) >= xmin ...
             & intercept_points(4,2) <= ymax & intercept_points(4,2) >= ymin;
num_in_bounds = sum(in_bounds)

if num_in_bounds<2
   disp(' Line out of bounds');
   return
end

if num_in_bounds == 2
   line_ends = zeros(2,3);
   l_index = 1;
   for l = 1:4
      if in_bounds(l) == 1
         line_ends(l_index,1:3) = intercept_points(l,1:3);
         l_index = l_index+1;
      end
   end
end

if num_in_bounds>2
   disp('Too many side intersections');
   return
end

image_out = zeros(num_rows, num_cols);
avg = (low_val+high_val)/2;
l_step = min(x_step,y_step)/4;
for t = 0:l_step:1
   p = line_ends(1,1:3) + t*(line_ends(2,1:3)-line_ends(1,1:3));
   r = floor(num_rows*((p(2)-ymin)/y_side))+1;
   c = floor(num_cols*((p(1)-xmin)/x_side))+1;
   if r <= num_rows & c <= num_cols
      image_out(r,c) = avg;
   end
end

q1 = line_ends(2,1:3) - line_ends(1,1:3);
for r = 1:num_rows
   y = h_lines(r,2) + y_step/2;
   for c = 1:num_cols
      x = v_lines(c,1) + x_step/2;
      if image_out(r,c) > 0
         % set actual value
         dummy = r;
      else
         q2 = [x y 0] - line_ends(1,1:3);
         v = cross(q1,q2);
         if v(3)<0
            image_out(r,c) = low_val;
         end
         if v(3)>0
            image_out(r,c) = high_val;
         end
      end
   end
end
avg_val = sum(sum(image_out))/(num_rows*num_cols);
