Pull to refresh

Comments 86

Не очень понятно — исходное изображение двумерное или все-таки это трехмерный массив точек и вы его просто визуализируете?
Я так понял что двухмерное, трёхмерность строится по глубине цвета. Райт?
шейдер явно лезет в texture3D
> Я так понял что двухмерное, трёхмерность строится по глубине цвета. Райт?
Райт. У меня, кстати, трёхмерный desktop 1600x900x32
Исходное изображение, судя по всему — не одно, а куча двумерных снимков одного объекта (слои/срезы).
на самом деле не важно, одно оно или на каждый слайс свой файл, или в файле по несколько слайсов. софту без разницы :) обработка от этого не зависит.
но суть медицинских изображений в этом, в том что есть некоторый набор двумерных слайсов. это идет от сути получения этих изображений. я ниже в комментарии чуть подробней описала
UFO just landed and posted this here
это массив точек. изображения в медицинской технике это по сути набор 2д слайсов X*Y, количеством Z. это обусловлено тем как именно получаются изображения, например ст-сканер делает двумерный снимок при определенных координатах (расположении в пространстве), после этого положение меняется и делается след снимок. а софт уже может этот набор снимков, который либо разделен по файлам (каждый слайс в отдельный файл), либо сразу объединённый в один файл, визуализирует его.
Там вроде ж чуток сложнее — снимки с камеры не параллельны друг другу, и содержат вид через всю толщину объекта, а не инфу об одном слое. А в параллельные поверхности воксельных слоёв их пересчитывает хитрая математика, типа «преобразования Радона» и т.п.
мы сейчас немного о разном говорим. во-первых это все отходит от темы собственно поста. и то, какую информацию содержат в себе снимки, не относится к данному делу.
если вам интересно, почитайте про DICOM формат — это стандарт, который должны поддерживать практически все медицинские аппараты такого плана.

в данном контексте говорилось только про один «тег» этого формата — непосредственно данные. ими мы и оперируем
Абсолютно случайно мне довелось иметь дело с DICOM. Так вот ничего особенного в этом формате нет — информация о пациенте в виде ключ-значение + набор изображений. Какой там набор изображений хранится — это уже на совести того, кто в этот формат сохраняет. Это могут быть совершенно не связанные друг с другом изображения, также могут быть изображения в таком виде как описал товарищ Nashev в предыдущем комментарии.
мне с ним приходилось иметь дело долгое время. и вот насчет несвязанных друг с другом изображений могу поспорить :) если взять отдельный обычный дайком файл, то raw данные которые там хранятся непосредственно описаны в предыдущих тегах. и они не могут быть никак не связанны. это не имеет смысла
так же это не имеет отношения к текущей статье :)
Полагаю, 3D облака в Crysis рисуются чем-то подобным
Если я ничего не путаю, то там для облаков и для лайт-шафт используются наслоения
На сколько я знаю, лайт-шафт (sunshafts) идут в постпроцессинге через z-buffer. А трёхмерность облаков отчётливо видно только в начале первой миссии и в предпоследней при полёте на VTOL — но там очень хорошая плавность по всем трём осям, если под наслоениями подразумевается куча наложенных 2D спрайтов — то их должно быть ОЧЕНЬ много. Хотя неоптимизированный raycasting через 3D поверхность является по сути тем же самым, разница только в порядке суммирования. Тем не менее, я всё же думаю, там raycasting, т.к. слышал, что алгоритм для облаков похож на применяемый в parallax occlusion mapping, а там точно raytracing.
По ходу — да, так. Вышел на этот же PDF с SIGGRAPH 2003
Вообще, как я понял, автор не новый алгоритм описывает, а скорее показывает применение известных техник в несколько необычной области.

Т. е. задача звучит так: по результатам МРТ/КТ или чего-то подобного отобразить проекцию костей, кровеносной системы или что другое (т. е. некоторой выборки вокселов по плотности) на экран.

Подход, соответственно, такой:
— собираем трехмерный массив плотностей
— трактуя плотность, как прозрачность воксела, применяя при этом пороговый фильтр, рендерим рейкастингом

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

