Об особенностях хранения 16 бит изображений в PNG формате

в 13:10, , рубрики: Без рубрики

Как обычно, в самом начале я постараюсь сэкономить время читателей — суть в том, что в PNG формате есть тег sBIT, который описывает "битность" (точнее, количество значащих бит) изображения, и хотя этот тег присутствует в libpng, которая является практически "стандартом де факто" для работы с PNG, далеко не все продукты, использующие эту библиотеку, учитывают его наличие, и в результате это может привести к неожиданному изменению интенсивности пикселей (хотя и не вызывает потерю данных). Ну а кому интересно более детально погрузиться, в том числе в приватные теги, могут продолжить чтение. Материал достаточно специфический, но возможно, кому-нибудь поможет сохранить немного времени в похожей ситуации. Я упомяну пару коммерческих продуктов типа LabVIEW и VDM, но рекламной нагрузки этот чисто технический пост не несёт.

Я работаю в основном с рентгеновскими изображениями, которые формально являются шестнадцатибитными. "Формально" в нашем контексте означает, что они могут быть и меньшей разрядности, например, 12 или 14 бит, но PNG в любом случае использует два байта для хранения одного пиксела (возможность хранения 16 бит — это часть стандарта PNG), при этом битовой упаковки там не происходит. Мы храним такие изображения в 16 бит PNG файлах для своих нужд. Почему именно PNG? Дело в том, что этот формат со сжатием без потери качества также позволяет хранить любые пользовательские данные, и мы пользуемся этим для хранения некоторой служебной информации прямо в файле изображения (например регионы интереса для сегментации, оверлеи, параметры рентгена и т.д.). Так можно делать не только с PNG, но библиотека, которой мы в основном пользуемся (NI Vision Development Toolkit), предлагает это "из коробки", и это очень удобно, значительно удобнее, чем писать в DICOM/DICONDE формат, общепринятый в радиографии. Оффтопик — одно время я пользовался известным онлайн заметочником, и в какой-то момент они сделали опцию хранения приаттаченных pdf платной опцией, так что я просто добавлял pdf как полезную нагрузку к PNG и хранил их как картинки. Это не стеганография, конечно, но под видом безобидной PNG может скрываться "волк в овечьей шкуре".

А началось всё с того, что в отдел техподдержки прилетело сообщение о том, что внутренние PNG изображения, которые хранят наши рентгеновские снимки, при открытии достигают небывалых величин интенсивности. Детектор там использовался 14 бит, так что больше чем 16383 градаций в теории быть не может, но заказчик утверждал, что там получается 60000 и выше. Я попросил прислать проблемную картинку, открыл её в Vision Assistant (довольно удобный инструмент, идёт в комплекте), не увидел вообще никаких аномалий, и написал, что не надо разгонять гистограмму на весь диапазон, они явно зачем-то домножают интенсивности на какой-нибудь множитель, пусть проверяют свой код, файл в порядке и задал встречный вопрос, чем они изображения открывают, на что заказчик написал о том, что они используют ImageJ. На сей раз я решил копнуть поглубже, вечерком пошёл в цех, засунул первый попавшийся под руку объект в систему, сделал "флюорографию", расчехлил Fiji (это клон ImageJ), открыл там картинку и вот оно, 912...12160 превратилось в 3648...48642:

Об особенностях хранения 16 бит изображений в PNG формате - 1

Чтобы не смущать вас кодом на LabVIEW — эта библиотека суть набор обычных DLL, все функции экспортированы, так что я могу воспользоваться и чистым Си, просто мне в LabVIEW быстрее код мышкой набрасывать:

#include <nivision.h>
int main (int argc, char *argv[])
{
	Image* image = imaqCreateImage (IMAQ_IMAGE_U16, 0); //создаём 16 бит картинку
	imaqReadFile2 (image, "Mouse.png", 0, NULL, NULL); //читаем файл
	HistogramReport* histo = imaqHistogram (image, 65536, 0, 65535, NULL);
	printf("Min = %f; Max = %fn", histo->min, histo->max);
	imaqDispose(histo);
	imaqDispose(image);
	return 0;
}

