Patent Publication Number: US-8983206-B2

Title: Method and system for increasing robustness of visual quality metrics using spatial shifting

Description:
RELATED APPLICATIONS 
     The present application claims benefit from the U.S. provisional application Ser. No. 61/482,472 filed on May 4, 2011, entire contents of which are incorporated herein by reference. 
    
    
     FIELD OF THE INVENTION 
     The present invention relates to the assessment of image quality, and in particular, to a method and system for determining an improved measure of quality or metric for images with enhanced perceptual performance. 
     BACKGROUND OF THE INVENTION 
     Image quality assessment plays an important role in the development and validation of various image and video applications, such as compression and enhancement. Objective quality models are usually classified based on the availability of reference images. If an undistorted reference image is available, the quality metric is considered as a full reference (FR) assessment method. Several image and video applications require that the quality be assessed in real-time. Consequently, the computational complexity associated with a quality assessment method becomes an important issue to address. Another important factor to consider is the accuracy of the quality assessment method, which indicates how well computed quality scores correspond with the opinions of human beings. 
     Accordingly, a method and system to improve the accuracy of a quality assessment method of low complexity is required. 
     SUMMARY OF THE INVENTION 
     There is an object of the present invention to provide a method and system for increasing robustness of visual quality metrics using spatial shifting. 
     According to one aspect of the invention, there is provided a method of computing a measure of quality for a distorted image Y, the measure characterizing a similarity between the image Y and an undistorted reference image X having the same number of rows and columns of pixels as the image Y, the method comprising employing at least one processor for:
         (a) transforming the distorted image Y and the undistorted reference image X by a transformation process to generate a transform of the distorted image Y and a transform of the undistorted reference image X respectively, the transformation process having a blindness property of not containing all information regarding an image to be transformed in a transform of the image, thus resulting in certain patterns of the image not to be recognized during the transformation process;   (b) spatially shifting pixels in the reference image X and the distorted image Y by a spatial shift, having the same predefined amount of pixels, to generate a spatially shifted version X 2  of the reference image X and a spatially shifted version Y 2  of the image Y respectively, the spatial shift being chosen so that the certain patterns become detectable during the transformation process applied to the spatially shifted versions X 2  and Y 2 ;   (c) transforming the spatially shifted version Y 2  and the spatially shifted version X 2  by the transformation process to generate a transform of the spatially shifted version Y 2  and a transform of the spatially shifted version X 2  respectively; and   (d) determining the measure of quality for the distorted image Y as a function of the transform of the distorted image Y, the transform of the undistorted reference image X; the transform of the spatially shifted version Y 2 , and the transform of the spatially shifted version X 2 .       

     The step (b) further comprises:
         dividing the pixels of the reference image X into blocks of W by W pixels, W being a smaller number than the number of rows and columns of pixels; and   dividing the pixels of the distorted image Y into corresponding blocks of W by W pixels.       

     The step (b) further comprises:
         spatially shifting horizontally by a first predefined amount; or   spatially shifting vertically by a second predefined amount; or   spatially shifting horizontally and vertically by first and second predefined amounts respectively.       

     Conveniently, the first predefined amount and the second predefined amount may be chosen to be the same. 
     In an embodiment of the invention, W is equal to 2^N, wherein an exponent N is a decomposition level of the transformation process. 
     For example, each of the first and second predefined amounts may be chosen as the larger of 1 and (W/2−1). 
     In an embodiment of the invention, the first and second predefined amounts are odd numbers less than W. 
     In the method described above, the transformation process is one of the following: a cosine transform; a sine transform; or a wavelet transform. 
     In the embodiment of the invention, the transformation process is a Haar transform. Alternatively, the transformation process can be a Daubechies transform. 
     In the method described above, the step (d) further comprises:
         computing a first approximation image quality metric IQMA using the transform of the distorted image Y and the transform of the undistorted reference image X;   computing a second approximation image quality metric IQM′A using the transform of the spatially shifted version Y 2  and the transform of the spatially shifted version X 2 ; and   processing the first and second approximation image quality metrics IQMA and IQM′A to generate the measure of quality for the distorted image Y.       

     In the embodiment of the invention, the first and second approximation quality metrics IQMA and IQM′A respectively are first and second approximation quality maps, and the step (d) further comprises:
         creating first and second contrast maps of the reference image X and of the spatially shifted version X 2  respectively;   performing weighted pooling of the first and second approximation quality maps using the first and second contrast maps respectively to produce an approximation quality score SA; and   determining the measure of quality using the approximation quality score SA.       

     The step (d) further comprises:
         (i) computing a first edge image quality metric IQME using the transform of the distorted image Y and the transform of the undistorted reference image X;   (ii) computing a second edge image quality metric IQM′E using the transform of the spatially shifted version Y 2  and the transform of the spatially shifted version X 2 ; and   (iii) determining the measure of quality for the distorted image Y as the function of IQMA, IQM′A, IQME and IQM′E.       

     The step (iii) further comprises:
         processing the first and second approximation image quality metrics IQMA and IQM′A to generate an approximation quality score SA;   processing the first and second edge quality metrics IQME and IQM′E to generate an edge quality score SE; and   determining the measure of quality for the distorted image Y as a function of the approximation quality score SA and the edge quality score SE.       

     In the embodiment of the invention, the measure of quality is computed as a final score SF=a*SA+(1−a)*SE, with SF value ranging from 0 to 1. 
     In the method described above, the first and second approximation quality metrics IQMA and IQM′A respectively are first and second approximation quality maps, and the first and second edge quality metrics IQME and IQM′E respectively are first and second edge quality maps, and the step (d) further comprises:
         (d1) creating first and second contrast maps of the reference image X and of the spatially shifted version X 2  respectively;   (d2) performing weighted pooling of the first and second approximation quality maps using the first and second contrast maps respectively to produce an approximation quality score SA;   (d3) performing weighted pooling of the first and second edge quality maps using the first and second contrast maps respectively to produce an edge quality score SE; and   (d4) processing the approximation quality score SA and the edge quality score SE to determine the measure of quality.       

     The step of creating first and second contrast maps further comprises assigning values to pixels of the first and second approximation quality maps, and the first and second edge maps of the reference image X and the spatially shifted version X 2  according to their respective importance to a human visual system. 
     The step (d2) further comprises one of averaging the first approximation quality map and the second approximation quality map, and of taking the minimum of the first approximation quality map and the second approximation quality map to produce the approximation quality score SA. 
     According to another aspect of the invention, there is provided a system for computing a measure of quality for a distorted image Y, the measure characterizing a similarity between the image Y and an undistorted reference image X having the same number of rows and columns of pixels as the image Y, the system comprising: 
     a processor; 
     a non-transitory computer readable storage medium having computer readable instructions stored thereon for execution by the processor, causing the processor to:
         (a) transform the distorted image Y and the undistorted reference image X by a transformation process to generate a transform of the distorted image Y and a transform of the undistorted reference image X respectively, the transformation process having a blindness property of not containing all information regarding an image to be transformed in a transform of the image, thus resulting in certain patterns of the image not to be recognized during the transformation process;   (b) spatially shift pixels in the reference image X and the distorted image Y by a spatial shift, having the same predefined amount of pixels, to generate a spatially shifted version X 2  of the reference image X and a spatially shifted version Y 2  of the image Y respectively, the spatial shift being chosen so that the certain patterns become detectable during the transformation process applied to the spatially shifted versions X 2  and Y 2 ;   (c) transform the spatially shifted version Y 2  and the spatially shifted version X 2  by the transformation process to generate a transform of the spatially shifted version Y 2  and a transform of the spatially shifted version X 2  respectively; and   (d) determine the measure of quality for the distorted image Y as a function of the transform of the distorted image Y, the transform of the undistorted reference image X, the transform of the spatially shifted version Y 2 , and the transform of the spatially shifted version X 2 .       

     In the system described above, the computer readable instructions further cause the processor to:
         divide the pixels of the reference image X into blocks of W by W pixels, W being a smaller number than the number of rows and columns of pixels; and   divide the pixels of the distorted image Y into corresponding blocks of W by W pixels.       

     The computer readable instructions further cause the processor to:
         spatially shift horizontally by a first predefined amount; or   spatially shift vertically by a second predefined amount; or   spatially shift horizontally and vertically by first and second predefined amounts respectively.       

     In the system described above, the first predefined amount and the second predefined amount are the same. 
     For example, W may be equal to 2^N, wherein an exponent N is a decomposition level of the transformation process; and each of the first and second predefined amounts may be chosen as the larger of 1 and (W/2−1). Alternatively, the first and second predefined amounts may be chosen as odd numbers less than W. 
     In the system described above, the transformation process is one of the following: a cosine transform; a sine transform; or a wavelet transform. 
     In the embodiments of the invention, the transformation process is a Haar transform. Alternatively, the transformation process may be a Daubechies transform. 
     In the system described above, the computer readable instructions further cause the processor to:
         compute a first approximation image quality metric IQMA using the transform of the distorted image Y and the transform of the undistorted reference image X;   compute a second approximation image quality metric IQM′A using the transform of the spatially shifted version Y 2  and the transform of the spatially shifted version X 2 ; and   process the first and second approximation image quality metrics IQMA and IQM′A to generate the measure of quality for the distorted image Y.       

     The first and second approximation quality metrics IQMA and IQM′A respectively are first and second approximation quality maps, and the computer readable instructions further cause the processor to:
         create first and second contrast maps of the reference image X and of the spatially shifted version X 2  respectively;   perform weighted pooling of the first and second approximation quality maps using the first and second contrast maps respectively to produce an approximation quality score SA; and   determine the measure of quality using the approximation quality score SA.       

     In the system described above, the computer readable instructions further cause the processor to:
         (i) compute a first edge image quality metric IQME using the transform of the distorted image Y and the transform of the undistorted reference image X;   (ii) compute a second edge image quality metric IQM′E using the transform of the spatially shifted version Y 2  and the transform of the spatially shifted version X 2 ; and   (iii) determine the measure of quality for the distorted image Y as the function of IQMA, IQM′A, IQME and IQM′E.       

     The computer readable instructions further cause the processor to:
         process the first and second approximation image quality metrics IQMA and IQM′A to generate an approximation quality score SA;   process the first and second edge quality metrics IQME and IQM′E to generate an edge quality score SE; and   determine the measure of quality for the distorted image Y as a function of the approximation quality score SA and the edge quality score SE.       

     In the system described above, the measure of quality is computed as a final score SF=a*SA+(1−a)*SE, with SF value ranging from 0 to 1. 
     In the embodiments of the invention, the first and second approximation quality metrics IQMA and IQM′A respectively are first and second approximation quality maps, the first and second edge quality metrics IQME and IQM′E respectively are first and second edge quality maps, and the computer readable instructions further cause the processor to:
         create first and second contrast maps of the reference image X and of the spatially shifted version X 2  respectively;   perform weighted pooling of the first and second approximation quality maps using the first and second contrast maps respectively to produce an approximation quality score SA;   perform weighted pooling of the first and second edge quality maps using the first and second contrast maps respectively to produce an edge quality score SE; and   process the approximation quality score SA and the edge quality score SE to determine the measure of quality.       

     The computer readable instructions further cause the processor to assign values to pixels of the first and second approximation quality maps, and the first and second edge maps of the reference image X and the spatially shifted version X 2  according to their respective importance to a human visual system. 
     The computer readable instructions further cause the processor to perform one of averaging the first approximation quality map and the second approximation quality map, and of taking the minimum of the first approximation quality map and the second approximation quality map to produce the approximation quality score SA. 
     According to yet another aspect of the invention, there is provided a system for computing a measure of quality for a distorted image Y, the measure characterizing a similarity between the image Y and an undistorted reference image X, the system comprising:
         a processor;   a memory having computer readable instructions stored thereon for execution by the processor, causing the processor to process the distorted image Y and an undistorted reference image X, comprising:   a Transform module transforming the distorted image Y and the undistorted reference image X by a transformation process to generate a transform of the distorted image Y and a transform of the undistorted reference image X respectively, the transformation process having a blindness property of not containing all information regarding an image to be transformed in a transform of the image, thus resulting in certain patterns of the image not to be recognized during the transformation process;   a Spatial Shift module for generating spatially shifted versions Y 2  and X 2  of the distorted image Y and the reference image X respectively;   the Transform module being further configured to transform the spatially shifted version Y 2  and the spatially shifted version X 2  by the transformation process to generate a transform of the spatially shifted version Y 2  and a transform of the spatially shifted version X 2  respectively   an Image Quality Metric (IQM) computation subsystem for determining the measure of quality for the distorted image Y as a function of the transform of the distorted image Y, the transform of the undistorted reference image X, the transform of the spatially shifted version Y 2 , and the transform of the spatially shifted version X 2 .       

     The IQM computation subsystem further comprises an Approximation Quality Map/Score Computation Module for:
         computing a first approximation image quality metric IQMA using the transform of the distorted image Y and the transform of the undistorted reference image X; and   computing a second approximation image quality metric IQM′A using the transform of the spatially shifted version Y 2  and the transform of the spatially shifted version X 2 .       

     The system further comprises a Final Score Computation module for combining the first and second approximation image quality metrics IQMA and IQM′A to generate the measure of quality for the distorted image Y. 
     The system further comprises an Edge Quality Map/Score Computation Module for:
         computing a first edge image quality metric IQME using the transform of the distorted image Y and the transform of the undistorted reference image X; and   computing a second edge image quality metric IQM′E using the transform of the spatially shifted version Y 2  and the transform of the spatially shifted version X 2 ;
 
