Наименьшие квадраты в общем виде для задач, Комп-Томографии

[B@R5uk]

Geen, ну получаются результаты с точностью до множителя, ну и что? Конечному пользователю главное относительное распределение плотности в объекте; он ещё потом её подкрутит, что бы контрастность в картинках слоёв нарезки была по-лучше.

-- 15.12.2021, 13:06 --

А вообще, я так понял, проблема ТС в гигантском объёме вычислений. Ну так бороться с ней! Подход к проблеме универсальный: всё, что раньше было неизвестными, заменяется какой-нибудь подходящей моделью (линейной разумеется) с меньшим числом неизвестных. Я не знаю как сейчас ТС делает (он постоянно во всём темнит, как будто его работа национальный секрет, а инфу приходится клещами тянуть): то ли у него простая кубическая/параллелепипедная сетка (на основе кубических вокселей), то ли что-то более серьёзное типа пространственной треугольной тетраэдрической сетки. Если первое, то стоит перейти ко второму. Такая сетка позволит сгущать точки разбиения (читай: искомые величины) там, где плотность объекта исследования испытывает резкие скачки, а там где плотности нет или она равномерна, сетку можно сделать значительно реже. Вычислений по интерполяции, трансформации и интегрированию, конечно, добавится. Тут трилинейной фильтрацией и простым суммированием не обойтись. Однако, это скорее изнурительная работа для программиста, а не для компьютера. Потому что усложнённые вычисления делаются один раз перед оптимизацией после каждого изменения сетки для расчёта матриц $Z_k$ минимизируемой величины

$$M=\sum_k(F^TZ_kG-h_k)^2$$

За счёт меньшего числа неизвестных ускоряется самая трудоёмкая часть вычислений, а точность в результате будет та же.

[Евгений Машеров]

Это если две точки "с точностью до множителя". А если больше - "до неизвестного преобразования". Мне отчего-то видится главной проблемой не объём счёта, а неопределённость, которая на практике выльется в численную неустойчивость.
А объём - опираясь на прочный фундамент незнания того, что, собственно, ТС надо, предложу как-то через Фурье считать...

[ilghiz]

Спасибо большое всем за интересные советы, комментарии и обсужение!

Я тоже не уверен, что только с точностью до множителя. Если до множителя - то задача поправима - нормируем например $f(x)$ и решаем.

Для простоты, положим КТ, объект - распределение плотности поглащения излучения внутри человека.

Неизвестная $g(x)$ - это то, на чем происходит поглащение.

Я бы с радостью взял бы что-то отличное от равномерной кубической дискретизациии, но об оъекте ничего заранее неизвестно. Известно, только, что внутри надо иметь разрешение хотя бы 1мм, а объект имеет размеры десятки литров, вот тут-то и вылезают десятки и сотни миллионов неизвестных.

Пока дискретизировал на кусочно-постоянных (для простоты). Планирую перейти на что-то более точное, но пока еще не решил, но, предполагаю, что сетку трогать нельзя.

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

Как я говорил, в идеале, $g(x)$ - должна быть в идеале бесконечно тонкая поверхность, и совсем в идеале - близкая к плоскости. Но это не так. Во-первых - она имеет форму куска эллипса, во-вторых, имеет ненулевую толщину. Если ее на той же сетке дискретизовать и посчитать только то, что имеет скажем, 10% от максимума, то там будет примерно в 100 раз меньше конечных элементов, чем если на той же сетке дискретизовать $f(x)$.

Хочется Фурье, или чего-то в этом роде, но нет понимания как.

[ilghiz]

B@R5uk в сообщении #1543021 писал(а):

(он постоянно во всём темнит, как будто его работа национальный секрет, а инфу приходится клещами тянуть)

если вдруг такое мнение сложилось, очень сильно прошу прощения, и реально пытаюсь рассказать все, что знаю, другое дело, не всегда во время. Железо, которое эти данные снимает - свое и самопально разработано и довольно хорошо протестированно. Ожидалось получить задачу вида КТ, но, судьба не сложилась, и приходится уточнять вот таким зверским методом.

[B@R5uk]

