Доля_ПГ_Радон МатЛаб (857495), страница 2
Текст из файла (страница 2)
Команда imagesc(x,y,A) отображает матрицу A в виденабора цветов (графического образа) на координатной сетке с диапазономопределяемым векторами x и y.Иногда желательно отобразить матрицу как поверхность над прямоугольнойобластью. Это можно сделать с помощью функции surface(Z), гдематрица Z содержит вещественные числа, которые будут интерпретироватьсякак высота поверхности над плоскостью XY. Но матрицы, прочитанные из7графического файла, содержали только целые числа. Здесь нам можетпотребоваться преобразование типа. Например, командаA = imread('Pogorelov32x32.bmp');C = double(A)/16;создаст вещественную матрицу double(A), каждый элемент которой будетвещественным числом, которые после деления на 16 создадут матрицу C сдействительными элементами в диапазоне от 0 до 1.C=double(A)/16;surface(C);% рисует поверхность матрицы с doubleэлементами.surface(C*10); colormap(hot);Аналогично работает команда surf(Z), которая создает 3 – х мерныйграфик поверхности по z координатам точек, заданных в матрице Zиспользуя x=1:n иy=1:m, где [m,n]=size(Z).
Высоты точекповерхности Z являются вещественными числами, а цвет поверхностипропорционален высоте точек.A = imread('Pogorelov32x32.bmp')C=double(A)/16;surfc(C);colormap(gray(17));Полезной может быть также функция построения каркасного изображенияповерхности. Один из наиболее часто используемых ее синтаксисовследующийmesh(X,Y,Z);где матрицы X,Y,Z должны быть одинакового размера. График поверхностистроится над прямоугольной областью плоскости XY, которая разбиваетсясеткой c размером таким же, как у матриц. В узлах сетки значения функцииравны соответствующим значениям элементов матрицы Z. Матрицы X и Yсодержат x и y координаты узлов сетки. В следующем примере мы строимединичную полусферу по уравнению⎧⎪ 1 − x 2 − y 2 , x 2 + y 2 ≤ 1,z=⎨22⎪⎩ 0,x + y >1где − 1 ≤ x ≤ 1 , − 1 ≤ y ≤ 1 .[X Y]=meshgrid(-1:0.1:1);XY=X.^2+Y.^2;8C=XY<=1;XYZ=zeros(size(XY));XYZ(C)=1-XY(C);Z=sqrt(XYZ);mesh(X,Y,Z);% нулевая матрица% логическое индексированиеЗдесь мы использовали логическое индексирование так же, как описановыше в этом параграфе.3.
Прямое преобразование Радона. Функция radon.Функция прямого преобразования Радона radon для своего выполнениятребует установки двух основных параметров – исходного изображения I иуглов theta. Она имеет следующий упрощенный синтаксис вызоваR = radon(I);R = radon(I,theta);Исходное изображение задается в форме матрицы I. Аргумент thetaзадает вектор углов, для которых должно выполняться преобразование. Еслиаргумент theta не указан, то он принимается равным0:179.Возвращаемое значение R представляет собой матрицу, в которой каждыйстолбец является преобразованием Радона для каждого из углов,содержащихся в векторе theta.
Матрица проекций R имеет формат данныхdouble.Пример 1. Создадим матрицу, состоящую из единиц и нулей.I = zeros(100,100); I(25:75, 25:75) = 1; imshow(I);Здесь единицы в матрице I соответствуют белым точкам рисунка, а нули –черным. Теперь получим матрицу проекций, т.е. выполним преобразованиеРадона.R=radon(I);Поскольку матрица R содержит вещественные элементы, а не типа uint8, топеред построением ее образа она должна быть масштабирована под текущуюпалитру.
Используем для этого функцию imagesc(R)imagesc(R); colormap(hot); % рисунок слева9Рисунок справа получен после выполнения командcolormap(gray(129)); colorbar;Команда colorbar рисует в графическом окне справа от образа текущуюпалитру цветов. Команды[R,xp]=radon(I); mesh(0:179,xp,R); colormap(gray);строят поверхность радоновского образа – матрицы R□Пример 2.
Построим ПР характеристической функции (ХФ) единичногокруга. Вначале создадим матрицу, значения которой представляют ХФединичного круга. Разберем методику на матрице небольшого размера[X,Y]=meshgrid(-1:0.5:1)XY=X.^2+Y.^2C=XY<=1% создаем логическую матрицуZ=ones(size(C))% единичная матрица вещественнаяZF=zeros(size(C))Все строки мы не завершаем точкой с запятой. Это приводит к выводуматриц в командное окно. Для краткости мы здесь не приводим этих матриц.Теперь посмотрим, как выделить нужные элементы.
Команда ZF(C)выделяет вектор – столбец элементов нулевой матрицы ZF, для которых влогической матрице C значения равны единице.ZF(C)'1111% отображаем транспонированный вектор - столбец111111111Элементам нулевой матрицы ZF, стоящим на таких местах, что в логическойматрице C стоят единицы, присваиваем значения матрицы Z из аналогичныхячеек (с теми же номерами строк и столбцов). Остальные значения матрицыZF не меняются. Это делает командаZF(C)=Z(C)ZF =0010001110111110111000100Отобразим матрицу ZF, например, командой10imagesc(ZF); colormap(gray);Теперь сделаем то же, но для матриц большего размера.[X,Y]=meshgrid(-1:0.1:1);XY=X.^2+Y.^2;C=XY<=1;Z=ones(size(C));ZF=zeros(size(C));ZF(C)=Z(C);imagesc(ZF); colormap(gray); axis image; % левый рисунокsurf(ZF*5); axis image; colormap(gray); % правый рисунокТеперь выполним ПР для характеристической функции единичного круга,используя матрицу ZF.R=radon(ZF,0:1:179);imagesc(R); colormap(gray); colorbar;По горизонтали у нас отложены углы в градусах, а по вертикали значенияномеров строк матрицы R (их количество n определено по умолчанию).
Есливыполнить командуmax(max(R))ans =20.9740то найдем максимальное значение в матрице R. Функция colorbarпоказала нам, что самым темным точкам образа соответствуют нулевыезначения матрицы R, а самым светлым – значения близкие к 20.Команда[R,xp] = radon(ZF,1:179);xp’-16-15-14-13...14151611показывает, что вычисления выполнены для 33 точек (значения s вдольпроекции при фиксированном угле) изменяется от -16 до +16. Кроме тогоясно, что начало координат XY изображения находится в центре матрицы.Мы начинали с ХФ единичного круга, которую представили матрицейнад прямоугольником [-1 1] x [-1 1]. Для приведения графика проекции креальным физическим координатам мы должны показать диапазонизменения угла [0, π ] , а радиального параметра [− 2 , 2 ] .RU=uint8(R);max(max(RU))ans =21image([0 pi],[-sqrt(2) sqrt(2)],RU);axis image; colormap(gray(21));□Рассмотрим подробнее синтаксис функции radon.
Он может быть одним изследующих:R=radon(I, theta)R=radon(I, theta, n)[R, xp]=radon(...)Функция R=radon(I, theta) выполняет преобразование Радонаполутонового изображения, задаваемого матрицей I, и помещает результат вматрицу проекций R. Преобразование Радона представляет собой вычислениепроекций изображения на прямые, задаваемые углами в градусах,отсчитываемых против часовой стрелки от горизонтали. Эти углыпередаются в виде вектора в параметре theta. Если theta – скаляр, то Rявляется вектор – столбцом, содержащим преобразование Радона для углаtheta.
Если theta – вектор, то R представляет собой матрицу, в которойкаждый столбец является преобразованием Радона для одного из углов,содержащихся в векторе theta. Если при вызове функции параметр thetaопущен, то в theta записываются значения углов от 0 до 179 с шагом в 1°.ВформатеR=radon(I,theta,n)функциявыполняетпреобразование Радона изображения I, значения каждой проекциивычисляется в n точках, и матрица R имеет n строк.
Если при вызовефункции параметр n опущен, то он устанавливается равным2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3.Если дополнительно определить выходной параметр xp[R, xp]=radon(...),то в него помещаются значения координат, в которых вычислялись значенияпроекции. Значения в k-й строке R соответствуют координате xp(k).12Матрица проекций R может рассматриваться как полутоновое изображение иимеет формат представления данных double.Напомним, что ПР двумерной функции f(x,у) на ось s представляетсобой линейный интегралR (s , ϕ ) =∞∫ f (s cosϕ − t sin ϕ , s sin ϕ + t cosϕ )dt−∞где оси s и t задаются поворотом на угол j против часовой стрелки⎡ s ⎤ ⎡ cosϕ sin ϕ ⎤ ⎡ x ⎤⎢t ⎥ = ⎢− sin ϕ cosϕ ⎥ ⎢ y ⎥⎣ ⎦ ⎣⎦⎣ ⎦Исходное полутоновое изображение I рассматривается как двумернаяфункция.
Начало координат в системе s t соответствует центральномупикселю изображения I в пиксельной системе координат. Центральныйпиксель I можно определить согласно выражению floor((size(I) +l)/2).Рассмотрим, как в преобразованиях Радона вычисляются проекции.Рассмотрим проекции под углом 0°, 15°, 30° и 45° квадрата из примера 1.I = zeros(100,100); I(25:75, 25:75) = 1; imshow(I);[R,xp] = radon(I,[0 15 30 45]);% Построение графиков для углов 0o, 15o, 30o и 45o.figure; plot(xp,R(:,1),'LineWidth',2);title('R_{45^o}(s)');grid on;figure; plot(xp,R(:,2),'LineWidth',2);title('R_{45^o}(s)');grid on;figure; plot(xp,R(:,3),'LineWidth',2);title('R_{45^o}(s)');grid on;figure; plot(xp,R(:,4),'LineWidth',2);title('R_{45^o}(s)');grid on;13Все графики можно показать в одном окнеplot(xp,R(:,1),xp,R(:,2),xp,R(:,3),xp,R(:,4),'LineWidth',2); grid on;Но это не всегда удобно.
Преобразования Радона при большом числе угловчасто строится в виде изображения матрицы. Мы это уже много раз делалиtheta = 0:179;[R,xp] = radon(I,theta);imagesc(theta,xp,R);title('R_{\psi} (S)');xlabel('\psi (degrees)'); ylabel('S');set(gca,'XTick',0:20:180);colormap(hot); colorbar;Пример. Cформируем небольшое тестовое изображение, состоящее из трехпрямоугольников различной яркости, и построим его ПР.% Пример демонстрирует преобразование Радона.% и формирование изображения.I=zeros(100, 100); I(20:80,20:80)=0.4;I(30:70,45:55)=0.6; I(45:55,45:55)=1;imshow(I);14[R, xp]=radon(I, 0:179); % Прямое преобразование Радона.% Вывод на экран проекций на оси X, Y и под углом в 45°.Проекции данного изображения на ось X, и на оси с углами под углом в 15°,30°, 45°, 75° и на ось Y приведены соответственно на следующих рисунках.По горизонтали на всех графиках отложен радиальный параметр s.figure, plot(xp,R(:,1),'LineWidth',2);title('0 (градусов)'); grid on;figure, plot(xp,R(:,16),'LineWidth',2);title('15 (градусов)'); grid on;figure, plot(xp,R(:,31),'LineWidth',2);title('30 (градусов)'); grid on;figure, plot(xp,R(:,46),'LineWidth',2);title('45 (градусов)'); grid on;figure, plot(xp,R(:,76),'LineWidth',2);title('75 (градусов)'); grid on;figure, plot(xp,R(:,91),'LineWidth',2);title('90 (градусов)'); grid on;figure, imshow( R,[],'XData',0:179, 'YData', xp);154.
Обратное преобразование Радона. Функция iradon.Функция прямого преобразования Радона radon для своего выполнениятребует установки двух основных параметров – исходного изображения I иуглов theta.R = radon(I,theta);Функция iradon реализует восстановление изображения I на основепроекционных данных, которые содержатся в массиве R.IR = iradon(R,theta);Параметр theta описывает углы (в градусах), под которыми полученакаждая проекция.Пример 1. Определи прообраз функции⎧⎪2 1 − s 2 , s ≤ 1R (s , ϕ ) = ⎨, s >1⎪⎩ 0Как он выглядит, мы уже знаем – это характеристическая функцияединичного круга22⎪⎧1, x0 + y0 ≤ 1f ( x, y ) = ⎨⎪⎩0 , x02 + y02 > 1Проверим это с помощью функции iradon. Создадим матрицу R,соответствующую функции R(s,ϕ ) .[Ph,S]=meshgrid( 0:180, -2:0.1:2);C=abs(S)<=1;RSPh=zeros(size(C));S2=S.^2;RSPh(C)=2.*sqrt(1-S2(C).^2);mesh(S,Ph,RSPh);Теперь выполним обратное ПР.I=iradon(RSPh,0:180);imagesc(I); colormap(gray); % образ слеваmesh(I);% график справа16Результат похож на ХФ круга, но точность недостаточная и, кроме того,требуется приведение к физическим координатам.Рассмотрим подробнее синтаксис функции iradon.
Он может быть однимиз следующих:I=iradon(P, theta)I=iradon(P, theta, interp, filter, d, n)[I, h]=iradon(...)Функция I=iradon(P,theta) выполняет реконструкцию изображения Iпо его проекционным данным, которые содержатся в массиве P. Строки Pявляются данными параллельно-лучевых проекций. В функции iradonцентр вращения является центральной точкой проекций и определяется повыражению ceil(size(P,1)/2).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