the Final Score Computation module being further configured to determine the measure of quality for the distorted image Y as the function of IQMA, IQMA′A, IQME and IQM′E.
       

     The system further comprises:
         a Common Pooling of Approximation Quality Map Module for combining the first and second approximation image quality metrics IQMA and IQM′A to generate an approximation quality score SA;   a Common Pooling of Edge Map Module for combining the first and second edge quality metrics IQME and IQM′E to generate an edge quality score SE; and   the Final Score Computation module being further configured to determine the measure of quality for the distorted image Y as a function of the approximation quality score SA and the edge quality score SE.       

     In the system described above, the function of the approximation quality score SA and the edge quality score SE is computed as a final score SF=a*SA+(1−a)*SE, with SF value ranging from 0 to 1. 
     In the embodiments of the invention, the first and second approximation quality metrics IQMA and IQM′A respectively are first and second approximation quality maps, and the system further comprises:
         a Contrast Map Generation Module for creating first and second contrast maps of the reference image X and of the spatially shifted version X 2  respectively, and   the Common Pooling of Approximation Quality Map Module being further configured to perform weighted pooling of the first and second approximation quality maps using the first and second contrast maps respectively to produce the approximation quality score SA.       

     In the system described above, the first and second approximation quality metrics IQMA and IQM′A respectively are first and second approximation quality maps, the first and second edge quality metrics IQME and IQM′E respectively are first and second edge quality maps, the system further comprising:
         a Contrast Map Generation Module for creating first and second contrast maps of the reference image X and of the spatially shifted version X 2  respectively;   the Common Pooling of Approximation Quality Map Module being further configured to perform weighted pooling of the first and second approximation quality maps using the first and second contrast maps respectively to produce the approximation quality score SA, and to perform weighted pooling of the first and second edge quality maps using the first and second contrast maps respectively to produce the edge quality score SE.       

     The Contrast Map Generation Module is further configured to create first and second contrast maps, including assigning corresponding values to pixels of the first and second approximation quality maps, and the first and second edge maps of the reference image X and the spatially shifted version X 2  according to their respective importance to a human visual system. 
     The Final Score Computation module is further configured to perform one of:
         averaging the approximation quality score SA and the edge quality score SE; and   taking the minimum of the approximation quality score SA and the edge quality score SE to determine the measure of quality.       

     The Transform Module is further configured to:
         divide the pixels of the reference image X into blocks of W by W pixels, W being a smaller number than the number of rows and columns of pixels; and   divide the pixels of the distorted image Y into corresponding blocks of W by W pixels.       

     The Spatial Shift module is further configured to:
         spatially shift horizontally by a first predefined amount; or   spatially shift vertically by a second predefined amount; or   spatially shift horizontally and vertically by first and second predefined amounts respectively.       

     Thus, an improved method and systems for increasing robustness of visual quality metrics using spatial shifting have been provided. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       Embodiments of the invention will now be described, by way of example, with reference to the accompanying drawings in which: 
         FIG. 1  shows a block diagram of a discrete wavelet domain image quality assessment framework  10  for determining a measure of quality; 
         FIG. 1A  shows a flowchart of a corresponding method  140  for determining the measure of quality; 
         FIG. 2  shows an example diagram  200  depicting wavelet sub-bands for a two-level decomposition of the reference image X (N=2) of  FIG. 1 ; 
         FIG. 3A  shows a typical original image; 
         FIG. 3B  demonstrates the resized contrast map obtained by eq. (9) for the original image of  FIG. 3A ; 
         FIG. 4  shows a Matlab code fragment  900  for generating and displaying the 4×4DCT basis used in H.264compression; 
         FIG. 5  shows a number array  1000  illustrating results of the 4×4DCT basis from  FIG. 4 , calculated at an arbitrary scale; 
         FIG. 6  illustrates a first quality metrics computation system  1100 , including two instances ( 1102 . 1  and  1102 . 2 ) of a Quality Assessment module  1102 , a Spatial Shift module  1104 , and a Final Score Computation module  1106 ; 
         FIG. 7  illustrates a second quality metrics computation system  1200  including two instances ( 1102 . 1  and  1102 . 2 ) of the Quality Assessment module  1102 , two instances ( 1104 . 1  and  1104 . 2 ) of the Spatial Shift module  1104 , and the Final Score Computation module  1106 ; 
         FIG. 8  illustrates a third quality metrics computation system  1300  including the Quality Assessment module  1102  and the Spatial Shift module  1104 ; 
         FIG. 9  shows a block diagram of an embodiment of a discrete wavelet domain image quality assessment system  1400 ; 
         FIG. 10  shows a flowchart of a method  1500  of computing a combined score which is a measure of quality of a distorted image by comparison with a corresponding original image, combined with a comparison of a spatially shifted version of the distorted image with an equally spatially shifted version of the corresponding original image, in a sequence chart associated with the system block diagram of  FIG. 9 ; 
         FIG. 10A  shows an expansion of the five-step group  1502  of  FIG. 10 ; 
         FIG. 10B  shows an expansion of the five-step group  1508  of  FIG. 10 ; 
         FIG. 11(   a ) shows a standard image “Image House” as an original reference image; 
         FIG. 11(   b ) shows a first distorted image created from the reference image of  FIG. 11(   a ) with 2 nd  row and 2 nd  column of the 4×4DCT transformed image set to zero; 
         FIG. 11(   c ) shows a second distorted image created from the reference image of  FIG. 11(   a ) with 3 rd  row and 3 rd  column of the 4×4DCT transformed image set to zero; 
         FIG. 11(   d ) shows a third distorted image created from the reference image of  FIG. 11(   a ) with 4 th  row and 4 th  column of 4×4DCT transformed image set to zero; and 
         FIG. 12  shows an exemplary hardware view  1700  of the discrete wavelet domain image quality assessment system  1400  including a Processor  1702  and a Computer Memory  1704  according to the embodiments of the invention. 
     
    
    
     DETAILED DESCRIPTION OF THE EMBODIMENTS OF THE INVENTION 
     Embodiments of the present invention provide a quality assessment method which aims to achieve a good compromise between good accuracy and low complexity. 
     Glossary 
     
       
         
           
               
               
             
               
                   
               
             
            
               
                 Approximation-map 
                 characterizes the main content of an image and is determined 
               
               
                 (or approximation quality map): 
                 as the approximation sub-band of the image obtained after 
               
               
                   
                 one-level (or higher level) discrete wavelet decomposition. 
               
               
                 Approximation sub-band: 
                 sub-band containing main content of an image, i.e. the 
               
               
                   
                 low frequency content. 
               
               
                 Approximation quality score: 
                 measure of similarity of the approximation sub-band of 
               
               
                   
                 image X and the approximation sub-band of image Y. 
               
               
                 Approximation quality measure: 
                 measure of image quality that provides an approximation 
               
               
                   
                 quality score when a map-based metric is used for quality 
               
               
                   
                 assessment (e.g., VIF) or an approximation quality map (also 
               
               
                   
                 called approximation-map) when a non map-based metric is 
               
               
                   
                 used for quality assessment (e.g., SSIM W&amp;B ). 
               
               
                 Contrast map: 
                 weighting function for automatically assigning weights to 
               
               
                   
                 pixels in different regions of an image based on their visual 
               
               
                   
                 importance. 
               
               
                 Correlation coefficient: 
                 a statistical measure computed between two image patches 
               
               
                   
                 based on Pearson product-moment correlation coefficient. 
               
               
                 DC: 
                 (Direct Current) indicates the zero-frequency component of a 
               
               
                   
                 transform, i.e. the average value. 
               
               
                 Detail sub-bands: 
                 sub-bands containing fine edges of an image. 
               
               
                 Diagonal detail sub-band: 
                 sub-band containing fine diagonal edges of an image. 
               
               
                 Discrete wavelet transform (DWT) 
                 transform applied to an image to separate its low frequency 
               
               
                   
                 components from its high frequency components. 
               
               
                 Discrete Cosine Transform (DCT) 
                 discrete cosine transform used for multiresolution 
               
               
                   
                 decomposition of images. 
               
               
                 Edge-map (or edge quality map): 
                 characterizes the fine edges of an image and is determined as 
               
               
                   
                 the square root of the mean square of horizontal, vertical, and 
               
               
                   
                 diagonal detail sub-bands of the image obtained after one-level 
               
               
                   
                 (or higher level) discrete wavelet decomposition. 
               
               
                 Edge quality score: 
                 measure of similarity between the detail sub-bands of two 
               
               
                   
                 images. 
               
               
                 Edge quality measure: 
                 measure of image quality that provides an edge quality score 
               
               
                   
                 when a map-based metric is used for quality assessment (e.g., 
               
               
                   
                 VIF) or an edge quality map (also called edge-map) when a 
               
               
                   
                 non map-based metric is used for quality assessment (e.g., 
               
               
                   
                 SSIM W&amp;B ). 
               
               
                 Gaussian sliding window: 
                 a set of N coefficients with unit sum and Gaussian probability 
               
               
                   
                 distribution. 
               
               
                 Horizontal sub-band: 
                 sub-band containing fine horizontal edges of an image. 
               
               
                 Multiresolution decomposition: 
                 a transform applied on digital images for generating sub- 
               
               
                   
                 bands, one of which contains the low frequency main content 
               
               
                   
                 of original image and other sub-bands contain fine edges of the 
               
               
                   
                 original image. 
               
               
                 S A  and S E : 
                 approximation and edge quality scores respectively. 
               
               
                 S f : 
                 Final score. 
               
               
                 Structural SIMilarity (SSIM): 
                 an accurate measure of image quality that captures similarity 
               
               
                   
                 between two images and can be computed with high, but still 
               
               
                   
                 acceptable, computational complexity compared to other 
               
               
                   
                 image quality metrics such as VIF. 
               
               
                 SSIM W&amp;B : 
                 SSIM method of quality assessment proposed by Z. Wang and 
               
               
                   
                 A. C. Bovik. 
               
               
                 Quality score: 
                 an intermediary similarity value, which is used for calculating 
               
               
                   
                 the final score. There are two kinds of quality scores: 
               
               
                   
                 approximation quality score and edge quality score. 
               
               
                 Quality map: 
                 an intermediary similarity value, which is used for calculating 
               
               
                   
                 the quality measure in map-based visual quality assessment 
               
               
                   
                 methods (e.g., when SSIM W&amp;B  is used for quality assessment). 
               
               
                   
                 There are two kinds of quality maps: approximation quality 
               
               
                   
                 maps (or approximation-maps) and edge quality maps (or 
               
               
                   
                 edge-maps). 
               
               
                 Visual Information Fidelity (VIF): 
                 a method of quality assessment for images proposed by H. R. 
               
               
                   
                 Sheikh and A. C. Bovik. 
               
               
                 Vertical sub-band: 
                 sub-band containing fine vertical edges of an image. 
               
               
                   
               
            
           
         
       
     
     The peak signal-to-noise ratio (PSNR) is among the oldest and the most widely used FR image quality evaluation measure, because it is simple, has clear physical meaning, is parameter-free, and performs superbly in an optimization context. However, conventional PSNR cannot adequately reflect the human perception of image fidelity, that is, a large PSNR gain may correspond to only a small improvement in visual quality. This has led researchers to develop a number of other quality measures. Tables I and II show the accuracy of PSNR against other quality metrics, see H. R. Sheikh, Z. Wang, L. Cormack, and A. C. Bovik, “LIVE Image Quality Assessment Database Release 2,” available at: http://live.ece.utexas.edu/research/quality. This release of the LIVE Image Quality Assessment Database consists of 779 distorted images derived from 29 original colour images using five types of distortion: JPEG compression, JPEG2000 compression, Gaussian white noise (GWN), Gaussian blurring (GBlur), and the Rayleigh fast fading (FF) channel model. The realigned subjective quality data for the database are used in all experiments, see H. R. Sheikh, Z. Wang, L. Cormack, and A. C. Bovik, “LIVE Image Quality Assessment Database Release 2”. The metric used in Table I is the Pearson correlation coefficient (LCC) between the Difference Mean Opinion Score (DMOS) and the objective model outputs after nonlinear regression. The correlation coefficient gives an evaluation of prediction accuracy. We use the five-parameter logistical function defined in H. R. Sheikh, M. F. Sabir, and A. C. Bovik, “A statistical evaluation of recent full reference image quality assessment algorithms,”  IEEE Trans. Image Process. , vol. 15, no. 11, pp. 3440-3451, November 2006 for nonlinear regression. The metric used in Table II is the Spearman rank correlation coefficient (SRCC), which provides a measure of prediction monotonicity. 
     It can be seen that PSNR (denoted PSNR spatial  in the tables) has a significantly lower correlation coefficient (LCC) and Spearman Rank Order (SRCC) compared to other visual quality metrics such as Structural SIMilarity (SSIM) spatial and autoscale as well as Visual Information Fidelity (VIF). Unfortunately, these more accurate metrics are obtained at the cost of very high computational complexity. 
     
       
         
           
               
             
               
                 TABLE I 
               
             
            
               
                   
               
               
                 LCC VALUES AFTER NONLINEAR REGRESSION 
               
               
                 FOR THE LIVE IMAGE DATABASE 
               
            
           
           
               
               
               
               
               
               
               
            
               
                   
                   
                   
                   
                   
                   
                 All 
               
               
                 Model 
                 JPEG 
                 JPEG2000 
                 GWN 
                 GBlur 
                 FF 
                 Data 
               
               
                   
               
            
           
           
               
               
               
               
               
               
               
            
               
                 SSIM spatial   
                 0.9504 
                 0.9413 
                 0.9747 
                 0.8743 
                 0.9449 
                 0.9038 
               
               
                 SSIM autoscale   
                 0.9778 
                 0.9669 
                 0.9808 
                 0.9483 
                 0.9545 
                 0.9446 
               
               
                 mean 
                 0.9762 
                 0.9699 
                 0.9645 
                 0.9548 
                 0.9625 
                 0.9412 
               
               
                 SSIM A   
               
               
                 S A   
                 0.9782 
                 0.9705 
                 0.9724 
                 0.9724 
                 0.9730 
                 0.9534 
               
               
                 SSIM DWT   
                 0.9835 
                 0.9747 
                 0.9791 
                 0.9690 
                 0.9735 
                 0.9556 
               
               
                 PSNR spatial   
                 0.8879 
                 0.8996 
                 0.9852 
                 0.7835 
                 0.8895 
                 0.8701 
               
               
                 wSNR 
                 0.9692 
                 0.9351 
                 0.9776 
                 0.9343 
                 0.8983 
                 0.9211 
               
               
                 PSNR A   
                 0.9793 
                 0.9542 
                 0.9806 
                 0.9241 
                 0.8868 
                 0.9288 
               
               
                 PSNR DWT   
                 0.9787 
                 0.9549 
                 0.9838 
                 0.9234 
                 0.8994 
                 0.9300 
               
               
                 AD A   
                 0.9817 
                 0.9587 
                 0.9637 
                 0.9307 
                 0.9005 
                 0.9350 
               
               
                 AD DWT   
                 0.9807 
                 0.9579 
                 0.9678 
                 0.9258 
                 0.9064 
                 0.9344 
               
               
                 VIF 
                 0.9864 
                 0.9773 
                 0.9901 
                 0.9742 
                 0.9677 
                 0.9593 
               
               
                 DWT-VIF A   
                 0.9856 
                 0.9735 
                 0.9904 
                 0.9615 
                 0.9611 
                 0.9639 
               
               
                 DWT-VIF 
                 0.9852 
                 0.9740 
                 0.9906 
                 0.9652 
                 0.9650 
                 0.9654 
               
               
                   
               
            
           
         
       
     
     
       
         
           
               
             
               
                 TABLE II 
               
             
            
               
                   
               
               
                 SRCC VALUES AFTER NONLINEAR REGRESSION 
               
               
                 FOR THE LIVE IMAGE DATABASE 
               
            
           
           
               
               
               
               
               
               
               
            
               
                   
                   
                   
                   
                   
                   
                 All 
               
               
                 Model 
                 JPEG 
                 JPEG2000 
                 GWN 
                 GBlur 
                 FF 
                 Data 
               
               
                   
               
            
           
           
               
               
               
               
               
               
               
            
               
                 SSIM spatial   
                 0.9449 
                 0.9355 
                 0.9629 
                 0.8944 
                 0.9413 
                 0.9104 
               
               
                 SSIM autoscale   
                 0.9764 
                 0.9614 
                 0.9694 
                 0.9517 
                 0.9556 
                 0.9479 
               
               
                 mean 
                 0.9738 
                 0.9634 
                 0.9490 
                 0.9620 
                 0.9622 
                 0.9441 
               
               
                 SSIM A   
               
               
                 S A   
                 0.9779 
                 0.9634 
                 0.9577 
                 0.9703 
                 0.9699 
                 0.9573 
               
               
                 SSIM DWT   
                 0.9819 
                 0.9678 
                 0.9683 
                 0.9707 
                 0.9708 
                 0.9603 
               
               
                 PSNR spatial   
                 0.8809 
                 0.8954 
                 0.9854 
                 0.7823 
                 0.8907 
                 0.8756 
               
               
                 wSNR 
                 0.9610 
                 0.9292 
                 0.9749 
                 0.9330 
                 0.8990 
                 0.9240 
               
               
                 PSNR A   
                 0.9647 
                 0.9499 
                 0.9777 
                 0.9219 
                 0.8853 
                 0.9307 
               
               
                 PSNR DWT   
                 0.9648 
                 0.9494 
                 0.9818 
                 0.9230 
                 0.9004 
                 0.9325 
               
               
                 AD A   
                 0.9666 
                 0.9553 
                 0.9805 
                 0.9335 
                 0.9067 
                 0.9421 
               
               
                 AD DWT   
                 0.9661 
                 0.9546 
                 0.9835 
                 0.9290 
                 0.9131 
                 0.9412 
               
               
                 VIF 
                 0.9845 
                 0.9696 
                 0.9858 
                 0.9726 
                 0.9649 
                 0.9635 
               
               
                 DWT-VIF A   
                 0.9837 
                 0.9669 
                 0.9848 
                 0.9618 
                 0.9597 
                 0.9663 
               
               
                 DWT-VIF 
                 0.9829 
                 0.9680 
                 0.9853 
                 0.9657 
                 0.9641 
                 0.9681 
               
               
                   
               
            
           
         
       
     
     In the previous patent applications Ser. Nos. 12/906,111 and 12/906,112 filed on Oct. 17, 2010 of the same inventors, a framework using the Haar transform is disclosed which provides a good compromise between computational complexity and accuracy. For instance PSNR A , S A , DWT_VIF A  are examples of such metrics developed using the framework. As shown in Tables I and II, PSNR A  has much higher accuracy than PSNR. Still, PSNR A  has a computational complexity comparable to PSNR. The “A” subscript in the above metrics indicate that only approximation sub-bands are used to assess the quality. To reduce computations, the 2×2 Haar transform has been used as a transformation process. With the 2×2 Haar transform, the approximation sub-band is obtained by computing the average of the pixel intensity within each W×W block (where the block dimension W=2^N with “N” being the decomposition level). Metrics developed from the averages of each W×W block in the reference and the distorted images, are then used in computing a quality score. Computing the quality is thus a block-based process. Typical choices for the decomposition level N are, for a 256×256 image, N=1, and for a 512×512 image, N=2. 
     However, it can be observed that this averaging process creates a “blindness” in these metrics due to the fact that any distortion (however bad it might be) in an image that does not affect the average over each W×W block remains hidden. As a result, the distortion becomes undetectable. In other words, when considering the approximation band, if the original and the distorted images have the same pixel intensity average for each W×W block, then no distortion will be detected although there could be severe image degradation. Note that in this document, we use the 2×2 Haar transform, but other sizes of the Haar transform could be used (e.g. 3×3, 4×4, etc.), also creating “blindness” problems. When the context is clear, we will use Haar transform to refer to the 2×2 Haar transform. A person skilled in the art of digital signal processing will know how to adapt the invention presented herein to other Haar transform sizes. 
     However, as tested over large image and video databases such as the LIVE image database used to generate Table I and II, for most common distortions such as inherent in JPEG and JPEG 2000 encoding, Gaussian blur, fast fading and H.264compression, this does not typically become a problem and the metrics perform very well. This is due to the fact that, usually, if an image is distorted, the average will be affected as well as the high frequencies in a W×W block. 
     Unfortunately, the proposed metrics, using the 2×2 Haar transform, are blind to certain types of distortions. For instance, consider the operation of modifying (e.g. quantizing or suppressing) the frequency coefficients of a 4×4DCT transformed block, which is a common operation in H.264video compression. For a decomposition level of N=2, leading to a Haar transform operating on 4×4blocks, altering any of these frequency coefficients, except the DC, will not have any impact on the measured quality since only the DC of the 4×4block is considered in the metric computation. This observation is a consequence of the definition of 4×4DCT, namely that none of the frequency coefficients except for the DC have any impact on the average. 
     Furthermore, using again a 4×4DCT transformed block, we can observe that, due to symmetry of the 4×4DCT, altering certain frequency coefficients in the block will not have any influence on the DC of the 2×2 Haar transform. Consequently, a metric based on a Haar transform with a decomposition level of 1 would be blind to changes of these other frequency coefficients or distortions. 
     A complete description of a DWT-based framework for computing a general purpose full reference (FR) image quality metric (IQM) may also be found in the appendix. 
       FIG. 1  shows a block diagtam of a discrete wavelet domain image quality assessment framework  10  for determining a measure of quality IQM DWT . The discrete wavelet domain image quality assessment framework  10  includes:
         an IQM computation subsystem  100  comprising:   first and second N-level DWT modules  102  and  104  respectively,   first and second Detail Sub-band Aggregation Modules  106  and  108 ,   a Contrast Map Generation Module  110 ,   an Approximation Quality Map/Score Computation Module  112 , and   an Edge Quality Map/Score Computation Module  114 ; as well as   a Pooling of Approximation Quality Map Module  116 ;   a Pooling of Edge Map Module  118 ; and   a Combination Module  120 .       

       FIG. 1A  shows a flowchart of a corresponding method  140  for determining the measure of quality IQM DWT , including steps: 
       150  “Compute IQMs” which is a grouping of five steps:
           152  “Perform N-level DWTs”;     154  “Calculate Approximation Quality Map (or Score)”;     156  “Form Estimate of Image Edges”;     158  “Compute Edge Quality Map (or Score)”; and     160  “Conditionally Generate Contrast Map”;       

       170  “Conditionally Perform weighted Pooling”; and 
       180  “Combine Approximation and Edge Scores”. 
     The five steps ( 152  to  160 ) of the grouping  150  “Compute IQMs” are performed by the IQM computation subsystem  100  of  FIG. 1 . 
     In the step  152  “Perform N-level DWTs”, N-level DWT operations are preformed on both the reference image X and the distorted image Y in the first and second N-level DWT modules  102  and  104  respectively. The N-level DWT is based on the Haar wavelet filter. With N-level decomposition, the approximation sub-bands X A     N    and Y A     N   , as well as first and second detail sub-bands, are obtained. 
     The Haar wavelet has been used previously in some quality assessment and compression methods reported in M. R. Bolin and G. W. Meyer, “A visual difference metric for realistic image synthesis,” in  Proc. SPIE Human Vision, Electron. Imag. , vol. 3644, San Jose, Calif., 1999, pp. 106-120, and in Y.-K. Lai and C.-C. J. Kuo, “A Haar wavelet approach to compressed image quality measurement,”  J. Visual Commun. Imag. Represent. , vol. 11, no. 1, pp. 17-40, March 2000. For our framework, we chose the Haar filter for its simplicity and good performance. The Haar wavelet has very low computational complexity compared to other wavelets. In addition, based on our simulations, it provides more accurate quality scores than other wavelet bases. The reason for this is that symmetric Haar filters have a generalized linear phase, so the perceptual image structures can be preserved. Also, Haar filters can avoid over-filtering the image, owing to their short filter length. The number of decomposition levels (N) selected for structural or information-theoretic strategies, such as SSIM or VIF, is equal to one. The reason for this is that, for more than one level of decomposition, the resolution of the approximation sub-band is reduced exponentially and it becomes very small. Consequently, a large number of important image structures or information will be lost in that sub-band. But, for error-based approaches, like PSNR or absolute difference (AD), we can formulate the required decomposition levels N as follows: when an image is viewed at distance d from a display of height h, we have (as reported in Y. Wang, J. Ostermann, and Y-Q. Zhang,  Video Processing and Communications . New Jersey: Prentice-Hall, 2002): 
                 N   ≥   0     ⇒   N     =     max   (     0   ,     round   (       log   2     (       min   ⁡     (     H   ,   W     )         (     344   /   k     )       )     )       )           