Можно ещё воспользоваться следующим приближённым методом. Сначала делается грубая оценка искомых величин (например, с помощью расчёта со значительно более крупными пространственными ячейками). Затем производится серия расчётов, но не во всём пространстве, а лишь в небольшой области. При этом вне этой области искомые величины принимаются известными и равными грубой оценке. По мере просчёта небольших областей результаты этих обсчётов обновляют содержимое искомых величин во всём пространстве. А обсчитываемые области одна за другой покрывают всё рабочее пространство с небольшим перекрытием. Таким образом всё искомое пространство проходится несколько раз всё лучше и лучше уточняя искомое решение.

Это должно сходится по той причине, что флуктуации в искомых величинах из соседних областей влияют на решение в текущей слабо, а из удалённых — пренебрежимо мало. Этот факт легко проверить по величине удалённых от диагонали членов обратной матрицы квадратичной формы, соответствующих перекрытию искомых величин из разных областей. Они должны быть много меньше (порядка желаемой точности) диагональных значений. Ускорение по времени этот подход даст по той причине, что методы счёта у вас явно сложнее, чем $O(N)$. Проблема как всегда будет в поиске $g(x)$, так как эта функция должна оцениваться по всему пространству измерений, а тут у будут небольшие его кусочки обсчитываться. Ну, надо будет что-нибудь придумать. Или воспользоваться моим предыдущим советом на счёт эталонного объекта в рабочем пространстве (чтобы сразу откалиброваться).

Более продвинутым вариантом этого подхода будет следующая модернизация. Рассматриваемая область обсчитывается на мелкой (желаемой) сетке, соседние области — на грубой сетке, а содержимое дальних областей считается известным. Результат сохраняется только для мелкой сетки. Остальное — всё как написано выше.

[ilghiz]

Спасибо большое, B@R5uk за Ваши советы!!!

Полностью с Вами согласен, что многосеточные методы и их разумное применение здесь должны дать существенное ускорение.

Пока успел поэкспериментировать с набором сеток с шагом дискретизации от 1см до 1мм. При 1см у меня имеется около 15 тысяч конечных элементов, а при 1мм - 12 миллионов. Чтобы надежнее сходится я начинаю решение на самой грубой сетке и далее мельчу ее. Если переносить решение с грубой сетки на мелкую, то время счета на каждой сетке слегка более, чем линейно зависит от числа конечных элементов, то есть на четке 1см решение получается за 20 секунд, а на сетке 1мм - за пол дня. Я перехожу с сетки на сетку деля шаг в 1.25 раз, то есть 1см, 8мм, 6.4мм, ...

Пока еще с локальным сгущением еще не игрался, но, повидимому, это основной и наиболее понятный способ куда можно двигаться.

По физике процесса - мне также "светит" уточнение решения. То есть вначале решаем задачу, когда и $g$ и $f$ - неизвестны, а потом начали собирать данные, чтобы уточнить $f$. То есть по сути получили какое-то решение для $f$, но юзер захотел уточнить какой-то регион и начал в нем снимать данные, и тогда там надо начинать уточнять решение. Тут конечно все в первого взгляда выглядит линейно, но если во время первой фазы, когда мы решали совместно $g$ и $f$ решение для $g$ было не точным, то вторая фаза не будет получаться.

[Евгений Машеров]

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

[ilghiz]

Спасибо большое, Евгений Машеров за вопрос!

Евгений Машеров в сообщении #1543137 писал(а):

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

В моем случае, это переносная МРТшка (не гиперфайновская тележка, а именно карманная переносная МРТшка). $g(x)$ - это поверхность с постоянным магнитным полем, а ее интенсивность определяется тем, как туда сфокусировали возбуждение. А отклик - просто интегральный, со всей поверхности, как получилось. Поверхность можно двигать, поворачивать, крутить, по-разному располагая в исследуемом объекте.

В моей постановке $f(x) \in {\mathbb R^3}$, а не $f(x) \in {\mathbb R}$, как я озвучил в исходной постановке, но, вроде это сильно сути не меняет, просто слегка усложняет постановку и решение.

Поверхность в переносном варианте плывет из-за того, что магниты не стабилизированы по температуре. Можно конечно затабулировать все температурные коэффициенты, но не хочу, так как тогда надо вкрячить с десяток термосенсоров и как-то все это затабулировать, а, как человек много и долго занимавшийся МРТ и ЯМРом, понимаю, что это будет не надежное решение, которое будет регулярно давать сбой, поэтому хочу решать с неизвестной $g(x)$, которую находить в процессе съема данных. Примерная форма $g(x)$ мне заранее известна, но, как я говорил, она может отклоняться на сантиметр или, совсем в военное время, на два сантиметра.