Для еще большей оптимизации скорости выполнения и уменьшения нагрузки на ЦП, реализуем это дело на шейдерах.
на самом деле в этой области данный алгоритм, рейкастинг, является самым обычным и используемым.
плотность при этом не трактуется как прозрачность, есть просто определенное значение плотности, которое мы ищем для того чтобы построить изоповерхность.
но если брать только цвет вокселя, получится равномерно окрашенный силует. Чтоб была поверхность надо же и наклон поверхности распознавать — как это делается? А то глубину отслеживать, для затенения более дальних мест…
опять таки.
в коде шейдера есть директива SHADER. именно она дает там затемнения и вообще реалистичность отображения объекта. подробнее об этом можно почитать тут — Shading, хотя вы скорей всего знаете этот алгоритм.

так вот, данная тема не относится к самой сути рендеринга изоповерхностей. это так сказать, общий алгоритма для «отрисовки» теней. Поэтому если shading не использоваться — у нас да, получится равномерно закрашенный силуэт. Я завтра, когда буду на работе, могу вставить пример картинки, пока же например можно посмотреть на изображение с той же ссылки на википедию, получится что-то типа такого

image
как-то так:
  1. sampler3D tex_;
  2. float sample(float3 texCoord) {
  3.   float v =  tex3D(tex_, texCoord).w;
  4.   return v;
  5. }
  6.  
  7. float gradMag_;
  8. float sizevol = 256.0;
  9.  
  10. float3 getNormal(float3 coord) {
  11.  
  12.   float dx = 1.0f/sizevol;
  13.   float dy = 1.0f/sizevol;
  14.   float dz = 1.0f/sizevol;
  15.   float vx = sample(coord+float3(dx,0,0));
  16.   float vy = sample(coord+float3(0,dy,0));
  17.   float vz = sample(coord+float3(0,0,dz));
  18.   float vxb = sample(coord-float3(dx,0,0));
  19.   float vyb = sample(coord-float3(0,dy,0));
  20.   float vzb = sample(coord-float3(0,0,dz));
  21.   float3 n = -float3(vx-vxb,vy-vyb,vz-vzb);
  22.   float3 normal;
  23.  
  24.   gradMag_ = length(n);
  25.   if(gradMag_<1e-5) {
  26.     normal = float3(0,0,0);
  27.   } else {
  28.     normal = normalize(n);
  29.   }
  30.  
  31.   return normal;
  32. }
  33.  

и далее в теле шейдера:
  1.       if(color.a>=0.001) {
  2.         float3 normalVec = getNormal(coord);
  3.         if(dot(normalVec,ray)>0) normalVec*=-1;
  4.  
  5.         float ndotl = dot(normalVec, lightVec);
  6.         float ndoth = dot(normalVec, halfVec);
  7.         float4 light = lit(ndotl, ndoth, shininess);
  8.         color.xyz = saturate(ambient*color.xyz + diffuse*color.xyz*light.y + specular*color.w*float3(1,1,1)*light.z);
  9.  
  10.         compositeColor = saturate((1-compositeColor.w)*color + compositeColor);
  11.       }
  12.  
немного избыточно в плане строк кода, но в целом вроде все верно )
Это какая-то уличная магия — отнимаем 2 соседние плотности из текстуры и получаем нормаль %)
Внимание вопрос! Что хранится в текстуре из которой вы делаете выборку? :) Это же не карта высот, где отняв две соседние высоты получается нормаль. У неё в текстуре плотность тканей, шифрованная в RGB %)
Я исхожу из предположения, что функция непрерывна. Вычисление градиента в данном случае позволяет получить нормаль к поверхности. Я вижу в этом математику, а никак не «уличную магию». Предложите свой алгоритм вычисления нормали, и мы продолжим беседу.
PS: если вы предлагаете добавить деление на (dx,dy,dz), то, возможно, я приму Вашу критику в общем случае. В данном же конкретном примере деление не играет роли, поскольку dx=dy=dz и после выполнения Normalize вектор становится единичным.
Причем здесь деление? Вы вытаскиваете из текстуры 2 значения, отнимаете их и считаете что это нормаль. Почему вдруг? Объясните.
float vx = sample(coord+float3(dx,0,0));
float vxb = sample(coord-float3(dx,0,0));
float3 n = -float3(vx-vxb,vy-vyb,vz-vzb);

Почему n — это направление нормали? Если vx и vxb — это плотность тканей в соседних вокселах? Объясните почему отняв плотность кости от плотности мяса получится нормаль?

