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