# fields, Tools for spatial data
# Copyright 2004-2007, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html

"Wtransform.image.cov" <-
function (ind1, ind2=ind1,Y, cov.obj) 
{

#
# define useful local function that does the H multiplcation
# taking advantage of the block diagonal form 
#  

    Mult.H <- function(u, H.obj) {
        IND <- H.obj$ind0
        u[IND] <- H.obj$H0 %*% c(u[IND])
        u * H.obj$H1
    }

         IND <- cov.obj$ind0 
         CUT<- cov.obj$cut.min         

      if( missing(ind1)){ 
      # do multiplcation for full matrix. 
         if( (nrow(Y)!=cov.obj$m) | (ncol(Y)!=cov.obj$m)){ 
             stop("bad dimensions for Y")}
         hold <- Wtransform.image(Y,cut.min=CUT,inv=TRUE,transpose = TRUE)
         hold <- Mult.H(hold, cov.obj)
         hold <- Mult.H(hold, cov.obj)
         return(
              Wtransform.image(hold, cut.min = CUT, inv = TRUE))
        }

        else{
        # multiplication for subset of indices

         temp <- matrix(0, nrow = cov.obj$m, ncol = cov.obj$n)
         temp[ind2] <- Y 
#
         hold <- Wtransform.image(temp,cut.min=CUT,inv = TRUE,transpose=TRUE)
         hold <- Mult.H(hold, cov.obj)
         hold <- Mult.H(hold, cov.obj)
         return( 
              Wtransform.image(hold, cut.min = CUT,inv =TRUE)[ind1] )
         }

}