function [XC,YC,R] = drawregions(IMAGE,TOL,NPTS,DR)
F = imread(IMAGE);
FG = rgb2gray(F);
[DY,DX] = size(FG);
XC = 0;
YC = 0;
K = 0;
for J = 1:DX
    for I = 1:DY
        if (F(I,J) <= TOL)
            XC = XC + J;
            YC = YC + I;
            K = K + 1;
        end;
    end;
end;
XC = uint8(XC/K);
YC = uint8(YC/K);
I = 0;
while F(YC+I,XC)-F(YC,XC) <= TOL
    I = I + 1;
end;
R = I;
imshow(F);
hold all;
THETA = linspace(5.0*pi/8.0,19.0*pi/8.0,NPTS);
for I = 1:5
    X = (R+I*DR)*cos(THETA)+double(XC);  
    Y = (R+I*DR)*sin(THETA)+double(YC);
    plot(X,Y,'y-');
end;
RHO = linspace(R+DR,R+5*DR,2);
X = RHO*cos(5.0*pi/8.0)+double(XC);
Y = RHO*sin(5.0*pi/8.0)+double(YC);
plot(X,Y,'y-');
X = RHO*cos(19.0*pi/8.0)+double(XC);
Y = RHO*sin(19.0*pi/8.0)+double(YC);
plot(X,Y,'y-');
THETA = linspace(3.0*pi/4.0,11.0*pi/8.0,NPTS);
for I = 1:3
    X = (R+(I+5)*DR)*cos(THETA)+double(XC);  
    Y = (R+(I+5)*DR)*sin(THETA)+double(YC);
    plot(X,Y,'y-');
end;
THETA = linspace(-3.0*pi/8.0,pi/4.0,NPTS);
for I = 1:3
    X = (R+(I+5)*DR)*cos(THETA)+double(XC);  
    Y = (R+(I+5)*DR)*sin(THETA)+double(YC);
    plot(X,Y,'y-');
end;
RHO = linspace(R+5*DR,R+11*DR,2);
X = RHO*cos(3.0*pi/4.0)+double(XC);
Y = RHO*sin(3.0*pi/4.0)+double(YC);
plot(X,Y,'y-');
X = RHO*cos(pi/4.0)+double(XC);
Y = RHO*sin(pi/4.0)+double(YC);
plot(X,Y,'y-');
RHO = linspace(R+5*DR,R+8*DR,2);
X = RHO*cos(11.0*pi/8.0)+double(XC);
Y = RHO*sin(11.0*pi/8.0)+double(YC);
plot(X,Y,'y-');
X = RHO*cos(-3.0*pi/8.0)+double(XC);
Y = RHO*sin(-3.0*pi/8.0)+double(YC);
plot(X,Y,'y-');
THETA = linspace(3.0*pi/4.0,5.0*pi/4.0,NPTS);
for I = 1:3
    X = (R+(I+8)*DR)*cos(THETA)+double(XC);  
    Y = (R+(I+8)*DR)*sin(THETA)+double(YC);
    plot(X,Y,'y-');
end;
THETA = linspace(-pi/4.0,pi/4.0,NPTS);
for I = 1:3
    X = (R+(I+8)*DR)*cos(THETA)+double(XC);  
    Y = (R+(I+8)*DR)*sin(THETA)+double(YC);
    plot(X,Y,'y-');
end;
RHO = linspace(R+8*DR,R+11*DR,2);
X = RHO*cos(5.0*pi/4.0)+double(XC);
Y = RHO*sin(5.0*pi/4.0)+double(YC);
plot(X,Y,'y-');
X = RHO*cos(-pi/4.0)+double(XC);
Y = RHO*sin(-pi/4.0)+double(YC);
plot(X,Y,'y-');
end