function [out2, pp] = drawball(fname,scale, threshold, func, shading_thres)
clf
colormap(gray(256));
%axis([0 4000 -2000 1500 0 1500]);
r=0.4;
[x,y,z] = sphere(10);
[nx,ny,nz] = surfnorm(x,y,z);
refl=diffuse(nx,ny,nz,[0.5 -1 1]);
out2=load_matrix(fname);
%axis([ min(arcHo(:,1))-r max(arcHo(:,1))+r
%      min(arcHo(:,2))-r max(arcHo(:,2))+r
%      min(arcHo(:,3))-r max(arcHo(:,3))+r ]);
out2(:,4) = (out2(:,4)>=threshold).*out2(:,4)/scale;
switch func
  case 1, out2(:,4) = log2(out2(:,4));
  case 2, out2(:,4) = log(out2(:,4));
  case 3, out2(:,4) = log10(out2(:,4));
  case {4,5,6,7,8,9,10}, out2(:,4) = (out2(:,4).^(1/(func-3)));
  % otherwize, out2(:,4) = out2(:,4);
end
%out2=[out2(:,1) out2(:,2) out2(:, 3) out2(:,4)/40.0];
[m,n]=size(out2);
balls_drawn = 0;
msg = sprintf('Total balls to be drawn: %d', sum(out2(:,4)>0));
disp(msg);
for i=1:m,   
    if out2(i,4)>0, 
      a=out2(i,1);b=out2(i,2);c=out2(i,3);    
      r=out2(i,4);
      p=surf(r*x+a,r*y+b,r*z+c,refl);
      if    out2(i,4) > shading_thres,
         shading interp;
      end;
      hold on;
      balls_drawn = balls_drawn +1;
      if mod(balls_drawn,100)==0,
         msg = sprintf('%d drawn already', balls_drawn);
         disp(msg);
      end;
    end;
end;
%rotate3d on;
axis equal;
box on;
%title('Circle Hough Transform for signal No. 1');
msg = sprintf('Spheres drawn %d', balls_drawn);
disp(msg);
pp=get(p,'Parent');
set(pp,'XTickLabel',[]);
set(pp,'YTickLabel',[]);
set(pp,'ZTickLabel',[]);
set(pp,'ZTick',[]);     
set(pp,'XTick',[]);
set(pp,'YTick',[]);
%set(pp,'XLabel', ['L^*']);
%set(pp,'YLabel',['u^*']);
%set(pp,'ZLabel',['v^*']);
axis('on');
set(pp,'XLim',[0 100]);
set(pp,'YLim',[-45 100]);
set(pp,'ZLim',[-53 70]);
p=get(pp,'Parent');
set(p,'Color',[1 1 1]);