(5)
 
     Where the image is of size H×W and the viewing distance is d=k×h, and k is a constant factor typically set to 3, as typical viewers are positioned at a distance of approximately 3 times the picture height. 
     For a derivation of “N”, please refer to an article by S. Rezazadeh and S. Coulombe, “A novel discrete wavelet transform framework for full reference image quality assessment” in Signal, Image and Video Processing, pp. 1-15, September 2011, which is incorporated herein by reference. 
     In the step  154  “Calculate Approximation Quality Map (or Score)”, which is performed in the Approximation Quality Map/Score Computation Module  112 , the quality map (or score) is calculated by applying IQM between the approximation sub-bands of X A     N    and Y A     N   , which will be called an approximation quality map (or score) IQM A . We have a quality map for map-based metrics such as SSIM. We have quality score otherwise (e.g. VIF). Examples of IQM computations applied to various quality metrics, such as SSIM and VIF, are presented further below. 
     In the step  156  “Form Estimate of Image Edges”, which is performed in the first and second Detail Sub-band Aggregation Modules  106  and  108 , first and second estimates X E  and Y E  of the image edges are formed for the reference image X and the distorted image Y respectively. The image using an aggregate of first detail sub-bands, where the aggregation of the first detail sub-bands is performed in the first Detail Sub-band Aggregation Modules  106 . If the N-level DWT is applied to the images, the edge map (estimate) of image X is defined as: 
                       X   E     ⁡     (     m   ,   n     )       =       ∑     L   =   1     N     ⁢       X     E   ,   L       ⁡     (     m   ,   n     )                 (   6   )               
