programs/source/draw_contigalign.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <errno.h>
   5
   6 #include <getopt.h>
   7
   8 #include <math.h>
   9
  10 #include <gd.h>
  11 #include <gdfontl.h>
  12
  13 extern gdFontPtr gdFontSmall;
  14 extern gdFontPtr gdFontTiny;
  15 extern gdFontPtr gdFontLarge;
  16 extern gdFontPtr gdFontMediumBold;
  17
  18 static gdFontPtr localFont;
  19
  20
  21 #define ARRAY_STEP (10)
  22
  23 /* The sizes of source_id and seq_id ought to be extracted from this
  24    structure, and also from any place that implicitly depends on
  25    agreement with these sizes (format string widths, mostly). */
  26 typedef struct {
  27   char source_id[64]; /* Warning: Marty changed this from 32,
  28                          2006-03-24. We don't print the source_id
  29                          verbatim, so all that matters is that the
  30                          structure is big enough to hold inputs (see
  31                          the sscanf below for input parsing). */
  32   char seq_id[32];
  33   char strand;
  34   int  start_loc;
  35   int  end_loc;
  36   int  start_trim;
  37   int  end_trim;
  38   int  highlight;
  39 } contig_align_t;
  40
  41 typedef struct {
  42   int position;
  43   int height;
  44 } height_data_t;
  45
  46 typedef struct height_list {
  47   int position;
  48   int height;
  49   struct height_list *next;
  50 } tmp_height_data_t;
  51
  52 static FILE *pngout, *mapfile;
  53 static char *image_name = NULL;
  54 static int use_thumbnail = 0;
  55
  56 static char *link_basename = "/cgi-bin/SGN/search/tomato_est_search_result.pl?esttigr=";
  57
  58 static void usage(char *argv[], int exit_code) {
  59
  60   fprintf(stderr,"\n  %s: Options\n\n  --link_basename=<url>\n  --imagefile=<filename.png>\n  --mapfile=<filename>\n  --thumbnail\n  --image_name=<string>\n\n  Alignment data is expected on STDIN, as tab delimited text\n\n  EST LABEL,  EST ID,  ASSEMBLY DIRECTION, START, END, START TRIM, END TRIM\n\n  Where EST ID is the ID to be used in the link_basename URL,\n  ASSEMBLY DIRECTION is either + or - for forward or reverse complement\n\n  IF START TRIM and END TIRM are positive, START to START + START TRIM will\n  be considered \"trimmed\" portions of the sequence, like wise with \n  END to END + END TRIM.\n",argv[0]);
  61   exit(-1);
  62 }
  63
  64 static void comline(int argc, char *argv[]) {
  65   static struct option long_options[] = {
  66     { "thumbnail", no_argument, NULL, 't' },
  67     { "imagefile", required_argument, NULL, 'o' },
  68     { "mapfile", required_argument, NULL, 'm' },
  69     { "image_name", required_argument, NULL, 'i' },
  70     { "link_basename", required_argument, NULL, 'b' }
  71   };
  72   const char optstring[] = "to:m:i:l:";
  73   int done;
  74   char *image_filename = NULL;
  75   char *imagemap_filename = NULL;
  76
  77   opterr = 0;
  78   done = 1;
  79   while(done) {
  80     switch(getopt_long(argc, argv, optstring, long_options, NULL)) {
  81     case -1:
  82       done = 0;
  83       break;
  84     case 't':
  85       use_thumbnail = 1;
  86       break;
  87     case 'o':
  88       image_filename = strdup(optarg);
  89       break;
  90     case 'm':
  91       imagemap_filename = strdup(optarg);
  92       break;
  93     case 'i':
  94       image_name = strdup(optarg);
  95       break;
  96     case 'b':
  97       link_basename = strdup(optarg);
  98       break;
  99     default:
 100       fprintf(stderr,"Unknown option %c\n",optopt);
 101       usage(argv, -1);
 102       break;
 103     }
 104   }
 105
 106   if (!image_filename) usage(argv, -1);
 107
 108   if (strcmp(image_filename,"-")==0) {
 109     pngout = stdout;
 110   } else {
 111     pngout = fopen(image_filename,"w");
 112     if (pngout == NULL) {
 113       fprintf(stderr,"Fatal Error: Can not open image output file \"%s\""
 114               "(%s)\n",image_filename, strerror(errno));
 115       usage(argv, -1);
 116     }
 117   }
 118   free(image_filename);
 119
 120   if (!use_thumbnail && imagemap_filename) {
 121     if (strcmp(imagemap_filename,"-")==0) {
 122       if (pngout == stdout) {
 123         fprintf(stderr,"Image output and imagemap output cannot both be"
 124                 " stdout\n");
 125         usage(argv, -1);
 126       } else {
 127         mapfile = stdout;
 128       }
 129     } else {
 130       mapfile = fopen(imagemap_filename,"w");
 131       if (mapfile == NULL) {
 132         fprintf(stderr,"Fatal Error: Can not open imagemap output file \"%s\""
 133                 " (%s)\n", imagemap_filename, strerror(errno));
 134         usage(argv, -1);
 135       }
 136     }
 137     free(imagemap_filename);
 138   } else {
 139     mapfile = NULL;
 140   }
 141
 142   if (image_name == NULL) {
 143     fprintf(stderr,"Warning: Image name not specified\n");
 144     image_name="(Unspecified)";
 145   }
 146
 147 }
 148
 149 static void process_arguments(int argc, char *argv[]) {
 150
 151   if (argc < 4 || argc > 6)  usage(argv,-1);
 152   if (strcmp(argv[1],argv[2])==0) {
 153     fprintf(stderr,"Fatal Error: Mapfile and Imagefile cannot be the same.\n");
 154     usage(argv, -1);
 155   }
 156   if (strcmp(argv[1],"-")==0) {
 157     pngout = stdout;
 158   } else {
 159     pngout = fopen(argv[1], "w");
 160     if (pngout == NULL) {
 161       fprintf(stderr,"Fatal Error: Can not open output file (%s)\n",
 162               strerror(errno));
 163       usage(argv, -1);
 164     }
 165   }
 166
 167   if (argc == 5) {
 168     if (strcmp(argv[4],"thumbnail")==0) {
 169       use_thumbnail = 1;
 170     } else {
 171       fprintf(stderr,"Unrecognized option %s\n",argv[4]);
 172       usage(argv,-1);
 173     }
 174   } else {
 175     use_thumbnail = 0;
 176   }
 177
 178   if (!use_thumbnail && strcmp(argv[2],"-")==0) {
 179     mapfile = stdout;
 180   } else {
 181     mapfile = fopen(argv[2], "w");
 182     if (mapfile == NULL) {
 183       fprintf(stderr,"Fatal Error: Can not open output file (%s)\n",
 184               strerror(errno));
 185       usage(argv,-1);
 186     }
 187   }
 188   strncpy(image_name, argv[3], 40);
 189 }
 190
 191 /* Function passed to qsort() to sort the height data array by
 192    position */
 193 static int height_data_compare(const void *a, const void *b) {
 194
 195   return ((height_data_t *)a)->position - ((height_data_t *)b)->position;
 196 }
 197
 198 /* This reads the array of contig alignment information and computes
 199    the points in the histogram where the hieght changes. Afterwards,
 200    we end up with "height data" which is a list of positions in the
 201    contig with corresponding height of the coverage histrogram. The
 202    histogram height is the same upto but not including the next
 203    height-position pair in the height data array */
 204 static void compute_heights(height_data_t **return_data, int *return_size,
 205                             contig_align_t *align_data, int n_sequences)  {
 206   int i, j, n_heights;
 207   height_data_t *height_data;
 208
 209
 210   n_heights = n_sequences*2;
 211   height_data = calloc(n_heights, sizeof(height_data_t));
 212   for(i=0;i<n_sequences;i++) {
 213     height_data[i].position = align_data[i].start_loc +
 214       align_data[i].start_trim;
 215     /* Add one to the end the sequence, so we have the property that the right
 216        end point is not included in the coverage range for the sequence, thus
 217        we can use this number as an "event" for decreasing the size of the
 218        histogram */
 219     height_data[i + n_sequences].position =
 220       align_data[i].end_loc - align_data[i].end_trim + 1;
 221   }
 222
 223   /* Sort and eliminate duplicates positions. */
 224   qsort(height_data, n_heights, sizeof(height_data_t),height_data_compare);
 225   for(i=1;i<n_heights;i++) {
 226     if (height_data[i].position == height_data[i-1].position) {
 227       int start, stop, difference;
 228       start = i;
 229       stop = start+1;
 230       while(stop<n_heights &&
 231             height_data[stop].position == height_data[start].position)
 232         stop++;
 233       difference = stop - start;
 234       for(;stop<n_heights;stop++,start++)
 235         height_data[start].position = height_data[stop].position;
 236       n_heights -= difference;
 237     }
 238   }
 239   height_data = realloc(height_data, sizeof(height_data_t)*n_heights);
 240
 241   /* For each sequence, increment the hieghts at positions that span the
 242      region covered by the sequence */
 243   for(i=0;i<n_sequences;i++) {
 244     int left, right;
 245     left = align_data[i].start_loc + align_data[i].start_trim;
 246     right = align_data[i].end_loc - align_data[i].end_trim;
 247     for(j=0;j<n_heights;j++) {
 248       if (height_data[j].position >= left && height_data[j].position < right)
 249         height_data[j].height++;
 250     }
 251   }
 252
 253   *return_data = height_data;
 254   *return_size = n_heights;
 255 }
 256
 257
 258 /* This function reads input from the STDIN which is assumed to be
 259    the contig alignment data read from the database. */
 260 static void read_input(contig_align_t **return_array, int *return_size) {
 261   char inputline[256];
 262   contig_align_t *array, *p;
 263   int array_size;
 264   char strand;
 265
 266   array = NULL;
 267   array_size = 0;
 268
 269   /* Put the fgets() at the end of the loop as it will need to be called
 270      on the last (END OF FILE) character in the file before feof(stdin) will
 271      be TRUE. */
 272   fgets(inputline, 256, stdin);
 273   while(!feof(stdin)) {
 274     int items;
 275
 276     if (array_size % ARRAY_STEP == 0) {
 277       array = realloc(array, (array_size + ARRAY_STEP)*sizeof(contig_align_t));
 278       if (array == NULL) {
 279         fprintf(stderr,"FATAL ERROR: Out of memory reading input.\n");
 280         exit(-1);
 281       }
 282     }
 283     p = &array[array_size++];
 284
 285     /* Set defaults to guard against any weird non-matching cases */
 286     p->highlight = 0;
 287     p->start_loc = p->end_loc = p->start_trim = p->end_trim = 0;
 288     p->source_id[0] = 0;
 289
 290     /* Warning: Marty changed the first field width in this sscanf to
 291        60, 2006-03-24. It was formerly 30, which, given that
 292        contig_align_t's source_id was 32, made him fear that Koni was
 293        doing something untoward with the 31st byte of that structure
 294        member, and for voodoo's sake, Marty provided for twice the
 295        untowardness in the doubly-large source_id member. */
 296     items = sscanf(inputline, "%60[^\t]\t%30[^\t]\t%c\t%d\t%d\t%d\t%d\t%d", p->source_id,
 297                    &p->seq_id, &p->strand, &p->start_loc, &p->end_loc,
 298                    &p->start_trim, &p->end_trim, &p->highlight);
 299     if (items != 8) {
 300       fprintf(stdout,"Only matched %d items at (stdin) line %d: Skipping "
 301               "line\n",items,array_size);
 302     }
 303
 304     fgets(inputline, 256, stdin);
 305   }
 306
 307   *return_array = array;
 308   *return_size = array_size;
 309 }
 310
 311 /* This function is passed to qsort() to sort the list of contig
 312    alignment data by starting position relative to the contig. */
 313 static int start_compare(const void *a, const void *b) {
 314   contig_align_t *c;
 315   contig_align_t *d;
 316
 317   c = (contig_align_t *) a;
 318   d = (contig_align_t *) b;
 319
 320   if ((c->start_loc+c->start_trim) != (d->start_loc+d->start_trim))
 321     return (c->start_loc+c->start_trim) - (d->start_loc+d->start_trim);
 322   else
 323     return (c->end_loc-c->end_trim) - (d->end_loc-d->end_trim);
 324 }
 325
 326 /* Find the maximum height of the histogram, round up. For layout in
 327    the image. Histogram must be at least 30 pixels high. */
 328 static int compute_histogram_height(height_data_t *height_data,
 329                                     int n_heights) {
 330   int i;
 331   int max_height;
 332
 333   max_height = 0;
 334   for(i=0;i<n_heights;i++)
 335     if (max_height < height_data[i].height)
 336       max_height = height_data[i].height;
 337
 338   if (max_height < 60) return 60;
 339
 340   max_height += (30 - (max_height % 30));
 341   return max_height;
 342 }
 343
 344 static int compute_ytic(int histogram_height) {
 345
 346   if (histogram_height <= 60) return 15;
 347   else return 20;
 348 }
 349
 350 static int compute_xmax(height_data_t *height_data, int n_heights) {
 351   int last_basepair, xmax;
 352
 353   last_basepair = height_data[n_heights-1].position;
 354   xmax = last_basepair + (100 - (last_basepair % 100));
 355
 356   return xmax;
 357 }
 358
 359 /* For the moment, we are going to hard code the number of x axis tics */
 360 #define N_XTICS (10)
 361 static int compute_xtic(int xmax) {
 362   return xmax / N_XTICS;
 363 }
 364
 365 /* Edges, counterclockwise from the top */
 366 typedef struct {
 367   int north;
 368   int west;
 369   int south;
 370   int east;
 371 } pane_t;
 372 #define BORDER_WIDTH (3)
 373 #define PAD (10)
 374 #define HIST_WIDTH (350)
 375
 376 enum { HIST_PANE = 0, ALIGN_PANE, INFO_PANE, LEGEND_PANE, N_PANES };
 377
 378 static int compute_panes(pane_t panes[], int n_sequences,
 379                          int histogram_height) {
 380   int font_height;
 381   int image_height;
 382
 383   /* Make height of font (for purpose of computing the space needed to write
 384      the sequence alignment and their names, and odd number so that the
 385      mid-point is an integer. */
 386   font_height = localFont->h;
 387   if (!(font_height & 0x1)) font_height += 1;
 388
 389   /* Space needed by the coverage histogram, plus the labeling */
 390   image_height = histogram_height + font_height*3;
 391
 392   /* Space needed by the sequence alignment pane */
 393   image_height += font_height*n_sequences;
 394
 395   /* Room for the bevel'd border */
 396   image_height += BORDER_WIDTH*2;
 397
 398   /* Padding from the bevel to the content area, plus each pane is padded */
 399   image_height += PAD*4;
 400
 401   panes[HIST_PANE].north = BORDER_WIDTH + PAD;
 402   panes[HIST_PANE].west = BORDER_WIDTH + PAD;
 403   panes[HIST_PANE].south = panes[HIST_PANE].north + histogram_height +
 404     font_height*3 + PAD;
 405   /* Account for ytics (up to 4 characters) and axis label, 1 line vertical */
 406   panes[HIST_PANE].east = panes[HIST_PANE].west + HIST_WIDTH +
 407     localFont->w*4 + font_height;
 408
 409   panes[ALIGN_PANE].north = panes[HIST_PANE].south + PAD;
 410   panes[ALIGN_PANE].west = panes[HIST_PANE].west;
 411   /* Allocate space for each sequence, but every 5 sequences we'll draw a rule
 412      so people can align the sequence with the information on the right */
 413   panes[ALIGN_PANE].south = panes[ALIGN_PANE].north + n_sequences*font_height;
 414   panes[ALIGN_PANE].east = panes[HIST_PANE].east;
 415
 416   panes[INFO_PANE].north = panes[ALIGN_PANE].north;
 417   panes[INFO_PANE].west = panes[ALIGN_PANE].east;
 418   panes[INFO_PANE].south = panes[ALIGN_PANE].south;
 419   /* Allow 42 characters maximum in the info PANE */
 420   panes[INFO_PANE].east = panes[INFO_PANE].west + localFont->w*48;
 421
 422   panes[LEGEND_PANE].north = panes[HIST_PANE].north;
 423   panes[LEGEND_PANE].west = panes[HIST_PANE].east;
 424   panes[LEGEND_PANE].east = panes[INFO_PANE].east;
 425   panes[LEGEND_PANE].south = panes[HIST_PANE].south;
 426
 427   return image_height;
 428 }
 429
 430 static void histogram_drawgrid(gdImagePtr im, pane_t pane, int xmax, int xtic,
 431                                int ymax, int ytic) {
 432
 433   float basis;
 434   int black, i;
 435   int x,y;
 436   char label[15];
 437
 438   /* draw rectange */
 439   black = gdImageColorResolve(im, 0, 0, 0);
 440   gdImageRectangle(im, pane.west, pane.north, pane.east,
 441                    pane.south, black);
 442
 443   /* Draw gridlines and tic marks */
 444   basis = (float) xmax/(float) (pane.east - pane.west);
 445   for(i=0;i<=N_XTICS;i++) {
 446     x = pane.west;
 447     x += (int) ((float) ((i*xtic)/basis + 0.5));
 448     if (i!=0 && i!=N_XTICS)
 449       for(y=pane.north;y<pane.south;y+=3)
 450         gdImageSetPixel(im, x, y, black);
 451     gdImageLine(im, x, pane.north, x, pane.north-3, black);
 452     gdImageLine(im, x, pane.south, x, pane.south+1, black);
 453
 454     sprintf(label,"%d",i*xtic);
 455     x = x - (strlen(label)/2.0)*localFont->w;
 456     y = pane.north - 4 - localFont->h;
 457     gdImageString(im, localFont, x, y, label, black);
 458   }
 459   strcpy(label,"position (bp)");
 460   x = pane.west + (pane.east - pane.west)/2 - (strlen(label)/2.0)*localFont->w;
 461   y = pane.north - 5 - localFont->h*2;
 462   gdImageString(im, localFont, x, y, label, black);
 463
 464   for(i=0;i<=ymax;i+=ytic) {
 465     y = i + pane.north;
 466     if (i!=0 && i!=ymax)
 467       for(x=pane.west;x<pane.east;x+=3)
 468         gdImageSetPixel(im, x, y, black);
 469     gdImageLine(im, pane.west, y, pane.west-3, y, black);
 470     gdImageLine(im, pane.east, y, pane.east+3, y, black);
 471
 472     sprintf(label,"%d",i);
 473     x = pane.west - strlen(label)*localFont->w - 4;
 474     y = y - localFont->h/2;
 475     gdImageString(im, localFont, x, y, label, black);
 476   }
 477   x = pane.west - localFont->w*4 - localFont->h - 1;
 478   y = pane.north + (pane.south - pane.north)/2.0 + localFont->w*2.5;
 479   gdImageStringUp(im, localFont, x, y, "Depth", black);
 480
 481 }
 482
 483 static void draw_histogram(gdImagePtr im, pane_t pane,
 484                            height_data_t *height_data, int n_heights,
 485                            int xmax, int xtic, int histogram_height,
 486                            int ytic) {
 487
 488   pane_t grid_pane;
 489   int font_height, i, mediumgray;
 490   float basis;
 491
 492   font_height = localFont->h;
 493   if ((font_height & 0x1) == 0) font_height++;
 494
 495   /* Establish the grid area where we actually draw the histogram. Other
 496      area is for axis label and tic marks */
 497   grid_pane.north = pane.north + font_height * 3;
 498   grid_pane.west = pane.west + localFont->h *1.5 + localFont->w * 4;
 499   grid_pane.east = pane.east - PAD;
 500   grid_pane.south = pane.south - PAD;
 501
 502   histogram_drawgrid(im, grid_pane, xmax, xtic, histogram_height, ytic);
 503   basis = (float) xmax/(grid_pane.east - grid_pane.west);
 504
 505   mediumgray = gdImageColorResolve(im, 0x80, 0x80, 0x80);
 506   for(i=0;i<n_heights-1;i++) {
 507     int left, right, top, bottom;
 508     left = grid_pane.west +
 509       (int) ((float) height_data[i].position/basis + 0.5);
 510     right = grid_pane.west +
 511       (int) ((float) height_data[i+1].position/basis + 0.5);
 512     top = grid_pane.north + 1;
 513     bottom = grid_pane.north + height_data[i].height + 1;
 514
 515     gdImageFilledRectangle(im, left, top, right, bottom, mediumgray);
 516   }
 517
 518 }
 519
 520 static void draw_leftarrow(gdImagePtr im, int x, int y, int color) {
 521
 522   gdImageLine(im,x,y,x-15,y,color);
 523   gdImageLine(im,x-15,y,x-5,y-2,color);
 524   gdImageLine(im,x-15,y,x-5,y+2,color);
 525 }
 526
 527 static void draw_rightarrow(gdImagePtr im, int x, int y, int color) {
 528
 529   gdImageLine(im,x,y,x+15,y,color);
 530   gdImageLine(im,x+15,y,x+5,y-2,color);
 531   gdImageLine(im,x+15,y,x+5,y+2,color);
 532 }
 533
 534 static void draw_dottedline(gdImagePtr im, int x, int y, int color) {
 535   int i;
 536
 537   for(i=0;i<18;i+=6)
 538     gdImageLine(im, x + i, y, x + i + 2, y, color);
 539 }
 540
 541 static void draw_alignments(gdImagePtr im, pane_t pane,
 542                             contig_align_t *align_data, int n_sequences,
 543                             int xmax, int xtic)  {
 544
 545   int i, x, y;
 546   float basis;
 547   int west_edge, east_edge;
 548   int black, red, blue;
 549   int font_height;
 550   int start, covered_start, covered_end, end;
 551   char label[15];
 552
 553   /* First, compute the west and east boundary that corresponds to the
 554      histogram image above our pane */
 555   west_edge = pane.west + localFont->h *1.5 + localFont->w * 4;
 556   east_edge = pane.east - PAD;
 557
 558   font_height = localFont->h;
 559   if ((font_height & 0x1) == 0) font_height++;
 560
 561   basis = (float) xmax/((float) east_edge - west_edge);
 562
 563   black = gdImageColorResolve(im, 0, 0, 0);
 564   red = gdImageColorResolve(im, 0xFF, 0, 0);
 565   blue = gdImageColorResolve(im, 0x33, 0x11, 0xFF);
 566
 567   for(i=0;i<=N_XTICS;i++) {
 568     x = west_edge + ((float) ((i*xtic)/basis + 0.5));
 569     for(y=pane.north;y<pane.south;y+=3) {
 570       gdImageSetPixel(im, x, y, black);
 571     }
 572   }
 573
 574   for(i=0;i<n_sequences;i++) {
 575     y = pane.north + (int) ((float) font_height*(i+0.5) + 0.5);
 576
 577     covered_start = align_data[i].start_loc + align_data[i].start_trim;
 578     covered_start = west_edge + (int) ((float)covered_start/basis + 0.5);
 579     covered_end = align_data[i].end_loc - align_data[i].end_trim;
 580     covered_end = west_edge + (int) ((float)covered_end/basis + 0.5);
 581     if (align_data[i].strand == '+') {
 582       gdImageLine(im, covered_start, y, covered_end, y, black);
 583     }
 584     else {
 585       gdImageLine(im, covered_start, y, covered_end, y, blue);
 586     }
 587
 588     /* draw leading red segment if trimmed */
 589     if (align_data[i].start_trim > 0) {
 590       /* Draw left arrow if untrimmed sequence starts before contig */
 591       if (align_data[i].start_loc < 0) {
 592         //draw_leftarrow(im, west_edge, y, red);
 593         draw_dottedline(im, west_edge - 18, y, red);
 594         sprintf(label,"%dbp",-1*align_data[i].start_loc);
 595         x = west_edge - (strlen(label)+1)*gdFontTiny->w;
 596         //x = west_edge - 18;
 597         gdImageString(im, gdFontTiny, x, y - gdFontTiny->h, label, red);
 598         start = west_edge;
 599       } else {
 600         start = west_edge + (int)((float) align_data[i].start_loc/basis + 0.5);
 601       }
 602       gdImageLine(im, start, y, covered_start-1, y, red);
 603     }
 604
 605     if (align_data[i].end_trim > 0) {
 606       if (align_data[i].end_loc > xmax) {
 607         draw_dottedline(im, east_edge, y, red);
 608         sprintf(label,"%dbp",align_data[i].end_loc - xmax);
 609         x = east_edge + gdFontTiny->w;
 610         gdImageString(im, gdFontTiny, x, y - gdFontTiny->h, label, red);
 611         end = east_edge;
 612       } else {
 613         end = west_edge + (int)((float) align_data[i].end_loc/basis + 0.5);
 614       }
 615       gdImageLine(im, covered_end + 1, y, end, y, red);
 616     }
 617
 618   }
 619 }
 620
 621 /* Allocate an image, draw the bevel border */
 622 static gdImagePtr new_image(int width, int height, int gray) {
 623   gdImagePtr im;
 624   int bgcolor, frontshade, backshade;
 625   int i;
 626
 627   im = gdImageCreate(width+1, height+1);
 628
 629   /* Interlacing allows browsers that support it to show it immedtiately
 630      upon loading the first portion of it, then update the quality as the rest
 631      of the image loads */
 632   //  gdImageInterlace(im, 1);
 633
 634   if (gray<30) gray = 30;
 635   bgcolor = gdImageColorAllocate(im, gray, gray, gray);
 636   frontshade = gdImageColorAllocate(im, gray-35, gray-35, gray-35);
 637   backshade = gdImageColorAllocate(im, gray-75, gray-75, gray-75);
 638
 639   /* Draw lighter bevel on north and west sides */
 640   for(i=0;i<BORDER_WIDTH;i++) {
 641     gdImageLine(im, i, i, width-(i+1), i, frontshade);
 642     gdImageLine(im, i, i, i, height-(i+1), frontshade);
 643   }
 644
 645   for(i=0;i<BORDER_WIDTH;i++) {
 646     gdImageLine(im, i, height-i, width-i, height-i, backshade);
 647     gdImageLine(im, width-i, i, width-i, height-i, backshade);
 648   }
 649
 650   return im;
 651 }
 652
 653 static void draw_seqinfo(gdImagePtr im, pane_t pane,
 654                          contig_align_t *align_data, int n_sequences) {
 655   int font_height;
 656   int black, blue, color;
 657   int west_edge;
 658   int y, i;
 659   char label[80];
 660   int used, total;
 661   int source_id_len;
 662   char short_source_id[32];
 663
 664   font_height = localFont->h;
 665   if ((font_height & 0x1) == 0) font_height++;
 666
 667   west_edge = pane.west + localFont->w*2;
 668
 669   black = gdImageColorResolve(im, 0, 0, 0);
 670   blue = gdImageColorResolve(im, 0x00, 0x00, 0xFF);
 671
 672   if (mapfile) fprintf(mapfile,"<map name=\"contigmap_%s\">", image_name);
 673
 674   for(i=0;i<n_sequences;i++) {
 675     y = pane.north + i * font_height;
 676
 677     total = align_data[i].end_loc - align_data[i].start_loc;
 678     used = total - (align_data[i].start_trim + align_data[i].end_trim);
 679
 680     /* If the source_id is too wide for the display, lop off the right end
 681        and replace it with dotdotdot.  Marty, 2006-03-24 */
 682     source_id_len = strlen(align_data[i].source_id);
 683     if (source_id_len > 29) { /* The width of the image was more or
 684                                  less fixed when I got here.
 685                                  -- Marty. */
 686       strncpy ((char *)short_source_id, align_data[i].source_id, 26);
 687       strcpy ((char *)(short_source_id+26), "... "); /* Mind the space! */
 688     } else
 689       strcpy (short_source_id, align_data[i].source_id);
 690
 691     sprintf(label,"%-30.30s %dbp (%1.0f%%)", short_source_id, total,
 692             used/(float) total * 100);
 693     if (align_data[i].strand == '+') color = black;
 694     else color = blue;
 695
 696     if (align_data[i].highlight) {
 697       /* This used to bold-face here, but it looks silly with the more recent gd
 698          versions -- (using gdFontMediumBold) -- leaving the structure here for
 699          experimentation later with other fonts that don't blow out as much */
 700       gdImageString(im, localFont, west_edge, y, label, color);
 701     } else {
 702       gdImageString(im, localFont, west_edge, y, label, color);
 703     }
 704
 705     if (mapfile)
 706       fprintf(mapfile,"<area coords=\"%d,%d,%d,%d\" "
 707             "href=\"%s%s\">\n",
 708             west_edge,y,west_edge+strlen(label)*localFont->w,y+localFont->h,
 709             link_basename, align_data[i].seq_id);
 710   }
 711   if (mapfile) {
 712     fprintf(mapfile,"</map>\n");
 713     fflush(mapfile);
 714     fclose(mapfile);
 715   }
 716 }
 717
 718 static void draw_legend(gdImagePtr im, pane_t pane, int contig_length,
 719                         int cluster, int contig) {
 720   int line;
 721   int north_edge, west_edge;
 722   int black, red, blue;
 723   char label[80];
 724
 725   north_edge = pane.north + PAD;
 726   west_edge = pane.west + PAD*5;
 727
 728   line = localFont->h;
 729
 730   black = gdImageColorResolve(im, 0, 0, 0);
 731   red = gdImageColorResolve(im, 0xCC, 0, 0);
 732   blue = gdImageColorResolve(im, 0, 0, 0xCC);
 733
 734   sprintf(label,"Alignment Image: %s",image_name);
 735   gdImageString(im, localFont, west_edge, north_edge, label, black);
 736
 737   strcpy(label,"Reverse Complement Strand");
 738   gdImageString(im, localFont, west_edge + 20, north_edge + line*1.5, label, blue);
 739   gdImageLine(im, west_edge, north_edge + line*2, west_edge+15, north_edge + line*2, blue);
 740
 741   strcpy(label,"Given Strand");
 742   gdImageString(im, localFont, west_edge + 20, north_edge + line*2.5, label, black);
 743   gdImageLine(im, west_edge, north_edge + line*3, west_edge + 15, north_edge + line*3, black);
 744
 745   strcpy(label,"Trimmed (non-matching) sequence");
 746   gdImageString(im, localFont, west_edge + 20, north_edge + line*3.5, label, red);
 747   gdImageLine(im, west_edge, north_edge + line*4, west_edge + 15, north_edge + line*4, red);
 748
 749   /* Actual contig length is one less the number passed in here, because that
 750      number is actually the position of the last trimmed EST member, and the
 751      right position is not inclusive by design. See the height building code. */
 752   sprintf(label,"Total Length: %d bp", contig_length-1);
 753   gdImageString(im, localFont, west_edge + 20, north_edge + line*5, label, black);
 754 }
 755
 756
 757 static int gray(int r, int g, int b) {
 758
 759   r |= 0x1F; g |= 0x1f; b |= 0x1f;
 760   if (((r<<1) - g - b) == 0) return 1;
 761
 762   return 0;
 763 }
 764
 765 static gdImagePtr build_thumbnail3(gdImagePtr im, pane_t pane, int w_factor,
 766                                    int h_factor) {
 767   int width, height, tn_width, tn_height;
 768   int x_cells, y_cells, x, y, source_x, source_y;
 769   int basis;
 770   int *average_red, *average_blue, *average_green;
 771   gdImagePtr tn_im;
 772
 773
 774   width = gdImageSX(im);
 775   height = gdImageSY(im);
 776
 777   width -= (width % w_factor);
 778   height -= (height % h_factor);
 779   basis = (w_factor) * (h_factor) * 4;
 780
 781   tn_width = width/w_factor;
 782   tn_height = height/h_factor;
 783
 784   average_red = malloc(tn_width*tn_height*sizeof(int));
 785   average_green = malloc(tn_width*tn_height*sizeof(int));
 786   average_blue = malloc(tn_width*tn_height*sizeof(int));
 787
 788   tn_im = gdImageCreate(tn_width, tn_height);
 789   gdImagePaletteCopy(tn_im, im);
 790
 791   x_cells = tn_width - 1;
 792   y_cells = tn_height - 1;
 793
 794   source_x = source_y = 0;
 795   for(x=0;x<x_cells;x++) {
 796     for(y=0;y<y_cells;y++) {
 797       int red, green, blue, index, i, j, pixel_color;
 798       red = green = blue = 0;
 799       source_x = x*w_factor;
 800       for(i=0;i<w_factor;i++) {
 801         source_x++;
 802         source_y = y*h_factor;
 803         for(j=0;j<h_factor;j++) {
 804           source_y++;
 805           pixel_color = im->pixels[source_y][source_x];
 806           red += gdImageRed(im, pixel_color);
 807           green += gdImageGreen(im, pixel_color);
 808           blue += gdImageBlue(im, pixel_color);
 809         }
 810       }
 811
 812       index = x*y_cells + y;
 813       average_red[index] = (red/basis);
 814       average_green[index] = (green/basis);
 815       average_blue[index] = (blue/basis);
 816     }
 817   }
 818
 819   for(x=0;x<x_cells;x++) {
 820     for(y=0;y<y_cells;y++) {
 821       int red, green, blue, color;
 822       int index;
 823       index = x*y_cells + y;
 824       red = (average_red[index] + average_red[index + y_cells] +
 825         average_red[index+1] + average_green[index + y_cells+1]);
 826       green = (average_green[index] + average_green[index + y_cells]
 827         +  average_green[index+1] + average_green[index + y_cells+1]);
 828       blue = (average_blue[index] + average_blue[index + y_cells] +
 829         average_blue[index+1] + average_blue[index + y_cells+1]);
 830       color = gdImageColorResolve(tn_im, red |0x7, green|0x7, blue|0x7);
 831       tn_im->pixels[y][x] = color;
 832     }
 833   }
 834
 835   free(average_red);
 836   free(average_blue);
 837   free(average_green);
 838   return tn_im;
 839 }
 840
 841 static void build_image(FILE *pngout, contig_align_t *align_data,
 842                         int n_sequences, height_data_t *height_data,
 843                         int n_heights) {
 844   int histogram_height, image_height, ytic;
 845   int xtic, xmax;
 846   int lightgray, i, font_height;
 847   gdImagePtr im;
 848   pane_t panes[N_PANES];
 849
 850   /* Compute parameters for the vertical aspect of the image */
 851   histogram_height = compute_histogram_height(height_data, n_heights);
 852   ytic = compute_ytic(histogram_height);
 853
 854   /* Compute parameters for the horizontal aspect of the image */
 855   xmax = compute_xmax(height_data, n_heights);
 856   xtic = compute_xtic(xmax);
 857
 858   /* This gives us four panes to draw stuff in, like this
 859
 860      ---------------------------------------
 861      |                              |      |
 862      |                              |      |
 863      | Coverage Histrogram          |LEGEND|
 864      |                              |      |
 865      |                              |      |
 866      ---------------------------------------
 867      |                              |      |
 868      |                              |      |
 869      |  Sequence Alignment          | Seq  |
 870      |                              | Info |
 871      |                              |      |
 872      |                              |      |
 873      |                              |      |
 874      --------------------------------------- */
 875   image_height = compute_panes(panes, n_sequences, histogram_height);
 876
 877   im = new_image(panes[LEGEND_PANE].east+PAD+BORDER_WIDTH, image_height, 0xFF);
 878
 879   /* Draw gray'd bars to disinquish sequences in alignment/info pane */
 880   lightgray = gdImageColorResolve(im, 0xDD, 0xDD, 0xDD);
 881   font_height = localFont->h;
 882   if (n_sequences > 6) {
 883     if (! (font_height & 0x1)) font_height ++;
 884     for(i=0;i<n_sequences-3;i+=6) {
 885       int y1, y2;
 886       y1 = panes[ALIGN_PANE].north + i*font_height;
 887       y2 = panes[ALIGN_PANE].north + (i+3)*font_height;
 888       gdImageFilledRectangle(im, panes[ALIGN_PANE].west, y1,
 889                              panes[INFO_PANE].east, y2, lightgray);
 890     }
 891     if (n_sequences % 6 > 0 && n_sequences % 6 <= 3) {
 892       int y1, y2;
 893       y1 = panes[ALIGN_PANE].north + i*font_height;
 894       y2 = panes[ALIGN_PANE].north + (n_sequences)*font_height;
 895       gdImageFilledRectangle(im, panes[ALIGN_PANE].west, y1,
 896                              panes[INFO_PANE].east, y2, lightgray);
 897     }
 898   }
 899   for(i=0;i<n_sequences;i++) {
 900     if (align_data[i].highlight) {
 901       int pink;
 902       int y1,y2;
 903       pink = gdImageColorAllocate(im,0xFF,0x00,0x00);
 904
 905       y1 = panes[ALIGN_PANE].north + i*font_height;
 906       y2 = panes[ALIGN_PANE].north + (i+1)*font_height;
 907       gdImageRectangle(im, panes[ALIGN_PANE].west, y1,
 908                        panes[INFO_PANE].east, y2, pink);
 909     }
 910   }
 911
 912   draw_histogram(im, panes[HIST_PANE], height_data, n_heights, xmax, xtic,
 913                  histogram_height, ytic);
 914
 915   draw_alignments(im, panes[ALIGN_PANE], align_data, n_sequences, xmax, xtic);
 916
 917   draw_seqinfo(im, panes[INFO_PANE], align_data, n_sequences);
 918
 919   draw_legend(im, panes[LEGEND_PANE], height_data[n_heights-1].position, 0, 0);
 920
 921   if (use_thumbnail) {
 922     gdImagePtr im_thumbnail;
 923
 924     im_thumbnail = build_thumbnail3(im, panes[HIST_PANE], 4, 4);
 925     gdImagePng(im_thumbnail, pngout);
 926     gdImageDestroy(im_thumbnail);
 927   } else {
 928     gdImagePng(im, pngout);
 929   }
 930   fflush(pngout);
 931   gdImageDestroy(im);
 932 }
 933
 934 int main(int argc, char *argv[]) {
 935   contig_align_t *align_data;
 936   height_data_t *height_data;
 937   int n_sequences, n_heights;
 938
 939   comline(argc, argv);
 940
 941   localFont = gdFontSmall;
 942
 943   read_input(&align_data, &n_sequences);
 944   if (n_sequences == 0) return -1;
 945
 946   qsort(align_data, n_sequences, sizeof(contig_align_t), start_compare);
 947
 948   compute_heights(&height_data, &n_heights, align_data, n_sequences);
 949
 950   build_image(pngout, align_data, n_sequences, height_data, n_heights);
 951
 952   fclose(pngout);
 953
 954   /*fclose(pngout);*/
 955   return 0;
 956 }
 957