[Евгений Машеров]

А как-то калибровать нельзя? Поместив в известные точки маркеры, дающие интенсивный ответ? (Где-то было описано, что в качестве таких маркеров использовались капсулы с витамином D, не ради витамина, а ради жирового растворителя). И по ним калибровать?

[ilghiz]

Евгений Машеров в сообщении #1543254 писал(а):

А как-то калибровать нельзя? Поместив в известные точки маркеры, дающие интенсивный ответ? (Где-то было описано, что в качестве таких маркеров использовались капсулы с витамином D, не ради витамина, а ради жирового растворителя). И по ним калибровать?

Это для МЧС и скорой помощи. Грубо говоря, у тебя есть 2 минуты на скан и ты с этой МРТшкой с парашута спрыгнул во время спасательной операции. Понятно про 1мм точность там речь не идет, но на 3-4мм точность и небольшой объект - и данных за минуту собрать можно и смартфоном обработать, а вот поле - да, не калибруется. По крайней мере при этих усовиях.

В стационарных ситуациях там все можно калибровать и прикрутить термоэлементы на магниты, чтобы поле не плыло. Там проблем такого рода нет, вернее конечно на стадии изготовления и заводской калибровки нужно решать похожую задачу, но однажды. Но тут все методы хороши и подложить заранее известную болванку, для которой есть $f(x)$ очень просто.

[Евгений Машеров]

Ну, я имел в виду нечто вроде - перед снятием изображения на голову надевается нечто вроде шапочки с "яркими" объектами, скажем, инион, назион, вертекс, левый и правый висок. Делается снимок и поправляется, исходя из того, куда ушли отметки этих точек.

[ilghiz]

Евгений Машеров в сообщении #1543328 писал(а):

Ну, я имел в виду нечто вроде - перед снятием изображения на голову надевается нечто вроде шапочки с "яркими" объектами, скажем, инион, назион, вертекс, левый и правый висок. Делается снимок и поправляется, исходя из того, куда ушли отметки этих точек.

Не, так не получится.
Во-первых, форма поля - кривая, и, чтобы ее восстановить надо иметь известный и довольно большой объект. Размер аппаратуры - планшет (это и магниты, и электроника, и катушки). Для калибровки нужен объект с чемоданчик, то есть в разы больше, чем этот планшет, и в нем по всему объему надо иметь много водородосодержащих вещест. Если навесить что-то на тело, то одной точкой-двумя-тремя точками решить сложно. Ну будет у меня в $f(x)$ что-то извесно, а остальное - нет. Мне это как мертвому припарка, проще решать сразу в общем виде.

Потом те же МЧСы хотят, чтобы на тело ничего не навешивали, тело и так от боли корчится, а если на него что-то надеть, то совсем плохо получится. Эта технология в общем замена рентгену, плюс переносная, в отличии от рентгена. Люди хотят именно быстро и без касания картинку на месте. И, если поле не уплыло, картинка за пару минут получается (пока на лаптопе, но на смартфон перенести тоже можно, современные процессоры смартфонов 50 гигафлоп/с на раз показывают).

То есть я к тому, что мы уже третью страницу обсуждаем как заранее $g(x)$ затабулировать, но, у нас довольно большая группа медиков, химиков и физиков не смогла найти удобное решение как сдалать такую табуляцию в полевых условиях, поэтому единственно к чему я и возвращаюсь - надо искать решение в общем виде.

[Евгений Машеров]

И сколько Тесла получается на такой компактной аппаратуре?

[ilghiz]

200-400 МГц, но не ядерного, а электронного резонанса. А в теслах - там очень мало получается - милитеслы. Для электронного резонанса больше и не надо, электромагнитное излучение на частотах существенно больше 500МГц уже почти внутрь тела не проникает.

[Евгений Машеров]

Ну, могу лишь пожелать хорошего снятия сигнала ЭПР. Критиковать не могу, не настолько в теме.
Но вот обработка - без наличия калибровки по известным точкам изображение будет заведомо плыть.