bgzf.c

   1 /*
   2  * The Broad Institute
   3  * SOFTWARE COPYRIGHT NOTICE AGREEMENT
   4  * This software and its documentation are copyright 2008 by the
   5  * Broad Institute/Massachusetts Institute of Technology. All rights are reserved.
   6  *
   7  * This software is supplied without any warranty or guaranteed support whatsoever.
   8  * Neither the Broad Institute nor MIT can be responsible for its use, misuse,
   9  * or functionality.
  10  */
  11
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #include <string.h>
  15 #include <unistd.h>
  16 #include <fcntl.h>
  17 #include <sys/types.h>
  18 #include <sys/stat.h>
  19 #include "bgzf.h"
  20
  21 extern off_t ftello(FILE *stream);
  22 extern int fseeko(FILE *stream, off_t offset, int whence);
  23
  24 typedef int8_t byte;
  25
  26 static const int DEFAULT_BLOCK_SIZE = 64 * 1024;
  27 static const int MAX_BLOCK_SIZE = 64 * 1024;
  28
  29 static const int BLOCK_HEADER_LENGTH = 18;
  30 static const int BLOCK_FOOTER_LENGTH = 8;
  31
  32 static const int GZIP_ID1 = 31;
  33 static const int GZIP_ID2 = 139;
  34 static const int CM_DEFLATE = 8;
  35 static const int FLG_FEXTRA = 4;
  36 static const int OS_UNKNOWN = 255;
  37 static const int BGZF_ID1 = 66; // 'B'
  38 static const int BGZF_ID2 = 67; // 'C'
  39 static const int BGZF_LEN = 2;
  40 static const int BGZF_XLEN = 6; // BGZF_LEN+4
  41
  42 static const int GZIP_WINDOW_BITS = -15; // no zlib header
  43 static const int Z_DEFAULT_MEM_LEVEL = 8;
  44
  45
  46 inline
  47 void
  48 packInt16(uint8_t* buffer, uint16_t value)
  49 {
  50     buffer[0] = value;
  51     buffer[1] = value >> 8;
  52 }
  53
  54 inline
  55 int
  56 unpackInt16(const uint8_t* buffer)
  57 {
  58     return (buffer[0] | (buffer[1] << 8));
  59 }
  60
  61 inline
  62 void
  63 packInt32(uint8_t* buffer, uint32_t value)
  64 {
  65     buffer[0] = value;
  66     buffer[1] = value >> 8;
  67     buffer[2] = value >> 16;
  68     buffer[3] = value >> 24;
  69 }
  70
  71 inline
  72 int
  73 min(int x, int y)
  74 {
  75     return (x < y) ? x : y;
  76 }
  77
  78 static
  79 void
  80 report_error(BGZF* fp, const char* message) {
  81     fp->error = message;
  82 }
  83
  84 static
  85 BGZF*
  86 open_read(int fd)
  87 {
  88     FILE* file = fdopen(fd, "r");
  89     BGZF* fp;
  90         if (file == 0) return 0;
  91         fp = malloc(sizeof(BGZF));
  92     fp->file_descriptor = fd;
  93     fp->open_mode = 'r';
  94     fp->owned_file = 0;
  95     fp->file = file;
  96     fp->uncompressed_block_size = MAX_BLOCK_SIZE;
  97     fp->uncompressed_block = malloc(MAX_BLOCK_SIZE);
  98     fp->compressed_block_size = MAX_BLOCK_SIZE;
  99     fp->compressed_block = malloc(MAX_BLOCK_SIZE);
 100     fp->block_address = 0;
 101     fp->block_offset = 0;
 102     fp->block_length = 0;
 103     fp->error = NULL;
 104     return fp;
 105 }
 106
 107 static
 108 BGZF*
 109 open_write(int fd)
 110 {
 111     FILE* file = fdopen(fd, "w");
 112     BGZF* fp;
 113         if (file == 0) return 0;
 114         fp = malloc(sizeof(BGZF));
 115     fp->file_descriptor = fd;
 116     fp->open_mode = 'w';
 117     fp->owned_file = 0;
 118     fp->file = file;
 119     fp->uncompressed_block_size = DEFAULT_BLOCK_SIZE;
 120     fp->uncompressed_block = NULL;
 121     fp->compressed_block_size = MAX_BLOCK_SIZE;
 122     fp->compressed_block = malloc(MAX_BLOCK_SIZE);
 123     fp->block_address = 0;
 124     fp->block_offset = 0;
 125     fp->block_length = 0;
 126     fp->error = NULL;
 127     return fp;
 128 }
 129
 130 BGZF*
 131 bgzf_open(const char* __restrict path, const char* __restrict mode)
 132 {
 133     BGZF* fp = NULL;
 134     if (strcasecmp(mode, "r") == 0) {
 135                 int oflag = O_RDONLY;
 136                 int fd = open(path, oflag);
 137                 if (fd == -1) return 0;
 138         fp = open_read(fd);
 139     } else if (strcasecmp(mode, "w") == 0) {
 140                 int oflag = O_WRONLY | O_CREAT | O_TRUNC;
 141                 int fd = open(path, oflag, 0644);
 142                 if (fd == -1) return 0;
 143         fp = open_write(fd);
 144     }
 145     if (fp != NULL) {
 146         fp->owned_file = 1;
 147     }
 148     return fp;
 149 }
 150
 151 BGZF*
 152 bgzf_fdopen(int fd, const char * __restrict mode)
 153 {
 154         if (fd == -1) return 0;
 155     if (strcasecmp(mode, "r") == 0) {
 156         return open_read(fd);
 157     } else if (strcasecmp(mode, "w") == 0) {
 158         return open_write(fd);
 159     } else {
 160         return NULL;
 161     }
 162 }
 163
 164 static
 165 int
 166 deflate_block(BGZF* fp, int block_length)
 167 {
 168     // Deflate the block in fp->uncompressed_block into fp->compressed_block.
 169     // Also adds an extra field that stores the compressed block length.
 170
 171     byte* buffer = fp->compressed_block;
 172     int buffer_size = fp->compressed_block_size;
 173
 174     // Init gzip header
 175     buffer[0] = GZIP_ID1;
 176     buffer[1] = GZIP_ID2;
 177     buffer[2] = CM_DEFLATE;
 178     buffer[3] = FLG_FEXTRA;
 179     buffer[4] = 0; // mtime
 180     buffer[5] = 0;
 181     buffer[6] = 0;
 182     buffer[7] = 0;
 183     buffer[8] = 0;
 184     buffer[9] = OS_UNKNOWN;
 185     buffer[10] = BGZF_XLEN;
 186     buffer[11] = 0;
 187     buffer[12] = BGZF_ID1;
 188     buffer[13] = BGZF_ID2;
 189     buffer[14] = BGZF_LEN;
 190     buffer[15] = 0;
 191     buffer[16] = 0; // placeholder for block length
 192     buffer[17] = 0;
 193
 194     // loop to retry for blocks that do not compress enough
 195     int input_length = block_length;
 196     int compressed_length = 0;
 197     while (1) {
 198
 199         z_stream zs;
 200         zs.zalloc = NULL;
 201         zs.zfree = NULL;
 202         zs.next_in = fp->uncompressed_block;
 203         zs.avail_in = input_length;
 204         zs.next_out = (void*)&buffer[BLOCK_HEADER_LENGTH];
 205         zs.avail_out = buffer_size - BLOCK_HEADER_LENGTH - BLOCK_FOOTER_LENGTH;
 206
 207         int status = deflateInit2(&zs, Z_DEFAULT_COMPRESSION, Z_DEFLATED,
 208                                   GZIP_WINDOW_BITS, Z_DEFAULT_MEM_LEVEL, Z_DEFAULT_STRATEGY);
 209         if (status != Z_OK) {
 210             report_error(fp, "deflate init failed");
 211             return -1;
 212         }
 213         status = deflate(&zs, Z_FINISH);
 214         if (status != Z_STREAM_END) {
 215             deflateEnd(&zs);
 216             if (status == Z_OK) {
 217                 // Not enough space in buffer.
 218                 // Can happen in the rare case the input doesn't compress enough.
 219                 // Reduce the amount of input until it fits.
 220                 input_length -= 1024;
 221                 if (input_length <= 0) {
 222                     // should never happen
 223                     report_error(fp, "input reduction failed");
 224                     return -1;
 225                 }
 226                 continue;
 227             }
 228             report_error(fp, "deflate failed");
 229             return -1;
 230         }
 231         status = deflateEnd(&zs);
 232         if (status != Z_OK) {
 233             report_error(fp, "deflate end failed");
 234             return -1;
 235         }
 236         compressed_length = zs.total_out;
 237         compressed_length += BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH;
 238         if (compressed_length > MAX_BLOCK_SIZE) {
 239             // should never happen
 240             report_error(fp, "deflate overflow");
 241             return -1;
 242         }
 243         break;
 244     }
 245
 246     packInt16((uint8_t*)&buffer[16], compressed_length-1);
 247     uint32_t crc = crc32(0L, NULL, 0L);
 248     crc = crc32(crc, fp->uncompressed_block, input_length);
 249     packInt32((uint8_t*)&buffer[compressed_length-8], crc);
 250     packInt32((uint8_t*)&buffer[compressed_length-4], input_length);
 251
 252     int remaining = block_length - input_length;
 253     if (remaining > 0) {
 254         if (remaining > input_length) {
 255             // should never happen (check so we can use memcpy)
 256             report_error(fp, "remainder too large");
 257             return -1;
 258         }
 259         memcpy(fp->uncompressed_block,
 260                fp->uncompressed_block + input_length,
 261                remaining);
 262     }
 263     fp->block_offset = remaining;
 264     return compressed_length;
 265 }
 266
 267 static
 268 int
 269 inflate_block(BGZF* fp, int block_length)
 270 {
 271     // Inflate the block in fp->compressed_block into fp->uncompressed_block
 272
 273     z_stream zs;
 274     zs.zalloc = NULL;
 275     zs.zfree = NULL;
 276     zs.next_in = fp->compressed_block + 18;
 277     zs.avail_in = block_length - 16;
 278     zs.next_out = fp->uncompressed_block;
 279     zs.avail_out = fp->uncompressed_block_size;
 280
 281     int status = inflateInit2(&zs, GZIP_WINDOW_BITS);
 282     if (status != Z_OK) {
 283         report_error(fp, "inflate init failed");
 284         return -1;
 285     }
 286     status = inflate(&zs, Z_FINISH);
 287     if (status != Z_STREAM_END) {
 288         inflateEnd(&zs);
 289         report_error(fp, "inflate failed");
 290         return -1;
 291     }
 292     status = inflateEnd(&zs);
 293     if (status != Z_OK) {
 294         report_error(fp, "inflate failed");
 295         return -1;
 296     }
 297     return zs.total_out;
 298 }
 299
 300 static
 301 int
 302 check_header(const byte* header)
 303 {
 304     return (header[0] == GZIP_ID1 &&
 305             header[1] == (byte) GZIP_ID2 &&
 306             header[2] == Z_DEFLATED &&
 307             (header[3] & FLG_FEXTRA) != 0 &&
 308             unpackInt16((uint8_t*)&header[10]) == BGZF_XLEN &&
 309             header[12] == BGZF_ID1 &&
 310             header[13] == BGZF_ID2 &&
 311             unpackInt16((uint8_t*)&header[14]) == BGZF_LEN);
 312 }
 313
 314 static
 315 int
 316 read_block(BGZF* fp)
 317 {
 318     byte header[BLOCK_HEADER_LENGTH];
 319     int64_t block_address = ftello(fp->file);
 320     int count = fread(header, 1, sizeof(header), fp->file);
 321     if (count == 0) {
 322         fp->block_length = 0;
 323         return 0;
 324     }
 325     if (count != sizeof(header)) {
 326         report_error(fp, "read failed");
 327         return -1;
 328     }
 329     if (!check_header(header)) {
 330         report_error(fp, "invalid block header");
 331         return -1;
 332     }
 333     int block_length = unpackInt16((uint8_t*)&header[16]) + 1;
 334     byte* compressed_block = (byte*) fp->compressed_block;
 335     memcpy(compressed_block, header, BLOCK_HEADER_LENGTH);
 336     int remaining = block_length - BLOCK_HEADER_LENGTH;
 337     count = fread(&compressed_block[BLOCK_HEADER_LENGTH], 1, remaining, fp->file);
 338     if (count != remaining) {
 339         report_error(fp, "read failed");
 340         return -1;
 341     }
 342     count = inflate_block(fp, block_length);
 343     if (count < 0) {
 344         return -1;
 345     }
 346     if (fp->block_length != 0) {
 347         // Do not reset offset if this read follows a seek.
 348         fp->block_offset = 0;
 349     }
 350     fp->block_address = block_address;
 351     fp->block_length = count;
 352     return 0;
 353 }
 354
 355 int
 356 bgzf_read(BGZF* fp, void* data, int length)
 357 {
 358     if (length <= 0) {
 359         return 0;
 360     }
 361     if (fp->open_mode != 'r') {
 362         report_error(fp, "file not open for reading");
 363         return -1;
 364     }
 365
 366     int bytes_read = 0;
 367     byte* output = data;
 368     while (bytes_read < length) {
 369         int available = fp->block_length - fp->block_offset;
 370         if (available <= 0) {
 371             if (read_block(fp) != 0) {
 372                 return -1;
 373             }
 374             available = fp->block_length - fp->block_offset;
 375             if (available <= 0) {
 376                 break;
 377             }
 378         }
 379         int copy_length = min(length-bytes_read, available);
 380         byte* buffer = fp->uncompressed_block;
 381         memcpy(output, buffer + fp->block_offset, copy_length);
 382         fp->block_offset += copy_length;
 383         output += copy_length;
 384         bytes_read += copy_length;
 385     }
 386     if (fp->block_offset == fp->block_length) {
 387         fp->block_address = ftello(fp->file);
 388         fp->block_offset = 0;
 389         fp->block_length = 0;
 390     }
 391     return bytes_read;
 392 }
 393
 394 static
 395 int
 396 flush_block(BGZF* fp)
 397 {
 398     while (fp->block_offset > 0) {
 399         int block_length = deflate_block(fp, fp->block_offset);
 400         if (block_length < 0) {
 401             return -1;
 402         }
 403         int count = fwrite(fp->compressed_block, 1, block_length, fp->file);
 404         if (count != block_length) {
 405             report_error(fp, "write failed");
 406             return -1;
 407         }
 408         fp->block_address += block_length;
 409     }
 410     return 0;
 411 }
 412
 413 int
 414 bgzf_write(BGZF* fp, const void* data, int length)
 415 {
 416     if (fp->open_mode != 'w') {
 417         report_error(fp, "file not open for writing");
 418         return -1;
 419     }
 420
 421     if (fp->uncompressed_block == NULL) {
 422         fp->uncompressed_block = malloc(fp->uncompressed_block_size);
 423     }
 424
 425     const byte* input = data;
 426     int block_length = fp->uncompressed_block_size;
 427     int bytes_written = 0;
 428     while (bytes_written < length) {
 429         int copy_length = min(block_length - fp->block_offset, length - bytes_written);
 430         byte* buffer = fp->uncompressed_block;
 431         memcpy(buffer + fp->block_offset, input, copy_length);
 432         fp->block_offset += copy_length;
 433         input += copy_length;
 434         bytes_written += copy_length;
 435         if (fp->block_offset == block_length) {
 436             if (flush_block(fp) != 0) {
 437                 break;
 438             }
 439         }
 440     }
 441     return bytes_written;
 442 }
 443
 444 int
 445 bgzf_close(BGZF* fp)
 446 {
 447     if (fp->open_mode == 'w') {
 448         if (flush_block(fp) != 0) {
 449             return -1;
 450         }
 451         if (fflush(fp->file) != 0) {
 452             report_error(fp, "flush failed");
 453             return -1;
 454         }
 455     }
 456     if (fp->owned_file) {
 457         if (fclose(fp->file) != 0) {
 458             return -1;
 459         }
 460     }
 461     free(fp->uncompressed_block);
 462     free(fp->compressed_block);
 463     free(fp);
 464     return 0;
 465 }
 466
 467 int64_t
 468 bgzf_tell(BGZF* fp)
 469 {
 470     return ((fp->block_address << 16) | (fp->block_offset & 0xFFFF));
 471 }
 472
 473 int64_t
 474 bgzf_seek(BGZF* fp, int64_t pos, int where)
 475 {
 476     if (fp->open_mode != 'r') {
 477         report_error(fp, "file not open for read");
 478         return -1;
 479     }
 480     if (where != SEEK_SET) {
 481         report_error(fp, "unimplemented seek option");
 482         return -1;
 483     }
 484     int block_offset = pos & 0xFFFF;
 485     int64_t block_address = (pos >> 16) & 0xFFFFFFFFFFFFLL;
 486     if (fseeko(fp->file, block_address, SEEK_SET) != 0) {
 487         report_error(fp, "seek failed");
 488         return -1;
 489     }
 490     fp->block_length = 0;  // indicates current block is not loaded
 491     fp->block_address = block_address;
 492     fp->block_offset = block_offset;
 493     return 0;
 494 }
 495