1 #include <omp.h>
    2 
    3 /*************************************************************************/
    4 /*                                                                       */
    5 /*  Copyright (c) 1994 Stanford University                               */
    6 /*                                                                       */
    7 /*  All rights reserved.                                                 */
    8 /*                                                                       */
    9 /*  Permission is given to use, copy, and modify this software for any   */
   10 /*  non-commercial purpose as long as this copyright notice is not       */
   11 /*  removed.  All other uses, including redistribution in whole or in    */
   12 /*  part, are forbidden without prior written permission.                */
   13 /*                                                                       */
   14 /*  This software is provided with absolutely no warranty and no         */
   15 /*  support.                                                             */
   16 /*                                                                       */
   17 /*************************************************************************/
   18 
   19 /*************************************************************************
   20 *                                                                        *
   21 *     main.c:  Starting point for rendering system.                      *
   22 *                                                                        *
   23       Usage:  VOLREND num_processes input_file [-a]
   24 
   25       where input_file is head for the head data set. i.e. the filename
   26           without a suffix.
   27       and the -a option enables adaptive sampling of pixels.
   28 
   29 *************************************************************************/
   30 
   31 #include "incl.h"
   32 #include <sys/time.h>
   33 #include <sys/resource.h>
   34 #include "tiffio.h"
   35 
   36 #define SH_MEM_AMT 60000000
   37 
   38 
   39 
   40 #include "anl.h"
   41 
   42 struct GlobalMemory *Global;
   43 
   44 int image_section[NI];
   45 int voxel_section[NM];
   46 
   47 int num_nodes,frame;
   48 int num_blocks,num_xblocks,num_yblocks;
   49 PIXEL *image_address;
   50 MPIXEL *mask_image_address;
   51 PIXEL *image_block,*mask_image_block;
   52 PIXEL *shd_address;
   53 BOOLEAN *sbit_address;
   54 long shd_length;
   55 short image_len[NI], mask_image_len[NI];
   56 long image_length, mask_image_length;
   57 //char *filename[FILENAME_STRING_SIZE];
   58 char filename[FILENAME_STRING_SIZE];
   59 void Render_Loop();
   60 
   61 mclock(stoptime,starttime,exectime)
   62      unsigned int stoptime,starttime,*exectime;
   63 {
   64   if (stoptime < starttime)
   65     *exectime = ((MAX_INT - starttime) + stoptime)/1000;
   66   else
   67     *exectime = (stoptime - starttime)/1000;
   68 }
   69 
   70 
   71 int WriteRawText(char *filename, int width, int height, int scanbytes, unsigned char *data);
   72 
   73 
   74 
   75 main(argc, argv)
   76      int argc;
   77      char *argv[];
   78 {
   79   if ((argc < 3) || (strncmp(argv[1],"-h",strlen("-h")) == 0) || (strncmp(argv[1],"-h",strlen("-H")) == 0)){
   80     printf("usage:  VOLREND num_processes input_file\n");
   81     exit(-1);
   82   }
   83 
   84   {;};
   85 
   86   num_nodes = atoi(argv[1]);
   87 
   88   strcpy(filename,argv[2]);
   89 
   90   if (argc == 4) {
   91     if (strncmp(argv[3],"-a",strlen("-a")) == 0)
   92       adaptive = YES;
   93     else {
   94       printf("usage:  VOLREND num_processes input_file [-a] \n");
   95       exit(-1);
   96     }
   97   }
   98 
   99   Frame();
  100 
  101   if (num_nodes > 1)
  102     {;};
  103   {exit(0);};
  104 }
  105 
  106 
  107 Frame()
  108 {
  109   unsigned int starttime,stoptime,exectime,i;
  110 
  111   Init_Options();
  112 
  113   printf("*****Entering init_decomposition with num_nodes = %d\n",num_nodes);
  114   fflush(stdout);
  115 
  116   Init_Decomposition();
  117 
  118   printf("*****Exited init_decomposition with num_nodes = %d\n",num_nodes);
  119   fflush(stdout);
  120 
  121 
  122 
  123   Global = (struct GlobalMemory *)malloc(sizeof(struct GlobalMemory));;
  124   {;};
  125   {;};
  126   {;};
  127   {;};
  128   {;};
  129 
  130   /* load dataset from file to each node */
  131 #ifndef RENDER_ONLY
  132   {long time(); (starttime) = time(0);};
  133   Load_Map(filename);
  134   {long time(); (stoptime) = time(0);};
  135   mclock(stoptime,starttime,&exectime);
  136   printf("wall clock execution time to load map:  %u ms\n", exectime);
  137 #endif
  138 
  139   {long time(); (starttime) = time(0);};
  140 #ifndef RENDER_ONLY
  141   Compute_Normal();
  142 #ifdef PREPROCESS
  143   Store_Normal(filename);
  144 #endif
  145 #else
  146   Load_Normal(filename);
  147 #endif
  148   {long time(); (stoptime) = time(0);};
  149   mclock(stoptime,starttime,&exectime);
  150   printf("wall clock execution time to compute normal:  %u ms\n", exectime);
  151 
  152   {long time(); (starttime) = time(0);};
  153 #ifndef RENDER_ONLY
  154   Compute_Opacity();
  155 #ifdef PREPROCESS
  156   Store_Opacity(filename);
  157 #endif
  158 #else
  159   Load_Opacity(filename);
  160 #endif
  161   {long time(); (stoptime) = time(0);};
  162   mclock(stoptime,starttime,&exectime);
  163   printf("wall clock execution time to compute opacity:  %u ms\n", exectime);
  164 
  165   Compute_Pre_View();
  166   shd_length = LOOKUP_SIZE;
  167   Allocate_Shading_Table(&shd_address,&sbit_address,shd_length);
  168   /* allocate space for image */
  169   image_len[X] = frust_len;
  170   image_len[Y] = frust_len;
  171   image_length = image_len[X] * image_len[Y];
  172   Allocate_Image(&image_address,image_length);
  173 
  174   if (num_nodes == 1) {
  175     block_xlen = image_len[X];
  176     block_ylen = image_len[Y];
  177     num_blocks = 1;
  178     num_xblocks = 1;
  179     num_yblocks = 1;
  180     image_block = image_address;
  181   }
  182   else {
  183     num_xblocks = ROUNDUP((float)image_len[X]/(float)block_xlen);
  184     num_yblocks = ROUNDUP((float)image_len[Y]/(float)block_ylen);
  185     num_blocks = num_xblocks * num_yblocks;
  186     Lallocate_Image(&image_block,block_xlen*block_ylen);
  187   }
  188 
  189   {long time(); (starttime) = time(0);};
  190 #ifndef RENDER_ONLY
  191   Compute_Octree();
  192 #ifdef PREPROCESS
  193   Store_Octree(filename);
  194 #endif
  195 #else
  196   Load_Octree(filename);
  197 #endif
  198   {long time(); (stoptime) = time(0);};
  199   mclock(stoptime,starttime,&exectime);
  200   printf("wall clock execution time to compute octree:  %u ms\n", exectime);
  201 
  202 #ifdef PREPROCESS
  203   return;
  204 #endif
  205 
  206   if (adaptive) {
  207     for (i=0; i<NI; i++) {
  208       mask_image_len[i] = image_len[i];
  209     }
  210     mask_image_length = image_length;
  211     Allocate_MImage(&mask_image_address, mask_image_length);
  212     if (num_nodes == 1)
  213       mask_image_block = (PIXEL *)mask_image_address;
  214     else
  215       Lallocate_Image(&mask_image_block, block_xlen*block_ylen);
  216   }
  217 
  218 #ifndef RENDER_ONLY
  219   Deallocate_Map(&map_address);
  220 #endif
  221 
  222   Global->Index = NODE0;
  223 
  224   printf("\nRendering...\n");
  225   printf("node\tframe\ttime\titime\trays\thrays\tsamples trilirped\n");
  226 
  227   for (i=1; i<num_nodes; i++) {fprintf(stderr, "No more processors -- this is a uniprocessor version!\n"); exit(-1);}
  228 
  229   Render_Loop();
  230 }
  231 
  232 
  233 void Render_Loop()
  234 {
  235   int step,i,dim,pid;
  236   PIXEL *local_image_address;
  237   MPIXEL *local_mask_image_address;
  238   PIXEL *local_shd_address;
  239   char outfile[FILENAME_STRING_SIZE],cmd[FILENAME_STRING_SIZE];
  240   int image_partition,mask_image_partition,shd_table_partition;
  241   float inv_num_nodes;
  242   int my_node;
  243   int zz;
  244   int zzz;
  245   PIXEL *p;
  246 
  247   {;};
  248   my_node = Global->Index++;
  249   {;};
  250   my_node = my_node%num_nodes;
  251 
  252 /*  POSSIBLE ENHANCEMENT:  Here's where one might bind the process to a
  253     processor, if one wanted to.
  254 */
  255 
  256   inv_num_nodes = 1.0/(float)num_nodes;
  257   image_partition = ROUNDUP(image_length*inv_num_nodes);
  258   mask_image_partition = ROUNDUP(mask_image_length*inv_num_nodes);
  259 
  260 #ifdef DIM
  261   for (dim=0; dim<NM; dim++) {
  262 #endif
  263 
  264    #pragma omp parallel for private (i, local_mask_image_address, local_image_address) schedule(static)
  265     for (step=0; step<ROTATE_STEPS; step++) { /* do rotation sequence */
  266 
  267 /*  POSSIBLE ENHANCEMENT:  Here is where one might reset statistics, if
  268     		one wanted to.
  269 */
  270 
  271       frame = step;
  272       /* initialize images here */
  273       local_image_address = image_address + image_partition * my_node;
  274       local_mask_image_address = mask_image_address +
  275 	mask_image_partition * my_node;
  276 
  277       {;};
  278 
  279       if (my_node == num_nodes-1) {
  280 	for (i=image_partition*my_node; i<image_length; i++)
  281 	  *local_image_address++ = background;
  282 	if (adaptive)
  283 	  for (i=mask_image_partition*my_node; i<mask_image_length; i++)
  284 	    *local_mask_image_address++ = NULL_PIXEL;
  285       }
  286       else {
  287 	for (i=0; i<image_partition; i++)
  288 	  *local_image_address++ = background;
  289 	if (adaptive)
  290 	  for (i=0; i<mask_image_partition; i++)
  291 	    *local_mask_image_address++ = NULL_PIXEL;
  292       }
  293 
  294       if (my_node == ROOT) {
  295 #ifdef DIM
  296 	Select_View((float)STEP_SIZE, dim);
  297 #else
  298         Select_View((float)STEP_SIZE, Y);
  299 #endif
  300 }
  301 
  302       {;};
  303 
  304       Global->Counter = num_nodes;
  305       Global->Queue[num_nodes][0] = num_nodes;
  306       Global->Queue[my_node][0] = 0;
  307 
  308       Render(my_node);
  309 
  310       if (my_node == ROOT) {
  311 
  312          int slashpos = -1; int i = 0;
  313           char * fname;
  314           for (; i < strlen(filename); ++i) {
  315              if (filename[i] == '/') { slashpos = i; }
  316           }
  317 
  318           fname = &filename[slashpos + 1];
  319 
  320 
  321          if (ROTATE_STEPS > 1) {
  322 
  323 
  324 
  325 #ifdef DIM
  326             sprintf(outfile, "%s_%d",fname, 1000+dim*ROTATE_STEPS+step);
  327 #else
  328             sprintf(outfile, "%s_%d.tiff", fname, 1000 + step);
  329 #endif
  330             /*	  Store_Image(outfile);
  331              p = image_address;
  332              for (zz = 0;zz < image_length;zz++) {
  333              tiff_image[zz] = (int) ((*p)*256*256*256 + (*p)*256*256 +
  334              (*p)*256 + (*p));
  335              p++;
  336              }
  337              tiff_save_rgba(outfile,tiff_image,image_len[X],image_len[Y]);  */
  338             WriteGrayscaleTIFF(outfile, image_len[X], image_len[Y],
  339                   image_len[X], image_address);
  340          } else {
  341             /*	  Store_Image(filename);
  342              p = image_address;
  343              for (zz = 0;zz < image_length;zz++) {
  344              tiff_image[zz] = (int) ((*p)*256*256*256 + (*p)*256*256 +
  345              (*p)*256 + (*p));
  346              p++;
  347              }
  348              tiff_save_rgba(filename,tiff_image,image_len[X],image_len[Y]);    */
  349             /*          strcat(filename,".tiff");
  350              WriteGrayscaleTIFF(filename, image_len[X],image_len[Y],image_len[X], image_address);
  351              */
  352 
  353             {
  354 
  355                strcat(filename, ".tiff");
  356                WriteGrayscaleTIFF(fname, image_len[X], image_len[Y],
  357                      image_len[X], image_address);
  358             }
  359 
  360 
  361 
  362          }
  363       }
  364     }
  365 
  366     WriteRawText("result.txt", image_len[X], image_len[Y],
  367           image_len[X], image_address);
  368 
  369 
  370 #ifdef DIM
  371   }
  372 #endif
  373 }
  374 
  375 
  376 int WriteRawText(char *filename, int width, int height, int scanbytes, unsigned char *data) {
  377    FILE * out = fopen(filename, "w");
  378    int h,w;
  379 
  380    for (h = 0; h < height; h++) {
  381       for (w = 0; w < width; w++) {
  382          int pixel = data[h*height + w];
  383          fprintf(out, "%d ", pixel);
  384       }
  385       fprintf(out, "\n");
  386    }
  387 
  388    fclose(out);
  389 
  390 }
  391 
  392 
  393 
  394 Error(string,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8)
  395      char string[], *arg1, *arg2, *arg3, *arg4, *arg5, *arg6, *arg7, *arg8;
  396 {
  397   fprintf(stderr,string,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8);
  398   exit(1);
  399 }
  400 
  401 
  402 Allocate_Image(address, length)
  403      PIXEL **address;
  404      long length;
  405 {
  406   unsigned int i,j,size,type_per_page,count,block;
  407   unsigned int p,numbytes;
  408 
  409   printf("    Allocating image of %ld bytes...\n",
  410 	 length*sizeof(PIXEL));
  411 
  412   *address = (PIXEL *)malloc(length*sizeof(PIXEL));;
  413 
  414   if (*address == NULL)
  415     Error("    No space available for image.\n");
  416 
  417   for (i=0; i<length; i++) *(*address+i) = 0;
  418 
  419 }
  420 
  421 
  422 Allocate_MImage(address, length)
  423      MPIXEL **address;
  424      long length;
  425 {
  426   unsigned int i,j,size,type_per_page,count,block;
  427   unsigned int p,numbytes;
  428 
  429   printf("    Allocating image of %ld bytes...\n",
  430 	 length*sizeof(MPIXEL));
  431 
  432   *address = (MPIXEL *)malloc(length*sizeof(MPIXEL));;
  433 
  434   if (*address == NULL)
  435     Error("    No space available for image.\n");
  436 
  437 /*  POSSIBLE ENHANCEMENT:  Here's where one might distribute the
  438     opacity map among physical memories if one wanted to.
  439 */
  440 
  441   for (i=0; i<length; i++) *(*address+i) = 0;
  442 
  443 }
  444 
  445 
  446 Lallocate_Image(address, length)
  447      PIXEL **address;
  448      long length;
  449 {
  450   /*  char *calloc();*/
  451   int i;
  452   unsigned int p,numbytes;
  453 
  454   printf("    Allocating image of %ld bytes...\n",
  455 	 length*sizeof(PIXEL));
  456   *address = (PIXEL *)calloc(length,sizeof(PIXEL));
  457   if (*address == NULL)
  458     Error("    No space available for image.\n");
  459 
  460 }
  461 
  462 
  463 Store_Image(filename)
  464      char filename[];
  465 {
  466   char local_filename[FILENAME_STRING_SIZE];
  467   int fd;
  468   short pix_version;
  469   short local_image_len[NI+1]; /* dimension larger than NI for backwards     */
  470                                /* compatibility of .pix file with no color   */
  471 
  472   local_image_len[X] = image_len[X];
  473   local_image_len[Y] = image_len[Y];
  474   local_image_len[2] = 1;
  475 
  476   pix_version = PIX_CUR_VERSION;
  477   strcpy(local_filename,filename);
  478   strcat(local_filename,".pix");
  479   fd = Create_File(local_filename);
  480   Write_Shorts(fd,&pix_version,(long)sizeof(pix_version));
  481 
  482   Write_Shorts(fd,local_image_len,(long)sizeof(local_image_len));
  483   Write_Longs(fd,&image_length,(long)sizeof(image_length));
  484 
  485   Write_Bytes(fd,image_address,(long)(image_length*sizeof(PIXEL)));
  486   Close_File(fd);
  487 }
  488 
  489 
  490 Allocate_Shading_Table(address1,address2,length)
  491      PIXEL **address1,**address2;
  492      long length;
  493 {
  494   unsigned int i,size,type_per_page;
  495   unsigned int p,numbytes;
  496 
  497   printf("    Allocating shade lookup table of %ld bytes...\n",
  498 	 length*sizeof(PIXEL));
  499 
  500 /*  POSSIBLE ENHANCEMENT:  If you want to replicate the shade table,
  501     replace the macro with a simple malloc in the line below */
  502 
  503   *address1 = (PIXEL *)malloc(length);;
  504 
  505   if (*address1 == NULL)
  506     Error("    No space available for table.\n");
  507 
  508 /*  POSSIBLE ENHANCEMENT:  Here's where one might distribute the
  509     shading table among physical memories if one wanted to.
  510 */
  511 
  512   for (i=0; i<length; i++) *(*address1+i) = 0;
  513 
  514 }
  515 
  516 
  517 Init_Decomposition()
  518 {
  519   int factors[MAX_NUMPROC];
  520   double processors,newfactor;
  521   int i,sq,cu,maxcu,count;
  522 
  523   /* figure out what to divide dimensions of image and volume by to */
  524   /* partition data and computation to processors                   */
  525   if (num_nodes == 1) {
  526     image_section[X] = 1;
  527     image_section[Y] = 1;
  528     voxel_section[X] = 1;
  529     voxel_section[Y] = 1;
  530     voxel_section[Z] = 1;
  531   }
  532   else {
  533     count = 1;
  534     processors = (double)num_nodes;
  535     sq = (int)sqrt(processors);
  536     cu = (int)pow(processors,1.0/3.0);
  537     factors[0] = 1;
  538 
  539     for (i=2; i<sq; i++) {
  540       if (FRACT(processors/(double)i) == 0.0) {
  541         factors[count++] = i;
  542         if (i <= cu) {
  543           maxcu = i;
  544           newfactor = (double)(num_nodes/i);
  545         }
  546       }
  547     }
  548     count--;
  549     image_section[X] = factors[count];
  550     image_section[Y] = num_nodes/factors[count];
  551 
  552     sq = (int)sqrt(newfactor);
  553     count = 1;
  554 
  555     for (i=2; i<sq; i++) {
  556       if (FRACT(newfactor/(double)i) == 0.0)
  557         factors[count++] = i;
  558     }
  559     count--;
  560     voxel_section[X] = MIN(factors[count],maxcu);
  561     voxel_section[Y] = MAX(factors[count],maxcu);
  562     voxel_section[Z] = (int)newfactor/factors[count];
  563   }
  564 }
  565 
  566 /*
  567  * WriteGrayscaleTIFF
  568  *
  569  * Create a grayscale TIFF image.  The input is a sequence of bytes in an
  570  * array (each byte representing one pixel).  This is converted into a
  571  * compressed TIFF image file.
  572  *
  573  * Return value is 1 for success, 0 for failure.
  574  */
  575 
  576 int
  577 WriteGrayscaleTIFF(filename, width, height, scanbytes, data)
  578 char *filename;         /* file to write to */
  579 int width, height;      /* size of image */
  580 int scanbytes;          /* length of each scanline in bytes */
  581 unsigned char *data;    /* image data */
  582 {
  583     unsigned long y;
  584     double factor;
  585     int c;
  586     unsigned short cmap[256];           /* output color map */
  587     TIFF *outimage;                     /* TIFF image handle */
  588 
  589     /* create a grayscale ramp for the output color map */
  590     factor = (double)((1 << 16) - 1) / (double)((1 << 8) - 1);
  591     for (c = 0; c < 256; c++)
  592         cmap[c] = (int)(c * factor);
  593 
  594     /* open and initialize output file */
  595     if ((outimage = TIFFOpen(filename, "w")) == NULL)
  596         return(0);
  597     TIFFSetField(outimage, TIFFTAG_IMAGEWIDTH, width);
  598     TIFFSetField(outimage, TIFFTAG_IMAGELENGTH, height);
  599     TIFFSetField(outimage, TIFFTAG_BITSPERSAMPLE, 8);
  600     TIFFSetField(outimage, TIFFTAG_SAMPLESPERPIXEL, 1);
  601     TIFFSetField(outimage, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
  602     TIFFSetField(outimage, TIFFTAG_COMPRESSION, COMPRESSION_LZW);
  603     TIFFSetField(outimage, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT);
  604     TIFFSetField(outimage, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_PALETTE);
  605     TIFFSetField(outimage, TIFFTAG_COLORMAP, cmap, cmap, cmap);
  606 
  607     /* write the image data */
  608     for (y = 0; y < height; y++) {
  609         if (!TIFFWriteScanline(outimage, data, y, 0)) {
  610             TIFFClose(outimage);
  611             return(0);
  612         }
  613         data += scanbytes;
  614     }
  615 
  616     /* close the file */
  617     TIFFClose(outimage);
  618     return(1);
  619 }
  620 
  621 
  622