home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / rad386 / radiosit / src / radiosit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-22  |  19.0 KB  |  659 lines
  1. /*
  2.     Computation of Form factor for Full Radiosity Solution.
  3. */
  4. #include <stdio.h>
  5. #include <math.h>
  6. #ifdef __WATCOMC__
  7. #include <string.h>
  8. #else
  9. #include <memory.h>
  10. #endif
  11. #include <malloc.h>
  12.  
  13. #include "GraphicsGems.h"
  14. #include "data_structure.h"
  15. #include "objects.h"
  16. #include "vox.h"
  17.  
  18. #include "raddecl.h"
  19.  
  20. extern int hemicube_resolution;
  21. extern int verbose_flag;
  22. int ProjectionByRaytracing=0;
  23.  
  24. static double **TopFace_delta_form_factor;
  25. static double **SideFace_delta_form_factor;
  26. struct Hemicube{
  27.         short int *topface;
  28.         short int *rightface;
  29.         short int *leftface;
  30.         short int *frontface;
  31.         short int *backface;
  32. };
  33. static double **F; /* Form Factors */
  34. static double **A; /* N x N Matrix */
  35. static double **I;  /* Computed Intensity */
  36. static double **deltaI;  /* Delta Intensity for progressive radiosity */
  37. static double *B;
  38. static struct Hemicube hemicube;
  39.  
  40. #define Iter_Threshold 0.0001
  41. #define VeryHigh 1.e5
  42.  
  43. static double **alloc_matrix(rows,cols)
  44. int rows,cols;
  45. {
  46.     char *s="Cannot Allocate Matrix.";
  47.     double **m;
  48.     int i;
  49.     m = (double **)malloc(rows*sizeof(double *));
  50.     if (m==NULL) error(s);
  51.     for (i=0;i<rows;i++) {
  52.         m[i] = (double *)malloc(cols*sizeof(double));
  53.         if (m[i] == NULL) error(s);
  54.     }
  55.     return(m);
  56. }
  57. static void allocate_hemicube()
  58. {
  59.     TopFace_delta_form_factor=
  60.         alloc_matrix(hemicube_resolution,hemicube_resolution);
  61.     SideFace_delta_form_factor=
  62.         alloc_matrix(hemicube_resolution/2,hemicube_resolution);
  63.     hemicube.topface=(short int *)malloc(
  64.         hemicube_resolution*hemicube_resolution
  65.         *sizeof(short int));
  66.     hemicube.rightface=(short int *)malloc(
  67.                 hemicube_resolution/2*hemicube_resolution
  68.                 *sizeof(short int));
  69.     hemicube.leftface=(short int *)malloc(
  70.                 hemicube_resolution/2*hemicube_resolution
  71.                 *sizeof(short int));
  72.     hemicube.frontface=(short int *)malloc(
  73.                 hemicube_resolution/2*hemicube_resolution
  74.                 *sizeof(short int));
  75.     hemicube.backface=(short int *)malloc(
  76.                 hemicube_resolution/2*hemicube_resolution
  77.                 *sizeof(short int));
  78. }
  79. static void free_matrix(m,rows,cols)
  80. double **m;
  81. int rows,cols;
  82. {
  83.     int i;
  84.     for (i=0;i<rows;i++) free(m[i]);
  85.     free(m);
  86. }
  87. static void deallocate_hemicube()
  88. {
  89.     free_matrix(TopFace_delta_form_factor,
  90.         hemicube_resolution,hemicube_resolution);
  91.     free_matrix(SideFace_delta_form_factor,
  92.         hemicube_resolution/2,hemicube_resolution);
  93.     free(hemicube.topface);
  94.     free(hemicube.rightface);
  95.     free(hemicube.leftface);
  96.     free(hemicube.frontface);
  97.     free(hemicube.backface);
  98. }
  99.  
  100. /*
  101.         compute_delta_form_factors : For the faces of the hemicube.
  102.                 Side faces are symmetric, so only one computation
  103.                 is enough.
  104. */
  105. static compute_delta_form_factors()
  106. {
  107.         int i,j;
  108.         double x, y, y2, delta;
  109.         double factor;
  110.         double z,zfactor,z2;
  111.         double halfdelta;
  112.         double delta_area_by_pi;
  113.  
  114.     delta  = 2.0/hemicube_resolution;
  115.     halfdelta =delta/2.0;
  116.     delta_area_by_pi=delta*delta/PI;
  117.  
  118.         for(i=0,y=(-1.0+halfdelta); i < hemicube_resolution; i++,y+=delta){
  119.                 y2 = y*y;
  120.                 for (j=0,x=(-1.0+halfdelta); j < hemicube_resolution; j++,x+=delta){
  121.                         factor=1.0/(x*x+y2+1.0);
  122.                         TopFace_delta_form_factor[i][j]=
  123.                                 delta_area_by_pi*factor*factor;
  124.                 }
  125.         }
  126.         for (i=0,z=halfdelta; i < hemicube_resolution/2; i++,z+=delta){
  127.                 z2 = z*z; zfactor = z*delta_area_by_pi;
  128.                 for(j=0,y=(-1.0+halfdelta); j < hemicube_resolution; j++,y+=delta){
  129.                         factor=1.0/(y*y+z2+1.0);
  130.                         SideFace_delta_form_factor[i][j]=zfactor*factor*factor;
  131.                 }
  132.         }
  133. }
  134.  
  135. static void increment_F(rows,raster,f,inc)
  136. int rows;
  137. short int *raster;
  138. double f[];
  139. double **inc;
  140. {
  141.     int i,j,k,patchnum;
  142.     for(i=k=0;i<rows;i++)
  143.         for(j=0;j<hemicube_resolution;j++,k++)
  144.             if (patchnum=raster[k])
  145.                 f[patchnum-1] += inc[i][j];
  146. }
  147.  
  148. static void extract_F_i(i)
  149. int i;
  150. {
  151.     int nn;
  152.            for(nn=0;nn<total_surface_grid_points;nn++)F[i][nn]=0.;
  153.     increment_F(hemicube_resolution,hemicube.topface,
  154.         F[i],TopFace_delta_form_factor);
  155.     increment_F(hemicube_resolution/2,hemicube.rightface,
  156.         F[i],SideFace_delta_form_factor);
  157.     increment_F(hemicube_resolution/2,hemicube.leftface,
  158.         F[i],SideFace_delta_form_factor);
  159.     increment_F(hemicube_resolution/2,hemicube.frontface,
  160.         F[i],SideFace_delta_form_factor);
  161.     increment_F(hemicube_resolution/2,hemicube.backface,
  162.         F[i],SideFace_delta_form_factor);
  163. }
  164.  
  165. static init_fullmatrix_radiosity()
  166. {
  167.     /* Allocate space for form factor */
  168.     F = alloc_matrix(total_surface_grid_points,total_surface_grid_points);
  169.     A = alloc_matrix(total_surface_grid_points,total_surface_grid_points);
  170.     I = alloc_matrix(color_channels,total_surface_grid_points);
  171.     B = (double *)malloc(total_surface_grid_points*sizeof(double));
  172.     allocate_hemicube();
  173.     compute_delta_form_factors();
  174. }
  175. static deinit_fullmatrix_radiosity()
  176. {
  177.     /* Free space allocated for FormFactor */
  178.     free_matrix(F,total_surface_grid_points,total_surface_grid_points);
  179.     free_matrix(A,total_surface_grid_points,total_surface_grid_points);
  180.     free_matrix(I,color_channels,total_surface_grid_points);
  181.     free(B);
  182.     deallocate_hemicube();
  183. }
  184.  
  185. static init_progressive_radiosity()
  186. {
  187.     /* Allocate space for form factor */
  188.     F = alloc_matrix(1,total_surface_grid_points);
  189.     I = alloc_matrix(color_channels,total_surface_grid_points);
  190.     deltaI = alloc_matrix(color_channels,total_surface_grid_points);
  191.     allocate_hemicube();
  192.     compute_delta_form_factors();
  193. }
  194. static deinit_progressive_radiosity()
  195. {
  196.     free_matrix(F,1,total_surface_grid_points);
  197.     free_matrix(I,color_channels,total_surface_grid_points);
  198.     free_matrix(deltaI,color_channels,total_surface_grid_points);
  199.     deallocate_hemicube();
  200. }
  201.  
  202. static void from_specular_surface(ray,objectnum,t,u,v,n)
  203. Ray *ray;
  204. int objectnum;
  205. double t,u,v;
  206. int *n;
  207. {
  208.     Vector3 N;
  209.         /* .... set the new ray */
  210.         point_on_line(ray->start,ray->direction,t,ray->start);
  211.         N = ofunc[object[objectnum].surface_geometry_type].get_surface_normal
  212.                 (object[objectnum].object_specific_structure,u,v);
  213.         mirror_reflection(&N,&(ray->direction),&(object[objectnum]));
  214.         (ray->number)++;
  215.     *n=get_diffuse_surface(ray);
  216. }
  217. static void from_diffuse_surface(ray,objectnum,t,u,v,n)
  218. Ray *ray;
  219. int objectnum;
  220. double t,u,v;
  221. int *n;
  222. {
  223.     int uindex=grid_u(objectnum,u);
  224.     int vindex=grid_v(objectnum,v);
  225.     *n=(object[objectnum].start_grid_index+
  226.         gridindex(objectnum,uindex,vindex)+1);
  227. }
  228. static void from_other_surface(ray,objectnum,t,u,v,n)
  229. Ray *ray;
  230. int objectnum;
  231. double t,u,v;
  232. int *n;
  233. {
  234.         error("Not Implemented.");
  235. }
  236. static void (*get_patchnum[MAX_REFLECTION_TYPE])()={
  237.         from_diffuse_surface,
  238.         from_specular_surface,
  239.         from_other_surface
  240. };
  241.  
  242. static int get_diffuse_surface(ray)
  243. Ray *ray;
  244. {
  245.     double t_entry,t_exit;
  246.     int HitBoundingBox();
  247.         if (HitBoundingBox(&(volume_grid.extent.min),&(volume_grid.extent.max),
  248.             &(ray->start),&(ray->direction),&t_entry,&t_exit)){
  249.                 int get_nearest_object_in_voxel();
  250.                 Vlist *create_vlist();
  251.                 Vlist *vlist,*templist;
  252.                 double t,u,v;
  253.             int nearest_object= UNDEFINED;
  254.                 vlist=templist=create_vlist(
  255.                         &(ray->start),&(ray->direction),t_entry,t_exit
  256.                 );
  257.                 while (templist!=NULL){
  258.                         Voxel *vox = volume_grid.voxels+templist->voxnum;
  259.                         if(vox->nobjects)
  260.                 nearest_object=get_nearest_object_in_voxel
  261.                 (ray,templist->t_near,templist->t_far,vox,&t,&u,&v);
  262.                         if (nearest_object!= UNDEFINED) break;
  263.                         templist=templist->next;
  264.                 }
  265.                 purge_vlist(vlist);
  266.                 if (nearest_object!= UNDEFINED){
  267.             int n;
  268.             get_patchnum[object[nearest_object].surface_reflection_type]
  269.             (ray,nearest_object,t,u,v,&n);
  270.             return(n);
  271.         }
  272.     }
  273.     return(0);
  274. }
  275. /*
  276. int dump_buffer(rows,cols,b)
  277. int rows,cols;
  278. char *b;
  279. {
  280.     int i,j,k;
  281.     fprintf(stderr,"\nDumping buffer\n");
  282.     for(i=k=0;i<rows;i++){
  283.         fprintf(stderr,"%02d.....",i);
  284.         for(j=0;j<cols;j++,k++)
  285.             fprintf(stderr,"%1d",b[k]);
  286.         fprintf(stderr,"\n");
  287.     }
  288. }
  289. int dump_raster(rows,cols,b)
  290. int rows,cols;
  291. short int *b;
  292. {
  293.     int i,j,k;
  294.     fprintf(stderr,"\nDumping raster\n");
  295.     for(i=k=0;i<rows;i++){
  296.         fprintf(stderr,"%02d.....",i);
  297.         for(j=0;j<cols;j++,k++)
  298.             fprintf(stderr,"%1d",b[k]);
  299.         fprintf(stderr,"\n");
  300.     }
  301. }
  302. */
  303. static long project(C,N,U,V,rows,raster)
  304. Point3 *C;
  305. Vector3 *N,*U,*V;
  306. int rows;
  307. short int *raster;
  308. {
  309.     static long RayNumber=0;
  310.     Ray ray;
  311.     Vector3 Dir,delta_U,delta_V,half_delta_U,half_delta_V;
  312.     int i,j,k;
  313.  
  314.     half_delta_U=delta_U = *U;
  315.     V3Scale(&delta_U,2./hemicube_resolution);
  316.     V3Scale(&half_delta_U,1./hemicube_resolution);
  317.     half_delta_V=delta_V = *V;
  318.     V3Scale(&delta_V,2./hemicube_resolution);
  319.     V3Scale(&half_delta_V,1./hemicube_resolution);
  320.     V3Add(N,V,&Dir); V3Sub(&Dir,U,&Dir);
  321.     V3Add(&Dir,&half_delta_U,&Dir); V3Sub(&Dir,&half_delta_V,&Dir);
  322.     ray.number=total_rays+RayNumber;
  323. #if defined (DEBUG)
  324.     ray.intersections=0;
  325. #endif
  326.     ray.path_length= -1;
  327.     if (ProjectionByRaytracing)
  328.     for(i=k=0;i<rows;i++){
  329.         Vector3 dir;
  330.         dir = Dir;
  331.         for(j=0;j<hemicube_resolution;j++,k++){
  332.             ray.start = *C;
  333.             ray.direction=dir;
  334.             V3Normalize(&(ray.direction));
  335.             raster[k]=get_diffuse_surface(&ray);
  336.             (ray.number)++;
  337.             V3Add(&dir,&delta_U,&dir);
  338.         }
  339.         V3Sub(&Dir,&delta_V,&Dir);
  340.     }
  341.     else {
  342.     double *zbuffer;
  343.     char *surface_buffer;
  344.     Matrix4 view_matrix;
  345.     int n=rows*hemicube_resolution;
  346.     extern void (*geometry_specific_scan_convert[MAX_GEOMETRY_TYPE])();
  347.     /* Clean up the raster */
  348.     memset((char *)raster,0,n*sizeof(short int));
  349.     /* Set up a Z-buffer of equal dimension */
  350.     zbuffer=(double *)malloc(n*sizeof(double));
  351.     if (zbuffer == NULL)
  352.       error("Cannot allocate Z-buffer memory for scan conversion.");
  353.     for(i=0;i<n;i++)zbuffer[i]=VeryHigh;
  354.     surface_buffer=(char *)calloc(n,1);
  355.     if (surface_buffer == NULL)
  356.       error("Cannot allocate surface-buffer memory for scan conversion.");
  357.     /* Compute View transform matrix */
  358.     view_transform(C,U,V,N,rows,hemicube_resolution,&view_matrix);
  359.     for(i=0;i<number_objects;i++){
  360.         /* Project each object on the raster */
  361.         geometry_specific_scan_convert[
  362.             object[i].surface_geometry_type
  363.         ](i,N,&view_matrix,rows,hemicube_resolution,
  364.             zbuffer,surface_buffer,raster);
  365.     }
  366.     free((char *)zbuffer);
  367. /*
  368. dump_raster(rows,hemicube_resolution,raster);
  369. dump_buffer(rows,hemicube_resolution,surface_buffer);
  370. */
  371.     for(i=k=0;i<rows;i++)
  372.     for(j=0;j<hemicube_resolution;j++,k++)
  373.         if ((surface_buffer[k]-1)>DIFFUSE){
  374.             ray.start = *C;
  375.             ray.direction.x=Dir.x+j*delta_U.x-i*delta_V.x;
  376.             ray.direction.y=Dir.y+j*delta_U.y-i*delta_V.y;
  377.             ray.direction.z=Dir.z+j*delta_U.z-i*delta_V.z;
  378.             V3Normalize(&(ray.direction));
  379.             raster[k]=get_diffuse_surface(&ray);
  380.             (ray.number)++;
  381.         }
  382.     free(surface_buffer);
  383.     }
  384.     return(RayNumber = (ray.number-total_rays));
  385. }
  386.  
  387. static long project_on_hemicube(c,m)
  388. Point3 *c;
  389. Matrix4 *m;
  390. {
  391.     Vector3 X,Y,Z,Xneg,Yneg,Zneg;
  392.     long n;
  393.  
  394.     X.x = 1.; X.y=X.z = 0.;    
  395.     X = transform_vector(&X,m); V3Normalize(&X); Xneg=X; V3Negate(&Xneg);
  396.     Y.y = 1.; Y.x=Y.z = 0.;    
  397.     Y = transform_vector(&Y,m); V3Normalize(&Y); Yneg=Y; V3Negate(&Yneg);
  398.     Z.z = 1.; Z.x=Z.y = 0.;    
  399.     Z = transform_vector(&Z,m); V3Normalize(&Z); Zneg=Z; V3Negate(&Zneg);
  400.     /* Project on Top Face */
  401.     project(c,&Z,&X,&Y,hemicube_resolution,hemicube.topface);
  402.     /* Project on Left Face */
  403.     project(c,&Xneg,&Y,&Z,hemicube_resolution/2,hemicube.leftface);
  404.     /* Project on Right Face */
  405.     project(c,&X,&Yneg,&Z,hemicube_resolution/2,hemicube.rightface);
  406.     /* Project on Front Face */
  407.     project(c,&Yneg,&Xneg,&Z,hemicube_resolution/2,hemicube.frontface);
  408.     /* Project on Back Face */
  409.     n=project(c,&Y,&X,&Z,hemicube_resolution/2,hemicube.backface);
  410.     return(n);
  411. }
  412.  
  413. static hemicube_method(i,j,k,n,nrays)
  414. int i;    /* Object num */
  415. int j;  /* Vertical Grid */
  416. int k;  /* Horizontal grid */
  417. int n;  /* Patch num */
  418. long *nrays;
  419. {
  420.     Matrix4 m;
  421.     Vector3 N;
  422.     Point3 centre;
  423.     if (object[i].surface_reflection_type==MIRROR)return;
  424.     centre = ofunc[object[i].surface_geometry_type].
  425.         compute_center_point_in_the_grid_cell(i,j,k);
  426.     N=ofunc[object[i].surface_geometry_type].
  427.         get_surface_normal
  428.         (object[i].object_specific_structure,
  429.         (k+0.5)/object[i].grid_h_reso,
  430.         (j+0.5)/object[i].grid_v_reso);
  431.     m  = ofunc[object[i].surface_geometry_type].
  432.         get_local_matrix(
  433.         &N,object[i].object_specific_structure);
  434.     *nrays=project_on_hemicube(¢re,&m);
  435. }
  436. static long compute_F_ij()
  437. {
  438.     int i,j,k,n;
  439.     long nrays;
  440.     if (verbose_flag)fprintf(stderr,"From Factor:");
  441.     for (i=n=0;i<number_objects; i++){
  442.         for(j=0;j<object[i].grid_v_reso;j++)
  443.         for(k=0;k<object[i].grid_h_reso;k++,n++){
  444.             hemicube_method(i,j,k,n,&nrays);
  445.             extract_F_i(n);    
  446.             if (verbose_flag)fprintf(stderr,".");
  447.         }
  448.     }
  449.     if (verbose_flag)fprintf(stderr,"\n");
  450.     return(nrays);
  451. }
  452.  
  453. static double get_emission_color(onum,channel)
  454. int onum,channel;
  455. {
  456.     int i;
  457.     for(i=0;i<light_sources;i++)
  458.         if (source[i].oindex == onum)
  459.             return(source[i].strength[channel]);
  460.     return(0.);
  461. }
  462.  
  463. static void solve(channel,color)
  464. int channel;
  465. double color[];
  466. {
  467.     int i,j;
  468.     int m,n;
  469.     int diagonally_dominant();
  470.  
  471.         for (i=m=0; i < number_objects; i++){
  472.         int grid_points = object[i].grid_v_reso*object[i].grid_h_reso;
  473.         double rho = object[i].reflectance[channel];
  474.         double c = get_emission_color(i,channel);
  475.         for (j=0; j < grid_points; j++,m++){
  476.             for(n=0;n<total_surface_grid_points;n++)
  477.                 A[m][n] = -F[m][n]*rho;
  478.             A[m][m] += 1.;
  479.             B[m] = color[m] = c/object[i].grid[j].area;
  480.         }
  481.     }
  482.     if (diagonally_dominant(total_surface_grid_points,A)){
  483.         fprintf(stderr,"Warning : Not Diagonally Dominant.\n");
  484.         fprintf(stderr,
  485.               "Check for the reflectance of the surfaces.\n");
  486.         fprintf(stderr,"Iteration may not converge.\n");
  487.         fprintf(stderr,"So a predefined number of iterations will be tried.\n");
  488.     }
  489.     gauss_iter(total_surface_grid_points,A,color,B,Iter_Threshold);
  490. }
  491.  
  492. static set_intensity(diff_flag,fl)
  493. int diff_flag;
  494. FILE *fl;
  495. {
  496.     int i,j,k,m;
  497.     double max_intensity=0.,max_cum_intensity=0.;
  498.         if((fl != NULL)&& verbose_output_flag){
  499.                 int rows,cols;
  500.                 for (i=m=0; i < number_objects; i++){
  501.                         fprintf(fl,"intensity output for Surface : %d\n",i);
  502.             for(rows=0;rows<object[i].grid_v_reso;rows++)
  503.             for(cols=0;cols<object[i].grid_h_reso;cols++,m++){
  504.                                 Point3 centre;
  505.                                 centre = ofunc[object[i].
  506.                     surface_geometry_type].
  507.                    compute_center_point_in_the_grid_cell
  508.                      (i,rows,cols);
  509.                 fprintf(fl,"%g\t%g\t%g",centre.x,
  510.                     centre.y,centre.z);
  511.                                 for(k=0; k < color_channels; k++)
  512.                                         fprintf(fl,"\t%g", I[k][m]);
  513.                                 fprintf(fl,"\n");
  514.                         }
  515.                 }
  516.         }
  517.     /* Normalize */
  518.     for(k=0;k<color_channels;k++)
  519.     for(i=0;i<total_surface_grid_points;i++)
  520.         if (I[k][i] > max_intensity) max_intensity=I[k][i];    
  521.     for(k=0;k<color_channels;k++)
  522.     for(i=0;i<total_surface_grid_points;i++)
  523.         I[k][i] /= max_intensity;
  524.  
  525.     for (i=m=0; i < number_objects; i++){
  526.         double cum_area=0.,cum_intensity[MAXCHANNELS];
  527.         int grid_points = object[i].grid_v_reso*object[i].grid_h_reso;
  528.         Surface_Grid_structure *g=object[i].grid;
  529.         for(k=0; k<color_channels;k++) cum_intensity[k] = 0.;
  530.         for (j=0; j < grid_points; j++,m++,g++){
  531.             cum_area += g->area;
  532.             for(k=0;k<color_channels;k++){
  533.                 cum_intensity[k] += I[k][m]*g->area;
  534.                 g->normalized_flux_density[k] = ((diff_flag)
  535.                     ? fabs(I[k][m]-g->normalized_flux_density[k])
  536.                     :I[k][m]);
  537.             }
  538.         }
  539.         for(k=0;k<color_channels;k++){
  540.             cum_intensity[k]/=cum_area;
  541.             object[i].normalized_flux_density[k] = cum_intensity[k];
  542.             if (cum_intensity[k] > max_cum_intensity )
  543.                 max_cum_intensity = cum_intensity[k];
  544.         }
  545.     }
  546.     for(i=0;i<number_objects; i++)
  547.         for(k=0;k<color_channels;k++)
  548.             object[i].normalized_flux_density[k]/=max_cum_intensity;
  549. }
  550.  
  551.  
  552. long fullmatrix_radiosity(diff_flag,fl)
  553. int diff_flag;
  554. FILE *fl;
  555. {
  556.     int i;
  557.     long number=0;
  558.     init_fullmatrix_radiosity();
  559.     number=compute_F_ij();
  560.     for(i=0;i<color_channels;i++)
  561.         solve(i,I[i]);
  562.     set_intensity(diff_flag,fl);
  563.     write_out_intensity(fl);
  564.     deinit_fullmatrix_radiosity();
  565.     return(number);
  566. }
  567.  
  568. long progressive_radiosity(maxpasses,diff_flag,fl)
  569. int maxpasses;
  570. int diff_flag;
  571. FILE *fl;
  572. {
  573.     int i,j,k,npasses=0;
  574.     long nrays;
  575.     int max_h_grid,max_v_grid,maxpatch,maxobject= UNDEFINED;
  576.     double maxdelta=0.;
  577.  
  578.     init_progressive_radiosity();
  579.     for(k=0;k<color_channels;k++)
  580.         for(i=0;i<total_surface_grid_points;i++)
  581.             deltaI[k][i]=I[k][i]=0.;
  582.     for(i=0;i<light_sources;i++){
  583.         int index=source[i].oindex,n,nn;
  584.         double cumarea=0.;
  585.         for(j=nn=0;j<object[index].grid_v_reso;j++)
  586.         for(k=0;k<object[index].grid_h_reso;k++,nn++)
  587.             cumarea += object[index].grid[nn].area;
  588.         n=object[index].start_grid_index;
  589.         for(j=nn=0;j<object[index].grid_v_reso;j++)
  590.         for(k=0;k<object[index].grid_h_reso;k++,nn++,n++){
  591.             int l;
  592.             double area=object[index].grid[nn].area;
  593.             for(l=0;l<color_channels;l++){
  594.                I[l][n]=source[i].strength[l]/cumarea;
  595.                deltaI[l][n] = I[l][n]*area;
  596.                if (maxdelta < deltaI[l][n]){
  597.                 maxdelta=deltaI[l][n];
  598.                 max_h_grid=k;
  599.                 max_v_grid=j;
  600.                 maxpatch=n;
  601.                 maxobject=index;
  602.                }
  603.             }
  604.         }
  605.     }
  606.     if (maxpasses < 0) maxpasses=total_surface_grid_points;
  607.     /* Distribute */
  608.     while (maxpatch!= UNDEFINED){
  609.         int l,n,nn;
  610.         double area,increment[MAXCHANNELS];
  611.         hemicube_method(maxobject,
  612.             max_v_grid,max_h_grid,
  613.             maxpatch,&nrays);
  614.         extract_F_i(0);    
  615.         area=object[maxobject].grid[max_v_grid*
  616.                 object[maxobject].grid_h_reso+
  617.                 max_h_grid].area;
  618.         for(l=0;l<color_channels;l++){
  619.            increment[l]=deltaI[l][maxpatch];
  620.            deltaI[l][maxpatch]=0.;
  621.         }
  622.         maxobject= UNDEFINED; maxdelta=0.;
  623.         for(i=n=0;i<number_objects;i++){
  624.             for(j=nn=0;j<object[i].grid_v_reso;j++)
  625.             for(k=0;k<object[i].grid_h_reso;k++,n++,nn++)
  626.             if (object[i].surface_reflection_type !=DIFFUSE)
  627.                 continue;
  628.             else{
  629.                 int l;
  630.                 for(l=0;l<color_channels;l++){
  631.                    deltaI[l][n]+=
  632.                     object[i].reflectance[l]*
  633.                     increment[l]*F[0][n];
  634.                    I[l][n]+=
  635.                     object[i].reflectance[l]*
  636.                     increment[l]*F[0][n]/
  637.                         object[i].grid[nn].area;
  638.                    if(maxdelta < (deltaI[l][n])){
  639.                     maxdelta=deltaI[l][n];
  640.                     max_h_grid=k;
  641.                     max_v_grid=j;
  642.                     maxpatch=n;
  643.                     maxobject=i;
  644.                    }
  645.                 }
  646.             }
  647.         }
  648.         npasses++;
  649.         if (verbose_flag)fprintf(stderr,".");
  650.         if (npasses >= maxpasses) break;
  651.     }
  652.     if (verbose_flag)
  653.         fprintf(stderr,"\n%d Distribution steps.\n",npasses);
  654.     set_intensity(diff_flag,fl);
  655.     write_out_intensity(fl);
  656.     deinit_progressive_radiosity();
  657.     return(nrays);
  658. }
  659.