Ваш алгоритм работает только для карты высот, где в пикселах хранится расстояние.
Мы считаем градиент функции плотности в указанной точке, который затем принимаем за нормаль.
Почитайте, что пишет trish в ответ на Ваш комментарий.
Если честно, мне уже становится не интересно что вы тут все пишете. Это не имеет ничего общего ни с математикой, ни с отрендеренными картинками из поста, ни просто со здравым смыслом.
в текстуре плотность тканей. она не шифрована в RGB, просто эту плотность можно потом дальнейшем как угодно интерпретировать и получать необходимые значения на выходе.
возьмем например плотность мяса и плотность кости — на границе будет большая разница, то есть большой градиент (его длина), внутри одной сущности, например мяса, градиент тоже будет из-за шумов, но весьма маленький. И этот градиент и есть наша нормаль

так же по коду можно заметить, что берутся не просто соседние воксели, а нечто усредненное из всех 9 вокселей, окружающих искомую точку
ой, не из 9 конечно же, а из 26
А вам самой не кажется странным перерасчет градиента плотности в нормаль? Берем 2 черепа одной и той же формы (т.е. нормаль во всех точках одинакова у них). У одного черепа рядом более плотные ткани чем сам череп, у другого менее плотные чем череп. В итоге у вас получится не просто разная нормаль, а направленная вообще в разные стороны (знак при отнимании плотности черепа от плотности соседних тканей будет разный). Ну бред ведь.
> на самом деле в этой области данный алгоритм, рейкастинг, является самым обычным и используемым.

Это не алгоритм необычен (он-то большинству аудитории наверняка известен). А вот область не характерна для хабра.

> плотность при этом не трактуется как прозрачность

Зато с прозрачностью понятнее ;-)
да, возможно не характерна для хабра, но весьма характерна для визуализации в медицине :)

с прозрачностью может и понятней, учитывая вашу фразу о пороговом фильтре, но все же не совсем верно

наверно стоило начинать все же с подробной статьи именно о ray casting в визуализации медицинских изображений :)
А почему LUT не называют «палитрой»?
назвать его в переводе можно как угодно. и палитрой в том числе.
«бантиком» я б его не назвал…
вы могли бы заметить, что я его вообще никак не называю на русском :)
именно это я и заметил, как вы могли бы заметить ;)
а как определяется «наклон» поверхности в точке пересечения с лучем? По локальному градиенту?
может быть по глубине соседних точек?
вас интересует наклон поверхности в плане шейдинга (shading)?
потому что если брать сам алгоритм отрисовки изоповерхности, который тут рассматривается — то здесь нам нужен вектор направления луча. у нас есть положение камеры — gl_ModelViewMatrixInverse и текущее положение луча, которое мы берем из пиксельного шейдера. на основании этого мы и вычисляем вектор направления, по которому движемся
Да, меня тоже именно это интересует. Понятно, что вы пытались упростить, но люди то тут не дураки — отрендеренные картинки не соответствуют опубликованному вами шейдеру. В вашем шейдере цвет пиксела всегда равен isoColor, а значит абсолютно все равно какой sampleRate, полоски на лбу черепа ну никак не получатся :)
первые картинки и не будут соответствовать, потому что были сделаны с уже кучей всего добавленного, там и не написано про это — это был просто пример изоповерхностей. но например вторые картинки (на черном фоне) — это именно те, которые сделаны с применением данного алгоритма. и суть данной статьи была не в шейдинге, а именно в отрисовке изоповерхности. шейдинг сам по себе достаточно сложная тема, которую нужно освещать отдельно и подробно.

да, здесь алгоритм упрощен, так как например в реальном алгоритме куча директив, которые выбирают — какой режим рендеринга использовать, дополнительные параметры — preintegration, alpha correlation и т.д. но для данной темы все это было ни к чему
> но например вторые картинки (на черном фоне) — это именно те, которые сделаны с применением данного алгоритма
Но НЕ с помощью шейдера из статьи. Я хорошо знаю GLSL и рендеринг в целом, не сомневайтесь :)
я и не сомневаюсь в вашем профессионализме :)
но как я уже написала — здесь не рассматривается шейдинг, который несомненно использовался для получения данных картинок, иначе это был бы просто контур, заполненный однородным цветом — как на картинке, что я вставила выше

#ifdef SHADER
// do shading < — вот сюда вставить функцию, которая будет реализовывать этот алгоритм
#endif

