How can I copy the assigned color from one point to another?

Hello Again!

Another Day Another Snag..

Background: I am attempting to project an image from a 2-D plane onto a 3-D unit sphere. My code generates a colored image and then iteratively goes through each pixel in the 2-D image and calulates its position on a 3-D sphere using a Stereographic Projection. The color from the point on the plane must match the corresponding point on the sphere.

The goal is to create interesting functions on the plane and see how they look on the sphere.

Problem: Right now the code only seems to project onto the top half of the sphere. But I need the entire sphere to be mapped.

Currently the colors on the sphere appear random and have no relation to the plane colors. I want the color of the point on the sphere to match its corresponding color on the plane.

Here is my code currently:

-----------------------------------------------------------------------

clear all
clc
clf
% Parameters
res = 0.05;               % Resolution
x = -2:res:2;            % X range
y = x';                  % Y range (transpose)
depth = 32;              % Iteration depth
grid = zeros(length(x), length(y));
map = zeros(length(x), length(y));
c = 0.30 + 0.5*1i;      % Complex parameter
% Generate grid of complex numbers
for i = 1:length(x)
    for j = 1:length(y)
        grid(i,j) = x(i) + y(j)*1i;
    end
end    
% Fractal generation using iteration
for i = 1:length(x)
    for j = 1:length(y)
        for n = 1:depth
            if abs(grid(i,j)) > 2
                map(i,j) = n;
                break;
            else
                %grid(i,j) = (grid(i,j) + c)/(grid(i,j) - 1i);
                grid(i,j) = grid(i,j); % More Complicated Functions Will Go Here
            end
        end
    end
end
map(map==0) = depth; % Set unassigned points to max depth
% Projection onto a 3D unit sphere
[x_sphere, y_sphere, z_sphere] = sphere(50); % Create a unit sphere 

% Prepare figure for fractal and sphere projection
%figure;
% Plotting the fractal image on a flat surface
subplot(1, 2, 1);
image(map);
axis image;
colormap(flipud(jet(depth)));
title('Fractal Image');
cmap = colormap(flipud(jet(depth)));
% Plotting the projected points on the sphere
subplot(1, 2, 2);
hold on;
% Map each point from the fractal to the sphere
for i = 1:length(x)
    for j = 1:length(y)
        if map(i,j) <= depth   % Only plot points that were iterated over
            r = sqrt(x(i)^2 + y(j)^2);
            theta = atan2(y(j), x(i));
            phi = acos(r / sqrt(r^2 + 1)); 
            % Calculate angle from vertical
            % Spherical coordinates conversion to Cartesian coordinates
            X = sin(phi) * cos(theta);
            Y = sin(phi) * sin(theta);
            Z = cos(phi);
            plot3(X, Y, Z,'.','MarkerFaceColor',[cmap(depth,1) cmap(depth,2) cmap(depth,3)]); %  Attempt to copy color from the plane           
        end
    end
end
% Add sphere surface for reference
surf(x_sphere, y_sphere, z_sphere, 'FaceAlpha', 0, 'linestyle','none'); 
% Transparent sphere    surface
axis equal;
title('Projection onto Unit Sphere');
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
view(3)
%grid on;
hold off;

---------------------------------------------

Any additional comments, code changes, or hints are greatly apprectiated.

P.S

You guys are amazing at your prompt and thorough responses! I always leave here feeling confident!

[Cris LaPierre]

FYI - You've named a variable grid so you can't use grid on

I recommend changing your variable name.

[Walter Roberson]

I actually commented out the "grid on" in the code, as otherwise the code does not produce any graphics !

[Cris LaPierre]

In case it wasn't clear, my comment was for the OP.

[Walter Roberson]

The original code that was posted had the "grid on", but I commented it out as it was interfering with producing graphics.

Answer 1

[Cris LaPierre]

The colors on the sphere are not random. They are following those specified in colororder.

The reason they are not honoring the color you specify is because the marker you are using does not have a face color properity. In only supports edge color. Use 'MarkerEdgeColor' instead.

Also note it is not necessary to plot an invisible sphere to get a plot3 command to appear as a sphere. Set the view and other properties on the plot3 axes.

% Parameters
res = 0.05;               % Resolution
x = -2:res:2;            % X range
y = x';                  % Y range (transpose)
depth = 32;              % Iteration depth
grd = zeros(length(x), length(y));
map = zeros(length(x), length(y));
c = 0.30 + 0.5*1i;      % Complex parameter
% Generate grd of complex numbers
for i = 1:length(x)
    for j = 1:length(y)
        grd(i,j) = x(i) + y(j)*1i;
    end
