0001 function img = HalleyFractal(P)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 NITER = 100;
0013 threshold = .001;
0014 pixelnum = 800;
0015 colorArr = [7,3];
0016
0017 colorArr(1,1) = 1;colorArr(1,2) = 0;colorArr(1,3) = 0;
0018
0019 colorArr(2,1) = 0;colorArr(2,2) = 1;colorArr(2,3) = 0;
0020
0021 colorArr(3,1) = 0;colorArr(3,2) = 0;colorArr(3,3) = 1;
0022
0023 colorArr(4,1) = 1;colorArr(4,2) = 1;colorArr(4,3) = 0;
0024
0025 colorArr(5,1) = 1;colorArr(5,2) = 0;colorArr(5,3) = 1;
0026
0027 colorArr(6,1) = 0;colorArr(5,2) = 1;colorArr(5,3) = 1;
0028
0029 colorArr(7,1) = 1;colorArr(6,2) = 1;colorArr(6,3) = 1;
0030
0031
0032 [xs,ys] = meshgrid(linspace(-1,1,pixelnum), linspace(-1,1,pixelnum));
0033
0034
0035 solutions = xs(:) + 1i*ys(:);
0036
0037 select = 1:numel(xs);
0038
0039
0040 niters = NITER*ones(numel(xs), 1);
0041
0042
0043 Proots = roots(P);
0044 if isempty(Proots)
0045 disp('Polynomial has no roots');
0046 return;
0047 end
0048
0049
0050 Pderivative = zeros(length(P) - 1,1);
0051 for it = 1:length(P)-1
0052 Pderivative(it)=(length(P)-it)*P(it);
0053 end
0054
0055 P2derivative = zeros(length(Pderivative) - 1,1);
0056 for it = 1:length(Pderivative)-1
0057 P2derivative(it)=(length(Pderivative)-it)*Pderivative(it);
0058 end
0059
0060 for iteration = 1:NITER
0061
0062
0063 oldi = solutions(select);
0064
0065 solutions(select) = oldi - 2 .* polyval(P,oldi) .* ...
0066 polyval(Pderivative,oldi) ./ (2.*polyval(Pderivative,oldi).^2 ...
0067 - polyval(P,oldi).*polyval(P2derivative,oldi));
0068
0069
0070 differ = (oldi - solutions(select));
0071
0072 converged = abs(differ) < threshold;
0073
0074 problematic = isnan(differ);
0075
0076
0077 niters(select(converged)) = iteration;
0078
0079
0080 niters(select(problematic)) = NITER+1;
0081
0082 select(converged | problematic) = [];
0083 end
0084
0085 Max = max(niters);
0086 niters = reshape(niters,size(xs));
0087 solutions = reshape(solutions,size(xs));
0088
0089 A = zeros(pixelnum, pixelnum);
0090 B = uint8(round(A * 255));
0091 RowCol = size(solutions);
0092 rows = RowCol(1);
0093 cols = RowCol(2);
0094 for i1 = 1:rows
0095 for i2 =1:cols
0096
0097 tmp = abs(repmat(solutions(i1,i2), size(Proots))-Proots);
0098 rootIndex = find(tmp<threshold);
0099 if ~isempty(rootIndex)
0100
0101 B(i1,i2,1)=colorArr(rootIndex,1) * ...
0102 (1-(niters(i1,i2)/(Max))) * 255;
0103 B(i1,i2,2)=colorArr(rootIndex,2) * ...
0104 (1-(niters(i1,i2)/(Max))) * 255;
0105 B(i1,i2,3)=colorArr(rootIndex,3) * ...
0106 (1-(niters(i1,i2)/(Max))) * 255;
0107 end
0108 end
0109 end
0110 img = image(B);
0111 axis off;
0112 end