function display_square(img)
%function display_square(img)
%
% DISPLAY_SQUARE -- allows the user to select a rectangular portion
%                   of an image and see the pixel values in that
%                   square.  Also creates a 3-d bar graph of those
%                   pixels.  The selected area must not contain
%                   more than 300 pixels.
%
% Input Parameters:
%  img  -- The image to be examined
%
% Ouptut Parameters:
%  none

% Written by Mark Bartsch, Winter 2002
% Modification History:
%    8/16/02  --  Added modification history (MB)


cf = gcf;
close(cf);
figure(cf);
imagesc(img); colorbar; colormap(gray); axis image
xlabel('Click and drag to select a region of pixels to display.');

pause(.05);

rect = round(getrect);
if prod(rect(3:4)) > 300
    error('Please select a rectangle with fewer than 300 pixels.');
end

z = img(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3));

h = gcf;
scr_size = get(0,'ScreenSize');
set(h,'Position',[5,scr_size(4)*.4,scr_size(3)-10,scr_size(4)*.6-70]);
    
subplot(1,2,1);
plot_array(z',rect(1),rect(2));
xlabel('X');
ylabel('Y');

subplot(1,2,2);
make_bar(rect(1),rect(2),0);
hold on
make_bar(rect(1),rect(2),0,255);

for r=rect(2):rect(2)+rect(4)
    for c=rect(1):rect(1)+rect(3);
%        stem3(r,c,img(r,c));
        make_bar(c,r,img(r,c));
    end
end

xlabel('X');
ylabel('Y');
axis ij;
axis tight;
colorbar;
hold off
rotate3d on


function make_bar(x,y,z,c)

% Plot a 3-D bar with height z at position x, y

if ~exist('c','var') | isempty('c')
    c = z;
end

xx = [ 0 0 1 1 0; 0 0 1 1 0 ];
yy = [ 0 1 1 0 0; 0 1 1 0 0 ];
zz = [ 0 0 0 0 0; 1 1 1 1 1 ];

xx = xx*.8 + .6;
yy = yy*.8 + .6;

surf(xx + x, yy + y, zz*z, c*ones(size(zz)));

xx = [ 0 1 ; 0 1 ];
yy = [ 0 0 ; 1 1 ];
zz = [ 1 1 ; 1 1 ];

xx = xx*.8 + .6;
yy = yy*.8 + .6;

surf(xx + x, yy + y, zz*z, c*ones(size(zz)));




function plot_array(f,x1,y1)

% Plot an array 'f' in in a figure at position x1, y1.
% Adapted from code by Dr. Jeffery Fessler

[nx,ny] = size(f);
c = [0 0];
form = '%g';

imagesc(x1+1:x1+size(f,1),y1+1:y1+size(f,2),f',[0 255]);

rectx = [0 0 1 1 0];
recty = [0 1 1 0 0];
for ix=1:nx
    for iy=1:ny
        t = sprintf(form, f(ix,iy));
        h = text(ix-c(1)+x1, iy-c(2)+y1, t);
        set(h, 'HorizontalAlignment', 'center')
        set(h, 'VerticalAlignment', 'middle')
        set(h, 'Color', [1 0 0]);
    end
end
axis square
axis ij
%axis([[x1 nx+x1]+1/2-c(1) [y1 ny+y1]+1/2-c(2)])
axis tight;
colorbar
hold off



function rect = getrect(varargin)
%GETRECT Select rectangle with mouse.

global GETRECT_FIG GETRECT_AX GETRECT_H1 GETRECT_H2
global GETRECT_PT1 GETRECT_TYPE

if ((nargin >= 1) & (isstr(varargin{1})))
    % Callback invocation: 'ButtonDown', 'ButtonMotion', or 'ButtonUp'
    feval(varargin{:});
    return;
end

if (nargin < 1)
    GETRECT_AX = gca;
    GETRECT_FIG = get(GETRECT_AX, 'Parent');
else
    if (~ishandle(varargin{1}))
        error('First argument is not a valid handle');
    end
    
    switch get(varargin{1}, 'Type')
    case 'figure'
        GETRECT_FIG = varargin{1};
        GETRECT_AX = get(GETRECT_FIG, 'CurrentAxes');
        if (isempty(GETRECT_AX))
            GETRECT_AX = axes('Parent', GETRECT_FIG);
        end

    case 'axes'
        GETRECT_AX = varargin{1};
        GETRECT_FIG = get(GETRECT_AX, 'Parent');

    otherwise
        error('First argument should be a figure or axes handle');

    end
end

% Remember initial figure state
old_db = get(GETRECT_FIG, 'DoubleBuffer');
state = uisuspend(GETRECT_FIG);

% Set up initial callbacks for initial stage
set(GETRECT_FIG, 'Pointer', 'crosshair', 'WindowButtonDownFcn', 'getrect(''ButtonDown'');', 'DoubleBuffer', 'on');

% Bring target figure forward
figure(GETRECT_FIG);

% Initialize the lines to be used for the drag
GETRECT_H1 = line('Parent', GETRECT_AX, 'XData', [0 0 0 0 0], 'YData', [0 0 0 0 0], 'Visible', 'off', 'Clipping', 'off', 'Color', 'k', 'LineStyle', '-');

GETRECT_H2 = line('Parent', GETRECT_AX, 'XData', [0 0 0 0 0], 'YData', [0 0 0 0 0], 'Visible', 'off', 'Clipping', 'off', 'Color', 'w', 'LineStyle', ':');


% We're ready; wait for the user to do the drag
% Wrap the waitfor call in try-catch so
% that if the user Ctrl-C's we get a chance to
% clean up the figure.
errCatch = 0;
try
    waitfor(GETRECT_H1, 'UserData', 'Completed');
catch
    errCatch = 1;
end

% After the waitfor, if GETRECT_H1 is still valid
% and its UserData is 'Completed', then the user
% completed the drag.  If not, the user interrupted
% the action somehow, perhaps by a Ctrl-C in the
% command window or by closing the figure.

if (errCatch == 1)
    errStatus = 'trap';
    
elseif (~ishandle(GETRECT_H1) | ~strcmp(get(GETRECT_H1, 'UserData'), 'Completed'))
    errStatus = 'unknown';
    
else
    errStatus = 'ok';
    x = get(GETRECT_H1, 'XData');
    y = get(GETRECT_H1, 'YData');
end

% Delete the animation objects
if (ishandle(GETRECT_H1))
    delete(GETRECT_H1);
end
if (ishandle(GETRECT_H2))
    delete(GETRECT_H2);
end

% Restore the figure state
if (ishandle(GETRECT_FIG))
   uirestore(state);
   set(GETRECT_FIG, 'DoubleBuffer', old_db);
end

% Clean up the global workspace
clear global GETRECT_FIG GETRECT_AX GETRECT_H1 GETRECT_H2
clear global GETRECT_PT1 GETRECT_TYPE

% Depending on the error status, return the answer or generate
% an error message.
switch errStatus
case 'ok'
    % Return the answer
    xmin = min(x);
    ymin = min(y);
    rect = [xmin ymin max(x)-xmin max(y)-ymin];
    
case 'trap'
    % An error was trapped during the waitfor
    error('Interruption during mouse selection.');
    
case 'unknown'
    % User did something to cause the rectangle drag to
    % terminate abnormally.  For example, we would get here
    % if the user closed the figure in the drag.
    error('Interruption during mouse selection.');
end

%--------------------------------------------------
% Subfunction ButtonDown
%--------------------------------------------------
function ButtonDown

global GETRECT_FIG GETRECT_AX GETRECT_H1 GETRECT_H2
global GETRECT_PT1 GETRECT_TYPE

set(GETRECT_FIG, 'Interruptible', 'off', 'BusyAction', 'cancel');

[x1, y1] = getcurpt(GETRECT_AX);
GETRECT_PT1 = [x1 y1];
GETRECT_TYPE = get(GETRECT_FIG, 'SelectionType');
x2 = x1;
y2 = y1;
xdata = [x1 x2 x2 x1 x1];
ydata = [y1 y1 y2 y2 y1];

set(GETRECT_H1, 'XData', xdata, 'YData', ydata, 'Visible', 'on');
set(GETRECT_H2, 'XData', xdata, 'YData', ydata, 'Visible', 'on');

% Let the motion functions take over.
set(GETRECT_FIG, 'WindowButtonMotionFcn', 'getrect(''ButtonMotion'');', 'WindowButtonUpFcn', 'getrect(''ButtonUp'');');


%-------------------------------------------------
% Subfunction ButtonMotion
%-------------------------------------------------
function ButtonMotion

global GETRECT_FIG GETRECT_AX GETRECT_H1 GETRECT_H2
global GETRECT_PT1 GETRECT_TYPE

[x2,y2] = getcurpt(GETRECT_AX);
x1 = GETRECT_PT1(1,1);
y1 = GETRECT_PT1(1,2);
xdata = [x1 x2 x2 x1 x1];
ydata = [y1 y1 y2 y2 y1];

if (~strcmp(GETRECT_TYPE, 'normal'))
    [xdata, ydata] = Constrain(xdata, ydata);
end

set(GETRECT_H1, 'XData', xdata, 'YData', ydata);
set(GETRECT_H2, 'XData', xdata, 'YData', ydata);

%--------------------------------------------------
% Subfunction ButtonUp
%--------------------------------------------------
function ButtonUp

global GETRECT_FIG GETRECT_AX GETRECT_H1 GETRECT_H2
global GETRECT_PT1 GETRECT_TYPE

% Kill the motion function and discard pending events
set(GETRECT_FIG, 'WindowButtonMotionFcn', '', 'Interruptible', 'off');

% Set final line data
[x2,y2] = getcurpt(GETRECT_AX);
x1 = GETRECT_PT1(1,1);
y1 = GETRECT_PT1(1,2);
xdata = [x1 x2 x2 x1 x1];
ydata = [y1 y1 y2 y2 y1];
if (~strcmp(GETRECT_TYPE, 'normal'))
    [xdata, ydata] = Constrain(xdata, ydata);
end

set(GETRECT_H1, 'XData', xdata, 'YData', ydata);
set(GETRECT_H2, 'XData', xdata, 'YData', ydata);

% Unblock execution of the main routine
set(GETRECT_H1, 'UserData', 'Completed');

%-----------------------------------------------
% Subfunction Constrain
% 
% constrain rectangle to be a square in
% axes coordinates
%-----------------------------------------------
function [xdata_out, ydata_out] = Constrain(xdata, ydata)

x1 = xdata(1);
x2 = xdata(2);
y1 = ydata(1);
y2 = ydata(3);
ydis = abs(y2 - y1);
xdis = abs(x2 - x1);

if (ydis > xdis)
   x2 = x1 + sign(x2 - x1) * ydis;
else
   y2 = y1 + sign(y2 - y1) * xdis;
end

xdata_out = [x1 x2 x2 x1 x1];
ydata_out = [y1 y1 y2 y2 y1];




function [x,y] = getcurpt(axHandle)
%GETCURPT Get current point.

pt = get(axHandle, 'CurrentPoint');
x = pt(1,1);
y = pt(1,2);