where X E  is the edge map of X; and X E,L  is the image edge map at decomposition level L, computed as defined in eq. (7). In eq. (7), x H     L   , x V     L   , and x D     L    denote the horizontal, vertical, and diagonal detail subbands obtained at the decomposition level L for image X respectively. x H     L     ,A     N-L   , x V     L     ,A     N-L   , and x D     L     ,A     N-L    are the wavelet packet approximation subbands obtained by applying an (N-L)-level DWT on x H     L   , x V     L   , and x D     L    respectively. The parameters μ, λ, and ψ are constant. As the HVS is more sensitive to the horizontal and vertical subbands and less sensitive to the diagonal one, greater weight is given to the horizontal and vertical subbands. We arbitrarily propose μ=λ=4.5ψ in this document, which results in μ=λ=0.45 and ψ=0.10 to satisfy eq. (8).
 
                       X     E   ,   L       ⁡     (     m   ,   n     )       =     {               μ   ·       (       X     H   L       ⁡     (     m   ,   n     )       )     2       +       λ   ⁡     (       X     V   L       ⁡     (     m   ,   n     )       )       2     +       ψ   ⁡     (       X     D   L       ⁡     (     m   ,   n     )       )       2             if         L   =   N                   μ   ·       (       X       H   L     ,     A     N   -   L           ⁡     (     m   ,   n     )       )     2       +       λ   ⁡     (       X       V   L     ,     A     N   -   L           ⁡     (     m   ,   n     )       )       2     +       ψ   ⁡     (       X       D   L     ,     A     N   -   L           ⁡     (     m   ,   n     )       )       2             if         L   &lt;   N                     (   7   )               μ+λ+ψ=1  (8)
 
     The edge map of Y is defined in a similar way for X. 
       FIG. 2  shows an example diagram  200  depicting wavelet sub-bands for a two-level decomposition of the reference image X (N=2). The six sub-bands involved in computing the edge map are shown in bold outline in this figure. It is notable that the edge map is intended to be an estimate of image edges. Thus, the most informative sub-bands are used in forming the edge map, rather than considering all of the sub-bands. In our method, we use only 3N edge bands. If we considered all the bands in our edge map, we would have to use 4 N −1 bands. When N is greater than or equal to 2, the value 4 N −1 is much greater than 3N. Thus, our proposed edge map helps save computation effort. According to our simulations, considering all the image sub-bands in calculating the edge map does not have a significant impact on increasing prediction accuracy. It is notable that the edge maps only reflect the fine-edge structures of images. 
     The step  158  “Compute Edge Quality Map (or Score)” is performed in the Edge Quality Map/Score Computation Module  114 . The IQM computation is applied between the edge maps X E  and Y E . The resulting quality map (or score) is called the edge quality map (or score), IQM E . 
     Some metrics, like AD, also generate an intermediate quality map which should be pooled to reach the final score. In the step  160  “Conditionally Generate Contrast Map”, when we have a map-based metric, a contrast map function is formed for pooling the approximation and edge quality maps. It is well known that the HVS is more sensitive to areas near the edges, as shown in Z. Wang and A. C. Bovik,  Modern Image Quality Assessment . USA: Morgan &amp; Claypool, 2006. Therefore, the pixels in the quality map near the edges should be given more importance. At the same time, high-energy (or high-variance) image regions are likely to contain more information to attract the HVS, as described in Z. Wang and X. Shang, “Spatial Pooling Strategies for Perceptual Image Quality Assessment,” in  Proc. IEEE Int. Conf. Image Process. , Atlanta, Ga., October 2006, pp. 2945-2948. Thus, the pixels of a quality map in high-energy regions must also receive higher weights (more importance). Based on these facts, we can combine our edge map with the computed variance to form a contrast map function. The contrast map is computed in the Contrast Map Generation Module  110  within a local Gaussian square window, which moves (pixel by pixel) over the entire edge maps X E  and Y E . As described in Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality assessment: From error visibility to structural similarity,”  IEEE Trans. Image Process. , vol. 13, no. 4, pp. 600-612, April 2004, we define a Gaussian sliding window W={w k |k=1, 2, . . . , K} with a standard deviation of 1.5 samples, normalized to unit sum. Here, we set the number of coefficients K to 16, that is, a 4×4window. This window size is not too large and can provide accurate local statistics. The contrast map is defined as follows: 
     Contrast 
     
       
         
           
             
               
                 
                   
                     ( 
                     
                       
                         x 
                         E 
                       
                       , 
                       
                         x 
                         
                           A 
                           N 
                         
                       
                     
                     ) 
                   
                   = 
                   
                     
                       ( 
                       
                         
                           μ 
                           
                             x 
                             E 
                           
                           2 
                         
                         ⁢ 
                         
                           σ 
                           
                             x 
                             
                               A 
                               N 
                             
                           
                           2 
                         
                       
                       ) 
                     
                     0.15 
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
                     σ     x     A   N       2     =       ∑     k   =   1     K     ⁢         w   k     (       x       A   N     ,   k       -     μ     x     A   N           )     2               (   10   )                   μ     x   E       =       ∑     k   =   1     K     ⁢       w   k     ⁢     x     E   ,   k             ,       μ     x     A   N         =       ∑     k   =   1     K     ⁢       w   k     ⁢     x       A   N     ,   k                     (   11   )               
where x E  and x A     N    denote image patches of X E  and X A     N    the sliding window. It is notable that the contrast map merely exploits the original image statistics to form the weighted function for quality map pooling.
 
       FIG. 3A  shows a typical original image, and  FIG. 3B  demonstrates the resized contrast map obtained by eq. (9) for the original image of  FIG. 3A . The sample values of the contrast map are scaled between [0, 255] for easy observation. As can be seen in  FIGS. 3A and 3B , the contrast map nicely shows the edges and important image structures to the HVS. Brighter (higher) sample values in the contrast map indicate image structures that are more important to the HVS and play an important role in judging image quality. 
     For map-based metrics only, the contrast map in (9) is used in the step  170  “Conditionally Perform weighted Pooling” for weighted pooling of the approximation quality map IQM A  and the edge quality map IQM E . 
                     S   A     =         ∑     j   =   1     M     ⁢       Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )       ·       IQM   A     ⁡     (       x       A   N     ,   j       ,     y       A   N     ,   j         )               ∑     j   =   1     M     ⁢     Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )                   (   12   )                 S   E     =         ∑     j   =   1     M     ⁢       Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )       ·       IQM   E     ⁡     (       x     E   ,   j       ,     y     E   ,   j         )               ∑     j   =   1     M     ⁢     Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )                   (   13   )               
