Детальніше про розробку софта рентгенівського томографа



Вчені з Томського державного університету створили микротомограф. Томограф дозволяє з точністю до мікрона дізнатися про внутрішній структурі різних матеріалів, наприклад, алмазів.

Але ж цікавіше в нього запхати муху.



Перед EDISON Software Developement поставили завдання написати софт для микротомографа. Про те, як вони успішно впоралися із завданням, була стаття на Хабре (Як за 5233 людино-години створити софт для микротомографа з описом алгоритмів, математичних методів, реалізації і налагодження.

Ненаситні читачі засипали нас запитаннями, на які ми, нарешті, сформулювали відповіді…

Томограф може просвітити матеріал з роздільною здатністю до мікрона. Це у 100 разів тонше людського волосся. Після сканування програма створює 3D-модель, де можна подивитися не тільки на зовнішню сторону деталі, але і дізнатися, що у неї всередині.

Відео результатів сканування
















Пристрій томографаСписок основних технічних характеристик приладу.

• Кількість дозвільних елементів детектора: 2048 x 2018 осередків при розмірі одного елемента не більш 13,3 x 13,3 мкм.
• Дозвіл: 13 мкм.
• Габаритні характеристики: 504 x 992,5 x 1504 мм.
• Маса: 450 кг.
• Поле зору: 1 мм.
• Робочий діапазон довжин хвиль: 0,3–2,3 А.
• Точність позиціонування електромеханічних модулів руху в системі позиціонування РМТ: ±1 мкм.
• Відповідність вимогам безпеки: ГОСТ 12.1.030-81.
• Захист від рентгенівського випромінювання: 1-3 мкЗв/год.


Рис. Зовнішній вигляд


Рис. Принцип роботи електро-механічної частини


Рис. Джерело рентгенівського випромінювання

• Напруга: 20-160 кВ.
• Сила струму: 0-250 мкА.
• Потужність: 10 Вт.
• Діаметр фокальної плями: 1-5 мкм.


Рис. Позиціонування

• Хід ротора: СЭМД — 360°, ЛЕМД — 100 мм.
• Точність позиціонування, не менше: ±0,5 мкм.
• Швидкість переміщення ротора: від 0,01 до 20 мм/с.
• Потужність: 0,7 кВт при напрузі 70 Ст.


Рис. Детектор рентгенівського випромінювання на базі ПЗС-матриці

• Чутлива область ПЗС: 2048 x 2048 пікселів.
• Геометричний розмір пікселя: 13 x 13 мкм.
• Геометричний розмір чутливої області ПЗС: 26,6 x27,6 мм.
• Вбудований двошвидкісний АЦП: 16 біт, 100 кГц і 16 біт 2МГц.


Питання №1

З задоволенням почитав про звернення перетворення Радону з точковим джерелом і про реалізацію цього всього з використанням CUDA Uranix

Розрахунок Радону для точкового джерела, якщо описати схематично, аналогічний розрахунку для паралельного режиму, але тільки напрямок променів не паралельно, що впливає на координати точок в кубі об'єкта, куди додаються щільності. При такому схематичному описі, з точки зору реконструкції, немає великої різниці, який алгоритм використовувати, змінюються тільки координати відрізків/променів та формули їх розрахунку.

Далі, коли на основі тіньової проекції ми отримали промені від джерела до кожної точки проекції в просторі координат микротомографа, знаючи положення обертового куба об'єкта, вздовж кожного з променів розраховуються воксели, що лежать на промені. Кожному вокселю присвоюється значення щільності пікселів тіньової проекції, в яку входить промінь від джерела, тип даних ваг перетворений до float. В результаті для кожної тіньової проекції повторюємо цю операцію проектування на воксели з урахуванням кута проекції. Після цього всі сумарні значення вокселей нормалізуються до розрядності, в якій вони будуть зберігатися в подальшому, звичайно це було 16-бітне беззнакове ціле. Даний опис максимально спрощено пояснює схему реконструкції та отримання кінцевої воксельної моделі, далі наведено більш детальний опис цього алгоритму.

Рентгенівська система генерує плоске тіньове зображення повної внутрішньої об'ємної структури зразка, але в окремій тіньової проекції глибинного розподілу структур зразка повністю змішуються. Тільки рентгенівська томографія дозволяє візуалізувати і виміряти просторові структури зразків без їх хімічної і механічної обробки.


Рис. 1. Геометрія паралельних пучків

Будь рентгенівське тіньове зображення є плоскою проекцію тривимірного об'єкта. У найбільш простому випадку ми можемо описати його, як зображення, отримане в паралельних рентгенівських променях (Рис. 1). В даному наближенні кожна точка тіньового зображення містить загальну інформацію по адсорбції конкретного рентгенівського пучка на всьому обсязі тривимірного об'єкта. Для паралельної геометрії рентгенівських пучків реконструкція об'ємного зображення зразка з двомірної тіньової проекції реалізується з допомогою реконструкцій серії двомірних зрізів зразка вздовж одновимірних тіньових ліній. Можливість такого роду реконструкцій демонструється на простому прикладі: розглянемо об'єкт з єдиною точкою з високою адсорбцією в невідомому місці.


Рис 1.1. Схематичне зображення трьох різних положень поглинаючої області та відповідна реконструкція з отриманих тіньових проекцій

У одномірної тіньової лінії буде спостерігатися зменшення інтенсивності внаслідок її поглинання на адсорбирующем об'єкті (Рис. 2). Тепер ми можемо змоделювати в комп'ютерній пам'яті порожній ряд пікселів (елементів зображення), що відповідає передбачуваному зміщення об'єкта. Природно, слід впевнитися, що всі частини реконструйованого об'єкта будуть знаходитися в полі зору. Оскільки ми маємо координати тіней від поглинаючих областей об'єкта, ми можемо виділити в реконструюється області в пам'яті комп'ютера всі можливі положення поглинаючих областей всередині об'єкта у вигляді ліній. Тепер почнемо обертати наш об'єкт і повторювати цю операцію. У кожному новому положенні об'єкта ми будемо додавати до реконструйованої області лінії можливих положень об'єкта у відповідності з положенням його тіньових проекцій. Ця операція називається зворотним проектування. Після кількох обертів ми можемо локалізувати положення поглинаючої області в межах обсягу реконструкції. Зі збільшенням числа тіньових проекцій з різних напрямків, ця локалізація стає все більш чіткої (Рис. 3).


Рис. 2. Реконструкція точкового об'єкта з використанням різного числа зсувів

У разі реконструкції на основі нескінченного числа проекцій виходить зображення з хорошою чіткістю визначення позиції області поглинання всередині досліджуваного об'єкта. У той же час точкове зображення буде супроводжуватися розмитою областю, оскільки воно було отримано в ході накладення ліній з усіма можливими відхиленнями. Тепер, оскільки нам відомо, що зображення утворено точковим об'єктом, ми можемо провести попередню корекцію початкової інформації в лініях сорбції щоб зробити кінцеве зображення більш відповідним реальному об'єкту. Ця корекція додає деяку кількість негативної адсорбції по зовнішній межі точки, щоб усунути розмиття, притаманну процесу зворотного проектування (Рис. 4). Цей алгоритм дає не тільки зображення перерізів окремих точкових структур, але і дозволяє досліджувати реальні об'єкти. Кожен матеріальний об'єкт може бути представлений як велика кількість окремих елементарних поглинаючих обсягів і лінійна адсорбція в кожному рентгенівському пучку відповідає сумарній адсорбції на всіх поглинаючих структурах, зустрінутих пучком.


Рис. 3. Фільтрація зворотної проекції

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


Рис. 4. Геометрія розходиться пучка

У віялової геометрії реконструйовані перерізу будуть мати деякі спотворення в областях, віддалених від оптичної осі пучка. Для вирішення цієї проблеми ми повинні використовувати алгоритм тривимірної реконструкції конічного пучка (Фельдкамп), щоб брати до уваги товщину об'єкта. Іншими словами, промені, що проходять крізь фронтальну і задню поверхні об'єкта, не будуть спроектовані на той самий ряд детектора (Рис. 6). Таким чином, найбільш швидко реконструюється перерізом є осьовий переріз, оскільки воно не вимагає корекції на товщину зразка.


Рис. 5. Геометрія конічного пучка

У разі рентгеноскопії зразка зображення містить інформацію про зменшення інтенсивності падаючого випромінювання всередині об'єкта. Оскільки поглинання рентгенівських променів описується експоненціальною залежністю (закон Ламберта-Біра), лінійна інформація про адсорбції випромінювання з тіньового зображення може бути розкрито з використанням логарифмічної функції. Оскільки логарифмування є нелінійною операцією, то найменший шум у сигналі призводить до значних помилок у реконструкціях. Для позбавлення від такого роду помилок може бути використано усереднення початкових даних. З іншого боку, ми могли б спробувати поліпшити відношення сигнал-шум в тіньовому зображенні за допомогою оптимізації часу витримки для накопичення найбільш репрезентативної інформації. Найбільш ефективним шляхом шумозаглушення в процесі реконструкції є правильний вибір функції корекції/фільтрації для конволюции перед зворотним проектування. У найпростішому випадку (описаному вище) функція корекції виробляє дві реакції негативної адсорбції навколо будь-якого піку сигналу або шуму в тіньовий лінії, і така поведінка стає дуже небезпечним при роботі з гучним сигналом. Вибір для проведення корекцій функції конволюции зі спектральними обмеженнями (вікно Хамминга) дозволяє вирішити цю проблему.

Опис роботи CUDA і реконструкції

Реконструкція даних за тіньовими проекціями
Максимальний розмір куба для реконструкції 1024x1024x1024 вокселя. Назвемо такий куб одиничним. Алгоритм дозволяє реконструювати менші куби, але при цьому розмірність по всіх напрямках повинна бути кратна 32. За один підхід реконструюється шар куба, обсягом не більше 128 МБ, тобто 1/8 одиничного куба. Якщо в системі встановлено декілька CUDA пристроїв, то реконструкція на них буде йти паралельно.

Завантаження в RAM тіньових проекцій проводиться при першому проході, в обсязі, достатньому для реконструкції куба. Тобто якщо для реконструкції куба потрібно 1200 рядків з тіньової проекції, а висота файлу 2000, то саме ці 1200 рядків з файлу завантажуються в RAM і кешуються. Оскільки при реконструкції наступного куба, швидше за все, сегменти проекції будуть перетинатися, можна буде використовувати вже завантажені дані повторно, догрузив лише нові рядки.

Виявлення серверів і розподіл завдань
Для виявлення клієнтів і серверів використовується трансляція UDP пакетів. При старті і кожні 5 секунд сервер розсилає таке повідомлення. Клієнт при старті розсилає широкомовний запит, який змушує оголосити про себе сервер негайно. Клієнт, отримавши повідомлення про присутність сервера, встановлює з ним TCP з'єднання, за яким відбувається відсилання команд управління, завдань і отримання повідомлень про завершення завдання. Всі клієнти, сервери, завдання володіють унікальним ідентифікатором.

При розриві з'єднання з одним сервером вважається, що сервер зупинився, і його завдання перерозподіляються серед інших. При підключенні нового сервера завдання від інших не відбираються і не перенаправляються на нього. Тому для рівномірного завантаження важливо виявити всі сервера до початку розподілу завдань реконструкції. По мірі готовності кубиків стає можливий запуск завдань об'єднання. Завдання спрямовуються на випадковий сервер. Задачі реконструкції та об'єднання виконуються серверами паралельно.

Завдання можуть надходити від кількох клієнтів. Сервер обробляє їх у порядку надходження (FIFO). Клієнт вважає за краще роздавати завдання реконструкції з одного шару одиничних кубиків на один сервер, щоб зменшити обсяг даних, необхідних сервера реконструкції. Об'єм оперативної пам'яті сервера визначає стратегію реконструкції та стиску кубиків. Якщо пам'яті багато (більше 20 ГБ), то стиснення попереднього кубика і реконструкція наступного здійснюються паралельно, інакше реконструкція наступного кубика не починається, поки куб не буде повністю оброблений, тобто стислий і збережено на диск.

Основні кроки реконструкції
Спочатку визначаються сегменти проекцій, які проектується реконструюється обсяг з урахуванням перспективи. На підставі цих даних визначається, які проекції потрібні для реконструкції. Якщо в RAM-кеші є сегменти файлів, що не входять в цей список, то вони видаляються з пам'яті.

Потік реконструкції звертається до кеш за сегментом проекції, блок пам'яті, що містить проекцію (з урахуванням розміру сторінки 4096 байт) замикається в пам'яті, що забороняє її вивантаження в своп, і покажчик передається в CUDA-процедуру. Всі подальші дії відбуваються в пам'яті прискорювача. Відбувається розпакування з 12 біт в float, далі всі обчислення ведуться в float. При розпакуванні сегмент обрізається ліворуч і праворуч до ширини реконструюється блоку плюс зазор, необхідний для згортки при фільтрації проекції. Відбуваються попередні обробки: інверсія, гамма-корекція, фільтрація усередненої проекцією, необхідну кількість ітерацій згладжування. Далі відбувається згортка, і сегмент ще раз обрізається ліворуч і праворуч до мінімальної ширини, необхідної при зворотному проектуванні. Так обробляються проекції, поки вистачає пам'яті на відеокарті. Тобто на відеокарті формується масив сегментів, готових для зворотного проектування.

Суть зворотного проекціювання — підсумувати для кожного вокселя зі всіх проекцій вага точки, з якої потрапляв промінь джерела рентгенівського випромінювання. Один виклик CUDA-ядер проводить розрахунок за 32 вокселям, розташованих один над одним. При цьому робиться цикл по всіх завантажених в пам'ять сегментами проекцій. Це дозволяє зберігати всі проміжні дані в регістрах і не робити проміжних записів в пам'ять. По закінченні розрахунків в пам'ять записується значення щільності вокселя. Ядро зворотного проекціювання має 2 реалізації, одна повільна з перевірками, що на проекції є точка, з якої брати дані, т. к. на кінцях проекцій не завжди є точка, звідки береться значення. Швидка реалізація виконує жодних перевірок. Потрібна реалізація вибирається автоматично. Вибірка значення проекції проводиться за алгоритмом «найближчого сусіда». Використання текстурної пам'яті на перших етапах показало різке зниження продуктивності, а на даному етапі при використанні масивів проекцій виявилося неможливим створення масивів текстур з-за обмежень CUDA 4.2.

Якщо пам'яті недостатньо, то для розрахунків може знадобитися кілька ітерацій завантаження проекцій та зворотного проекціювання.

По завершенні розрахунків робиться оцінка отриманих вокселей. Вважається максимальне і мінімальне значення щільності, будується гістограма. Побудова гістограми та копіювання результатів в RAM проводиться паралельно.

При наявності декількох CUDA пристроїв сегменти реконструюється куба обробляються паралельно. По завершенні реконструкції куба стисненням куба починає займатися інший потік, а потік реконструкції відразу береться за завдання, за умови достатньої кількості оперативної пам'яті. Як виявилося, виділення і звільнення великих сегментів пам'яті — дуже витратна за часом операція, з цього всі отримані великі блоки пам'яті використовуються повторно. Так на CUDA робиться виділення пам'яті під блок вокселей, кілька проміжних буферів, а вся інша пам'ять виділятися під проекції, і управління розміщенням проекцій в ній виконується вже вручну.

Стиснення
Стиснення в октодерево проводиться в кілька потоків. Кількість потоків залежить від кількості процесорних ядер. Кожен потік обробляє свій октант. Далі відбувається збереження у файл, який був вказаний клієнтом реконструкції в задачі, і відправляється повідомлення, що обробка кубика завершена.

Клієнт
Клієнт веде базу даних кубиків у вигляді файлу SQLite. У ньому зберігаються характеристики кубів, час реконструкції і т. д. По завершенні реконструкції програма клієнта закривається, таким чином можна організовувати пакетні завдання на реконструкцію.

Приклад коду на CUDA
Ядро (measureModel.cu), яке робить перший прохід при розрахунку гістограми. Файл наведено, щоб показати загальний вигляд коду на CUDA в реальному робочому проекті.

#include < iostream>
#include <cuda_runtime.h>
#include <cuda.h>
#include <device_functions.h>

#include "CudaUtils.h"
#include "measureModel.h"

#define nullptr NULL

struct MinAndMaxValue
{
float minValue;
float maxValue;
};
__global__ void
cudaMeasureCubeKernel(float* src, MinAndMaxValue* dst, int ySize, int sliceSize, int blocksPerGridX)
{
int x = blockDim.x * blockIdx.x + threadIdx.x;
int z = blockDim.y * blockIdx.y + threadIdx.y;
int indx = z * blocksPerGridX * blockDim.x + x;
MinAndMaxValue v;
v.minValue = src[indx];
v.maxValue = v.minValue;
//indx += sliceSize;
for (int i = 0; i < ySize; ++i)
{
float val = src[indx];
indx += sliceSize;
if (v.minValue > val)
v.minValue = val;
if (v.maxValue < val)
v.maxValue = val;
}

//поміщаємо дані по нашій рядку в поділювану пам'ять
__shared__ MinAndMaxValue tmp[16][16];
tmp[threadIdx.y][threadIdx.x] = v;
//всі нитки нашого блоку порахували свої значення
__syncthreads();

if (threadIdx.y == 0)
{//першої рядком потоків робимо згортання тимчасових даних до 16 значень
for (int i = 0; i < blockDim.y; i++)
{
if (v.minValue > tmp[i][threadIdx.x].minValue)
v.minValue = tmp[i][threadIdx.x].minValue;
if (v.maxValue < tmp[i][threadIdx.x].maxValue)
v.maxValue = tmp[i][threadIdx.x].maxValue;
}
tmp[0][threadIdx.x] = v;
}
__syncthreads();
///робимо звірку до одного значення
if (threadIdx.x == 0 && threadIdx.y == 0)
{
for (int i = 0; i < blockDim.y; i++)
{
if (v.minValue > tmp[0][i].minValue)
v.minValue = tmp[0][i].minValue;
if (v.maxValue < tmp[0][i].maxValue)
v.maxValue = tmp[0][i].maxValue;
}
dst[blockIdx.y * blocksPerGridX + blockIdx.x] = v;
}
}

void Xrmt::Cuda::CudaMeasureCube(MeasureCubeParams ¶ms)
{
//будемо робити розрахунки квадаратными блоками по 256 потоку
int threadsPerBlockX = 16;
int threadsPerBlockZ = 16;
int blocksPerGridX = (params.xSize + threadsPerBlockX - 1) / threadsPerBlockX;
int blocksPerGridZ = (params.zSize + threadsPerBlockZ - 1) / threadsPerBlockZ;
int blockCount = blocksPerGridX * blocksPerGridZ;
dim3 threadDim(threadsPerBlockX, threadsPerBlockZ);
//розмір гріду
dim3 gridDim(blocksPerGridX, blocksPerGridZ);

MinAndMaxValue *d_result = nullptr;
size_t size = blockCount * sizeof(MinAndMaxValue);
checkCudaErrors(cudaMalloc( (void**) &d_result, size));

cudaMeasureCubeKernel<<<gridDim, threadDim>>>(
params.d_src,
d_result,
params.ySize,
params.xSize * params.zSize,
blocksPerGridX
);

checkCudaErrors(cudaDeviceSynchronize());

MinAndMaxValue *result = new MinAndMaxValue[blockCount];
checkCudaErrors(cudaMemcpy(result, d_result, size, cudaMemcpyDeviceToHost));

params.maxValue = result[0].maxValue;
params.minValue = result[0].minValue;
for (int i = 1; i < blockCount; i++)
{
params.maxValue = std::max(params.maxValue, result[i].maxValue);
params.minValue = std::min(params.minValue, result[i].minValue);
}

checkCudaErrors(cudaFree(d_result));
}

__global__ void
cudaCropDensityKernel(float* src, int xSize, int ySize, int sliceSize, float minDensity, float maxDensity)
{
int x = blockDim.x * blockIdx.x + threadIdx.x;
int z = blockDim.y * blockIdx.y + threadIdx.y;
int indx = z * xSize + x;
for (int i = 0; i < ySize; ++i)
{
float val = src[indx];
if (val < minDensity)
src[indx] = 0;
else if (val > maxDensity)
src[indx] = maxDensity;
indx += sliceSize;
}
}

void Xrmt::Cuda::CudaPostProcessCube(PostProcessCubeParams ¶ms)
{
//будемо робити розрахунки квадаратными блоками по 256 потоку
int threadsPerBlockX = 16;
int threadsPerBlockZ = 16;
int blocksPerGridX = (params.xSize + threadsPerBlockX - 1) / threadsPerBlockX;
int blocksPerGridZ = (params.zSize + threadsPerBlockZ - 1) / threadsPerBlockZ;
dim3 threadDim(threadsPerBlockX, threadsPerBlockZ);
//розмір гріду
dim3 gridDim(blocksPerGridX, blocksPerGridZ);

cudaCropDensityKernel<<<gridDim, threadDim>>>(
params.d_src,
params.xSize,
params.ySize,
params.xSize * params.zSize,
params.minDensity,
params.maxDensity
);

checkCudaErrors(cudaDeviceSynchronize());
}

//на скільки частин ділимо одиницю
#define DivideCount 8
#define ValuesPerThread 32
#define SplitCount 16
#define SplitCount2 64

__global__ void
cudaBuildGistogramKernel(
float *src, 
int count, 
int *d_values,
int minValue,
int barCount
)
{
//індекс значення, починаючи з якого вважаємо
int barFrom = threadIdx.x * ValuesPerThread;
int blockIndx = SplitCount * blockIdx.x + threadIdx.y;
int blockLen = count / (SplitCount * SplitCount2);

unsigned int values[ValuesPerThread];
for (int i = 0; i < ValuesPerThread; i++)
{
values[i] = 0;
}

src += blockLen * blockIndx;
for(int i = blockLen - 1; i >= 0; i--)
{
float v = src[i];
int indx = (int)((v (float)minValue) * DivideCount) - barFrom;
if (indx >= 0 && indx < ValuesPerThread)
values[indx]++;
}

for (int i = 0; i < ValuesPerThread; i++)
{
d_values[blockIndx * barCount + barFrom + i] = values[i];
}
}

__global__ void
cudaSumGistogramKernel(
int *d_values,
int barCount
)
{
//індекс стовпчика, який підсумовуємо
int barIndx = blockIdx.x * ValuesPerThread + threadIdx.x;

unsigned int sum = 0;
d_values += barIndx;
for (int i = 0 ; i < SplitCount * SplitCount2; i++)
{
sum += d_values[i * barCount];
}

d_values[0] = sum;
}

void Xrmt::Cuda::CudaBuildGistogram(BuildGistogramParams ¶ms)
{
checkCudaErrors(cudaMemcpyAsync(
params.h_dst, params.d_src, params.count * sizeof(float), cudaMemcpyDeviceToHost, 0
));

//визначаємо мінімальне значення в гістограмі
int minValue = (int)(params.minDensity - 1);// для від'ємних значень дробова частина відкидається, а цього не достатньо!
minValue = std::max(minValue, -200);
//максимальне значення у гістограмі
int maxValue = (int)params.maxDensity;
maxValue = std::min(maxValue, 200);
//кількість стовпчиків у гістограмі (діапазон 1 ділимо на DivideCount частин)
int barCount = (maxValue - minValue + 1) * DivideCount;
///робимо кількість стовпчиків кратним ValuesPerThread
int x32 = barCount % ValuesPerThread;
barCount += (x32 == 0) ? 0 : ValuesPerThread - x32;

///налагодження потоків:
/// x - кожен потік вважає 32 стовпчика, тобто x пропорційний ширині гістограми
/// y - ділимо блок даних на SplitCount частин
dim3 threadDim(barCount / 32, SplitCount);
///розмір гріду, ділимо всі дані на SplitCount2 частин
dim3 gridDim(SplitCount2);
///таким чином всі дані ділимо по блокам на SplitCount2 * SplitCount частин
///за кожну таку частину буде обробляти barCount / 32 потоків

///резервуємо пам'ять для зберігання проміжних результатів
int *d_values = nullptr;
///для кожного блоку свої значення
size_t blockValuesCount = SplitCount2 * SplitCount * barCount;
size_t d_size = blockValuesCount * sizeof(int);
checkCudaErrors(cudaMalloc( (void**) &d_values, d_size));

///виконуємо обробку по блокам
cudaBuildGistogramKernel<<<gridDim, threadDim>>>(
params.d_src,
params.count,
d_values,
minValue,
barCount
);
//підсумовуємо блоки на GPU
cudaSumGistogramKernel<<<dim3(barCount / ValuesPerThread), dim3(ValuesPerThread)>>>(
d_values,
barCount
);
checkCudaErrors(cudaDeviceSynchronize());
//отримуємо значення гістограми
params.histogramValues.resize(barCount);
checkCudaErrors(cudaMemcpy(
¶ms.histogramValues[0], 
d_values, 
barCount * sizeof(int), 
cudaMemcpyDeviceToHost
));
checkCudaErrors(cudaFree(d_values)); 
params.histogramStep = 1.0 / DivideCount;
params.histogramStart = (float)minValue;
}

Питання №2

Камера, яку ви використовували (PIXIS-XF, треба думати), віддає картинки 2048х2048, а у статті ви пишете «до 8000х8000». Це ви отримали? Ви переміщаєте також зразок або камеру по вертикалі і горизонталі і склеює потім картинки? AndreyDmitriev

Так, у мікро томографі використовувалася камера PIXIS-XF:2048B. Розмір реконструкції по ТЗ 8000x8000x20000, що по ширині і глибині 8000 означає, що камера бачить тільки частину зображення, але реконструюється більше, це дає падіння якості реконструкції по краях. По висоті просто переміщується столик, робиться нова реконструкція, і результати склеюються.

Питання №3

Демо-зображення, які відео в кінці статті, всі отримані з 360 проекцій? Якщо так, то це добре, адже 360 проекцій з кроком у градус — досить мало, зазвичай йдуть з кроком третину або чверть градуса, інакше будуть артефакти реконструкції. Начебто є формула оптимальної кількості проекцій для заданого дозволу, але ось забув. AndreyDmitriev

У нашому випадку використовувалося 180-360 проекцій. Різні алгоритми можуть давати різну якість реконструкцій в залежності від кількості проекцій, але в загальному випадку виконується правило, що із збільшенням кількості проекцій якість зростає. Для наших завдань було достатньо використовувати 180-360 проекцій для отримання хорошої якості.

Алгоритм нормально реагує на пропуски кадрів ціною невеликого зниження дозволу в напрямку, де є пропуски. Крім того, камера мала особливість перегріватися, і у неї плавала «яскравість», тому частина кадрів відкидалась як браковані, інші нормалізувалися за середньої яскравості.

Питання №4

Ще я не дуже зрозумів про частоту камери. По специфікації вона двочастотна на 100 кілогерц або два мегагерца. Якщо у неї чотири мегапікселу, це означає, що вона віддає кадр кожні дві секунди на максимальній частоті? AndreyDmitriev

Так, ви праві.
Згідно документації камери, характеристики частоти наступні: ADC speed/bits 100 kHz/16-bit and 2 MHz/16-bit

Те, тобто якщо одна точка зображення має розрядність 16біт, то отримуємо, що частота для всього зображення буде: 2000000 / (2048*2048) = 0,47 герц

Питання №5

Ви переміщаєте маніпулятор покроково або безперервно? Скільки часу займає сканування типових зразків, що представлені на відео? AndreyDmitriev

Маніпулятор (столик) пересувається шагово. Повернули на кут, зупинили, відзняли, знову повернули і т. д.
Середній час сканування — 11 хвилин на 378 кадрів.

Питання №6

Ну і про 12 біт — дуже цікаво. Те, що відрізати чотири біта і упакувати кожні два пікселя в три байта можна для зберігання і передачі — це зрозуміло. Але для реконструювання вам же доведеться знову розгорнути кожен піксел як мінімум в два байти? Або у вас вся математика на 12-ти бітах? У такому разі як ви вирішили проблему того, що пікселі займають півтора байта і не вирівнюються на кордон? AndreyDmitriev

Формат файлів 12 бітних бічних проекцій.
  1. Пікселі побітно упаковані і розташовуються в рядки, так що дві точки займають 3 байти, при цьому передбачається, що ширина зображення проекції кратна 2, за формулою: lineWidthInBytes = width() * 3 / 2.
  2. Зазвичай відразу весь обсяг даних не завантажується, а використовуються «вікна», які підвантажується тільки необхідна порція зображення, при цьому передбачається вирівнювання лівої і правої меж кратно чотирьом пікселях.
  3. Після завантаження вікна» подальша робота ведеться за допомогою функцій повного імпорту/експорту цього вікна в масиви типу float, або quint16. В процесі імпорту/експорту і відбувається перетворення з 12 біт в потрібну розрядність.
  4. У нашому випадку виграш від скорочення розмірів файлів проекцій був більше, ніж втрата продуктивності при цьому. Крім того, так як завантажувалися файли не повністю, а тільки «вікна», призначені для обробки на конкретній машині кластера, то втрати продуктивності ще менше за рахунок розпаралелювання завдань.

Об'ємний (воксельный) візуалізація в реальному часі

Так само цікавою завданням було відображення моделі з використанням вокселей в реальному часі. Крім вокселей самої моделі потрібно було відображати у 3D такі інструменти, як обмежуючий паралелепіпед, січну площину, приховувати половину об'єкта вздовж січної площини, а також приховувати об'єкт за межами паралелепіпеда. На січній площині відображається елемент інтерфейсу користувача у вигляді «важеля» (який в той же час відповідає вектору нормалі), що дозволяє переміщати вікно зрізу проекції по січній площині і саму площину зміщувати вздовж вектора нормалі. Вікно зрізу проекції можна розтягувати за його межі і обертати. На ограничивающем паралелепіпеді крім його ребер відображаються такі ж «важелі» як на площині, які так само дозволяють змінювати його розміри. Всі ці маніпуляції виконуються у реальному часі та обсязі простору моделі.

Достатня частота кадрів при роботі забезпечувалася за рахунок динамічного рівня деталізації, так як підвантаження даних відбувається при зменшенні обмежує прямокутника автоматично.

Для маніпуляцій елементами інтерфейсу 3D були реалізовані спеціальні набори математичних функцій. Первинне створення січної площини якщо описати її в загальних поняттях без використання формул, виглядає приблизно так: користувач двічі клацає мишкою в будь-яку точку моделі, після чого ця точка фіксується, і з неї при подальшому русі мишею відображається вектор нормалі до створюваної площині. При цьому фактично площину вже створена, і вона лежить на промені, що йде з екранних координат курсору миші вглиб простору, не паралельно лініям перспективи, так що проектується промінь, як точка на екран, і орієнтація площини змінюється при переміщенні миші динамічно, і користувач бачить це відразу на екрані. Для завершення будівництва потрібно ще раз клацнути мишкою, щоб зафіксувати орієнтацію площини.

Крім елементів маніпуляції моделлю в обсязі є інструменти для управління колірними спектрами щільності і рівнями прозорості.

Безпосередньо сам рендеринг вокселей і обрізка січною площиною і обмежує параллелепипедом були реалізовані з допомогою бібліотеки VTK.

Перелік математичних функцій

Математичні класи та функції, використовувані для геометричних обчислень у проекті.
1. Клас плоскости.
  • Опис площині по нормалі і точці.
  • Опис по трьом точкам.
  • Кут повороту 2D-системи координат всередині площині навколо вектора нормалі.
  • Відступи кордонів 2D — прямокутної області всередині площині з точкою відліку на вектор нормалі.
  • Опції переходу від 2D-точки всередині площині оберненої навколо нормалі до реальних 3D-координатами і зворотне перетворення.
2. Набір математичних функцій для геометричних обчислень.
3. Математичні функції, які використовуються при реконструкції.




Більше проектів:
Як за 5233 людино-години створити софт для микротомографа
SDK для підтримки впровадження електронних книг у форматі FB2
«Сфера»: як моніторити мільярди кіловат-годин
Управління доступом до електронних документів. Від DefView до Vivaldi
Розробка простого плагіна для JIRA для роботи з базою даних
На допомогу DevOps: складальник прошивок для мережевих пристроїв на Debian за 1008 годин
Автооновлення служби Windows через AWS для бідних
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.