function disco3d(sector1)

if (~exist('sector1','var'))
    sector1 = -90;
end

[fname folder] = uigetfile('*.jpg;*.tiff;*.tif;*.bmp','Choose an image file');
if (isnumeric(fname))
    disp('Aborted!')
    return
end

in=imread(fullfile(folder,fname));
xmax = size(in,1);
ymax = size(in,2);

figure
imagesc(in)
title('Choose region of interest or <ENTER> for full image')
axis equal

ROI=ginput(2);
if (isempty(ROI))
    ROI=[1 1;ymax xmax];
end
    
ROI = int16(ROI);
x_in = ROI(:,1);
y_in = ROI(:,2);

data = in(min(y_in):max(y_in),min(x_in):max(x_in),:);

figure
imagesc(data)
title('Cropped image - choose three colors or <ENTER> for RGB')
axis equal
colors=ginput(3);

if (isempty(colors))
    color_vec = [255 0 0;0 255 0;0 0 255];
else
    colors = int16(colors);
    color_vec = [data(colors(1,2),colors(1,1),1),data(colors(1,2),colors(1,1),2),data(colors(1,2),colors(1,1),3);data(colors(2,2),colors(2,1),1),data(colors(2,2),colors(2,1),2),data(colors(2,2),colors(2,1),3);data(colors(3,2),colors(3,1),1),data(colors(3,2),colors(3,1),2),data(colors(3,2),colors(3,1),3)];
end

figure
subplot(2,3,1)
image(data(:,:,1))
axis equal tight
%caxis([0 255]);
colormap gray
title('{\color{red}Red}')

subplot(2,3,2)
image(data(:,:,2))
axis equal tight
title('{\color{green}Green}')

subplot(2,3,3)
image(data(:,:,3))
axis equal tight
title('{\color{blue}Blue}') 

color_vec=single(color_vec)/255
image1 = zeros(size(data,1),size(data,2));
image2 = zeros(size(image1));
image3 = zeros(size(image1));

for i=1:size(data,1)
    for j=1:size(data,2)
        image1(i,j) = color_vec(1,:)*single(squeeze(data(i,j,:)));
        image2(i,j) = color_vec(2,:)*single(squeeze(data(i,j,:)));
        image3(i,j) = color_vec(3,:)*single(squeeze(data(i,j,:)));
    end
end

subplot(2,3,4)
image(image1)
axis equal tight
%caxis([0 255]);
colormap gray
str=strcat('\color[rgb]{', num2str(color_vec(1,:)), '}Channel 1');
title(str)

subplot(2,3,5)
image(image2)
axis equal tight
str=strcat('\color[rgb]{', num2str(color_vec(2,:)), '}Channel 2');
title(str)

subplot(2,3,6)
image(image3)
axis equal tight
str=strcat('\color[rgb]{', num2str(color_vec(3,:)), '}Channel 3');
title(str)

% I have defined sector 1 to be the -y illumination sector, which was
% yellow in this case, so if the channels are RGY, and the sectors are 231,
% we then define the mapping vector this way.
mapping = [2 3 1]; 
sectors=[sector1;sector1+120;sector1+240];
%sectors(mapping)

A=[1 1 1.5]; % strength of each channel

image1=image1-mean(mean(image1));
image2=image2-mean(mean(image2));
image3=image3-mean(mean(image3));

c1=A(1)*cos(sectors(mapping(1))*pi/180)*image1;
c2=A(2)*cos(sectors(mapping(2))*pi/180)*image2;
c3=A(3)*cos(sectors(mapping(3))*pi/180)*image3;
s1=A(1)*sin(sectors(mapping(1))*pi/180)*image1;
s2=A(2)*sin(sectors(mapping(2))*pi/180)*image2;
s3=A(3)*sin(sectors(mapping(3))*pi/180)*image3;

sl_x = (c1+c2+c3)/3;
sl_y = -(s1+s2+s3)/3;

figure
subplot(2,3,1)
image(c1)
axis equal tight
title('Cosine Sector 1')
subplot(2,3,2)
image(c2)
axis equal tight
title('Cosine Sector 2')
subplot(2,3,3)
image(c3)
axis equal tight
title('Cosine Sector 3')
subplot(2,3,4)
image(s1)
axis equal tight
title('Sine Sector 1')
subplot(2,3,5)
image(s2)
axis equal tight
title('Sine Sector 2')
subplot(2,3,6)
image(s3)
axis equal tight
title('Sine Sector 3')

figure
imagesc(sl_x)
title('X Slope')
caxis([-60 60]);
colorbar

figure
imagesc(sl_y)
title('Y Slope')
caxis([-60 60]);
colorbar

figure
imagesc(sqrt(sl_x.^2+sl_y.^2))
title('Slope Magnitude')
colorbar

map = integrateImage(sl_x,sl_y,0);
map1 = real(map);
map2 = imag(map);
map3 = sqrt(map.*conj(map));

figure
imagesc(map1)
title('Surface Map')
axis equal

figure
imagesc(map2)
title('Surface Map')
axis equal

figure
imagesc(map3)
title('Surface Map')
axis equal

figure
surf(flipud(map3))
shading interp
colormap hot
% [X,Y] = meshgrid(1:size(map,2),1:size(map,1));
% flatmap=reshape(map,[1,size(map,2)*size(map,1)]);
% figure
% plot3(reshape(X,[1,size(map,2)*size(map,1)]),reshape(Y,[1,size(map,2)*size(map,1)]),flatmap,'b.');
% title('Surface Map')

return