Как забавно — я больше двадцати лет работаю с этой библиотекой, но так и не удосужился открыть изображения в сторонних просмотрщиках. Вечер переставал быть томным.

Итак, в моей картинке интенсивности укладываются в 912...12160, а ImageJ показывает 3648...48642. И если 3648 ещё можно получить домножением 912 на 4, то вот как 12160 превращается в 48642 — ещё предстоит объяснить, поскольку тут ещё +2 градации, немного похоже на мистику, но любое наблюдаемое и воспроизводимое поведение обычно имеет рациональное объяснение.

Прежде всего, надо было понять — это только ImageJ так чудит, или проблема глубже. Нет проблем, навскидку Питон наше всё:

import cv2
image = cv2.imread('Mouse.png', cv2.IMREAD_UNCHANGED)
min_val,max_val, min_indx, max_indx=cv2.minMaxLoc(image)
print(min_val, max_val)

Кто бы сомневался:

>python Mouse.py
3648.0 48642.0

Поскольку тут OpenCV используется, то на Си с OpenCV я даже проверять не стал.

А что нам скажет ImageMagick 7.1.1?

magick.exe identify -verbose Mouse.png

И вот, здесь тоже 3648...48642:

Image:
  Filename: Mouse.png
  Format: PNG (Portable Network Graphics)
  Mime type: image/png
  Class: DirectClass
  Geometry: 1024x1024+0+0
  Colorspace: Gray
  Type: Grayscale
  Depth: 16-bit
  Channels: 1.0
  Channel depth:
    Gray: 16-bit
  Channel statistics:
    Pixels: 1048576
    Gray:
      min: 3648  (0.0556649)
      max: 48642 (0.742229)

Ситуация начинает напоминать известный анекдот, когда жена звонит мужу:

Дорогой, будь осторожен на дороге, тут по радио передали, что какой-то идиот едет по автобану против движения...

Один?! Да тут ВСЕ едут против движения!

Тут в самый раз сделать простую тестовую картинку, я брошу в неё пару пикселей — 1000 и 10000 и сохраню, для меня самый простой способ вот так, на LabVIEW. Один из возможных претендентов, вызывающих такое непонятное поведение — сжатие (и фильтр), так что я его выключу для чистоты эксперимента (это делается установкой Image Quality (compression) в 1000):

Об особенностях хранения 16 бит изображений в PNG формате - 2

На Си, что б понятнее:

#include <nivision.h>
int main (int argc, char *argv[])
{
	const unsigned short Pixels[2] = {1000, 10000};
	Image* image = imaqCreateImage (IMAQ_IMAGE_U16, 0);
	imaqArrayToImage(image, Pixels, 2, 1);
	imaqWritePNGFile2 (image, "Test-1000-10000.png", 1000, NULL, FALSE);
	imaqDispose(image);
	return 0;
}

Смотрим: записав 1000 и 10000 в ответ получаем в ImageJ 4000...40002. То есть видим снова домножение на 4 и две градации. А ну ка, если сохранить 0 и 1024 — тогда будет 0...32784. Тут 1024 уже на 32 домножилось, это 32768, и прибавилось 16 градаций. Всё чудесатее и чудесатее.

Теперь настало самое время заглянуть во внутренности PNG и прежде всего проверить, что же там реально сохраняется. К счастью, мне не надо расписывать внутреннее устройство этого формата, так как на Хабре об этом уже писали — PNG — not GIF! и Разбираем самый маленький PNG в мире.

Вкратце, кому лень читать — PNG разбит на секции-"чанки", каждый со своим тегом, в которых хранятся как сами пикселы, так и служебная информация, у них четырёхбуквенные имена. Нас интересует секция IDAT, она там одна (в большом файле их несколько). Смотреть удобнее всего hex редактором (либо любым онлайн просмотрщиком), я могу порекомендовать 010 Нех, он сразу разбирает структуру. Вот как хранится файл с двумя пикселами 1000 и 10000:

Об особенностях хранения 16 бит изображений в PNG формате - 3

Байты наши упакованы при помощи Deflate, использующим комбинацию алгоритмов LZ77 и Хаффмана, так что вначале их надо распаковать.

