Моделирование процесса получения гартманнограмм деформированного зеркала
Для проверки точности того или иного метода восстановления профиля зеркала по
гартманнограммам удобно воспользоваться моделированием гартманнограмм для
эталонной модели искаженного профиля зеркала.
Для моделирования процесса регистрации пред- и зафокальных гартманнограмм на телескопе БТА использован следующий алгоритм:
Моделирование поверхности деформированного зеркала
Поверхность зеркала моделируется в виде распределения координат поверхности зеркала и нормалей к ней по равномерной сетке. Параметры зеркала (фокусное расстояние и диаметр) задаются пользователем. Также задается размер итоговой квадратной матрицы, в которую будут занесены координаты точек на поверхности зеркала и векторы нормалей.
Нулевая поверхность зеркала строится на основе уравнения параболы: z=x²+y²/(4F), где x,y,z - декартовы координаты, F - фокусное расстояние зеркала. Первоначальное приближение строится на сетке 255x255 с таким масштабом, чтобы итоговый размер модели совпадал с диаметром зеркала.
Искажения формы зеркала (в микрометрах) задаются заранее предопределенной матрицей произвольного размера. Для коррекции нулевой поверхности зеркала на искажения матрица искажений масштабируется к матрице зеркала посредством интерполяции сплайнами. Далее результирующие две матрицы складываются, в результате чего получается грубое приближение к форме искаженного зеркала.
Масштабирование полученной матрицы координат поверхности зеркала на заданную пользователем сетку также осуществляется посредством интерполяции сплайнами. После получения матрицы координат поверхности зеркала рассчитывается матрица нормалей к поверхности зеркала по той же координатной сетке. Нормали рассчитываются простейшей линейной интерполяцией. В результате двойной интерполяции моделируемая поверхность зеркала получается сглаженной, наиболее приближенной к реальным условиям.
Моделирование процесса отражения
Отражение поверхностью деформированного зеркала падающего на него плоского волнового фронта будем моделировать также посредством дискретного приближения.
Рассмотрим сформированную на предыдущем этапе поверхность зеркала как набор бесконечно малых плоских зеркал, расположенных в точках (x,y,z) с нормалями n. При падении луча света с направляющим вектором f (см. рисунок выше) на площадку AB с нормалью n направляющий вектор отраженного луча r антисимметричен падающему. Из чертежа видно, что для того, чтобы получить значение r, достаточно к f прибавить удвоенное значение вектора k, являющегося произведением вектора нормали на абсолютную величину скалярного произведения нормали и направляющего вектора падающего луча: k = n·|f·n|. Так как вектор нормали и вектор падающего луча направлены противоположно, последнее выражение можно записать так: k = n·(-f·n).
Итак, направляющий вектор отраженного луча вычисляется по формуле r = f - 2n·(f·n). В дальнейшем для трассировки отраженных лучей через каждую точку сетки на поверхности зеркала строится прямая в соответствии с направляющими векторами.
Моделирование диафрагмы
Моделирование диафрагмы Гартманна производится в соответствии с параметрами реальной диафрагмы телескопа БТА. Диафрагма представляет собой набор отверстий, расположенных на пересечении концентрических окружностей с радиусами 175, 247, 295, 340, 379, 414, 448 и 478 миллиметров с 32 лучами, проведенными из центра диафрагмы с шагом π/16, начиная с угла π/32. Диаметр каждого отверстия - 15мм. Для отождествления ориентации гартманнограммы на снимке, помимо основных 256 отверстий на диафрагме присутствует два отверстия-маркера.
На рисунке сверху приведен негатив маски, соответствующей модели диафрагмы Гартманна. Модель представляет собой квадратную матрицу, размер которой задается пользователем. Матрица строится в соответствии с реальной гартманнограммой, диаметр которой, 1.2м, масштабируется на всю матрицу.
Соответствующее маске изображение квантуется с учетом масштаба на матрице, в ячейки матрицы, совпадающие с отверстиями, заносятся единицы, остальные ячейки заполняются нулями. В результате моделирования форма отверстий будет тем ближе к окружности, чем большим будет размер матрицы маски.
Для упрощения построения матрицы маски в модели этот процесс происходит в два этапа: сначала единицами помечаются центры окружностей, затем полученная матрица сворачивается с маской, характеризующей округлые отверстия на данной матрице.
Погрешности положения и формы отверстий маски, возникающие вследствие ошибок квантования, позволяют приблизить условия реалистичности моделирования (отверстия на маске Гартманна расположены с определенной точностью, кроме того, их положение не постоянно, т.к. положение створок стакана БТА, на которые и нанесена диафрагма, не фиксируется).
Моделирование прохождения отраженного света сквозь маску и регистрации его светоприемником
Моделирование процесса регистрации гартманнограммы заключается в трассировке каждого луча, отраженного элементом зеркала с координатами (x,y), на плоскость изображения с учетом формы маски и квантования на элементах светоприемника.
Координаты точки пересечения лучом, исходящим из точки R₀ и имеющим направляющий вектор r, плоскости z=z₁ будут равны R₁ = R₀ + r(z₁-z₀)(rz), где zi - координата z точки Ri, rz - z-проекция вектора r.
Вычисляя координаты пересечения каждым лучом плоскости, в которой расположена диафрагма, и квантуя координаты по ячейкам маски диафрагмы, получим соответствующие индексы для матрицы маски:
Mi=round(Si/2+i/s)
(если адресовать матрицу с единицы, а не ноля), где i - координата x или y, s - масштаб на маске (м/пикс), S - размер матрицы маски диафрагмы, round(x) - округление x. Из найденных индексов Mi удаляются выходящие за пределы матрицы маски (т.е. индексы Mi< и Mi>Si). Удаление индексов может производиться путем замены их на 0.
Аналогичным образом находятся индексы для матрицы, соответствующей светоприемнику. В результате для матрицы направляющих векторов отраженных лучей будут получены две матрицы индексов. Луч, для которого существует индекс в обеих матрицах, попал бы на светоприемник, если бы диафрагма была полностью прозрачной. Оценить виньетирование при полностью открытой маске диафрагмы (маска имеет прямоугольную форму!) можно, разделив количество попавших на светоприемник лучей к общему количеству отраженных лучей.
Проходя в цикле по всем прошедшим лучам, будем инкрементировать содержимое тех ячеек матрицы изображения, для которых соответствующее значение матрицы маски равно единице.
При моделировании приходится идти на компромисс между затратами ресурсов (времени вычисления и памяти) о точностью модели. Чем меньше будут размеры всех матриц в модели, тем меньшим будет потребление ресурсов, однако тем хуже будет точность. В то же время очень больших точностей достичь невозможно ввиду слишком высоких требований к объему оперативной памяти и слишком сильно возрастающим временем вычислений.
Следует учесть, что при рассогласовании размеров матрицы «зеркала» и матрицы «светоприемника» изображения пятен имеют разрывы (как на рис. сверху слева, где матрица «зеркала» имела размер 4000x4000, а матрица «светоприемника» - 3000x3000). Чтобы избавиться от этих разрывов, необходимо либо увеличить размер матрицы «зеркала» (однако, при этом очень сильно возрастает ресурсоемкость), либо уменьшить размер матрицы «светоприемника» (что приведет уменьшению точности модели), либо сгладить изображение при помощи простого гауссова фильтра (рис. сверху справа) или иного фильтра.
Результаты моделирования представлены на рис. снизу. Слева на рисунке изображена полученная модель предфокальной гартманнограммы при положении светоприемника в точке z=23.9м с использованием гауссова сглаживания для компенсации шумов дискретизации. Справа на рисунке - модель отклонений поверхности зеркала (в микрометрах), для которой получена данная модель. Размеры матриц модели были следующими: матрица зеркала имела размер 4000x4000, матрица-маска диафрагмы - 1200x1200, матрица светоприемника - 3000x3000 с размером пикселя 12.5мкм (для наилучшего приближения к реально используемому светоприемнику).
Оценку влияния отклонения поверхности зеркала от идеальной на результирующее изображение можно по приведенным ниже изображениям для моделей «плохого» (отклонения поверхности до 16нм) и «очень плохого» (отклонения поверхности до 75нм) зеркала можно по приведенным ниже рисункам.
Искажения «плохого» зеркала. Изолинии на карте поверхности проведены через 1нм.
Предфокальная гартманограмма для «плохого» зеркала
Искажения «очень плохого» зеркала. Изолинии на карте поверхности проведены через 10нм.
Предфокальная гартманограмма для «очень плохого» зеркала Исходники
function [Pre Post zn] = prepare_all() % % [Pre Post zn] = prepare_all() % формирует массивы Pre и Post - гартманнограммы, которые получились бы на БТА % при деформации зеркала zn (из build_mirror) % при положении светоприемника в точках z=24.148м и z=23.9м % % строим зеркало: матрица 4Kx4K [nx ny nz x y z zn] = build_mirror(4000,6.05,24.024); falling = single([0 0 -1]); % свет падает нормально [rx ry rz] = reflection(nx,ny,nz, falling); % получаем отраженные лучи % строим маску, 1200x1200 пикселей [mask scale] = make_mask(1200); % получаем "гартманнограммы", 3Kx3K Post = build_image(mask, 20.017, scale, 24.148, 12.5e-6, 3000, rx,ry,rz,x,y,z); Pre = build_image(mask, 20.017, scale, 23.9, 12.5e-6, 3000, rx,ry,rz,x,y,z); % выполняем фильтрацию для сглаживания неоднородностей f = fspecial('gaussian', 9, 3); Post = imfilter(Post, f); Pre = imfilter(Pre, f); endfunction
function [nx ny nz x y z znoice] = build_mirror(S, D, F, zmask) % build_mirror(S, R, F) % S - размер выходной квадратной матрицы % D - диаметр зеркала в метрах % F - фокусное расстояние зеркала % zmask - матрица отклонений поверхности зеркала от параболы (в метрах), если пропущена - берется значение по умолчанию % на выходе имеем матрицы % nx, ny,nz - матрицы SxS проекций нормалей к поверхности зеркала в точках (xi,yi) % z - такая же матрица высот (от центра зеркала) % x,y - матрицы координат % (x-xc)^2+(y-yc)^2 = 4fz % R = D/2; % радиус зеркала scale0 = D / 255; % масштаб для "заготовки" зеркала (256х256) printf('Making mirror prototype\n'); fflush(stdout); [x0, y0] = meshgrid([-127.5:127.5]*scale0); % сетка для заготовки z = (x0.^2 + y0.^2)/4/F; % подготавливаем заготовку if(nargin == 3) zmask = [ 0 1 2 3 2 1 0; 1 2 2 3 3 2 1; 2 3 2 2 2 3 2; 2 3 1 0 2 3 3; 2 3 2 3 3 1 0; 1 2 3 1 1 0 0; 0 1 2 0 1 1 0; ] * 2e-6; % отклонения в м endif tmp = size(zmask,1); A = tmp / 2 - 0.5; [x1, y1] = meshgrid([-A:A]*D/(tmp-1)); % сетка для маски отклонений printf('Interpolating deformation mask\n'); fflush(stdout); znoice = interp2(x1,y1,zmask,x0,y0,'spline'); % сглаживаем маску отклонений z = z + znoice; % добавляем искривление поверхности znoice(find((x0.^2 + y0.^2)>R^2+scale0^2)) = nan; scale = D / (S-1); % шаг по координате A = S / 2 - 0.5; % амплитуда (в пикселях) [x, y] = meshgrid([-A:A]*scale); % сетка координат на зеркале x = single(x); y = single(y); printf('Interpolating mirror matrix to %dx%d\n',S,S); fflush(stdout); z = single(interp2(x0,y0,z,x,y,'spline')); % интерполяцией получаем новые значения printf('Compute normals\n'); fflush(stdout); % вычисляем нормали: printf('nx ... '); fflush(stdout); tmp = single(zeros(size(z))); tmp1= single(zeros(size(z))); tmp(:,[2:end]) = z(:,[1:end-1]); % получаем сдвинутые матрицы (для nx) tmp1(:,[1:end-1]) = z(:,[2:end]); nx = (tmp - tmp1) / 2 / scale; % координата x ненормированной нормали printf('ny ... '); fflush(stdout); tmp = single(zeros(size(z))); tmp1= single(zeros(size(z))); tmp([2:end],:) = z([1:end-1],:); % получаем сдвинутую вниз матрицу (для ny) tmp1([1:end-1],:) = z([2:end],:); ny = (tmp - tmp1) / 2 / scale; % координата y ненормированной нормали % нормируем нормали printf('norm ... '); fflush(stdout); amp = sqrt(nx.^2+ny.^2+1); % длина нормали nx = nx ./ amp; ny = ny ./ amp; nz = ones(size(nx)) ./ amp; printf('remove outern ... '); fflush(stdout); idx = find((x.^2+y.^2)>R^2+scale^2); % находим ячейки, где зеркала уже нет nx(idx) = nan; ny(idx) = nan; % и удаляем несуществующие nz = sqrt(1-nx.^2-ny.^2); z(idx) = nan; % выбрасываем значения z в этих ячейках printf('OK!\n'); endfunction
function [rx ry rz] = reflection(nx, ny, nz, f) % reflection(nx, ny, nz, fx, fy, fz) % получение координат вектора отраженного луча % nx,ny,nz - нормали к отражающей поверхности % f == [fx fy fz] - вектор, характеризующий падающий пучок света % rx,ry,rz - векторное поле отраженных лучей с нормалями nx,ny,nz % ВСЕ ВЕКТОРЫ - НОРМИРОВАННЫЕ!!! fx = f(1); fy = f(2); fz = f(3); % знак "-" обусловлен противоположной направленностью падающего луча и нормали koeff = single(-2*(fx.*nx + fy.*ny + fz.*nz)); % удвоенная проекция падающего луча на нормаль % при отражении к падающему вектору надо прибавить удвоенное произведение % нормали на модуль проекции падающего вектора на нормаль, чтобы получить % отраженный вектор rx = koeff.*nx + fx; ry = koeff.*ny + fy; rz = koeff.*nz + fz; endfunction
function [mask scale] = make_mask(SS) % построение гартмановской маски % SS - размер маски R = [175 247 295 340 379 414 448 478] * 1e-3; % радиусы колец на гартманограмме HoleR = 7.5e-3; % радиус отверстий - 7.5мм R0 = .6; % радиус самой гартманограммы alpha0 = pi/32; % смещение лучей относительно секущих горизонтали/вертикали mask = logical(zeros(SS)); % маска гартманограммы scale = 0.6 / (SS / 2); % масштаб на маске (м на пиксель) Angles = [0:31] * 2 * alpha0 + alpha0; % углы, по которым располагаются лучи % для того, чтобы разместить на маске окружности, создадим маску % окружности: zeros(15) с единицами там, где должна быть дырка. Затем % пометим единицами в mask те точки, куда должен попадать левый верхний % край маски с дыркой и сделаем свертку (conv2 или xcorr2) HS = ceil(15e-3/scale); % размер матрицы "дырок" hole = logical(zeros(HS)); % !!! КООРДИНАТА Y ПЕРЕВЕРНУТА "ВВЕРХ НОГАМИ" !!! [X Y] = meshgrid(([1:HS]-ceil(HS/2))*scale); % координаты в матрице "дырок" dX = -X(1,1); % смещение относительно центра для верхнего левого угла матрицы дырок hole(find(X.^2+Y.^2 <= HoleR^2)) = 1; %[X Y] = meshgrid(([1:SS]-ceil(SS/2))*scale); % координаты в матрице гартманограммы for i = [1 : size(R,2)] % цикл по кольцам x = round(R(i) * cos(Angles)/scale + SS/2); y = round(-R(i) * sin(Angles)/scale + SS/2); mask((y-1)*SS+x) = 1; % помечаем положения матриц дырок на матрице гартманограммы end % помечаем маркеры x = round(R([4 8]) .* cos([-2 -7]* 2 * alpha0)/scale + SS/2); y = round(-R([4 8]) .* sin([-2 -7]* 2 * alpha0)/scale + SS/2); mask((y-1)*SS+x) = 1; tmp = xcorr2(mask, hole); mask = logical(tmp([1:SS], [1:SS])); endfunction
function ima = build_image(mask, maskz, maskS, imagez, imageS, imageN, rx,ry,rz, x0,y0,z0) % build_image(mask, maskz, maskS, imagez, imageN, rx,ry,rz, x0,y0,z0) % построение изображения через маску % ima - результирующее изображение % mask - булева двумерная матрица маски (1 - свет проходит, 0 - нет) % maskz - координата Z маски (от вершины зеркала) % maskS - масштаб в плоскости маски (м/пикс) % imagez - координата Z изображения % imageS - масштаб в плоскости изображения (м/пикс) % imageN - длина стороны изображения (пикс) % rx,ry,rz - координаты направляющих векторов отраженных от зеркала лучей % x0,y0,z0 - координаты точек на поверхности зеркала, соотв. направл. векторам % ima = int16(zeros(imageN)); k = (maskz - z0) ./ rz; % коэффициент увеличения вектора M_0 ==> M в плоскости маски x = x0 + rx.*k; % координаты лучей в плоскости маски y = y0 + ry.*k; idxX = int16(round(size(mask,2)/2+x./maskS)); % индексы на маске idxY = int16(round(size(mask,1)/2+y./maskS)); % удаляем индексы, выхоодящие за пределы размера маски % (для int nan переходит в 0) tmp = find(idxX < 1 | idxY < 1 | idxX > size(mask, 2) | idxY > size(mask, 1)); idxX(tmp) = 0; % для idxY не делаем обнуления, т.к. проверять будем лишь idxX k = (imagez - z0) ./ rz; % координаты лучей в плоскости изображения => пикс. x = int16(round(imageN/2+(x0 + rx.*k)./imageS)); y = int16(round(imageN/2+(y0 + ry.*k)./imageS)); % чистим выходящее за пределы tmp = find(x < 1 | y < 1 | x > imageN | y > imageN); x(tmp) = 0; % ищем области, через которые лучи пройдут tmp = find(idxX & x); NptsTotal = size(find(~isnan(rx)),1); printf('\n%g%% of rays falls onto the image (within transparent mask)\n', 100*size(tmp,1)/NptsTotal); printf('Start checking %d points\n', size(tmp,1)); fflush(stdout); for j = 1:size(tmp,1) i = tmp(j); if(mask(idxY(i), idxX(i)) == 1) ++ima(y(i), x(i)); % отмечаем на изображении прошедший луч end end printf('\n%g%% of rays falls onto the image\n', 100*sum(sum(ima))/NptsTotal); endfunction
Для моделирования процесса регистрации пред- и зафокальных гартманнограмм на телескопе БТА использован следующий алгоритм:
- моделирование поверхности главного зеркала БТА с предопределенными искажениями;
- моделирование процесса отражения деформированным зеркалом коллимированного пучка света;
- построение модели маски Гартманна;
- моделирование прохождения отраженного деформированным зеркалом света через маску Гартманна с последующей регистрацией его идеальным светоприемником.
Моделирование поверхности деформированного зеркала
Поверхность зеркала моделируется в виде распределения координат поверхности зеркала и нормалей к ней по равномерной сетке. Параметры зеркала (фокусное расстояние и диаметр) задаются пользователем. Также задается размер итоговой квадратной матрицы, в которую будут занесены координаты точек на поверхности зеркала и векторы нормалей.
Нулевая поверхность зеркала строится на основе уравнения параболы: z=x²+y²/(4F), где x,y,z - декартовы координаты, F - фокусное расстояние зеркала. Первоначальное приближение строится на сетке 255x255 с таким масштабом, чтобы итоговый размер модели совпадал с диаметром зеркала.
Искажения формы зеркала (в микрометрах) задаются заранее предопределенной матрицей произвольного размера. Для коррекции нулевой поверхности зеркала на искажения матрица искажений масштабируется к матрице зеркала посредством интерполяции сплайнами. Далее результирующие две матрицы складываются, в результате чего получается грубое приближение к форме искаженного зеркала.
Масштабирование полученной матрицы координат поверхности зеркала на заданную пользователем сетку также осуществляется посредством интерполяции сплайнами. После получения матрицы координат поверхности зеркала рассчитывается матрица нормалей к поверхности зеркала по той же координатной сетке. Нормали рассчитываются простейшей линейной интерполяцией. В результате двойной интерполяции моделируемая поверхность зеркала получается сглаженной, наиболее приближенной к реальным условиям.
Моделирование процесса отражения
Отражение поверхностью деформированного зеркала падающего на него плоского волнового фронта будем моделировать также посредством дискретного приближения.
Рассмотрим сформированную на предыдущем этапе поверхность зеркала как набор бесконечно малых плоских зеркал, расположенных в точках (x,y,z) с нормалями n. При падении луча света с направляющим вектором f (см. рисунок выше) на площадку AB с нормалью n направляющий вектор отраженного луча r антисимметричен падающему. Из чертежа видно, что для того, чтобы получить значение r, достаточно к f прибавить удвоенное значение вектора k, являющегося произведением вектора нормали на абсолютную величину скалярного произведения нормали и направляющего вектора падающего луча: k = n·|f·n|. Так как вектор нормали и вектор падающего луча направлены противоположно, последнее выражение можно записать так: k = n·(-f·n).
Итак, направляющий вектор отраженного луча вычисляется по формуле r = f - 2n·(f·n). В дальнейшем для трассировки отраженных лучей через каждую точку сетки на поверхности зеркала строится прямая в соответствии с направляющими векторами.
Моделирование диафрагмы
Моделирование диафрагмы Гартманна производится в соответствии с параметрами реальной диафрагмы телескопа БТА. Диафрагма представляет собой набор отверстий, расположенных на пересечении концентрических окружностей с радиусами 175, 247, 295, 340, 379, 414, 448 и 478 миллиметров с 32 лучами, проведенными из центра диафрагмы с шагом π/16, начиная с угла π/32. Диаметр каждого отверстия - 15мм. Для отождествления ориентации гартманнограммы на снимке, помимо основных 256 отверстий на диафрагме присутствует два отверстия-маркера.
На рисунке сверху приведен негатив маски, соответствующей модели диафрагмы Гартманна. Модель представляет собой квадратную матрицу, размер которой задается пользователем. Матрица строится в соответствии с реальной гартманнограммой, диаметр которой, 1.2м, масштабируется на всю матрицу.
Соответствующее маске изображение квантуется с учетом масштаба на матрице, в ячейки матрицы, совпадающие с отверстиями, заносятся единицы, остальные ячейки заполняются нулями. В результате моделирования форма отверстий будет тем ближе к окружности, чем большим будет размер матрицы маски.
Для упрощения построения матрицы маски в модели этот процесс происходит в два этапа: сначала единицами помечаются центры окружностей, затем полученная матрица сворачивается с маской, характеризующей округлые отверстия на данной матрице.
Погрешности положения и формы отверстий маски, возникающие вследствие ошибок квантования, позволяют приблизить условия реалистичности моделирования (отверстия на маске Гартманна расположены с определенной точностью, кроме того, их положение не постоянно, т.к. положение створок стакана БТА, на которые и нанесена диафрагма, не фиксируется).
Моделирование прохождения отраженного света сквозь маску и регистрации его светоприемником
Моделирование процесса регистрации гартманнограммы заключается в трассировке каждого луча, отраженного элементом зеркала с координатами (x,y), на плоскость изображения с учетом формы маски и квантования на элементах светоприемника.
Координаты точки пересечения лучом, исходящим из точки R₀ и имеющим направляющий вектор r, плоскости z=z₁ будут равны R₁ = R₀ + r(z₁-z₀)(rz), где zi - координата z точки Ri, rz - z-проекция вектора r.
Вычисляя координаты пересечения каждым лучом плоскости, в которой расположена диафрагма, и квантуя координаты по ячейкам маски диафрагмы, получим соответствующие индексы для матрицы маски:
Mi=round(Si/2+i/s)
(если адресовать матрицу с единицы, а не ноля), где i - координата x или y, s - масштаб на маске (м/пикс), S - размер матрицы маски диафрагмы, round(x) - округление x. Из найденных индексов Mi удаляются выходящие за пределы матрицы маски (т.е. индексы Mi< и Mi>Si). Удаление индексов может производиться путем замены их на 0.
Аналогичным образом находятся индексы для матрицы, соответствующей светоприемнику. В результате для матрицы направляющих векторов отраженных лучей будут получены две матрицы индексов. Луч, для которого существует индекс в обеих матрицах, попал бы на светоприемник, если бы диафрагма была полностью прозрачной. Оценить виньетирование при полностью открытой маске диафрагмы (маска имеет прямоугольную форму!) можно, разделив количество попавших на светоприемник лучей к общему количеству отраженных лучей.
Проходя в цикле по всем прошедшим лучам, будем инкрементировать содержимое тех ячеек матрицы изображения, для которых соответствующее значение матрицы маски равно единице.
При моделировании приходится идти на компромисс между затратами ресурсов (времени вычисления и памяти) о точностью модели. Чем меньше будут размеры всех матриц в модели, тем меньшим будет потребление ресурсов, однако тем хуже будет точность. В то же время очень больших точностей достичь невозможно ввиду слишком высоких требований к объему оперативной памяти и слишком сильно возрастающим временем вычислений.
Следует учесть, что при рассогласовании размеров матрицы «зеркала» и матрицы «светоприемника» изображения пятен имеют разрывы (как на рис. сверху слева, где матрица «зеркала» имела размер 4000x4000, а матрица «светоприемника» - 3000x3000). Чтобы избавиться от этих разрывов, необходимо либо увеличить размер матрицы «зеркала» (однако, при этом очень сильно возрастает ресурсоемкость), либо уменьшить размер матрицы «светоприемника» (что приведет уменьшению точности модели), либо сгладить изображение при помощи простого гауссова фильтра (рис. сверху справа) или иного фильтра.
Результаты моделирования представлены на рис. снизу. Слева на рисунке изображена полученная модель предфокальной гартманнограммы при положении светоприемника в точке z=23.9м с использованием гауссова сглаживания для компенсации шумов дискретизации. Справа на рисунке - модель отклонений поверхности зеркала (в микрометрах), для которой получена данная модель. Размеры матриц модели были следующими: матрица зеркала имела размер 4000x4000, матрица-маска диафрагмы - 1200x1200, матрица светоприемника - 3000x3000 с размером пикселя 12.5мкм (для наилучшего приближения к реально используемому светоприемнику).
Оценку влияния отклонения поверхности зеркала от идеальной на результирующее изображение можно по приведенным ниже изображениям для моделей «плохого» (отклонения поверхности до 16нм) и «очень плохого» (отклонения поверхности до 75нм) зеркала можно по приведенным ниже рисункам.
Искажения «плохого» зеркала. Изолинии на карте поверхности проведены через 1нм.
Предфокальная гартманограмма для «плохого» зеркала
Искажения «очень плохого» зеркала. Изолинии на карте поверхности проведены через 10нм.
Предфокальная гартманограмма для «очень плохого» зеркала Исходники
function [Pre Post zn] = prepare_all() % % [Pre Post zn] = prepare_all() % формирует массивы Pre и Post - гартманнограммы, которые получились бы на БТА % при деформации зеркала zn (из build_mirror) % при положении светоприемника в точках z=24.148м и z=23.9м % % строим зеркало: матрица 4Kx4K [nx ny nz x y z zn] = build_mirror(4000,6.05,24.024); falling = single([0 0 -1]); % свет падает нормально [rx ry rz] = reflection(nx,ny,nz, falling); % получаем отраженные лучи % строим маску, 1200x1200 пикселей [mask scale] = make_mask(1200); % получаем "гартманнограммы", 3Kx3K Post = build_image(mask, 20.017, scale, 24.148, 12.5e-6, 3000, rx,ry,rz,x,y,z); Pre = build_image(mask, 20.017, scale, 23.9, 12.5e-6, 3000, rx,ry,rz,x,y,z); % выполняем фильтрацию для сглаживания неоднородностей f = fspecial('gaussian', 9, 3); Post = imfilter(Post, f); Pre = imfilter(Pre, f); endfunction
function [nx ny nz x y z znoice] = build_mirror(S, D, F, zmask) % build_mirror(S, R, F) % S - размер выходной квадратной матрицы % D - диаметр зеркала в метрах % F - фокусное расстояние зеркала % zmask - матрица отклонений поверхности зеркала от параболы (в метрах), если пропущена - берется значение по умолчанию % на выходе имеем матрицы % nx, ny,nz - матрицы SxS проекций нормалей к поверхности зеркала в точках (xi,yi) % z - такая же матрица высот (от центра зеркала) % x,y - матрицы координат % (x-xc)^2+(y-yc)^2 = 4fz % R = D/2; % радиус зеркала scale0 = D / 255; % масштаб для "заготовки" зеркала (256х256) printf('Making mirror prototype\n'); fflush(stdout); [x0, y0] = meshgrid([-127.5:127.5]*scale0); % сетка для заготовки z = (x0.^2 + y0.^2)/4/F; % подготавливаем заготовку if(nargin == 3) zmask = [ 0 1 2 3 2 1 0; 1 2 2 3 3 2 1; 2 3 2 2 2 3 2; 2 3 1 0 2 3 3; 2 3 2 3 3 1 0; 1 2 3 1 1 0 0; 0 1 2 0 1 1 0; ] * 2e-6; % отклонения в м endif tmp = size(zmask,1); A = tmp / 2 - 0.5; [x1, y1] = meshgrid([-A:A]*D/(tmp-1)); % сетка для маски отклонений printf('Interpolating deformation mask\n'); fflush(stdout); znoice = interp2(x1,y1,zmask,x0,y0,'spline'); % сглаживаем маску отклонений z = z + znoice; % добавляем искривление поверхности znoice(find((x0.^2 + y0.^2)>R^2+scale0^2)) = nan; scale = D / (S-1); % шаг по координате A = S / 2 - 0.5; % амплитуда (в пикселях) [x, y] = meshgrid([-A:A]*scale); % сетка координат на зеркале x = single(x); y = single(y); printf('Interpolating mirror matrix to %dx%d\n',S,S); fflush(stdout); z = single(interp2(x0,y0,z,x,y,'spline')); % интерполяцией получаем новые значения printf('Compute normals\n'); fflush(stdout); % вычисляем нормали: printf('nx ... '); fflush(stdout); tmp = single(zeros(size(z))); tmp1= single(zeros(size(z))); tmp(:,[2:end]) = z(:,[1:end-1]); % получаем сдвинутые матрицы (для nx) tmp1(:,[1:end-1]) = z(:,[2:end]); nx = (tmp - tmp1) / 2 / scale; % координата x ненормированной нормали printf('ny ... '); fflush(stdout); tmp = single(zeros(size(z))); tmp1= single(zeros(size(z))); tmp([2:end],:) = z([1:end-1],:); % получаем сдвинутую вниз матрицу (для ny) tmp1([1:end-1],:) = z([2:end],:); ny = (tmp - tmp1) / 2 / scale; % координата y ненормированной нормали % нормируем нормали printf('norm ... '); fflush(stdout); amp = sqrt(nx.^2+ny.^2+1); % длина нормали nx = nx ./ amp; ny = ny ./ amp; nz = ones(size(nx)) ./ amp; printf('remove outern ... '); fflush(stdout); idx = find((x.^2+y.^2)>R^2+scale^2); % находим ячейки, где зеркала уже нет nx(idx) = nan; ny(idx) = nan; % и удаляем несуществующие nz = sqrt(1-nx.^2-ny.^2); z(idx) = nan; % выбрасываем значения z в этих ячейках printf('OK!\n'); endfunction
function [rx ry rz] = reflection(nx, ny, nz, f) % reflection(nx, ny, nz, fx, fy, fz) % получение координат вектора отраженного луча % nx,ny,nz - нормали к отражающей поверхности % f == [fx fy fz] - вектор, характеризующий падающий пучок света % rx,ry,rz - векторное поле отраженных лучей с нормалями nx,ny,nz % ВСЕ ВЕКТОРЫ - НОРМИРОВАННЫЕ!!! fx = f(1); fy = f(2); fz = f(3); % знак "-" обусловлен противоположной направленностью падающего луча и нормали koeff = single(-2*(fx.*nx + fy.*ny + fz.*nz)); % удвоенная проекция падающего луча на нормаль % при отражении к падающему вектору надо прибавить удвоенное произведение % нормали на модуль проекции падающего вектора на нормаль, чтобы получить % отраженный вектор rx = koeff.*nx + fx; ry = koeff.*ny + fy; rz = koeff.*nz + fz; endfunction
function [mask scale] = make_mask(SS) % построение гартмановской маски % SS - размер маски R = [175 247 295 340 379 414 448 478] * 1e-3; % радиусы колец на гартманограмме HoleR = 7.5e-3; % радиус отверстий - 7.5мм R0 = .6; % радиус самой гартманограммы alpha0 = pi/32; % смещение лучей относительно секущих горизонтали/вертикали mask = logical(zeros(SS)); % маска гартманограммы scale = 0.6 / (SS / 2); % масштаб на маске (м на пиксель) Angles = [0:31] * 2 * alpha0 + alpha0; % углы, по которым располагаются лучи % для того, чтобы разместить на маске окружности, создадим маску % окружности: zeros(15) с единицами там, где должна быть дырка. Затем % пометим единицами в mask те точки, куда должен попадать левый верхний % край маски с дыркой и сделаем свертку (conv2 или xcorr2) HS = ceil(15e-3/scale); % размер матрицы "дырок" hole = logical(zeros(HS)); % !!! КООРДИНАТА Y ПЕРЕВЕРНУТА "ВВЕРХ НОГАМИ" !!! [X Y] = meshgrid(([1:HS]-ceil(HS/2))*scale); % координаты в матрице "дырок" dX = -X(1,1); % смещение относительно центра для верхнего левого угла матрицы дырок hole(find(X.^2+Y.^2 <= HoleR^2)) = 1; %[X Y] = meshgrid(([1:SS]-ceil(SS/2))*scale); % координаты в матрице гартманограммы for i = [1 : size(R,2)] % цикл по кольцам x = round(R(i) * cos(Angles)/scale + SS/2); y = round(-R(i) * sin(Angles)/scale + SS/2); mask((y-1)*SS+x) = 1; % помечаем положения матриц дырок на матрице гартманограммы end % помечаем маркеры x = round(R([4 8]) .* cos([-2 -7]* 2 * alpha0)/scale + SS/2); y = round(-R([4 8]) .* sin([-2 -7]* 2 * alpha0)/scale + SS/2); mask((y-1)*SS+x) = 1; tmp = xcorr2(mask, hole); mask = logical(tmp([1:SS], [1:SS])); endfunction
function ima = build_image(mask, maskz, maskS, imagez, imageS, imageN, rx,ry,rz, x0,y0,z0) % build_image(mask, maskz, maskS, imagez, imageN, rx,ry,rz, x0,y0,z0) % построение изображения через маску % ima - результирующее изображение % mask - булева двумерная матрица маски (1 - свет проходит, 0 - нет) % maskz - координата Z маски (от вершины зеркала) % maskS - масштаб в плоскости маски (м/пикс) % imagez - координата Z изображения % imageS - масштаб в плоскости изображения (м/пикс) % imageN - длина стороны изображения (пикс) % rx,ry,rz - координаты направляющих векторов отраженных от зеркала лучей % x0,y0,z0 - координаты точек на поверхности зеркала, соотв. направл. векторам % ima = int16(zeros(imageN)); k = (maskz - z0) ./ rz; % коэффициент увеличения вектора M_0 ==> M в плоскости маски x = x0 + rx.*k; % координаты лучей в плоскости маски y = y0 + ry.*k; idxX = int16(round(size(mask,2)/2+x./maskS)); % индексы на маске idxY = int16(round(size(mask,1)/2+y./maskS)); % удаляем индексы, выхоодящие за пределы размера маски % (для int nan переходит в 0) tmp = find(idxX < 1 | idxY < 1 | idxX > size(mask, 2) | idxY > size(mask, 1)); idxX(tmp) = 0; % для idxY не делаем обнуления, т.к. проверять будем лишь idxX k = (imagez - z0) ./ rz; % координаты лучей в плоскости изображения => пикс. x = int16(round(imageN/2+(x0 + rx.*k)./imageS)); y = int16(round(imageN/2+(y0 + ry.*k)./imageS)); % чистим выходящее за пределы tmp = find(x < 1 | y < 1 | x > imageN | y > imageN); x(tmp) = 0; % ищем области, через которые лучи пройдут tmp = find(idxX & x); NptsTotal = size(find(~isnan(rx)),1); printf('\n%g%% of rays falls onto the image (within transparent mask)\n', 100*size(tmp,1)/NptsTotal); printf('Start checking %d points\n', size(tmp,1)); fflush(stdout); for j = 1:size(tmp,1) i = tmp(j); if(mask(idxY(i), idxX(i)) == 1) ++ima(y(i), x(i)); % отмечаем на изображении прошедший луч end end printf('\n%g%% of rays falls onto the image\n', 100*sum(sum(ima))/NptsTotal); endfunction