эту дополнительную функцию пришлось бы отдельно рассматривать и объяснять, а она не маленькая. и в общем не входит в формат именно данной статьи )
Меня не инетересуют алгоритмы освещения, как я уже писал ниже :)
Я спрашиваю откуда вы данные для этого алгоритма берете. Томограф же вам дает лишь 3D текстуру, без нормалей. Значит вы их строите каким-то образом сами. Это как бы не относится к алгоритмам шейдинга — в любом алгоритме используются нормали как данное ;)
Если честно — я до сих пор не очень понял формат статьи. Она содержит лоскутное поверхностное ведение в воксельную графику т медицины и огрызочек примитивного шейдера. И оставляет ощущение «что это было?». Но направление, из которого приведённые лоскутки, само по себе интересно — потому и в комментах обсуждается, и плюсуется.

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

я в свое время, 2 года назад когда занималась именно конкретно этой темой, столкнулась с небольшими проблемами именно в рендеринге изоповерхностей. когда разобралась — все оказалось легким. но я тогда подумала — а может кто-то еще с этим столкнется. вот и решила написать. но руки только сейчас дошли :)
Могли бы не заморачиваться :)
Судя по самым первым комментариям к статье — всё равно вышло непонятно для тех кто не знаком с рендерингом.
:))) ну ничего, зато хоть интересно
Даешь вторую статью с описанием как получаются картинки как в начале статьи — с освещением. Нормали в каждой точке откуда берутся? В вашем шейдере этого нет :)
это не относится к данной теме. для этого нужно просто раcсмотреть любой алгоритм Shading, я например знаю и здесь использовала Phong shading, который можно использовать для рендеринга любых поверхностей, не только изо.
Если народ интересуется — я могу написать попозже статью и об этом алгоритме
Алгоритмы освещения меня не интересуют, это я вам и сам расскажу. Но в любом алгоритме освещения используется нормаль к поверхности, иначе вы никак не посчитаете отраженный свет от источника. Откуда вы берете нормали к поверхностям — вот в чем вопрос. Это для меня не очевидно.
Думаю что тем же алгоритмом, что и расчет normalmap'ы из heightmap'ы — по соседним точкам смотрим градиент, считаем углы, суммируем, нормалиуем… и вот уже вроде как нормаль )
К тому же я бы рекомендовал Blinn–Phong — работает быстрее, более аккуратно моделирует BRDF.
Тогда нужно рисовать в 2 прохода. Автор же уверяет, что в этом же шейдере делается освещение.
да, по соседним точкам строится градиент, считается наклон, нормализуется и рассчитывается освещение.
соседние точки находятся как текущее положение луча+-смещение, которое передается извне и равно соответственно (1/x,1/y,1/z). если интересует могу подробней описать ))
Да, опишите, пожалуйста.
Потому что для нахождения соседних точек вам нужно первым проходом сохранить в текстуру «положение луча» (глубину), а вторым проходом уже из этой текстуры делать выборку по соседним точкам. Как у вас выходит это в одном шейдере — загадка.
так у меня есть положение луча. у меня есть vec3 rayPosition, которое я и нахожу в цикле. у меня есть смещения (задаются извне), и я просто вызываю функцию, которая по заданным смещениям берет из текстуры соседние от rayPosition точки
У меня создается впечатление, что вы сами не понимаете о чем говорите. Ладно, подождем вторую статью.
как-то некорректно вы выражаетесь.
возможно мы просто говорим с вами о разных вещах или об одном и том же, но с разных позиций
>как-то некорректно вы выражаетесь.
Есть такая беда. Постараюсь быть вежливее.
Я так понял, что алгоритм для каждого луча, упёршегося в воксел нужной плотности, смотрит непосредственно в 3D-текстуре на вокселы, соседние к найденному (которых в плоскости получилось бы 6, а в объёме — выходит втрое больше), и ищет какие из них с той же плотностью, т.е. типа на той же изо-поверхности. И дальше яркость точки считает опираясь на них.

trish, угадал?
О, это уже хоть что-то приближенное к реальности. Но точность такого метода просто ужасна. Фактически у вас нормаль с точностью в 45 градусов будет. Их бы потом сгладить с соседними нормалями…
В целом картинка достаточно плавно освещена.
да, так и есть. shading хорошо работает на границах материалов, а внутри мы получаем в принципе одно направление, которое тем не менее достаточно случайное, так как зависит от шума еще
в плоскости их бы получилось бы 8. в объеме соответственно 26. но по сути вы правы.
да, с плоскостью считал-считал да обсчитался. А в объёме я имел ввиду что вокселы непосредственно над и под центральным вроде бы не интересуют — потому получается ровно втрое.

