function outmap = encode_layer(in,model)
if model.numlayer==1 && model.type_zerolayer==3
   % the zeroth layer is the gradient map
   num_orients=model.smap;
   if num_orients == 1
      fprintf('size of feature maps should be greater than one\n');
      return;
   end
   [DY DX]=gradient(in);
   rho=sqrt(DY.^2+DX.^2);
   ind=find(rho(:));
   DY(ind) = DY(ind) ./  rho(ind);
   DX(ind) = DX(ind) ./  rho(ind);
   % theta is between [-pi pi]
   outmap=zeros(size(in,1),size(in,2),num_orients);
   subsampling_theta=2*pi/num_orients;
   sigma=sqrt((1-cos(subsampling_theta))^2 + sin(subsampling_theta)^2);
   tabtheta=0:subsampling_theta:2*pi-subsampling_theta;
   for ii=1:num_orients
      outmap(:,:,ii) = rho .* exp(-(1/(sigma*sigma)) * ((DX - cos(tabtheta(ii))).^2 + (DY - sin(tabtheta(ii))).^2));
   end
else
   X=mexExtractPatches(in,model.npatch,1);
   sx=size(in,1)-model.npatch+1;
   nchannels=size(in,3);
   % contrast_normalization of input maps
   if model.centering
      X=centering(X,nchannels);
   end
   nrm=sqrt(sum(X.^2));
   if isfield(model,'median_contrast_normalization') && model.median_contrast_normalization
      X=contrast_normalize_median(X);
   else
      X=contrast_normalize(X);
   end

   % convolution with filters 
   outmap=bsxfun(@plus, model.W' * X, model.b);
   % nonlinearities
   outmap=exp(outmap);
   nchannels=size(outmap,1);
   outmap=bsxfun(@times,outmap,nrm);
   outmap=outmap';
   outmap = reshape(outmap,[sx sx nchannels]);
end

% subsampling with Gaussian smoothing
if model.subsampling > 1
   outmap=subsampling(outmap,model.subsampling);
end