where x E,j  and x A     N     ,j  in the contrast map function denote image patches in the j-th local window; x A     N     ,j , y A     N     ,j , and y E,j  in the quality map (or score) terms are image patches (or pixels) in the j-th local window position; M is the number of samples (pixels) in the quality map; and S A  and S E  represent the approximation and edge quality scores respectively.
 
     It is noted that, for non map-based metrics like PSNR, the step  170  “Conditionally Perform weighted Pooling” is simply a copy function, so that S A =IQM A  and S E =IQM E . In the step  180  “Combine Approximation and Edge Scores”, which is executed in the Combination Module  120 , the approximation and edge quality scores are combined linearly, as defined in eq. (14), to obtain the overall quality score IQM DWT  between images X and Y: 
                         IQM   DWT     ⁡     (     X   ,   Y     )       =       β   ·     S   A       +       (     1   -   β     )     ·     S   E           ⁢     
     ⁢     0   &lt;   β   ≤   1             (   14   )               
where IQM DWT  gives the final quality score between the images; and β is a constant. As the approximation sub-band contains the main image contents, β should be close to 1 to give the approximation quality score (S A ) much greater importance. We set β to 0.85 in our simulations, which means the approximation quality score constitutes 85% of the final quality score and only 15% is made up of the edge quality score.
 
     Examples of Framework Applications 
     In this section, the application of the framework to different quality assessment methods is described. The PSNR approach is discussed in the error-based (bottom-up) category, where the AD metric is an example of a map-based metric, and PSNR is an example of a non map-based metric. 
     PSNR 
     The conventional PSNR and the mean square error (MSE) are defined as in eqs. (32) and (33): 
                     PSNR   ⁡     (     X   ,   Y     )       =     10   ·       log   10     ⁡     (       X   max   2       MSE   ⁡     (     X   ,   Y     )         )                 (   32   )                 MSE   ⁡     (     X   ,   Y     )       =       1     N   P       ·       ∑     m   ,   n       ⁢       (       X   ⁡     (     m   ,   n     )       -     Y   ⁡     (     m   ,   n     )         )     2                 (   33   )               
where X and Y denote the reference and distorted images respectively; X max  is the maximum possible pixel value of the reference image X (the minimum pixel value is assumed to be zero); and N P  is the number of pixels in each image. Although the PSNR is still popular because of its ability to easily compute quality in decibels (dB), it cannot adequately reflect the human perception of image fidelity. Other error-based techniques, such as wSNR[6] and NQM both proposed in N. Damera-Venkata, T. D. Kite, W. S. Geisler, B. L. Evans, and A. C. Bovik, “Image quality assessment based on a degradation model,”  IEEE Trans. Image Process. , vol. 9, no. 4, pp. 636-650, April 2000, and VSNR proposed in D. M. Chandler and S. S. Hemami, “VSNR: A wavelet-based visual signal-to-noise ratio for natural images,”  IEEE Trans. Image Process. , vol. 16, no. 9, pp. 2284-2298, September 2007, are more complex to use, as they follow sophisticated procedures to compute the human visual system (HVS) parameters. In this subsection, we explain how to calculate PSNR-based quality accurately in the discrete wavelet domain using the proposed framework. The first step is to determine the right number of decomposition levels (N) required to calculate the PSNR DWT  value. This number can be calculated using eq. (5). Eq. (5) is used to compute the appropriate number of decomposition levels for the IVC database. For that database, k is equal to 6. Thus,
 
 N   IVC =max(0, round(log 2 (512/57.33)))=3   (34)
 
     Please note that the LCC has low sensitivity to small variations in β, that is, the proposed β=0.85 does not drastically affect PSNR DWT  performance compared with the optimum β value for the quality prediction across different image databases. 
     In the second step, the edge-map functions of images X and Y are computed by eq. (6). Then the approximation quality score PSNR A  and the edge quality score PSNR E  are calculated using eq. (32), as defined in eq. (35) and eq. (36):
 
PSNR A =PSNR( X   A     N     , Y   A     N   )   (35)
 
PSNR E =PSNR( X   E   , Y   E )   (36)
 
     Finally, the overall quality score PSNR DWT  is computed by combining approximation and edge quality scores according to eq. (14):
 
PSNR DWT ( X, Y )=β·PSNR A +(1−β)·PSNR E , 0&lt;β≦1   (37)
 
     where PSNR DWT  gives the final quality score, or measure of quality, of the images in dB. 
     Absolute Difference (AD) 
     To verify the performance of the proposed framework of  FIG. 1  more generally, it is appropriate to investigate how it works if the AD of the images is considered as the IQM. As in previous cases, it is necessary to know the required number of decomposition levels in order to calculate the AD DWT  value. 
     In the second step of this map-based method, an approximation AD map, AD A  is calculated, between the approximation sub-bands of X and Y.
 
AD A ( m, n )=| X   A     N   ( m, n )− Y   A     N   ( m, n )|  (38)
 
     where (m,n) shows a sample position in the approximation subband. 
     In the third step, the edge-map function images X and Y are defined in eqs. (6), (7), and the edge AD map, AD E , is calculated between the edge maps X E  and Y E  in the next step.
 
AD E ( m,n )=| X   E ( m, n )− Y   E ( m, n )|  (39)
 
     In the fourth step, the contrast map is obtained using eq. (9), and then AD A  and AD E  are pooled using the contrast map to calculate the approximation and edge quality scores S A  and S E . 
     
       
         
           
             
               
                 
                   
                     S 
                     A 
                   
                   = 
                   
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           1 
                         
                         M 
                       
                       ⁢ 
                       
                         
                           Contrast 
                           ⁡ 
                           
                             ( 
                             
                               m 
                               , 
                               n 
                             
                             ) 
                           
                         
                         · 
                         
                           
                             AD 
                             A 
                           
                           ⁡ 
                           
                             ( 
                             
                               m 
                               , 
                               n 
                             
                             ) 
                           
                         
                       
                     
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           1 
                         
                         M 
                       
                       ⁢ 
                       
                         Contrast 
                         ⁡ 
                         
                           ( 
                           
                             m 
                             , 
                             n 
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   40 
                   ) 
                 
               
             
             
               
                 
                   
                     S 
                     E 
                   
                   = 
                   
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           1 
                         
                         M 
                       
                       ⁢ 
                       
                         
                           Contrast 
                           ⁡ 
                           
                             ( 
                             
                               m 
                               , 
                               n 
                             
                             ) 
                           
                         
                         · 
                         
                           
                             AD 
                             E 
                           
                           ⁡ 
                           
                             ( 
                             
                               m 
                               , 
                               n 
                             
                             ) 
                           
                         
                       
                     
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           1 
                         
                         M 
                       
                       ⁢ 
                       
                         Contrast 
                         ⁡ 
                         
                           ( 
                           
                             m 
                             , 
                             n 
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   41 
                   ) 
                 
               
             
           
         
       
     
     The final quality score, AD DWT , is calculated using eq. (42).
 