Кстати, для прикидки насчёт нормали вроде бы вообще можно обойтись и вовсе диагональными, угловыми — а таких всего 8 для куба.
и как раз таки из-за того, что этот rayPosition зависит от скорости прохода, при большой скорости и большом шаге мы может пропустить точное положение изоповерхности и получаются такие вот круги, о которых я в конце упомянула. и которые надеюсь еще рассмотреть в следующей статье
т.е. Вы на скорости, типа, рисуете каждый второй пиксел, но не точкой а кружочком?
нет :))
при большом шаге мы можем проскочить точное положение изоповерхности. рисуется все одинаково.
по соседним точкам смотрим градиент, считаем углы, суммируем, нормализуем…

Совершенно верно, я добавил пример для наглядной демонстрации.
Не буду придираться ни к стилю изложения, ни к коду шейдера — по ним и так уже прошлись товарищи выше.

Мне вот другое непонятно: почему рендеринг изоповерхностей? Если сканер (MRI/CAT/что еще) на выходе дает массив фотографий срезов, их ведь надо как-то проанализировать, чтобы построить непрерывную поверхность. Есть в планах рассказ об этом?
Зачем её строить, если нужно лишь рисовать? Рисовать и так можно, и именно об этом рассказ.

А если очень хочется строить — то тут сильно всё зависит от того, что считать «построеной поверхностью» — триангулированную сетку, или ещё что.
Ну, я надеялся, что простой визуализаций дело не ограничится :-)

Входные данные ведь не являются ни изо-, ни поверхностью. Облако точек — не более.
Дык построишь поверхность, а чего дальше-то делать? Задача-то «посмотреть»! Она и на имеющихся вокселях решается. Более того, физику — столкновения — тоже на вокселях считать вообще милое дело.

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

Но вот смотреть хочется не только изоповерхности, но (в основном!) срезы, а по изоповерхностям срезы строить — убийственно. А вот по вокселям отрисовать срез — раз плюнуть. Даже легче, чем изоповерхность отрисовать.
режим рендеринга изоповерхностей лишь один из режимов визуализации, причем даже не основной
и вы не построите векторную графику в данной области, здесь физическая суть другая )
Как это не построю?? Векторные изоповерхности — легко.

Другое дело, что их как правило совершенно не достаточно для практической работы. Но вот устроить из них что-то типа персонального google body было б занятно.
В симуляции жидкостей (CFD), очень хорошие результаты получаются при моделировании на решетке с последующим построением поверхности для рендеринга. Это позволяет сгладить артефакты, вызванные низким разрешением модели, не увеличивая время симуляции.
Тут при крупном зуме тоже придётся один воксел много раз рисовать, и высчитывать субвоксельные сдвиги… Векторизация происходит непосредственно при отрисовке пиксела экрана, весьма локальная и не сохраняемая никуда.

Зачем в CFD придумано векторизовать решётку заранее?
1. Кеширование. Рейкастинг ведь отнюдь не бесплатный и выполнять его несколько раз на одних и тех же данных бессмысленно.

2. Интерполяция. Возьмите картинку с undersampled черепом из статьи. Если построить поверхность, то артефакт исчезнет.
Наверное затем, что обработав данные 1 раз, построив триангулированную поверхность и посчитав нормали, всё будет рисоваться на порядки быстрее и намного красивее. Все таки рейкастинг, да еще куча выборок соседних вокселей из текстуры для освещения, не располагают к скорости. Автор сама пишет про 13 fps для несчастной головы. И это еще не известно на каком железе. Судя по тому, что этот монстряковый шейдер не превысил количество допустимых инструкций, видяшки у них хорошие :)
ну это да — такой подход действительно позволит покрутить выбранную поверхность с хорошим fps.

С другой стороны можно сдублировать 3D-текстуру, и пересчитать в копии плотность каждого вокселя в нормаль соответствующей изоплоскости (благо размерность идентична). Шейдеру работы поубавится, кажись в разы… Если памяти хватит.

Нукась — оно разрешением порядка пары точек на миллиметр нынче, к сожалению не больше — получается для одной только головы, которая приблизительно влазит в кубик 20х20х20 см, нужно 4003*4 байта с выравниванием = 256 мегабайт… А для всего туловища, которое по пропорциям в два-три раза шире, и в 8 выше?..
Sign up to leave a comment.

Articles