function [noisyField] = constructNoisyField(field, lats2D, longs2D, ...
                                            borderCoords, noiseLevel)
%constructNoisyFiled()    Function to construct a noisy field image,
% usually nightlights
%
% PROCEDURE
%    1. We consider the field (e.g., nighlights), and scale them so max 
%       field = 100
%    2. We zero-pad the field to be dyadic.
%    3. We wavelet transform it. The results is a dyadic matrix
%    4. We find the median wavelet magnitude and add a gaussian noise, of 
%       std = noiseLevel * median to all wavelets
%    5. we do an inverse wavelet transform and crop out the 27 x 64 matrix,
%       which is the reconstructed, noisy field
%    6. We zero out non-US areas (oceans)
%
%  INPUTS :
%   field       : a field, usually upscaled nightlights, at 1deg X 1deg resoution, 
%                 in the domain (corners) [51N, 127W]  [24, 63W]. This is a
%                 27 x 64 matrix.
%   lats2D      : 2D matrix, containg the lats of the cell-centers of the 2D
%   longs2D     : ditto, but longitudes
%   borderCoords: [long, lat] of points along US bndry, lower 48 states
%   levelsToZero: zero out the wavelets at the finest 'levelsToZero' levels.
%
%  OUTPUTS
%    noisyField : a noisy version of the field
%
% Jaideep Ray, 5/17/2012
% =========================================================================

% ---- Definition of zero
small = 1.0e-9 ;

% ---- Do a wavelet transform of it. But first make the matrix dyadic
dyadicSize1 = ceil( log(size(field, 1)) / log(2) ) ;
dyadicSize2 = ceil( log(size(field, 2)) / log(2) ) ;
dyadicSize = max( [dyadicSize1 dyadicSize2] ) ;
if ( (2^dyadicSize1 ~= size(field, 1)) || ...
     (2^dyadicSize2 ~= size(field, 2)) )
    paddedDataSize = [2^dyadicSize 2^dyadicSize] ;
else
    paddedDataSize = size(field) ;
end
D = zeros(paddedDataSize) ;   % make a dynadic array
D(1:size(field, 1), 1:size(field,2)) = field ; % and fill with data

% ---- Do the wavelet transform
coarseLevel = 0 ;
qmf = MakeONFilter('Haar', 2) ;
wc = FWT2_PO(D, coarseLevel, qmf) ;

% ---- The bulk of the wavelets in wc are zero. Find the max and scale up
% to 100. Then find the non-zero wavelets, find median and add noise
maxwc = max( wc(:) ) ; factor = 100.0 / maxwc ; wc = factor * wc ;
ii    =  abs(wc(:)) > small  ;
med   = median( abs( wc(ii) ) ) ;
fprintf('The median wavelet weight, from %d wavelets is %f\n', nnz(ii), ...
         med) ;
fprintf('Noise fraction = %f\n', noiseLevel);     
noise = noiseLevel * med * randn( size(wc, 1), size(wc, 2) ) ;
wc    = wc + noise ;

% ---- Get back the original nightlight image
D = IWT2_PO(wc, coarseLevel, qmf) ;
noisyField = D(1:size(field, 1), 1:size(field,2)) ;

% ---- Zero out non-US nightlights
inOrOut = inpolygon(longs2D(:), lats2D(:), borderCoords(:, 1), ...
                    borderCoords(:, 2)) ;
noisyField(~inOrOut) = 0 ;
end