AD DWT ( X, Y )=β· S   A   30  (1−β)· S   E    (42)
 
     Detailed Analysis of Problem to be Solved 
     Let&#39;s consider an image segmented into 4×4blocks. Each block is subject to 4×4DCT transform and alteration of its frequency coefficients. These alterations can include quantization, suppression (set to zero), etc. Such operation is commonplace in H.264compression. 
       FIG. 4  shows a Matlab code fragment  900  for generating and displaying the 4×4DCT basis used in H.264compression. 
       FIG. 5  shows a number array  1000  illustrating results of the 4×4DCT basis from  FIG. 4 , calculated at an arbitrary scale, and showing numeric patterns in each of the 4×4heavily outlined squares. The square at the top-left position (reference numeral  1002 ) corresponds to the DC basis, representing a constant pixel level.  FIG. 5  shows various combinations of horizontal and vertical frequency patterns comprising 16 values in each 4×4square. 
     Each pattern square comprises 16 numerical values which, by definition, average to zero, except the top left pattern square  1002  corresponding to the DC component. Therefore the proposed metrics would not be able to detect any distortions by summing any of these patterns when using a Haar transform with a decomposition level of 2, which is operating on 4×4blocks as the approximation band for a 2×2 Haar transform with level of 2 becomes the average over 4×4blocks. Furthermore, computing the average of each of these patterns on each of their 2×2 sub-blocks, we obtain, using the same Matlab code, a mask showing which coefficients have zero average for each of their 2×2 sub-blocks, shown with a “1” in the following array denoted “sym”. 
     
       
         
           
             sym 
             = 
             
               
                 
                   0 
                 
                 
                   0 
                 
                 
                   1 
                 
                 
                   0 
                 
               
               
                 
                   0 
                 
                 
                   0 
                 
                 
                   1 
                 
                 
                   0 
                 
               
               
                 
                   1 
                 
                 
                   1 
                 
                 
                   1 
                 
                 
                   1 
                 
               
               
                 
                   0 
                 
                 
                   0 
                 
                 
                   1 
                 
                 
                   0 
                 
               
             
           
         
       
     
     This may also be determined directly in  FIG. 5  which shows that individual 2×2 sub-blocks in the third column (reference numeral  1004 ) as well in the third row (reference numeral  1006 ) sum to zero. 
     It can be seen that a metric, based on using a 2×2 Haar transform, would be blind to any modification affecting frequency coefficients in the 3rd row or in the 3rd column. 
     Besides, the approximation band of the 2×2 Haar transform being blind to such image patterns (i.e. the pixel patterns associated with the DCT basis frequency coefficients in the 3rd row or in the 3rd column), other block operations belonging to a transform may be blind to certain image patterns (in the previous example, in the approximation band, the block operation computing the average of the pixels in the block). For a block operation T, the set of block patterns of the image that cannot be recognized, i.e. for which T has a blindness, is defined as the set of blocks B for which T(B) is zero, or substantially zero. We denote this set “Null(T)” as null space due to its similarity to the concept of null space when dealing with vector spaces. Indeed, a null space is a subspace of a vector space having vectors that under a given linear transformation are mapped onto zero, or substantially zero. Therefore, we define: Null(T)={B|T(B)=0}. 
     Let C be a block composed of B and its horizontal, vertical and diagonal repetitions. Mathematically, we have: 
     
       
         
           
             C 
             = 
             
               [ 
               
                 
                   
                     B 
                   
                   
                     B 
                   
                 
                 
                   
                     B 
                   
                   
                     B 
                   
                 
               
               ] 
             
           
         
       
     
     For instance, if: 
     
       
         
           
             B 
             = 
             
               [ 
               
                 
                   
                     a 
                   
                   
                     b 
                   
                 
                 
                   
                     c 
                   
                   
                     d 
                   
                 
               
               ] 
             
           
         
       
     
     Then 
     
       
         
           
             C 
             = 
             
               [ 
               
                 
                   
                     a 
                   
                   
                     b 
                   
                   
                     
                         
                     
                   
                   
                     a 
                   
                   
                     b 
                   
                 
                 
                   
                     c 
                   
                   
                     d 
                   
                   
                     
                         
                     
                   
                   
                     c 
                   
                   
                     d 
                   
                 
                 
                   
                     
                         
                     
                   
                   
                     
                         
                     
                   
                   
                     
                         
                     
                   
                   
                     
                         
                     
                   
                   
                     
                         
                     
                   
                 
                 
                   
                     a 
                   
                   
                     b 
                   
                   
                     
                         
                     
                   
                   
                     a 
                   
                   
                     b 
                   
                 
                 
                   
                     c 
                   
                   
                     d 
                   
                   
                     
                         
                     
                   
                   
                     c 
                   
                   
                     d 
                   
                 
               
               ] 
             
           
         
       
     
     Let S=(sx,sy) be a spatial shift and B′ be a block extracted from C starting at location S and with same dimension as B (i.e. both of size W×W). We have: 
     B′(i,j)=C(i+sx,j+sy), for i,j both between 0 and W-1. 
     where both sx and sy are between 1 and W-1. 
     In the previous example, with S=(1,1), we have: 
     
       
         
           
             
               B 
               ′ 
             
             = 
             
               [ 
               
                 
                   
                     d 
                   
                   
                     c 
                   
                 
                 
                   
                     b 
                   
                   
                     a 
                   
                 
               
               ] 
             
           
         
       
     
     We denote this spatial shift operation over the repetitions of B as R(B,S) where B is the block and S the spatial shift. We have B′=R(B,S). 
     What we want is to select a spatial shift S for which as many elements of Null(T) as possible become recognized under the spatial shift (i.e. the block transform applied to spatially shifted images is not blind anymore). Ideally, we would like that the spatial shift S be selected such that T(R(B,S)) is not zero for each B element of Null(T) (except when B is zero everywhere). 
     In the 4×4DCT basis example above, if we apply a spatial shift of S=(1,1), then the approximation band of the 2×2 Haar transform will not be blind to the 3rd row or 3rd column of frequency coefficients. For other block operations, spatial shift values can be determined which would reduce and often minimize the blindness of the transformation process. 
     Generally, the transformation process has a blindness property of not containing all information regarding an image to be transformed in a transform of the image, thus resulting in certain patterns of the image not to be recognized during the transformation process, and the spatial shift being chosen so that the above mentioned certain patterns become detectable during the transformation process applied to the spatially shifted versions X 2  and Y 2 . 
     Spatial Shifting 
     According to embodiments of the invention, a spatial shift horizontally and vertically is introduced in the block alignment of the W×W pixel blocks in images for computing Haar-based visual quality metrics, in order to reduce “blindness”, thereby increasing the accuracy and robustness of the visual quality metrics calculation. 
     Spatial shifting may also be used beneficially for deriving robust quality metrics which are based on using averages of square pixel blocks. 
       FIG. 6  illustrates a first quality metrics computation system  1100 , including two instances ( 1102 . 1  and  1102 . 2 ) of a Quality Assessment module  1102 , a Spatial Shift module  1104 , and a Final Score Computation module  1106 . 
       FIG. 7  illustrates a second quality metrics computation system  1200 , which is a different arrangement of the same modules used in the first quality metrics computation system  1100 . The second quality metrics computation system  1200  includes two instances ( 1102 . 1  and  1102 . 2 ) of the Quality Assessment module  1102 , two instances ( 1104 . 1  and  1104 . 2 ) of the Spatial Shift module  1104 , and the Final Score Computation module  1106 . 
       FIG. 8  illustrates a third quality metrics computation system  1300 , which is yet another different arrangement of some of the same modules used in the first quality metrics computation system  1100 . The third quality metrics computation system  1300  includes the Quality Assessment module  1102  and the Spatial Shift module  1104 . 
     All three quality metrics computation systems  1100 ,  1200 , and  1300  have in common that a final quality score S f  is computed by processing of the two images, the original image X and the distorted image Y. The Spatial Shift module  1104  has the general ability to shift the pixels of both the original image X and the distorted image Y by predefined horizontal and vertical amounts, or number of pixels, relative to the image frame. The Quality Assessment module  1102  has the ability to compare the original image X with the distorted image Y, or shifted versions of the 2 images, and generate a Quality Metric Q. The Final Score Computation module  1106  is used to combine the Quality Metrics from two Quality Assessment modules  1102 . 
     The effect of spatial shifting is to cause different pixels from the image to be selected for the transformation blocks. By combining the Quality Metrics computed with differently shifted source images, the aforementioned “blindness” of Haar-based metrics may be overcome. 
     According to one embodiment of the invention, spatial shifts are simply applied horizontally and vertically such that the metric is not aligned with the 4×4block alignments used in H.264as illustrated in the third quality metrics computation system  1300  of  FIG. 8 . 
     In another embodiment, designed for use in cases where it is not known where the 4×4block alignment is located, the metric is computed at two locations separated with a spatial shift, which is preferably an odd number of pixels both horizontally and vertically as shown in  FIGS. 6 and 7 . In  FIG. 6 , one of the two locations is at the origin (no spatial shift corresponding to that location) while both spatial shifts may be arbitrarily selected in  FIG. 7 . 
     The preferred embodiment is illustrated in  FIG. 6 , for W×W blocks, the spatial shift is by first and second predefined amounts (numbers of pixels), both selected to be the larger of 1 and (W/2-1), horizontally and vertically respectively. This corresponds to a near quincunx pattern for the quality metric application. In the first quality assessment, blocks are shifted by 0. In the second quality assessment, the blocks are shifted by (x, y). Quincunx is a pattern for which the pixel content of the blocks of the shifted version by (x, y) is maximally different from the blocks shifted by 0. If block size is W, this maximum differentiation would happen when (x, y)=(W/2, W/2). We describe the spatially shifted pattern as near-quincunx because we want to have this maximum differentiation property by having (x, y) near (W/2, W/2) but also want to have x and y odd to avoid symmetries which could be causing problems with the Haar transform. 
     The results of the two metrics may be either averaged or the maximum distortion value of each individual application of the metric may be selected, to provide the final result. 
     Simulation results have shown that these alternatives increase the metric robustness without affecting overall accuracy, see experimental results below. 
     The quality metrics computation system  1300  of  FIG. 8  above has the same computational complexity as the original metric. 
     The quality metrics computation systems  1100  and  1200  of  FIG. 6  and  FIG. 7  respectively possess double the computational complexity compared to the quality metrics computation system  1300 , but its complexity is still much lower compared to that of other accurate visual quality metrics such as SSIM and VIF. 
     The quality metrics computation system  1100  computes, using the quality assessment module  1102 . 1 , a quality score Q 1  between the original image X and the distorted image Y. In addition, it computes, using the quality assessment module  1102 . 2 , a quality score Q 2  between shifted images X 2  and Y 2 , which are spatially shifted versions of the original image X and the distorted image Y respectively, generated by the spatial shift module  1104 . Formally:
 
 X   2 ( m,n )= X ( m+s   2,x   ,n+s   2,y )
 
 Y   2 ( m,n )= Y ( m+s   2,x   ,n+s   2,y )
 
     where (m,n) is the pixel position in the image, s 2,x  is the horizontal shift and s 2,y  is the vertical shift of the spatial shift module  1104 .
 
 Q   i =VQM( X, Y )
 
 Q   2 =VQM( X   2   , Y   2 )
 
     where VQM(X, Y) is a visual quality metric function performing quality assessment, in the quality assessment module  1102 . 1 , from images X and Y. VQM(X 2 , Y 2 ) is a visual quality metric function performing quality assessment, in quality assessment module, from the shifted images X 2  and Y 2 . VQM may be any visual quality metric such as the SSIM, PSNR A , etc. 
     The final quality score, or the measure of quality, S f  is computed in the final score computation module  1106  based on Q 1  and Q 2 . The final score can be computed using a weighted sum of Q 1  and Q 2 , taking the minimum or other appropriate operations depending on the context. 
     The Weighted sum may be computed as follows:
 
 S   f   =a ( Q   1 )+(1 −a )( Q   2 ) with 0 &lt;=a &lt;=1.
 
     As a special case, one would obtain the mean by setting a=0.5. 
     The weighted sum is useful in assigning different levels of importance to each quality score. 
     The Minimum may be computed as follows:
 
 S   f =Min( Q   1   , Q   2 ).
 
     The minimum is useful to take the worst case of the two quality scores. 
     More complex operations: 
     Sometimes, the final score is advantageously a more complex function of the quality scores. For instance, for PSNR A =10 log 10 ((2 L (255)) 2 /mse), the quality scores to combine would be the mse. And one would have:
         Q 1 =mse 1      Q 2 =mse 2          

     Then several options are possible:
         mse=Max(Q 1 , Q 2 )   mse=a(Q 1 )+(1−a)(Q 2 ) with 0&lt;=a&lt;=1
           such as mse=(0.5)(Q 1 )+(0.5)(Q 2 ) when a=0.5   
               

     S f =10 log 10 ((2 L (255)) 2 /mse), where L is the number of decomposition levels in PSNR A  and mse represents the mean squared error and is defined as: 
             mse   =       1   MN     ⁢       ∑       i   =   m     ,     j   =   n         M   ,   N       ⁢       (     〚       X   ⁡     (     m   ,   n     )       -     Y   ⁡     (     m   ,   n     )         )     〛     2               
