Approximating the integral over hemisphere domain with sample data

The integral part of the rendering equation performs an integral over a hemisphere range (theta = 0~PI/2 and phi = 0~2PI). I can generate samples that carry the value of corresponding integrand with this code:

n = 10; rho_s = 0.5; rho_d = 0.5;

light_phi = degtorad(30); light_theta = degtorad(60); [lx ly lz] = sph2cart(light_phi, light_theta, 1);

refDir = cat(3, -lx, -ly, lz); normal = cat(3, 0, 0, 1);

[sample_phi sample_theta] = meshgrid(0:12:360, 0:3:90); [sample_x sample_y sample_z] = sph2cart(degtorad(sample_phi), degtorad(sample_theta), 1); viewDir = cat(3, sample_x, sample_y, sample_z);

dotRL = sum(bsxfun(@times, viewDir, refDir), 3); dotRL = max(dotRL, 0);

brdf_s = rho_s*(n+2)*(dotRL.^n)/(2*pi); brdf_d = rho_d/pi; brdf = brdf_s + brdf_d;

x = sample_x.*brdf; y = sample_y.*brdf; z = sample_z.*brdf; mesh(x, y, z); line([0 lx],[0 ly],[0 lz]);

axis equal; axis vis3d; xlim([-1 1]); ylim([-1 1]); zlim([0 1]);

The question I face now is how to do the integral or approximation with these samples point ?

Thanks.