% fourieroperator2dat.m                            Steve Mann 1992 Jan 20
%
% creates the fourieroperator W = exp(j2pi|f><t|)
%
% visualization of the continuous W; math image_inagtion
%
% use:  fourierop2dat(200);  % creates 200 by 200 in datfile
%       fourierop2dat(M,N);  % creates  M  by  N  operator
%       fourierop2dat(512,512,.4) % scales down freq. range
%                                         (to about same as 200 by 200)
%       fourierop2dat(M,N,frac); % frac=1 (default): fracfreq -1/2 to 1/2
%
%       matrix converted to 3 real color matrices, prodeces 3ch (color) datfile
%       on Munsell system (Munsell, 1929)
%       fix value = 0 (black = 0; white =10)
%       vary chroma by magnitude;  vary hue by phase
%       uses green, yellow, red, purple, blue, on equal 72 deg space
%
% References:
% 1  Principles of Color Photography; Evans, Hanson, Brewer; John Wiley and Sons
%    1953 December, page 78   thode lib, tr 510 e92p            (Munsell space)
% 2  The Chirplet Transform; Steve Mann; World Scientific Press (Fourierop)

function dummy(M,N,frac);   % returns nothing

if nargin<2
  N = M;
end%if

if nargin<3
  frac = 1;   % default is -1/2 to 1/2 fracfreqlim
end%if

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%             INITIALIZE DATFILE HEADERS:  3 BINARY DATFILES                %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t = fix(clock);

!rm -r deletemer
!rm -r deletemer
!mkdir deletemer
!rm -r deletemeg
!rm -r deletemeg
!mkdir deletemeg
!rm -r deletemeb
!rm -r deletemeb
!mkdir deletemeb
!rm deletemer.tmp
!rm deletemeg.tmp
!rm deletemeb.tmp

fprintf('deletemer/descriptor','(_data_sets 1)\n')
fprintf('deletemer/descriptor','(_channels 1)\n')
fprintf('deletemer/descriptor','(_dimensions %g %g )\n',M,N)
fprintf('deletemer/descriptor','(_data_type unsigned_1)\n')
fprintf('deletemer/descriptor','(_history from fourieroperator2dat\n')
fprintf('deletemer/descriptor',')\n')  % end bracket for history
fprintf('deletemer/descriptor','(_date %g %g %g ',t(1),t(2),t(3))
fprintf('deletemer/descriptor','%g %g %g)\n',t(4),t(5),t(6)) % rest o date

fprintf('deletemeg/descriptor','(_data_sets 1)\n')
fprintf('deletemeg/descriptor','(_channels 1)\n')
fprintf('deletemeg/descriptor','(_dimensions %g %g )\n',M,N)
fprintf('deletemeg/descriptor','(_data_type unsigned_1)\n')
fprintf('deletemeg/descriptor','(_history from fourieroperator2dat\n')
fprintf('deletemeg/descriptor',')\n')  % end bracket for history
fprintf('deletemeg/descriptor','(_date %g %g %g ',t(1),t(2),t(3))
fprintf('deletemeg/descriptor','%g %g %g)\n',t(4),t(5),t(6)) % rest o date

fprintf('deletemeb/descriptor','(_data_sets 1)\n')
fprintf('deletemeb/descriptor','(_channels 1)\n')
fprintf('deletemeb/descriptor','(_dimensions %g %g )\n',M,N)
fprintf('deletemeb/descriptor','(_data_type unsigned_1)\n')
fprintf('deletemeb/descriptor','(_history from fourieroperator2dat\n')
fprintf('deletemeb/descriptor',')\n')  % end bracket for history
fprintf('deletemeb/descriptor','(_date %g %g %g ',t(1),t(2),t(3))
fprintf('deletemeb/descriptor','%g %g %g)\n',t(4),t(5),t(6)) % rest o date

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%   WRITE FOURIEROPERATOR INTO DATFILES IN REASONABLY SIZED BLOCKS          %
%              (each block is integer number of lines;                      %
%          blocks may differ slightly in size due to roundoff)              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%t = (0:N-1)/N;                     % time span 1 second:  [0,1) second
%t = linspace(-1/2,1/2,N);          % time span 1 second:  [-.5,.5] second
f = (-(M-1)/2:(M-1)/2)*frac;        % time span 1 second:  (-.5,.5) second
%f = linspace(-M/2*frac,M/2*frac,M);% count by ones, nyquist at M/2 Hertz
%f = (-(M-1)/2:(M-1)/2)/M);         % or M/2 cycles in M pts