Мне проще всего воспользоваться опять же LabVIEW, там есть подходящая функция в OpenG тулките:

Об особенностях хранения 16 бит изображений в PNG формате - 4

Ну или вот так, если угодно:

import zlib
import binascii
data = bytes([0x78, 0x01, 0x01, 0x05, 0x00, 0xFA, 0xFF, 0x00, 0x0F, 0xA0, 0x9C, 0x42, 0x03, 0x9B, 0x01, 0x8E])
out = bytes(zlib.decompress(data))
print (binascii.hexlify(out))

Вывод:

b'000fa09c42

Таким образом у нас есть пять байтов. Первый байт — это фильтр (о нём можно отдельно написать, вкратце — PNG может хранить не только абсолютные значения интенсивности пикселей, но и относительные изменения по горизонтали или вертикали, или среднее, что в ряде случаев позволяет повысить степень компрессии, там пять возможных вариантов), ну а оставшиеся четыре байта — наши два пиксела по два байта, и тут всё встаёт на свои места: 0x0FA0 — это действительно 4000, а 0x9C42 — это 40002. Таким образом ImageJ, равно как и Питон честно показывают то, что реально записано в файле, чудес не бывает.

Минуточку, но ведь тот же OpenCV никогда не был замечен в таком преобразовании, это легко проверить:

#include <print>
#include "include/opencv2/opencv.hpp"
#pragma comment( lib, "lib/opencv_world490" )
using namespace cv;
int main()
{
    const unsigned short Pixels[2] = { 1000, 10000 };
    Mat src = Mat(1, 2, CV_16U, (unsigned short*)Pixels);
    imwrite("1000-10000-OpenCV.png", src);
    Mat dst = imread("1000-10000-OpenCV.png", IMREAD_UNCHANGED);
    double min, max;
    minMaxLoc(dst, &min, &max);
    std::println("Min = {}; Max = {}", min, max);
}

Логичный вопрос, который возникает — как же библиотека NI возвращает значения обратно? Собственно, когда у нас в руках два файла, сгенерированных разными библиотеками, это легко проверить, сравнив их. Для разнообразия можно воспользоваться tweakpng.

Об особенностях хранения 16 бит изображений в PNG формате - 5

Протерев очки, и немного поэкспериментировав, я и обнаружил различие — вот они наши 14 бит — сидят в теге sBIT! Также мы видим пару тегов zTXt и tEXt, хранящих кое-какую служебную информацию (один сжатый, другой нет), к примеру NI Image Type 7 означает беззнаковый 16 бит.

Таким образом, согласно документации, при сохранении 14 бит изображения в 16 бит контейнер происходит сдвиг на два бита, и нам нужно задвинуть эти данные обратно на разницу между 16 и тем, что стоит в sBIT. В случае, если данные в файле укладываются, скажем, в диапазон 0...1023, то реально будет записано 0 и 65535, а sBIT будет 0х0A, что суть 10, так что сдвигать обратно надо будет уже на шесть бит, а не на два. Ничего нестандартного здесь нет, просто по умолчанию все этот тег игнорируют (в силу редкости его использования, надо полагать). Таким образом, если вы сохраните библиотекой NI Vision, скажем однопиксельный PNG с пикселом 1023 или 32767, то результат, записанный в IDAT, будет абсолютно идентичен, единственное отличие будет в теге sBIT — он будет либо 0A, либо 0F (ну и контрольная сумма этого тега изменится, поэтому отличие в пяти байтах):

>fc/b 1023.png 32767.png
Comparing files 1023.png and 32767.PNG
00000029: 0A 0F
0000002A: 08 78
0000002B: 04 6E
0000002C: 3A CE
0000002D: B5 3A

Зачем так сделано? Стандарт пишет, что это позволяет упростить декодирование, и изначально введён как раз для случая медицинских изображений, которые зачастую обладают 12- или 14-бит глубиной. На самом деле это хорошо и правильно, поскольку в некоторых случаях нам действительно надо знать "истинную" глубину наших данных. Если указано, что данные наши — 14 бит, а интенсивность пикселей лежит в районе 1000, то это не значит, что мои данные 10 бит, а значит то, что весь динамический диапазон не используется полностью. Это можно учитывать, скажем при нормировке данных, когда мы переходим в формат float, и нормируем наши данные на диапазон 0...1, тогда максимум можно брать исходя из максимально возможного числа бит, а не динамически из максимума конкретного изображения.