end    
% Fractal generation using iteration
for i = 1:length(x)
    for j = 1:length(y)
        for n = 1:depth
            if abs(grd(i,j)) > 2
                map(i,j) = n;
                break;
            else
                %grd(i,j) = (grd(i,j) + c)/(grd(i,j) - 1i);
                grd(i,j) = grd(i,j); % More Complicated Functions Will Go Here
            end
        end
    end
end
map(map==0) = depth; % Set unassigned points to max depth
% Projection onto a 3D unit sphere
[x_sphere, y_sphere, z_sphere] = sphere(50); % Create a unit sphere 

% Prepare figure for fractal and sphere projection
figure;
% Plotting the fractal image on a flat surface
subplot(1, 2, 1);
image(map);
axis image;
colormap(flipud(jet(depth)));
title('Fractal Image');
cmap = colormap(flipud(jet(depth)));
% Plotting the projected points on the sphere
subplot(1, 2, 2);
hold on;
% Map each point from the fractal to the sphere
for i = 1:length(x)
    for j = 1:length(y)
        if map(i,j) <= depth   % Only plot points that were iterated over
            r = sqrt(x(i)^2 + y(j)^2);
            theta = atan2(y(j), x(i));
            phi = acos(r / sqrt(r^2 + 1)); 
            % Calculate angle from vertical
            % Spherical coordinates conversion to Cartesian coordinates
            X = sin(phi) * cos(theta);
            Y = sin(phi) * sin(theta);
            Z = cos(phi);
            plot3(X, Y, Z,'.','MarkerEdgeColor',[cmap(depth,1) cmap(depth,2) cmap(depth,3)]); %  Attempt to copy color from the plane           
        end
    end
end

axis equal;
title('Projection onto Unit Sphere');
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
view(3)
grid on;
hold off;
axis([-1 1 -1 1 -1 1])

●

[Christopher Scott]

Thanks for the reply!

However the issue is still two-fold

The colors are not being compied onto the sphere from the plane. As per the figure:

There should definitley be some color variations on the sphere corresponding to the image in the plane. (The equation is grd(i,j) = grd(i,j)^2 -0.791 +0.15i;)

And I am trying to populate the entire sphere instead of just the top hemisphere. Image the sphere sitting on the plane as in the figure:

It is still the unit sphere R = 1 but it rests ontop of the plane instead of being bisected by the plane.

[Walter Roberson]

            plot3(X, Y, Z,'.','MarkerEdgeColor',[cmap(depth,1) cmap(depth,2) cmap(depth,3)]); %  Attempt to copy color from the plane           

You are not changing depth in your loops, so you are always using the same MarkerEdgeColor

By the way, it would be more efficient to use

            plot3(X, Y, Z,'.','MarkerEdgeColor', cmap(depth,:)); %  Attempt to copy color from the plane           

[Cris LaPierre]

I'd also add the code you have shared does not create the fractal image you show above. The projection onto a sphere image is also the same, whether you add a sphere surface for reference or not.

[Christopher Scott]

Okay I fixed the color problem and now the code looks like this:

---------------------------------------------------------------------------------

% Parameters
res = 0.05;               % Resolution
x = -2:res:2;            % X range
y = x';                  % Y range (transpose)
depth = 32;              % Iteration depth
grd = zeros(length(x), length(y));
map = zeros(length(x), length(y));
c = 0.30 + 0.5*1i;      % Complex parameter
% Generate grd of complex numbers
for i = 1:length(x)
    for j = 1:length(y)
        grd(i,j) = x(i) + y(j)*1i;
    end
end    
% Fractal generation using iteration
for i = 1:length(x)
    for j = 1:length(y)
        for n = 1:depth
            if abs(grd(i,j)) > 2
                map(i,j) = n;
                break;
            else
                %grd(i,j) = (grd(i,j) + c)/(grd(i,j) - 1i);
                grd(i,j) = grd(i,j)^2 -0.791 +0.15i; % More Complicated Functions Will Go Here
            end
        end
    end
