Как обычно, в самом начале я постараюсь сэкономить время читателей — суть в том, что в 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:
Чтобы не смущать вас кодом на 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 = %f\n", 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
Permissions: rw-rw-rw-
Format: PNG (Portable Network Graphics)
Mime type: image/png
Class: DirectClass
Geometry: 1024x1024+0+0
Units: Undefined
Colorspace: Gray
Type: Grayscale Endianness: Undefined
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):
На Си, что б понятнее:
#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:
Байты наши упакованы при помощи Deflate, использующим комбинацию алгоритмов LZ77 и Хаффмана, так что вначале их надо распаковать.
Мне проще всего воспользоваться опять же LabVIEW, там есть подходящая функция в OpenG тулките:
Ну или вот так, если угодно:
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.
Протерев очки, и немного поэкспериментировав, я и обнаружил различие — вот они наши 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, хотя и там и сям интенсивность пикселей формально одинакова:
Теоретически (вероятно) можно попросить 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:
Очевидно, что для открытия такого файла придётся проверять наличие этого тега и соответствующим образом пересчитывать пиксели, что библиотека и делает. Кстати, если вы обратили внимание, то отдельные буквы в наименовании тегов находятся в разных регистрах (IDAT, sBIT, scAl) — это не случайно, а часть стандарта. Глядя на тег scAl, можно сразу сказать, что он нестандартный.
Сохранение 16 бит изображения без тега sBIT
Очевидно, что если тег sBIT будет отсутствовать, то у открывающей программы нет шансов определить "истинный" диапазон данных и в этом случае мы получим ровно то, что записано в IDAT. Если в нашем случае мы всё-таки запишем туда значения 1000 и 10000, то такими они там и будут, именно так по умолчанию и поступает OpenCV. К счастью, библиотека NI позволяет записать данные и так, хотя это сделано несколько неочевидным образом, для этого надо насильно установить глубину в 16 при помощи IMAQ Image Bit Depth и установить флаг при сохранении "использовать эту глубину", вот тогда в файл будут записаны "честные" значения:
Ну а в случае 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)
Несмотря на то, что количество пользователей Хабра, использующих эту библиотеку, исчезающе мало, оно всё же отлично от нуля, так что надеюсь, что этот материал может кому-нибудь пригодиться.
Всем добра!