Как приятный бонус такого открытия PNG мы получаем полный динамический диапазон в превьюшках Эксплорера. Он вполне понимает 16 бит, но PNG всегда выглядит светлее и контрастнее, по сравнению с TIFF, хотя и там и сям интенсивность пикселей формально одинакова:

Об особенностях хранения 16 бит изображений в PNG формате - 6

Теоретически (вероятно) можно попросить libpng сразу отдавать правильные значения, там сдвиг отрабатывается и о sBIT эта библиотека знает, либо получить его программно через png_get_sBIT() и обрабатывать извне.

Паззл в общем-то уже сложился, причина понятна, и то, как получить обратно оригинальные значения — тоже, но остаётся три небольших вопроса. Во-первых, каким образом библиотека позволяет хранить в PNG представление со знаком, во-вторых, можно ли записать файл так, чтобы "внешней коррекции" не потребовалось и, наконец, откуда берутся две дополнительные градации?

Знаковое представление

Наряду с "классическим" беззнаковым представлением шестнадцатибитных данных в диапазоне 0...65535 библиотека NI позволяет оперировать и со знаковым представлением signed short, когда наши данные находятся в диапазоне -32768...32767. Кстати, в ранних версиях этой библиотеки такое представление 16 бит было единственным, что добавляло некоторое количество головной боли разработчику, беззнаковый тип появился позже, где-то в 2007 году. Как бы то ни было, библиотека честно сохраняет и читает правильные значения и в этом случае. Надо отметить, что данное представление не является стандартом PNG, но мне стало любопытно, как же разработчики выкрутились, и я заглянул в PNG ещё раз. Ларчик открывался просто — был добавлен приватный тег scAl, который описывает сдвиг в отрицательную область (полагаю, его название происходит от Scale). Ну то есть, если данные содержат минимальное значение, скажем, -500, то в этот тег пишется FE0C, что означает 65036, и вычитая 65536, получаем как раз -500. Если же, например, данные имеют минимальное значение -1000, то будет записано 0хFC18, и это даст 64536-65536=-1000:

Об особенностях хранения 16 бит изображений в PNG формате - 7

Очевидно, что для открытия такого файла придётся проверять наличие этого тега и соответствующим образом пересчитывать пиксели, что библиотека и делает. Кстати, если вы обратили внимание, то отдельные буквы в наименовании тегов находятся в разных регистрах (IDAT, sBIT, scAl) — это не случайно, а часть стандарта. Глядя на тег scAl, можно сразу сказать, что он нестандартный.

Сохранение 16 бит изображения без тега sBIT

Очевидно, что если тег sBIT будет отсутствовать, то у открывающей программы нет шансов определить "истинный" диапазон данных и в этом случае мы получим ровно то, что записано в IDAT. Если в нашем случае мы всё-таки запишем туда значения 1000 и 10000, то такими они там и будут, именно так по умолчанию и поступает OpenCV. К счастью, библиотека NI позволяет записать данные и так, хотя это сделано несколько неочевидным образом, для этого надо насильно установить глубину в 16 при помощи IMAQ Image Bit Depth и установить флаг при сохранении "использовать эту глубину", вот тогда в файл будут записаны "честные" значения:

Об особенностях хранения 16 бит изображений в PNG формате - 8

Ну а в случае 14 бит именно здесь можно указать "истинную" глубину, и использовать именно её. По умолчанию же библиотека вычисляет sBIT динамически при сохранении исходя из максимальной интенсивности в сохраняемом файле.

Две дополнительные градации