end
map(map==0) = depth; % Set unassigned points to max depth
% Projection onto a 3D unit sphere
[x_sphere, y_sphere, z_sphere] = sphere(50); % Create a unit sphere 
% Prepare figure for fractal and sphere projection
figure;
% Plotting the fractal image on a flat surface
subplot(1, 2, 1);
image(map);
axis image;
colormap(flipud(jet(depth)));
title('Fractal Image');
cmap = colormap(flipud(jet(depth)));
% Plotting the projected points on the sphere
subplot(1, 2, 2);
hold on;
% Map each point from the fractal to the sphere
for i = 1:length(x)
    for j = 1:length(y)
        % Only plot points that were iterated over
        r = sqrt(x(i)^2 + y(j)^2);
        theta = atan2(y(j), x(i));
        phi = acos(r / sqrt(r^2 + 1)); 
        % Calculate angle from vertical
        % Spherical coordinates conversion to Cartesian coordinates
        X = sin(phi) * cos(theta);
        Y = sin(phi) * sin(theta);
        Z = cos(phi);
        color = map(i,j);
        plot3(X, Y, Z,'.','MarkerEdgeColor',[cmap(color,:)]); %  Attempt to copy color from the plane           
        
    end
end
axis equal;
title('Projection onto Unit Sphere');
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
view(3)
grid on;
hold off;
axis([-1 1 -1 1 -1 1])

●

-------------------------------------------------------------

Now the problem is how to map over the entire sphere instead of just the top hemisphere. The stereographic projection I need has the unit sphere sitting ON the plane so that z = -1 on the sphere coincides with (x,y) = o,o on the plane.

You guys have any ideas?

[Walter Roberson]

        plot3(X, Y, Z,'.','MarkerEdgeColor',[cmap(color,:)]); %  Attempt to copy color from the plane           
        hold on
        plot3(X, Y, -Z,'.','MarkerEdgeColor',[cmap(color,:)]); %  Attempt to copy color from the plane           
        hold off

[Cris LaPierre]

The graphic you shared also only shows a 2D image mapping to a sphere.(look at the solid mesh lines in the plane and on the sphere).

[Cris LaPierre]

Given your description of how the sphere sits on the plane, the mapping of the fractal image to the sphere seemed incorrect. The image should be centered at one of the poles, but is instead split across the equator.

It looks like you are trying to perform an inverse stereographic projection. I did some digging in and couldn't resolve your equation for phi. If I use a solution I find online, I get the result I would expect.

After cleaning up your code, here is what my solution looks like.

• I've removed all but one for loop
• I swapped how you defined x and y so that the fractal image axes align with the sphere x and y axes
• I've displayed the fractal image using the actual x and y scaling
• I use scatter3 instead of plot3 since it allows me to specify a unique color for each point in one call
• I use inverse stereographic projection instead of spherical coordinates

% Parameters
res = 0.05;               % Resolution
[x,y] = meshgrid(-2:res:2);         % X,Y range
depth = 32;              % Iteration depth
map = ones(size(x))*depth;
c = 0.30 + 0.5*1i;      % Complex parameter

% Generate grd of complex numbers
grd = complex(x,y);

% Fractal generation using iteration
for n = 1:depth
    map(abs(grd) > 2 & map==depth ) = n;
    grd(abs(grd) <= 2) = grd(abs(grd) <= 2).^2 - 0.791 +0.15i; % More Complicated Functions Will Go Here
end

% Prepare figure for fractal and sphere projection
figure;
% Plotting the fractal image on a flat surface
subplot(1, 2, 1);
image(map,'XData',[-2,2],'YData',[-2,2]);
axis xy equal tight
colormap(flipud(jet(depth)));
title('Fractal Image');

% Inverse stereographic projection
X = 4*x./(x.^2 + y.^2 + 4);
Y = 4*y./(x.^2 + y.^2 + 4);
Z = (x.^2 + y.^2 - 4)./(x.^2 + y.^2 + 4);
C = map(:);

% Plotting the projected points on the sphere
subplot(1, 2, 2);
scatter3(X(:), Y(:), Z(:),[],C,"."); %  Attempt to copy color from the plane

axis equal;
title('Projection onto Unit Sphere');
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
grid on;
axis([-1 1 -1 1 -1 1])

You can remove the corners by setting those values of Z to nan.

The projection can be mapped to the top half of the sphere by negating Z.

Z(map==1) = nan;
scatter3(X(:), Y(:), -Z(:),[],C,".");

axis equal;
title('Projection onto Unit Sphere');
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
grid on;
axis([-1 1 -1 1 -1 1])

[Christopher Scott]

Wow! Thanks for cleaning up the code and digging into the stereographic projection.

