pro draw_grids, style, draw_cm, grid_info, slice_ori, slice_coord, xy_sl_size, image=image

;if style eq 2  then  erase

; determine grids in slice
if style eq 2 then grid_flag = grid_in_projection(grid_info, slice_ori, slice_coord, [0.,1.]) $
else grid_flag = grid_in_slice(grid_info, slice_ori, slice_coord) 
if total(grid_info[0].dim[*] eq 1) then grid_flag = grid_in_projection(grid_info, slice_ori, slice_coord, [0.,1.]) ; 2D data?

;print, grid_flag 

axes = (slice_ori gt 1.e-30)
const_sub = where(axes gt 1e-10)
other_subs = where(axes eq 0)
slice_size = [ xy_sl_size,  xy_sl_size]

min_left  = [slice_coord(0),slice_coord(1)]
max_right = [slice_coord(2),slice_coord(3)]
scale_up    = min(1.D/(max_right-min_left))

image = fltarr(xy_sl_size, xy_sl_size)
ind = where(grid_flag gt 0)
if ind[0] eq -1 then return

u_g = grid_info[ind]
index_range = u_g.End_index-u_g.Start_index+1
delta_distance = (u_g.Right_edge - u_g.Left_edge)/(index_range)

dim_grid_info = size(u_g)
num_of_grids  = dim_grid_info(1)


for i=0L, num_of_grids-1 DO BEGIN
      delta_dist   = delta_distance(*,i)
      data_points  = index_range(*,i)  
      slice_ori_l                = DBLarr(3)
      slice_ori_l(const_sub)     = slice_ori(const_sub)
      slice_ori_l(other_subs)    = [slice_coord(0),slice_coord(1)]
      slice_ori_r                = slice_ori_l 
      slice_ori_r(other_subs)    = [slice_coord(2),slice_coord(3)]

      i_s = ROUND((slice_ori_l(const_sub) - $
                   u_g(i).Left_edge(const_sub))/delta_dist(const_sub) >0 $
        < data_points(const_sub)-1)
      i_l = CEIL( (slice_ori_l(other_subs) - $
                   u_g(i).Left_edge(other_subs)) $
                    /delta_dist(other_subs) ) > 0
      i_r   = (FLOOR((slice_ori_r(other_subs) - $
                  u_g(i).Left_edge(other_subs))/ $
                    delta_dist(other_subs) - 1) < $
        (data_points(other_subs)-1) ) > 0

      L_p = (FIX(DOUBLE(slice_size) *scale_up * $
                         ((u_g(i).Left_edge(other_subs)+ $
                           (i_l)*delta_dist(other_subs)-min_left)))  < (xy_sl_size-1) ) > 0
      R_p= ((L_p + $
        Round(DOUBLE(slice_size) *scale_up * $
              DOUBLE(i_r-i_l+1)*delta_dist(other_subs))) < (xy_sl_size-1) ) > 0
  
    box = [[L_p(0), L_p(1)], $
           [L_p(0), R_p(1)], $
           [R_p(0), R_p(1)], $
           [R_p(0), L_p(1)], $
           [L_p(0), L_p(1)]]

    CASE style OF 
        1 : BEGIN
            plots, box, /device, LINESTYLE = 1
            xyouts, R_p(0), R_p(1), STRCOMPRESS(STRING(u_g(i).level)), $
              CHARSIZE = 0.7, /DEVICE, ALIGNMENT = 1.0
        END
        ELSE: BEGIN 
           if (R_P[0] ge L_p[0]) and (R_P[1] ge L_p[1]) then $
              image[L_p[0]:R_P[0],L_p[1]:R_P[1]] = u_g[i].level 
;           colors = indgen(45)*3+20
 ;           polyfill, box+DOUBLE([-1.,-1.,-1.,1.,1.,1.,1.,-1.,-1.,-1.]),$
 ;             /device, col = colors[u_g(i).level]
 ;           next_level = 0
;            if i lt num_of_grids-1 then next_level = u_g(i+1).level $
;            else next_level =200
;            current_level = u_g(i).level 
;            if next_level gt current_level then BEGIN
;    xyouts, R_p(0), R_p(1), STRCOMPRESS(STRING(u_g(i).num)), $
;                xyouts, R_p(0), R_p(1), STRCOMPRESS(STRING(u_g(i).level)), $
;                  CHARSIZE = 0.7, /DEVICE, ALIGNMENT = 1.0
;            ENDIF
        END
    END
ENDFOR

; for the drawe chart erase colormap window 
IF style ne 1 THEN BEGIN
    old_wind = !d.window
    widget_control, get_value = win_index, draw_cm
    wset, win_index
    erase
    wset, old_wind
ENDIF
END