Ну а что касается двух "лишних" градаций — то я уже был готов написать баг репорт, но ещё раз внимательно почитал стандарт, действительно, так положено делать, это рекомендованный метод для энкодеров. Суть в том, что в "избыточные" младшие биты мы должны поместить (продублировать) старшие, а не устанавливать их в нули. Кстати, дело не ограничивается 16 битами. Документация приводит пример, когда наши данные изначально пятибитные, и мы хотим хранить их в восьмибитном контейнере. Например, значение пиксела — 27 (десятичное). Битовое представление — 11011. Теперь мы сдвигаем наши данные на три бита влево, а в освободившиеся три младших пишем три старших 110, и в результирующий файл будет записано не 11011000 (216), а 11011110 (222):

4 3 2 1 0
---------
1 1 0 1 1

Сдвиг и дублирование даст нам значение 222:

7 6 5 4 3  2 1 0
----------------
1 1 0 1 1  1 1 0
|=======|  |===|
    |      Старшие биты повторяются для заполнения "открытых" битов
    |
Оригинальные биты

Теперь всё встало на свои места. Оригинальное значение 912 исходно представлено в 14 битах как 00 0011 1001 0000 (два старших бита нулевые) и после сдвига на два бита превращается в 16 бит 0000 1110 0100 0000, что даёт 3648, а вот оригинальное значение 12160 — 10 1111 1000 0000, и после сдвига два старших 10 дописываются в конец, таким образом в файл сохраняется 1011 1110 0000 0010, вот откуда берутся две дополнительные градации:

|=| - старшие биты
1 0 1 1 1 1 1 0 0 0 0 0 0 0 (12160)
1 0 1 1 1 1 1 0 0 0 0 0 0 0  1 0 (48642)
|=========================|  |=|
    |                       Старшие биты повторяются
    |
Оригинальные биты

Так что это не ошибка вовсе, а педантичное следование стандарту.

Ну и в конце я оставлю небольшой Питон скрипт, который смотрит в теги sBIT и scAl и правильно открывает такие изображения:

Я не питонист, я только учусь
import cv2
import struct
import numpy as np

def __getSBit(bytes):
  bytes = bytes[8:]
  sBit = 0
  while bytes:
      length = struct.unpack('>I', bytes[:4])[0]
      bytes = bytes[4:]
      chunk_type = bytes[:4]
      bytes = bytes[4:]
      chunk_data = bytes[:length]
      bytes = bytes[length:]
      if chunk_type == b'sBIT':
          sBit = int.from_bytes(chunk_data, "big")
          break
      bytes = bytes[4:]
  return sBit

def __getOffset(bytes):
  bytes = bytes[8:]
  Offset = 0
  sBit = 0
  while bytes:
      length = struct.unpack('>I', bytes[:4])[0]
      bytes = bytes[4:]
      chunk_type = bytes[:4]
      bytes = bytes[4:]
      chunk_data = bytes[:length]
      bytes = bytes[length:]
      if chunk_type == b'scAl':
          chunk_data_part = chunk_data[:2]
          Offset = 65536 - int.from_bytes(chunk_data_part, "big")
          break
      bytes = bytes[4:]
  return Offset

def getSigBits(filename):
  with open(filename, 'rb') as f:
      bytes = f.read()
  return __getSBit(bytes)

def getOffset(filename):
  with open(filename, 'rb') as f:
      bytes = f.read()
  return __getOffset(bytes)

def shift_offset(image_src, shift, offset):
  height, width = image_src.shape
  counter=0;
  if shift<16:
      temp=np.zeros(image_src.shape,dtype=np.int16)
      temp=(image_src>> shift).astype(np.int16)
      temp=temp-offset
      return temp

def load_lv_i16(file_path):
  offset=getOffset(file_path)
  sigbits=getSigBits(file_path)
  image_src=cv2.imread(file_path, cv2.IMREAD_UNCHANGED)
  shift=16-sigbits
  image=shift_offset(image_src,shift,offset)
  return image

if __name__ == '__main__':
  file_path="your_image.png"
  image=load_lv_i16(file_path)
  min_val,max_val, min_indx, max_indx=cv2.minMaxLoc(image)
  print(min_val, max_val)

Несмотря на то, что количество пользователей Хабра, использующих эту библиотеку, исчезающе мало, оно всё же отлично от нуля, так что надеюсь, что этот материал может кому-нибудь пригодиться.

Всем добра!

Автор: Андрей Дмитриев

Источник

* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js