where M and N are the image dimensions.
 
     Please note that specific values of spatial shifts s 2,x , and S 2,y  can be beneficial to make a metric robust to distortions. For instance, the 4×4discrete cosine transform (DCT) is used in H.264video compression. Often, the high frequency coefficients of a DCT-transforrned block of pixels are heavily quantized to increase compression in H.264. But this introduces distortions. Unfortunately, as seen in the previous section, the PSNR A  metric, when using a 2×2 Haar transform with a decomposition level of 2, which is operating on 4×4blocks aligned with the 4×4DCTs of H.264 compression (as is the case by default), would not be able to detect these distortions since they do not alter the average pixel value within each 4×4block. 
     A solution to the problem is to use an odd spatial shift horizontally and vertically so that not only the 2×2 Haar transform with a decomposition level of 2 would not be aligned with the H.264 4×4 DCT blocks but also that this misalignment would then not permit 4×4DCT coefficients having an even symmetry to remain undetected for any number of levels of the Haar transform (i.e. any average of W×W blocks with W being a power of 2). 
     Therefore, for PSNR A , it is suggested to set the shift values to max(1, W/2-1) for block sizes of W for the Haar transform. Computing the quality metric twice but with the second one being applied from a well-chosen shifted position, as described here, will avoid the situation in which both metrics miss distortions having specific symmetry patterns. This is especially useful when the alignment of a potential distortion is unknown, for example, if the actual alignment of the 4×4DCT transforms is not known, for instance in the case where the image was encoded and decoded with H.264and then cropped. 
     The quality metrics computation system  1200  computes, again using the quality assessment module  1102 . 1 , a quality score Q 1s  between shifted images X 1  and Y 1 , which are spatially shifted versions, generated by the spatial shift module  1104 . 1 , of the original image X and the distorted image Y respectively. In addition, the quality metrics computation system  1200  computes, using the quality assessment module  1102 . 2 , the quality score Q 2s  between the shifted images X 2  and Y 2 , which are spatially shifted versions, generated by the spatial shift module  1102 . 2  of the original image X and the distorted image Y respectively. Formally:
 
 X   1 ( m,n )= X ( m+s   1,x   ,n+s   1,y )
 
 Y   1 ( m,n )= Y ( m+s   1,x   ,n+s   1,y )
 
 X   2 ( m,n )= X ( m+s   2,x   ,n+s   2,y )
 
 Y   2 ( m,n )= Y ( m+s   2,x ,n+s 2,y )
 
     where (m,n) is the pixel position in the image, s 1,x  is the horizontal shift and s 1,y  is the vertical shift of the spatial shift module  1102 . 1 , and s 2,x  is the horizontal shift and s 2,y  is the vertical shift of the spatial shift module  1102 . 2 .
 
 Q   1s =VQM( X   1   , Y   1 )
 
 Q   2s =VQM( X   2   , Y   2 )
 
     where VQM(X, Y) is as described previously. 
     The final quality score S f  is computed in the final score computation module  1106  based on Q 1s  and Q 2s . As described previously. 
     The same comment made with regard to  FIG. 6  regarding the spatial shift applies here as well. It may further noted that specific values of spatial shifts s 1,x , s 1,y , s 2,x , and s 2,y  can be very beneficial to make the metric robust to distortions. For PSNR A , it is suggested that a difference between the shift values s 1,x −s 2,x , be an odd number, and the difference between the shift values s 1,y −s 2,y  also be an odd number. It is further recommended that s 1,x −s 2,x =s 1,y −s 2,y =max(1, W/2-1) for block sizes of W for the Haar transform. 
     The quality metrics computation system  1300  computes, using the quality assessment module  1102 , the quality score Q 1s  between image X 1  and Y 1 , which are spatially shifted versions, generated by the spatial shift module  1104 , of the original image X and the distorted image Y respectively. 
     Formally:
 
 X   1 ( m,n )= X ( m+s   1,x   ,n+s   1,y )
 
 Y   1 ( m,n )= Y ( m+s   1,x   ,n+s   1,y )
 
     where (m,n) is the pixel position in the image, s 1,x  is the horizontal shift value and s 1,y  is the vertical shift value used in the spatial shift module  1104 .
 
 Q   1s =VQM( X, Y )
 
     where VQM(X, Y) is as described previously. 
     The final quality score, or measure of quality, S f =Q 1s . 
     The same comment made with regard to  FIG. 6  regarding the spatial shift applies here as well. Here, the metric is computed only once but a spatial shift is applied in order to avoid alignment with distortions. This may be especially useful when the alignment of a distortion is known, for example, if the alignment of the 4×4DCT transforms is known to be on even locations, it is advantageous to apply an odd shift. For PSNR A , shift value equal to 1 is suggested. 
     Additional embodiments of the invention may be based on an analogous use of more combinations of locations, generated by spatial shift of the original and distorted images X and Y, for the computation of PSNR A . The same method, involving images and shifted versions of the images in combination, may also be applied to other wavelet based metrics besides Haar, and to other Haar transform sizes. Using more combinations of locations, i.e. shifted positions may, depending on each situation, increase robustness at the expense of increased computational complexity. A person skilled in the art of digital signal processing will know how to adapt the invention presented herein to other Haar transform sizes, and wavelets depending on the application of interest. 
       FIG. 9  shows a block diagram of an embodiment of a discrete wavelet domain image quality assessment system  1400 , as a preferred implementation of the first quality metrics computation system  1100  of  FIG. 6 . The discrete wavelet domain image quality assessment system  1400  comprises: first and second IQM computation subsystems  100 . 1  and  100 . 2 , each corresponding to an instance of the IQM computation subsystem  100  of  FIG. 1 , the Spatial Shift module  1104  of  FIG. 6 , an Common Pooling of Approximation Quality Map Module  1402 , an Common Pooling of Edge Quality Map Module  1404 , and the Final Score Computation module  1106  of  FIG. 6  which is equivalent to the Combination Module  120  of  FIG. 1 . 
     The inputs to the discrete wavelet domain image quality assessment system  1400  are the Reference (original) Image X and the Distorted Image Y, and the final output of the system  1400  is the final quality score S f  as in  FIG. 6 . 
     The Spatial Shift module  1104  includes a Shift X module  1406  for generating the spatially shifted version X 2  of the Original image X, and a Shift Y module  1408  for generating the spatially shifted version Y 2  of the distorted image Y. 
     The first IQM computation subsystems  100 . 1  computes a first approximation quality map (or score) IQM A , a first Edge Quality Map (Score) IQM E , and optionally a first contrast map labeled “Contrast”, based on a comparison of the images X and Y. Similarly, the second IQM computation subsystems  100 . 2  computes a second approximation quality map (or score) IQM′ A , a second Edge Quality Map (Score) IQM′ E , and optionally a second contrast map labeled “Contrast&#39; ”, based on the spatially shifted versions X 2  and Y 2  of the images X and Y respectively. The first and second approximation quality maps (scores) IQM A  and IQM′ A  for the approximation sub-bands are combined in the Common Pooling of Approximation Quality Map module  1402 . The first and second Edge Quality Maps (scores) IQM E  and IQM′ E  for the detail sub-bands are combined in the Common Pooling of Edge Quality Map module  1404 . 
     Both, the Common Pooling of Approximation Quality Map module  1402  and the Common Pooling of Edge Quality Map module  1404  also optionally receive the contrast maps “Contrast” and “Contrast&#39; ”, and output the approximation quality score S A  and the edge quality score S E  respectively to the Final Score Computation module  1106  which combines them into the Final Score S f , also referred to as a final measure of quality. 
     The inputs to the second IQM computation subsystem  100 . 2 , after the spatial shift operation performed in the Spatial Shift module  1104 , are:
 
 X   2 ( m,n )= X ( m+s   x   ,n+s   y )
 
 Y   2 ( m,n )= Y ( m+s   x ,n+s y )
 
       FIG. 10  shows a flowchart of a method  1500  of computing a combined score which is a measure of quality of a distorted image by comparison with a corresponding original image, combined with a comparison of a spatially shifted version of the distorted image with an equally spatially shifted version of the corresponding original image, in a sequence chart associated with the system block diagram of  FIG. 9 . 
     The method  1500  comprises steps: 
       1502  “Perform five steps on original image and distorted image”; 
       1504  “Apply spatial shift to original image”; 
       1506  “Apply spatial shift to distorted image”; 
       1508  “ Perform five steps on spatially shifted images”; 
       1510  “Conditionally Perform common weighted Pooling”; and 
       1512  “Compute Final Quality Score”. 
     Except for the initial spatial shift which is performed in the Spatial Shift module  1104  to provide the shifted images X 2  and Y 2  to the second IQM computation subsystem  100 . 2 , each subsystem is the same as the IQM computation subsystem  100  of  FIG. 1  described earlier, except the final step  1512  “Compute Final Quality Score”. Accordingly, the method  1500  includes a first step group  1502  and a second step group  1508 , each of essentially the same five steps as the steps  152  to  160  of  FIG. 1A . 
     The first step group  1502 , which corresponds to the actions of the first IQM computation subsystems  100 . 1  to perform the steps  152  to  160  on the original image and the distorted image, is followed by the steps  1504  “Apply spatial shift to original image” and  1506  “Apply spatial shift to distorted image”, in which spatially shifted versions of the original and the distorted images respectively are generated. 
     The next steps are in the step group  1508 , which corresponds to the actions of the second IQM computation subsystems  100 . 2  to perform the steps  152  to  160  on the spatially shifted versions of the original image and the distorted image. 
     For clarity, the two five-step groups  1502  and  1508  are expanded to individual steps, and shown in  FIGS. 10A and 10B  respectively, indicating the arguments of each step. 
     However, the step  170  “Conditionally Perform weighted Pooling” of  FIG. 1A  is replaced in  FIG. 10  with the new step  1510  “Conditionally Perform common weighted Pooling” which is performed in the Common Pooling of Approximation Quality Map module  1402  and the Common Pooling of Edge Quality Map module  1404  as described in more detail below. 
     In the final step  1512  “Compute Final Quality Score” of combining approximation and edge scores, the approximation quality score S A  and the edge quality score S E  are computed. 
     The step  1512  is performed in the Final Score Computation module  1106 , and corresponds exactly to the step  180  of  FIG. 1A . 
     The step  1510  “Conditionally Perform common weighted Pooling” uses the Contrast and Contrast&#39; maps for weighted pooling in the case of map-based metrics but can be omitted when other types of metrics are used. 
     For map-based metrics, the contrast map of equation (9) is used for weighted pooling of the approximation quality maps IQM A  and IQM′ A  and the edge quality maps IQM E  and IQM′ E . 
                     S   A     =       1   2     [           ∑     j   =   1     M     ⁢       Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )       ·       IQM   A     ⁡     (       x       A   N     ,   j       ,     y       A   N     ,   j         )               ∑     j   =   1     M     ⁢     Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )           +         ∑     j   =   1     M     ⁢         Contrast   ′     ⁡     (       x     E   ,   j     ′     ,     x       A   N     ,   j     ′       )       ·       IQM   A   ′     ⁡     (       x       A   N     ,   j     ′     ,     y       A   N     ,   j     ′       )               ∑     j   =   1     M     ⁢       Contrast   ′     ⁡     (       x     E   ,   j     ′     ,     x       A   N     ,   j     ′       )             ]             (   44   )                 S   E     =       1   2     [           ∑     j   =   1     M     ⁢       Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )       ·       IQM   E     ⁡     (       x     E   ,   j       ,     y     E   ,   j         )               ∑     j   =   1     M     ⁢     Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )           +         ∑     j   =   1     M     ⁢         Contrast   ′     ⁡     (       x     E   ,   j     ′     ,     x       A   N     ,   j     ′       )       ·       IQM   E   ′     ⁡     (       x     E   ,   j     ′     ,     y     E   ,   j     ′       )               ∑     j   =   1     M     ⁢       Contrast   ′     (       x     E   ,   j     ′     ,     x       A   N     ,   j     ′       )           ]             (   45   )               
where x E,j  and x A     N     ,j  in the contrast map function denote image patches in the j-th local window related to image X; x′ E,j  and x′ A     N     ,j  in the contrast map function denote image patches in the j-th local window related to image X′; x A     N     ,j , y A     N     ,j , x E,j , and y E,j  in the quality map (or score) terms are image patches (or pixels) in the j-th local window position related to images X and Y; x′ A     N     ,j , y′ A     N     ,j , x′ E,j , and y′ E,j  in the quality map (or score) terms are image patches (or pixels) in the j-th local window position related to images X′ and Y′; M is the number of samples (pixels) in the quality map; and S A  and S E  represent the approximation and edge quality scores respectively.
 
     As for the previous section, we can combine the quality scores differently. For instance: 
                       V   A     =     [           ∑     j   =   1     M     ⁢       Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )       ·       IQM   A     ⁡     (       x       A   N     ,   j       ,     y       A   N     ,   j         )               ∑     j   =   1     M     ⁢     Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )           ,         ∑     j   =   1     M     ⁢         Contrast   ′     ⁡     (       x     E   ,   j     ′     ,     x       A   N     ,   j     ′       )       ·       IQM   A   ′     ⁡     (       x       A   N     ,   j     ′     ,     y       A   N     ,   j     ′       )               ∑     j   =   1     M     ⁢       Contrast   ′     ⁡     (       x     E   ,   j     ′     ,     x       A   N     ,   j     ′       )             ]       ⁢     
     ⁢     index   =         arg   ⁢           ⁢   min     i     ⁢     (       V   A     ⁡     [   i   ]       )         ⁢     
     ⁢       S   A     =       V   A     ⁡     [   index   ]                 (   46   )                   V   E     =     [           ∑     j   =   1     M     ⁢       Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )       ·       IQM   E     ⁡     (       x     E   ,   j       ,     y     E   ,   j         )               ∑     j   =   1     M     ⁢     Contrast   ⁡     (       x     E   ,   j       ,     x       A   N     ,   j         )           ,         ∑     j   =   1     M     ⁢         Contrast   ′     ⁡     (       x     E   ,   j     ′     ,     x       A   N     ,   j     ′       )       ·       IQM   E   ′     ⁡     (       x     E   ,   j     ′     ,     y     E   ,   j     ′       )               ∑     j   =   1     M     ⁢       Contrast   ′     (       x     E   ,   j     ′     ,     x       A   N     ,   j     ′       )           ]       ⁢     
     ⁢       S   E     =       V   E     ⁡     [   index   ]                 (   47   )               
