Рендеринг изоповерхностей с использованием алгоритма рейкастинга

В данной статье хочу рассказать вам про Isosurface rendering или рендеринг изоповерхностей.

По долгу работы я постоянно использую документацию лишь на английском или немецком языке. Поэтому в тексте я буду часто использовать для всяких терминов английский вариант слова.

Итак, что такое Isosurface — это, как говорит нам Википедия, 3х-мерный вариант изолинии, которая в свою очередь «представляет собой линию, в каждой точке которой измеряемая величина сохраняет одинаковое значение». На словах это может быть не совсем понятно — посмотрим лучше на картинки.



На данных изображениях показаны поверхности с разным значением «измеряемое величины» — плотности (density). На первом значение меньше, чем на втором.

Для начала краткая теория.

Я занимаюсь обработкой медицинских изображений. Данные, получаемые из медицинских аппаратов в основном содержат так называемое raw или grayscale изображение — от 8 бит (256 оттенков) до 16 бит (65536 оттенков серого цвета). Так как на деле 16 бит — это очень много, часто используется сокращенный вариант — 12 бит, то есть используется часть от старшего байта, это дает нам 4096 оттенков серого цвета.
Далее для преобразования серого изображения в цветное можно использовать так называемые LUT — look up table. Эта таблица представляет собой набор строк, каждая из которых содержит RGBA значение. Для 12 бит эта таблица содержит 4096 строк. И когда мы встречаем воксель со значением, например, 542, мы соответственно берем из 542 строки этой таблицы указанные там значения RGBA.
Данные таблицы составляются квалифицированными медиками/инженерами, которые знают как коррелировать значение плотности с RGBA. Использование альфа-канала, то есть прозрачности, позволяет создать таблицы, применение которых позволит нам отобразить только кости, или только вены, или артерии — остальное будет прозрачным. То есть, меняя на лету LUT, мы получаем каждый раз новое отображение одного и того же объекта.

Основным способом рендеринга данных изображений является Volume ray casting (что-то похожее на русском — Объемный рендеринг). Слишком подробно его я тут описывать не буду. Но понимание его работы нам важно для понимания процесса рендеринга изоповерхностей.

Алгоритм ray casting состоит в том, что мы лучами пронизываем наш объект. Луч исходит из нашего глаза (камеры), проходит через каждую точку экрана (каждый пиксель), и пересекается с нашим объектом в определенном вокселе (если есть пересечение). На этом луч не останавливается, а идет дальше, пересекая дальнейшие воксели и определенным образом аккумулируя информацию из каждой точки. Критериев остановки луча может быть несколько, наиболее распространенный — когда альфа аккумулируемого значения близко к 1 (на практике используется значение a>0.95), либо, например, если мы вышли за границы изображения. То есть по сути во время трейсинга мы отбрасываем прозрачные воксели и определенным образом акумулируем значения полупрозрачных, пока не дойдем до цельного объекта, который дальше не пропускает наш луч в силу своей непрозрачности. Полученное в результате значение и используется для отрисовки на экране.

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

Данный алгоритм можно реализовать на CPU (и в принципе первые реализации были именно такие), но работать он будет очень медленно. Намного быстрее обрабатывается все на GPU с использованием фрагментного шейдера (GLSL). Необходимые параметры передаются извне — граничное значение, цвет, скорость прохода и т.д. Ниже представлен код самого простого фрагментного шейдера, который рендерит изоповерхность. Так как для изоповерхности не важны все предыдущие или последующие воксели, здесь происходит не аккумулирование значений, а просто используется только одно граничное значение цвета — isoColor.

1.  
2. uniform sampler3D volume;
3. uniform float sampleRate;
4. uniform float isoValue;
5. uniform vec4 isoColor;
6.  
7. varying vec4 texCoord0;
8.  
9. void main()
10. {
11.   floast steps = 1.0 / sampleRate;
12.   vec4 outputColor = vec4(0.,0.,0.,0.);    
13.   float isoThr;
14.  
15.   // get ray position and ray direction
16.   vec4 vPosition = gl_ModelViewMatrixInverse[3];
17.   vec3 rayPosition = texCoord0.xyz;
18.   vec3 vecDif = rayPosition - vPosition.xyz;
19.   vec3 rayDirection = sampleRate * normalize(vecDif);
20.  
21.   // for all samples along the ray
22.   while (steps)
23.   {
24.     steps--;
25.     // get trilinear interpolated value from 3d texture
26.     float value  = texture3D(volume, rayPosition);
27.     isoThr = value-isoValue;
28.    
29.     // check if we get isosurface line
30.     if (isoThr < 0)
31.     {
32.       // march to the next sample
33.       rayPosition = rayPosition + rayDirection;
34.        
35.       // get next density
36.       continue;
37.     }
38.    
39.     // else we do color transformation
40.     outputColor = isoColor;
41.     outputColor.rgb = outputColor.rgb * outputColor.a;
42.  
43. #ifdef SHADER
44.     // do shading
45. #endif
46.  
47.     break;
48.   }
49.   gl_FragColor.rgba = outputColor;
50. }
51.  

______________________

Текст подготовлен в Редакторе Блогов от © SoftCoder.ru

Одним из главных параметров ray casting является шаг трассировки — с какой скоростью луч проходит через наш объект. Если у нас рендерится куб размерностью X*Y*Z, то идеальной скоростью будет 1/MAX(X,Y,Z). То есть прирост должен быть не больше 1 вокселя для того, чтобы ничего не пропустить. Но с другой стороны часто не нужна такая степень детализации, которая к тому же влияет на производительность.



На представленных картинках мы видим изображения, отрендеренные с разным sampleRate. В первом случае использовался как раз таки идеальный 1/MAX(X,Y,Z) — изображение получилось плавным, без видимых переходов, так как мы не пропустили ни одного вокселя, НО! FPS при этом равняется всего лишь 13 кадров в секунду.
Для второго варианта использовался sampleRate в два раза больше — то есть мы обрабатывали каждый второй воксель, и при этом уже видны круги, которые образуются из-за погрешности в вычислении нахождения первого вокселя изоповерхности. Но FPS при этом повысилась до 37 кадров.
Как увеличить скорость прохода алгоритма и при этом не потерять в качестве я постараюсь рассмотреть в следующей статье.