The results you have do not appear to be consistent with the projection I am trying to achieve. This is really a comblex analysis problem as technically I am trying to project the extended complex plane onto a Riemann Sphere. Eventually I want to create some Fractional Linear Mappings (FLM), project them onto a Riemann Sphere, and try to find patterns. (The overall context is the Bloch Sphere used in Modeling a Quantum Computing Q-Bit).

Conceptually the projection is like the image:

Here is a copy of my code now (before your gracefully cleaned up bits!):

----------------------------------------------------------------------------------------------------------

clear all
clc
clf
% Parameters
res = 0.05;               % Resolution
x = -2:res:2;            % X range
y = x';                  % Y range (transpose)
depth = 32;              % Iteration depth
grid = zeros(length(x), length(y));
map = zeros(length(x), length(y));
c = 0.30 + 0.5*1i;      % Complex parameter
% Generate grid of complex numbers
for i = 1:length(x)
    for j = 1:length(y)
        grid(i,j) = x(i) + y(j)*1i;
    end
end    
% Fractal generation using iteration
for i = 1:length(x)
    for j = 1:length(y)
        for n = 1:depth
            if abs(grid(i,j)) > 2
                map(i,j) = n;
                break;
            else
                %grid(i,j) = (grid(i,j) + c)/(grid(i,j) - 1i);
                grid(i,j) = grid(i,j); % More Complicated Functions Will Go Here
            end
        end
    end
end
map(map==0) = depth; % Set unassigned points to max depth
% Projection onto a 3D unit sphere
[x_sphere, y_sphere, z_sphere] = sphere(50); % Create a unit sphere 

% Prepare figure for fractal and sphere projection
%figure;
% Plotting the fractal image on a flat surface
subplot(1, 2, 1);
image(map);
axis image;
colormap(flipud(jet(depth)));
title('Fractal Image');
cmap = colormap(flipud(jet(depth)));
% Plotting the projected points on the sphere
subplot(1, 2, 2);
hold on;
% Map each point from the fractal to the sphere
for i = 1:length(x)
    for j = 1:length(y)
        if map(i,j) <= depth   % Only plot points that were iterated over
            r = sqrt(x(i)^2 + y(j)^2);
            theta = atan2(y(j), x(i));
            phi = acos(r / sqrt(r^2 + 1)); 
            % Calculate angle from vertical
            % Spherical coordinates conversion to Cartesian coordinates
            X = sin(phi) * cos(theta);
            Y = sin(phi) * sin(theta);
            Z = cos(phi);
            plot3(X, Y, Z,'.','MarkerFaceColor',[cmap(depth,1) cmap(depth,2) cmap(depth,3)]); %  Attempt to copy color from the plane           
        end
    end
end
% Add sphere surface for reference
surf(x_sphere, y_sphere, z_sphere, 'FaceAlpha', 0, 'linestyle','none'); 
% Transparent sphere    surface
axis equal;
title('Projection onto Unit Sphere');
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
grid on;
hold off;

----------------------------------------------------------------------

(On the next code passover I will implement your code effiency suggestions)

The colors are being copied as expected but only the top hemisphere of the sphere is being projected onto instead of the entire sphere (as in the image above). It should look different than a literal stretching of the image onto a sphere.

My difficulty is in reconciling equations from "Wikipedia: Stereographic Projection" to the equations in my code (not by me). I have tried the equations from wikipedia but they do not produce a meaningful result.

So thats where I am right now. The last step to finish the model before the theory phase.

I wouldnt be anywhere near where I am without your Guys help!

Always Appreciated

[Cris LaPierre]

I'm not sure I can offer much more help, then. I'm just googling the same as you to try to find an answer.

One parting thought. My understanding of stereographic projections is that points are projected to the 2D plane by a vector that passes through the North pole and the point on the sphere. When that is the case, then a plane that sits at z=-1 and is defined with x/y limits of +/-2r (where r is the radius of the sphere) can only capture the projection of the lower hemisphere of the sphere.

r=1;
% Create a unit sphere with a projection to -2r
theta=linspace(0,2*pi,200); 
x=r*cos(theta); 
y=r*sin(theta); 
plot(x,y,[-2 -1 -1 0 0],[-1 -1 0 0 1],'k^--',[-2 -1 0],[-1 0 1],'.r-');
yline(-1,'c')
axis equal
legend(["Riemann sphere","Radius","Projection","Complex plane"],'Location','best')