diff --git a/src/init.h b/src/init.h
index ae4a1dd..20a65e4 100644
--- a/src/init.h
+++ b/src/init.h
@@ -1,2 +1,9 @@
+#include <gsl/gsl_interp.h>
+#include <gsl/gsl_spline.h>
 double chi(double cotpsi);
 int init(double rmax, double rela_err_limit, int *size, double **x, double **y);
+
+typedef struct {
+  gsl_spline *spline;
+  gsl_interp_accel *acc;
+} Spline_data ;
diff --git a/src/render.c b/src/render.c
index ebae401..cde9d38 100644
--- a/src/render.c
+++ b/src/render.c
@@ -1,8 +1,16 @@
 #include <gsl/gsl_interp.h>
+#include <gsl/gsl_monte.h>
+#include <gsl/gsl_rng.h>
 #include <gsl/gsl_spline.h>
+#include <gsl/gsl_monte_miser.h>
 #include "common.h"
 #include "init.h"
 
+typedef struct {
+  const System *system;
+  const Spline_data *spline_data;
+} Integrand_params;
+
 inline double xy_to_b(double x, double y) {
   double tanpsi2 = x*x+y*y;
   double cotpsi2 = 1/tanpsi2;
@@ -10,19 +18,79 @@ inline double xy_to_b(double x, double y) {
   return b;
 }
 
-int MC_pixel_render(const System *system, int i, int j, const gsl_spline *spline, double *rgb, int pixel_render_max, double pixel_render_err){
+double integrand_r(double *xy, size_t dim, void *integrand_params_void) {
+  double cotpsi = 1.0/hypot(xy[0], xy[1]);
+  Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
+  double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+  double phi = atan2(xy[1], xy[0]);
+  double rgb[3] = {0};
+  double theta_phi[2] = {theta, phi};
+  integrand_params->system->angle_to_pixel(theta_phi, rgb);
+  return rgb[0]*gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+}
+
+double integrand_g(double *xy, size_t dim, void *integrand_params_void) {
+  double cotpsi = 1.0/hypot(xy[0], xy[1]);
+  Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
+  double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+  double phi = atan2(xy[1], xy[0]);
+  double rgb[3] = {0};
+  double theta_phi[2] = {theta, phi};
+  integrand_params->system->angle_to_pixel(theta_phi, rgb);
+  return rgb[1]*gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+}
+
+double integrand_b(double *xy, size_t dim, void *integrand_params_void) {
+  double cotpsi = 1.0/hypot(xy[0], xy[1]);
+  Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
+  double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+  double phi = atan2(xy[1], xy[0]);
+  double rgb[3] = {0};
+  double theta_phi[2] = {theta, phi};
+  integrand_params->system->angle_to_pixel(theta_phi, rgb);
+  return rgb[2]*gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+}
+
+int MC_pixel_render(const System *system, int i, int j, const Spline_data spline_data, double *rgb, int pixel_render_max, double pixel_render_err, gsl_rng *r){
   double dw = (system->w)/(system->W); // width and hight of one pixel
   double x_lu = i*dw - (system->w)/2;
   double y_lu = -j*dw + (system->H)*dw/2;
+  double xl[2] = {x_lu, y_lu - dw};
+  double xu[2] = {x_lu + dw, y_lu};
 
   //TODO: fill the MC intergral
-  // linear rgb = ∫ rgb(χ, ϕ) sin^2(χ) χ'(ψ) dx dy
+  // linear rgb = ∫ rgb(χ, ϕ) sin^2(χ) |χ'(ψ)| dx dy
+  // χ'(ψ) = χ'(cot(ψ)) * cot'(ψ) = - χ'(cot(ψ))/\sin^2(ψ)
+  gsl_monte_miser_state *miser_state = gsl_monte_miser_alloc(2);
+
+  double color;
+  double err;
+  Integrand_params integrand_params = {system, &spline_data};
+  gsl_monte_function F;
+  F.params = &integrand_params;
+  F.dim = 2;
+  for (int c = 0; c < 3; c++) {
+    switch (c) {
+      case 0:
+        F.f = integrand_r;
+        break;
+      case 1:
+        F.f = integrand_g;
+        break;
+      case 2:
+        F.f = integrand_b;
+        break;
+    }
+    gsl_monte_miser_init(miser_state);
+    gsl_monte_miser_integrate(&F, xl, xu, 2, 100, r, miser_state, &color, &err);
+    rgb[c] = color;
+  }
 
   return 0;
 }
 
-int pixel_render(const System *system, int i, int j, const gsl_spline *spline, double *rgb, int pixel_render_max, double pixel_render_err){
-  return MC_pixel_render(system, i, j, spline, rgb, pixel_render_max, pixel_render_err);
+int pixel_render(const System *system, int i, int j, const Spline_data spline_data, double *rgb, int pixel_render_max, double pixel_render_err, gsl_rng *r){
+  return MC_pixel_render(system, i, j, spline_data, rgb, pixel_render_max, pixel_render_err, r);
 }
 
 int render(System *system, double **buffer, int pixel_render_max, double pixel_render_err, double chi_rela_err) {
@@ -31,6 +99,7 @@ int render(System *system, double **buffer, int pixel_render_max, double pixel_r
   double cotpsi2 = 1/tanpsi2;
   //double bmax = R0/sqrt(f(R0)*(1+cotpsi2));
   double rmax = sqrt(tanpsi2);
+  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  
   double *x=NULL;
   double *y = NULL;
@@ -40,11 +109,12 @@ int render(System *system, double **buffer, int pixel_render_max, double pixel_r
 
   gsl_interp_accel *acc = gsl_interp_accel_alloc();
   gsl_spline *spline = gsl_spline_alloc(gsl_interp_steffen, size);
+  Spline_data spline_data = {spline, acc};
   gsl_spline_init(spline, x, y, size);
 
   for(int j = 0; j < system->H; j++) {
     for (int i = 0; i < system->W; i++) {
-      pixel_render(system, i, j, spline, buffer[j*(system->W) + i], pixel_render_max, pixel_render_err);
+      pixel_render(system, i, j, spline_data, buffer[j*(system->W) + i], pixel_render_max, pixel_render_err, r);
     }
   }