function hw2(picture_file)
% Created by Prof. Braun at Plymouth Institute of Neuroscience
% Modified by Nao Tsuchiya Apr 15 03
% 
% Project a picture of your face onto visual cortex.  Assume that the
% picture subtends 20 by 20 deg of visual angle.  Use the cortical
% magnification formulas
%
% (epsilon[deg], alpha[rad]) -> (X[mm], Y[mm]) 
%
% for -pi/2 <= alpha < pi/2
%
% X      = -lambda * ln( 1 + epsilon / epsilon0 )         
% deltaX = -lambda * deltaepsilon/ ( epsilon0 + epsilon )
%
% Y      = -lambda * epsilon * alpha / ( epsilon0 + epsilon )
% deltaY = -lambda * epsilon * deltaalpha / ( epsilon0 + epsilon)
%

% read image file and scale to approx 20 b 20 deg

if nargin < 1
    A = imread( 'amandamug2.jpg' );
else
    A = imread(picture_file);
end
[height, width, depth] = size( A );

scale = sqrt( 800 ) / sqrt( width*width + height*height );

horiz = scale * ( [1:width] - 0.5*width );
verti = scale * ( [1:height] - 0.5*height );


figure(1);
image( horiz, verti, A );
axis equal;
xlabel( 'deg visual angle' );
ylabel( 'deg visual angle' );
title( 'World' );

% print hw3world -dpsc;
% close(1);

% compute cortical location and size of image pixels

global EPSILON0 LAMBDA;

EPSILON0 = 1;
LAMBDA = 12;

[H, V] = meshgrid( horiz, verti );

epsilon      = sqrt( H.*H + V.*V );                % in degrees
alpha        = atan2( V, H );                      % in radians

deltaepsilon = scale * ones(size(epsilon));
deltaalpha   = deltaepsilon ./ epsilon;

[X,Y,DX,DY] = world2cortex( epsilon, alpha, deltaepsilon, deltaalpha );

% paint pixels onto cortex

XI = X-0.5*DX;
XF = X+0.5*DX;

YI = Y-0.5*DY;
YF = Y+0.5*DY;

clear X Y DX DY H V;

figure(2);
hold on;

for i=1:height
	for j=1:width
		image( [XI(i,j) XF(i,j)], [YI(i,j) YF(i,j)], A(i,j,:) );
	end
end

hold off;
xlabel( 'mm cortex' );
ylabel( 'mm cortex' );
title( 'Cortex' );

clear all;

%print hw3cortex -dpsc;


function [X, Y, DX, DY] = world2cortex( epsilon, alpha, deltaepsilon, deltaalpha )
% Visual coordinates (epsilon[deg], alpha[rad]) 
% -> 
% Cortical coordinates (X[mm], Y[mm]) 
%
% for -pi/2 <= alpha < pi/2
%
% X      = -lambda * ln( 1 + epsilon / epsilon0 )         
% deltaX = -lambda * deltaepsilon/ ( epsilon0 + epsilon )
%
% Y      = -lambda * epsilon * alpha / ( epsilon0 + epsilon )
% deltaY = -lambda * epsilon * deltaalpha / ( epsilon0 + epsilon)
%
% for alpha < -pi/2
%
% X      =  lambda * ln( 1 + epsilon / epsilon0 )         
% deltaX =  lambda * deltaepsilon/ ( epsilon0 + epsilon )
%
% Y      = -lambda * epsilon * (-pi - alpha) / ( epsilon0 + epsilon )
% deltaY =  lambda * epsilon * deltaalpha / ( epsilon0 + epsilon)
%
% for alpha >= pi/2
%
% X      =  lambda * ln( 1 + epsilon / epsilon0 )         
% deltaX =  lambda * deltaepsilon/ ( epsilon0 + epsilon )
%
% Y      = -lambda * epsilon * ( pi - alpha) / ( epsilon0 + epsilon )
% deltaY =  lambda * epsilon * deltaalpha / ( epsilon0 + epsilon)
%

global LAMBDA EPSILON0;

if    size( epsilon ) ~= size( alpha ) ...
	| size( deltaepsilon ) ~= size( deltaalpha ) ...
	| size( alpha ) ~= size( deltaalpha )

    [' Oops, incompatible arguments of world2cortex!' ]
	size(epsilon)
	size(alpha)
	size(deltaepsilon)
	size(deltaalpha)
	return;
	
end

X  = zeros( size( epsilon ) );
Y  = zeros( size( epsilon ) );
DX = zeros( size( epsilon ) );
DY = zeros( size( epsilon ) );

M  = LAMBDA * log( 1 + epsilon / EPSILON0 );
DM = LAMBDA * ones(size(epsilon)) ./ ( epsilon + EPSILON0 );

ii = find( alpha > pi/2 );
	
    X(ii)    =  M(ii);
    Y(ii)    =  DM(ii) .* epsilon(ii) .* ( pi - alpha(ii) );
    DX(ii)   =  DM(ii) .* deltaepsilon(ii);
    DY(ii)   =  DM(ii) .* epsilon(ii) .* deltaalpha(ii);
	
jj = find( pi/2 >= alpha & alpha > -pi/2 );
	
    X(jj)    = -M(jj);
    Y(jj)    =  DM(jj) .* epsilon(jj) .* alpha(jj);
    DX(jj)   = -DM(jj) .* deltaepsilon(jj);
    DY(jj)   = -DM(jj) .* epsilon(jj) .* deltaalpha(jj);
	
kk = find( -pi/2 >= alpha );
	
    X(kk)    =  M(kk);
    Y(kk)    =  DM(kk) .* epsilon(kk) .* ( -pi - alpha(kk) );
    DX(kk)   =  DM(kk) .* deltaepsilon(kk);	
    DY(kk)   =  DM(kk) .* epsilon(kk) .* deltaalpha(kk);