where V A  is an array containing the pooled IQM A  and IQM′ A  and V E  is an array containing the pooled IQM E  and IQM′ E . In that case, we select the approximation quality map with minimal quality and its corresponding edge quality map to compute S A  and S E  respectively. Other more complex combining of the edge maps is also possible. Also the system can easily be extended to combinations of more shifted images at input (i.e. arbitrary number of shifts instead on a single one).
 
     For non map-based metrics like AD, we can compute as follows:
 
 S   A   =a (IQM A )+(1 −a )(IQM′ A ) and  S   E   =a (IQM E )+(1 −a )(IQM′ E ). with 0 &lt;=a &lt;=1.
         For instance, we can set a=0.5 to compute the average.
 
 V   A =[IQM A ,IQM′ A ], index=argmin i ( V   A   [i ]),  S   A   =V   A [index],  S   E   =V   E [index].
       

     For PSNR DWT  we compute:
 
mse A   =a (IQM A )+(1 −a )(IQM′ A ) and mse E   =a (IQM E )+(1 −a )(IQM′ E ). with 0 &lt;=a &lt;=1.
         For instance, we can set a=0.5 to compute the average.
 
 V   A =[IQM A ,IQM′ A ], index=argmax i ( V   A   [i]), mse   A   =V   A [index], mse E   =V   E [index].
 
 S   A =10 log 10 ((2 L (255)) 2 /mse A ) and  S   E =10 log 10 ((2 L (255)) 2 /mse E ).
       

     where L is the number of decomposition levels in PSNR DWT  and mse represents the mean squared error. 
     Other more complex combining of the edge maps is also possible. The system can also be extended to combinations of more shifted images (i.e. an arbitrary number of shifts instead of combining image quality metrics (IQMs) from only one shifted image pair X 2  and Y 2  and the un-shifted pair of images X and Y). Furthermore, the invention applies also even if the detailed sub-bands are ignored in the metric computation such as is the case of PSNR A . 
     The systems illustrated in FIGS.  1  and  6 - 9  of the embodiments of the present invention include a general purpose or specialized computer having a CPU and a computer readable storage medium e.g., memory, DVD, CD-ROM, floppy, magnetic tape or other storage medium, having computer readable instructions stored thereon, for execution by the CPU, forming modules of the above noted systems as described above Alternatively, the systems of FIGS.  1  and  6 - 9  can include specialized dedicated hardware or a combination of firmware and a specialized computer having a computer readable storage medium having computer readable instructions stored thereon for execution by the CPU for forming modules of these systems as described above. Each module of the systems of FIGS.  1  and  6 - 9  include a firmware or, alternatively, computer readable instructions stored in a computer readable storage medium such as memory device for execution by a processor. 
     The discrete wavelet domain image quality assessment system  1400  is realized on a computer equipped with software programs stored in the memory of the computer for enabling the features of the embodiments of the invention. Some or all of the functions of the system  1400  may also be implemented with an Application Specific Integrated Circuit (ASIC) or a number of ASICs. 
       FIG. 12  shows an exemplary hardware view  1700  of the discrete wavelet domain image quality assessment system  1400  including a Processor  1702  and a Computer Memory  1704  according to the embodiments of the invention. 
     The Processor  1702  may be any commercial processor capable of executing programs under an operating system such as, but not limited to, Linux, Microsoft Windows, or Mac Operating System  10  (OSX) for example, comprising a Central Processing Unit (CPU)  1706 , a Network Input/Output (I/O) system  1708  and a Command Interface  1710 . 
     The CPU  1706  may be any of a number of commercial processor implementations, including a multi-processor. The Network Input/Output (I/O) system  1708  provides an interface facility for receiving images for which quality assessment is to be performed. The Command Interface  1710  may include a keyboard, mouse, and visual interface, as well as removable storage devices, or any other hardware suitable for controlling a software configuration as well as the operations of discrete wavelet domain image quality assessment system  1400 . 
     The Processor  1702  is connected to the Computer Memory  1704  which is preferably a non-transitory memory device such as dynamic memory (DRAM) capable of storing software programs and related data. Software modules, which include computer executable instructions for execution by the Processor  1702  are stored in the Computer Memory  1704 , comprise: the IQM computation subsystem  100  which includes the first and second N-level DWT modules  102  and  104 , to be also referred to as Transform modules  102  and  104  respectively, the first and second Detail Sub-band Aggregation Modules  106  and  108 , the Contrast Map Generation Module  110 , the Approximation Quality Map/Score Computation Module  112 , the Edge Quality Map/Score Computation Module  114 , the Pooling of Approximation Quality Map Module  116  the Pooling of Edge Map Module  118  and the Combination Module  120 . Other software modules stored in the Computer memory  1704  include the Final Score Computation module  1106 , the Approximation Quality Map Module  1402 , the Common Pooling of Edge Quality Map Module  1404 , and the Spatial Shift module  1104  which includes the Shift X module  1406  and the Shift Y module  1408 . 
     The Computer Memory  1704  may also include a Cache  1712  for temporarily storing data such as the original and distorted images X and Y and the shifted versions X 2  and Y 2  which may be loaded from image databases through the Network I/O system  1708  as well as temporary results such as first and second contrast maps  1714  and  1716  generated by the Contrast Map Generation Module  110  for example. 
     Experimental Results 
     Table III shows accuracy results of various metric alternatives derived from PSNR A  compared to PSNR A  and PSNR on the LIVE image database. The description of each method is as follows: QPSNR A  applies PSNR A  at two spatial locations to form a near quincunx pattern. It averages the error values (mean squared error) obtained at these two spatial positions to compute the final dB score. There are two spatial positions: one at (0, 0) and one with a spatial shift of max(1, N/2-1) pixels horizontally and vertically. 
     QMPSNR A  applies PSNR A  at two spatial locations to form a near quincunx pattern. It keeps the largest error value (mean squared error) obtained at these two spatial positions to compute the final dB score. There are two spatial positions: one at (0, 0) and one with a spatial shift of max(1, N/2-1) pixels horizontally and vertically. 
     EPSNR A  applies PSNR A  at three spatial locations. It averages the error values (mean squared error) obtained at these three spatial positions to compute the final dB score. There are three spatial positions: one at (0, 0), one with a spatial shift of max(1, N/2-1) pixels horizontally only, and one with a spatial shift of max(1, N/2-1)) pixels horizontally and vertically. 
     SPSNR A  performs the PSNR A  but with a spatial shift of (1, 1) pixels. 
     
       
         
           
               
             
               
                 TABLE III 
               
             
            
               
                   
               
               
                 LCC, SRCC AND RMSE VALUES AFTER NONLINEAR 
               
               
                 REGRESSION FOR THE LIVE IMAGE DATABASE 
               
            
           
           
               
               
               
               
               
            
               
                   
                 Metric 
                 LCC 
                 SRCC 
                 RMSE 
               
               
                   
               
               
                   
                 PSNR 
                 0.8701 
                 0.8756 
                 13.4685 
               
               
                   
                 PSNR A  (1 location) 
                 0.9282 
                 0.9302 
                 10.1635 
               
               
                   
                 QPSNR A  (2 locations) 
                 0.9281 
                 0.9308 
                 10.1710 
               
               
                   
                 QMPSNR A  (2 locations) 
                 0.9279 
                 0.9299 
                 10.1883 
               
               
                   
                 EPSNR A  (3 locations) 
                 0.9287 
                 0.9314 
                 10.1331 
               
               
                   
                 SPSNR A  (1 location) 
                 0.9269 
                 0.9299 
                 10.2526 
               
               
                   
               
            
           
         
       
     
       FIGS. 11(   a ) to  11 ( d ) show four versions of a standard image “Image House” (image number 4.1.05) with a resolution of 256×256 available at http://sipi.usc.edu/database/database.php?volume=misc. 
       FIG. 11(   a ) shows the original (reference) image. 
       FIG. 11(   b ) shows a first distorted image created from the reference image with 2 nd  row and 2 nd  column of the 4×4DCT transformed image set to zero. 
       FIG. 11(   c ) shows a second distorted image created from the reference image with 3 rd  row and 3 rd  column of the 4×4DCT transformed image set to zero. 
       FIG. 11(   d ) shows a third distorted image created from the reference image with 4 th  row and 4 th  column of 4×4DCT transformed image set to zero. 
     Table IV shows the final score of various metric alternatives derived from shifting PSNR A  compared to PSNR A  and PSNR on the reference image of  FIG. 11(   a ) and the three distorted images below. The first distorted image of  FIG. 11(   b ) is obtained by segmenting the image of  FIG. 11(   a ) into 4×4 blocks, apply a 4×4 forward DCT transform, as defined in the H.264 standard, set to zero all frequency coefficients on the 2 nd  row and 2 nd  column of the forward DCT, and apply the 4×4 inverse DCT transform. The second distorted image of  FIG. 11(   c ) is obtained in a similar manner but 3 rd  row and 3 rd  column are set to zero instead. The third distorted image of  FIG. 11(   d ) is obtained in a similar manner but 4 th  row and 4 th  column are set to zero instead. We can visually observe that  FIG. 11(   b ) has worse quality than  FIG. 11(   c ) and that  FIG. 11(   c ) has worse quality than  FIG. 11(   d ). However, in Table VI, PSNR A  gives a much higher score to  FIG. 11(   c ) than  FIG. 11(   d ) which is highly unacceptable. We can see that the metrics having a spatial shift do not have this problem and rank the image qualities properly (i.e. they rank the images in the right order from best quality to worst quality). This shows that the proposed metrics have increased robustness. 
     
       
         
           
               
             
               
                 TABLE IV 
               
             
            
               
                   
               
               
                 METRIC VALUES FOR IMAGES OF FIG. 11(a) COMPARED 
               
               
                 TO FIGS. 11(b), 11(c), AND 11(d). 
               
            
           
           
               
               
               
               
            
               
                   
                 FIG. 11(a) 
                 FIG. 11(a) 
                 FIG. 11(a) 
               
               
                   
                 compared to 11(b) 
                 compared to 11(c) 
                 compared to 11(c) 
               
               
                   
                 (i.e. 2 nd  row and 
                 (i.e. 3 rd  row and 
                 (i.e. 4 th  row and 
               
               
                 Quality 
                 2 nd  column 
                 3 rd  column 
                 4 th  column 
               
               
                 metric 
                 removed) 
                 removed) 
                 removed) 
               
               
                   
               
               
                 PSNR 
                 26.04 
                 31.76 
                 36.80 
               
               
                 PSNR A   
                 26.93 
                 64.20 
                 45.42 
               
               
                 (1 location) 
                   
                   
                   
               
               
                 QPSNR A   
                 28.02 
                 36.50 
                 46.75 
               
               
                 (2 locations) 
                   
                   
                   
               
               
                 QMPSNR A   
                 26.93 
                 33.50 
                 45.42 
               
               
                 (2 locations) 
                   
                   
                   
               
               
                 EPSNR A   
                 28.16 
                 35.79 
                 47.12 
               
               
                 (3 locations) 
                   
                   
                   
               
               
                 SPSNR A   
                 29.47 
                 33.50 
                 48.67 
               
               
                 (1 location) 
               
               
                   
               
            
           
         
       
     
     In the embodiments described above, the transformation process is a Haar transform. It is understood that another transformation process may be also used, for example, a cosine transform, a sine transform, or a wavelet transform. Instead of Haar transform, alternatively, the transformation process can be a Daubechies transform. 
     Although specific embodiments of the invention have been described in detail, it should be understood that the described embodiments are illustrative and not restrictive. Various changes and modifications of the embodiments may be made within the scope of the following claims without departing from the scope of the invention in its broader aspect.