% divide image into blocks of size Less or Eq to pixelsperblock
% compute starting line for each block 
%pixelsperblock = 300*300;  %  300by300 was a good size for a sparc 2
pixelsperblock = 1024*1024;  %  on HP 735, we should have at least 8 meg RAM
linesperblock  = floor(pixelsperblock/N)
if ~linesperblock
  disp('even one line is too long to fit into a block')
end%if
% blocks go down cols of matrix, across rows of datfile
for Nstartofblock=1:linesperblock:N  % do in blocks: matlab mem <= 4meg ??
  disp(sprintf('Nstartofblock = %g',Nstartofblock))
  n = (Nstartofblock : min(Nstartofblock + linesperblock - 1,N));
  t = (n - (N+1)/2)/N;  % freq range for the given block
  %%TF = f(:)*t(:).'; % the outer product |f><t|
  W  = exp(j*2*pi*f(:)*t(:).');  % faster as a one-liner

  mag = abs(W);  %  change rowvec to colvec
  o = pi/180;   % the degree symbol
  p   = 180 - angle(W)/o;  % angles per degree shifted to all positive

  r = ones(size(mag))*NaN;  % init., save malloc calls
  g = ones(size(mag))*NaN;
  b = ones(size(mag))*NaN;

  g_to_r_indices = find (p < 144 &  p >= 0);  % red to green crossfade
  r(g_to_r_indices) = p(g_to_r_indices)/144;        % RED
  g(g_to_r_indices) = (144-p(g_to_r_indices))/144;  % GREEN
  b(g_to_r_indices) = zeros(size(g_to_r_indices));        % BLUE

  r2b_ind = find (p >= 144 & p < 288);        % r to b crossfade
  r(r2b_ind) = (288-p(r2b_ind))/144;
  g(r2b_ind) = zeros(size(r2b_ind));
  b(r2b_ind) = (p(r2b_ind) - 144)/144;

  b2g_ind = find ( p >= 288 & p < 360 );      % b to g crossfade
  r(b2g_ind) = zeros(size(b2g_ind));
  g(b2g_ind) = (p(b2g_ind) - 288)/72;
  b(b2g_ind) = (360-p(b2g_ind))/72;

  R = floor(r.*mag*255.99);   % scale by magnitudes; from black to red
  G = floor(g.*mag*255.99);   % scale by magnitudes; from black to green
  B = floor(b.*mag*255.99);   % scale by magnitudes; from black to blue

  %save deletemer R
  %!/mas/vision/bigu/steve/src/reduce/dbl2byt deletemer.mat 22 >> deletemer/data
  %save deletemeg G
  %!/mas/vision/bigu/steve/src/reduce/dbl2byt deletemeg.mat 22 >> deletemeg/data
  %save deletemeb B
  %!/mas/vision/bigu/steve/src/reduce/dbl2byt deletemeb.mat 22 >> deletemeb/data
  saveraw('deletemer.tmp',R)
  saveraw('deletemeg.tmp',G)
  saveraw('deletemeb.tmp',B)
  !cat deletemer.tmp >> deletemer/data
  !cat deletemeg.tmp >> deletemeg/data
  !cat deletemeb.tmp >> deletemeb/data

end%for

disp('created the 3 datfiles (r,g,b); now concatenating to one colorfile')
!rm -r deleteme
!rm -r deleteme
!cp -r deletemer deleteme
!cat deletemeg/data >> deleteme/data
!cat deletemeb/data >> deleteme/data
!setkey deleteme _channels 3

disp('done, produced Fourieroperator in deleteme/')
