Changes done in npr module
authorsiddharth <siddharthkherada27@gmail.com>
Thu, 17 Oct 2013 15:44:31 +0000 (21:14 +0530)
committersiddharth <siddharthkherada27@gmail.com>
Wed, 4 Dec 2013 13:48:00 +0000 (19:18 +0530)
modules/photo/src/npr.cpp
modules/photo/src/npr.hpp
modules/photo/src/seamless_cloning.cpp
modules/photo/src/seamless_cloning.hpp

index 08c7753..a3b38d4 100644 (file)
@@ -79,34 +79,24 @@ void cv::detailEnhance(InputArray _src, OutputArray _dst, float sigma_s, float s
 
     int h = I.size().height;
     int w = I.size().width;
-    int channel = I.channels();
     float factor = 3.0;
 
     Mat img = Mat(I.size(),CV_32FC3);
     I.convertTo(img,CV_32FC3,1.0/255.0);
 
-    Mat res = Mat(h,w,CV_32FC3);
+    Mat res = Mat(h,w,CV_32FC1);
     dst.convertTo(res,CV_32FC3,1.0/255.0);
 
     Mat result = Mat(img.size(),CV_32FC3);
     Mat lab = Mat(img.size(),CV_32FC3);
-    Mat l_channel = Mat(img.size(),CV_32FC1);
-    Mat a_channel = Mat(img.size(),CV_32FC1);
-    Mat b_channel = Mat(img.size(),CV_32FC1);
+    vector <Mat> lab_channel;
 
     cvtColor(img,lab,COLOR_BGR2Lab);
-
-    for(int i = 0; i < h; i++)
-        for(int j = 0; j < w; j++)
-        {
-            l_channel.at<float>(i,j) = lab.at<float>(i,j*channel+0);
-            a_channel.at<float>(i,j) = lab.at<float>(i,j*channel+1);
-            b_channel.at<float>(i,j) = lab.at<float>(i,j*channel+2);
-        }
+    split(lab,lab_channel);
 
     Mat L = Mat(img.size(),CV_32FC1);
 
-    l_channel.convertTo(L,CV_32FC1,1.0/255.0);
+    lab_channel[0].convertTo(L,CV_32FC1,1.0/255.0);
 
     Domain_Filter obj;
 
@@ -114,23 +104,13 @@ void cv::detailEnhance(InputArray _src, OutputArray _dst, float sigma_s, float s
 
     Mat detail = Mat(h,w,CV_32FC1);
 
-    for(int i = 0; i < h; i++)
-        for(int j = 0; j < w; j++)
-            detail.at<float>(i,j) = L.at<float>(i,j) - res.at<float>(i,j);
-
-    for(int i = 0; i < h; i++)
-        for(int j = 0; j < w; j++)
-            L.at<float>(i,j) = res.at<float>(i,j) + factor*detail.at<float>(i,j);
+    detail = L - res;
+    multiply(detail,factor,detail);
+    L = res + detail;
 
-    L.convertTo(l_channel,CV_32FC1,255);
+    L.convertTo(lab_channel[0],CV_32FC1,255);
 
-    for(int i = 0; i < h; i++)
-        for(int j = 0; j < w; j++)
-        {
-            lab.at<float>(i,j*channel+0) = l_channel.at<float>(i,j);
-            lab.at<float>(i,j*channel+1) = a_channel.at<float>(i,j);
-            lab.at<float>(i,j*channel+2) = b_channel.at<float>(i,j);
-        }
+    merge(lab_channel,lab);
 
     cvtColor(lab,result,COLOR_Lab2BGR);
     result.convertTo(dst,CV_8UC3,255);
@@ -174,64 +154,21 @@ void cv::stylization(InputArray _src, OutputArray _dst, float sigma_s, float sig
     int channel = img.channels();
 
     Mat res = Mat(h,w,CV_32FC3);
+    Mat magnitude = Mat(h,w,CV_32FC1);
 
     Domain_Filter obj;
     obj.filter(img, res, sigma_s, sigma_r, NORMCONV_FILTER);
 
-    vector <Mat> planes;
-    split(res, planes);
-
-    Mat magXR = Mat(h, w, CV_32FC1);
-    Mat magYR = Mat(h, w, CV_32FC1);
-
-    Mat magXG = Mat(h, w, CV_32FC1);
-    Mat magYG = Mat(h, w, CV_32FC1);
-
-    Mat magXB = Mat(h, w, CV_32FC1);
-    Mat magYB = Mat(h, w, CV_32FC1);
-
-    Sobel(planes[0], magXR, CV_32FC1, 1, 0, 3);
-    Sobel(planes[0], magYR, CV_32FC1, 0, 1, 3);
-
-    Sobel(planes[1], magXG, CV_32FC1, 1, 0, 3);
-    Sobel(planes[1], magYG, CV_32FC1, 0, 1, 3);
-
-    Sobel(planes[2], magXB, CV_32FC1, 1, 0, 3);
-    Sobel(planes[2], magYB, CV_32FC1, 0, 1, 3);
-
-    Mat magx = Mat(h,w,CV_32FC1);
-    Mat magy = Mat(h,w,CV_32FC1);
-
-    Mat mag1 = Mat(h,w,CV_32FC1);
-    Mat mag2 = Mat(h,w,CV_32FC1);
-    Mat mag3 = Mat(h,w,CV_32FC1);
-
-    magnitude(magXR,magYR,mag1);
-    magnitude(magXG,magYG,mag2);
-    magnitude(magXB,magYB,mag3);
-
-    Mat magnitude = Mat(h,w,CV_32FC1);
-
-    for(int i =0;i < h;i++)
-        for(int j=0;j<w;j++)
-        {
-            magnitude.at<float>(i,j) = mag1.at<float>(i,j) + mag2.at<float>(i,j) + mag3.at<float>(i,j);
-        }
-
-    for(int i =0;i < h;i++)
-        for(int j=0;j<w;j++)
-        {
-            magnitude.at<float>(i,j) = 1.0f -  magnitude.at<float>(i,j);
-        }
+    obj.find_magnitude(res,magnitude,2);
 
     Mat stylized = Mat(h,w,CV_32FC3);
 
-    for(int i =0;i < h;i++)
-        for(int j=0;j<w;j++)
-            for(int c=0;c<channel;c++)
-            {
-                stylized.at<float>(i,j*channel + c) = res.at<float>(i,j*channel + c) * magnitude.at<float>(i,j);
-            }
+    vector <Mat> temp;
+    split(res,temp);
+    multiply(temp[0],magnitude,temp[0]);
+    multiply(temp[1],magnitude,temp[1]);
+    multiply(temp[2],magnitude,temp[2]);
+    merge(temp,stylized);
 
     stylized.convertTo(dst,CV_8UC3,255);
 }
@@ -259,7 +196,7 @@ void cv::edgeEnhance(InputArray _src, OutputArray _dst, float sigma_s, float sig
 
     obj.filter(img, res, sigma_s, sigma_r, NORMCONV_FILTER);
 
-    obj.find_magnitude(res,magnitude);
+    obj.find_magnitude(res,magnitude,1);
 
     magnitude.convertTo(dst,CV_8UC1,255);
 }
index 3f5519d..733d418 100644 (file)
@@ -61,7 +61,7 @@ class Domain_Filter
         void getGradienty( const Mat &img, Mat &gy);
         void diffx(const Mat &img, Mat &temp);
         void diffy(const Mat &img, Mat &temp);
-        void find_magnitude(Mat &img, Mat &mag);
+        void find_magnitude(Mat &img, Mat &mag, int flags);
         void compute_boxfilter(Mat &output, Mat &hz, Mat &psketch, float radius);
         void compute_Rfilter(Mat &O, Mat &horiz, float sigma_h);
         void compute_NCfilter(Mat &O, Mat &horiz, Mat &psketch, float radius);
@@ -131,9 +131,8 @@ void Domain_Filter::getGradienty( const Mat &img, Mat &gy)
             }
 }
 
-void Domain_Filter::find_magnitude(Mat &img, Mat &mag)
+void Domain_Filter::find_magnitude(Mat &img, Mat &mag, int flags)
 {
-
     int h = img.rows;
     int w = img.cols;
 
@@ -149,17 +148,28 @@ void Domain_Filter::find_magnitude(Mat &img, Mat &mag)
     Mat magXB = Mat(h, w, CV_32FC1);
     Mat magYB = Mat(h, w, CV_32FC1);
 
-    getGradientx(planes[0], magXR);
-    getGradienty(planes[0], magYR);
+    if(flags == 1)
+    {
+        getGradientx(planes[0], magXR);
+        getGradienty(planes[0], magYR);
+
+        getGradientx(planes[1], magXG);
+        getGradienty(planes[1], magYG);
 
-    getGradientx(planes[1], magXG);
-    getGradienty(planes[1], magYG);
+        getGradientx(planes[2], magXR);
+        getGradienty(planes[2], magYR);
+    }
+    else if(flags == 2)
+    {
+        Sobel(planes[0], magXR, CV_32FC1, 1, 0, 3);
+        Sobel(planes[0], magYR, CV_32FC1, 0, 1, 3);
 
-    getGradientx(planes[2], magXR);
-    getGradienty(planes[2], magYR);
+        Sobel(planes[1], magXG, CV_32FC1, 1, 0, 3);
+        Sobel(planes[1], magYG, CV_32FC1, 0, 1, 3);
 
-    Mat magx = Mat(h,w,CV_32FC1);
-    Mat magy = Mat(h,w,CV_32FC1);
+        Sobel(planes[2], magXB, CV_32FC1, 1, 0, 3);
+        Sobel(planes[2], magYB, CV_32FC1, 0, 1, 3);
+    }
 
     Mat mag1 = Mat(h,w,CV_32FC1);
     Mat mag2 = Mat(h,w,CV_32FC1);
@@ -169,39 +179,21 @@ void Domain_Filter::find_magnitude(Mat &img, Mat &mag)
     magnitude(magXG,magYG,mag2);
     magnitude(magXB,magYB,mag3);
 
-    for(int i =0;i < h;i++)
-        for(int j=0;j<w;j++)
-        {
-            mag.at<float>(i,j) = mag1.at<float>(i,j) + mag2.at<float>(i,j) + mag3.at<float>(i,j);
-        }
-
-
-    for(int i =0;i < h;i++)
-        for(int j=0;j<w;j++)
-        {
-            mag.at<float>(i,j) = 1.0f -  mag.at<float>(i,j);
-        }
-
+    mag = mag1 + mag2 + mag3;
+    mag = 1.0f - mag;
 }
 
 void Domain_Filter::compute_Rfilter(Mat &output, Mat &hz, float sigma_h)
 {
-
-    float a;
-
     int h = output.rows;
     int w = output.cols;
     int channel = output.channels();
 
-    a = (float) exp((-1.0 * sqrt(2.0)) / sigma_h);
+    float a = (float) exp((-1.0 * sqrt(2.0)) / sigma_h);
 
     Mat temp = Mat(h,w,CV_32FC3);
 
-    for(int i =0; i < h;i++)
-        for(int j=0;j<w;j++)
-            for(int c=0;c<channel;c++)
-                temp.at<float>(i,j*channel+c) = output.at<float>(i,j*channel+c);
-
+    output.copyTo(temp);
     Mat V = Mat(h,w,CV_32FC1);
 
     for(int i=0;i<h;i++)
@@ -232,14 +224,7 @@ void Domain_Filter::compute_Rfilter(Mat &output, Mat &hz, float sigma_h)
         }
     }
 
-
-    for(int i =0; i < h;i++)
-        for(int j=0;j<w;j++)
-            for(int c=0;c<channel;c++)
-                output.at<float>(i,j*channel+c) = temp.at<float>(i,j*channel+c);
-
-    temp.release();
-    V.release();
+    temp.copyTo(output);
 }
 
 void Domain_Filter::compute_boxfilter(Mat &output, Mat &hz, Mat &psketch, float radius)
@@ -249,12 +234,8 @@ void Domain_Filter::compute_boxfilter(Mat &output, Mat &hz, Mat &psketch, float
     Mat lower_pos = Mat(h,w,CV_32FC1);
     Mat upper_pos = Mat(h,w,CV_32FC1);
 
-    for(int i=0;i<h;i++)
-        for(int j=0;j<w;j++)
-        {
-            lower_pos.at<float>(i,j) = hz.at<float>(i,j) - radius;
-            upper_pos.at<float>(i,j) = hz.at<float>(i,j) + radius;
-        }
+    lower_pos = hz - radius;
+    upper_pos = hz + radius;
 
     lower_idx = Mat::zeros(h,w,CV_32FC1);
     upper_idx = Mat::zeros(h,w,CV_32FC1);
@@ -334,19 +315,11 @@ void Domain_Filter::compute_boxfilter(Mat &output, Mat &hz, Mat &psketch, float
             upper_idx.at<float>(i,j) = temp_upper_idx.at<float>(0,j) + 1;
         }
 
-        lower_pos_row.release();
-        upper_pos_row.release();
-        temp_lower_idx.release();
-        temp_upper_idx.release();
     }
-    for(int i=0;i<h;i++)
-        for(int j=0;j<w;j++)
-            psketch.at<float>(i,j) = upper_idx.at<float>(i,j) - lower_idx.at<float>(i,j);
-
+    psketch = upper_idx - lower_idx;
 }
 void Domain_Filter::compute_NCfilter(Mat &output, Mat &hz, Mat &psketch, float radius)
 {
-
     int h = output.rows;
     int w = output.cols;
     int channel = output.channels();
@@ -377,24 +350,25 @@ void Domain_Filter::compute_NCfilter(Mat &output, Mat &hz, Mat &psketch, float r
     Mat a = Mat::zeros(h,w,CV_32FC1);
     Mat b = Mat::zeros(h,w,CV_32FC1);
 
+    // Compute the box filter using a summed area table.
     for(int c=0;c<channel;c++)
     {
         Mat flag = Mat::ones(h,w,CV_32FC1);
-        for(int i=0;i<h;i++)
-            for(int j=0;j<w;j++)
-                flag.at<float>(i,j) = (c+1)*flag.at<float>(i,j);
+        multiply(flag,c+1,flag);
 
-        for(int i=0;i<h;i++)
-            for(int j=0;j<w;j++)
-            {
-                a.at<float>(i,j) = (flag.at<float>(i,j) - 1) * h * (w+1) + (lower_idx.at<float>(i,j) - 1) * h + indices.at<float>(i,j);
-                b.at<float>(i,j) = (flag.at<float>(i,j) - 1) * h * (w+1) + (upper_idx.at<float>(i,j) - 1) * h + indices.at<float>(i,j);
+        Mat temp1, temp2;
+        multiply(flag - 1,h*(w+1),temp1);
+        multiply(lower_idx - 1,h,temp2);
+        a = temp1 + temp2 + indices;
 
-            }
+        multiply(flag - 1,h*(w+1),temp1);
+        multiply(upper_idx - 1,h,temp2);
+        b = temp1 + temp2 + indices;
 
         int p,q,r,rem;
         int p1,q1,r1,rem1;
 
+        // Calculating indices
         for(int i=0;i<h;i++)
         {
             for(int j=0;j<w;j++)
@@ -416,7 +390,6 @@ void Domain_Filter::compute_NCfilter(Mat &output, Mat &hz, Mat &psketch, float r
                     q=q-1;
                 }
 
-
                 r1 = (int) a.at<float>(i,j)/(h*(w+1));
                 rem1 = (int) a.at<float>(i,j) - r1*h*(w+1);
                 q1 = rem1/h;
@@ -427,18 +400,13 @@ void Domain_Filter::compute_NCfilter(Mat &output, Mat &hz, Mat &psketch, float r
                     q1=q1-1;
                 }
 
-
                 final.at<float>(i,j*channel+2-c) = (box_filter.at<float>(p-1,q*channel+(2-r)) - box_filter.at<float>(p1-1,q1*channel+(2-r1)))
                     /(upper_idx.at<float>(i,j) - lower_idx.at<float>(i,j));
             }
         }
     }
 
-    for(int i=0;i<h;i++)
-        for(int j=0;j<w;j++)
-            for(int c=0;c<channel;c++)
-                output.at<float>(i,j*channel+c) = final.at<float>(i,j*channel+c);
-
+    final.copyTo(output);
 }
 void Domain_Filter::init(const Mat &img, int flags, float sigma_s, float sigma_r)
 {
@@ -482,20 +450,15 @@ void Domain_Filter::init(const Mat &img, int flags, float sigma_s, float sigma_r
 
     Mat final = Mat(h,w,CV_32FC3);
 
-    for(int i = 0; i < h; i++)
-        for(int j = 0; j < w; j++)
-        {
-            horiz.at<float>(i,j) = (float) 1.0 + (sigma_s/sigma_r) * distx.at<float>(i,j);
-            vert.at<float>(i,j) = (float) 1.0 + (sigma_s/sigma_r) * disty.at<float>(i,j);
-        }
+    Mat tempx,tempy;
+    multiply(distx,sigma_s/sigma_r,tempx);
+    multiply(disty,sigma_s/sigma_r,tempy);
 
+    horiz = 1.0f + tempx;
+    vert = 1.0f + tempy;
 
     O = Mat(h,w,CV_32FC3);
-
-    for(int  i =0;i<h;i++)
-        for(int  j =0;j<w;j++)
-            for(int c=0;c<channel;c++)
-                O.at<float>(i,j*channel+c) =  img.at<float>(i,j*channel+c);
+    img.copyTo(O);
 
     O_t = Mat(w,h,CV_32FC3);
 
@@ -588,17 +551,14 @@ void Domain_Filter::pencil_sketch(const Mat &img, Mat &sketch, Mat &color_res, f
     init(img,2,sigma_s,sigma_r);
     int h = img.size().height;
     int w = img.size().width;
-    int channel = img.channels();
 
     /////////////////////// convert to YCBCR model for color pencil drawing //////////////////////////////////////////////////////
 
     Mat color_sketch = Mat(h,w,CV_32FC3);
-    Mat Y_channel = Mat(h,w,CV_32FC1);
-    Mat U_channel = Mat(h,w,CV_32FC1);
-    Mat V_channel = Mat(h,w,CV_32FC1);
 
     cvtColor(img,color_sketch,COLOR_BGR2YCrCb);
 
+    vector <Mat> YUV_channel;
     Mat vert_t = ct_V.t();
 
     float sigma_h = sigma_s;
@@ -635,33 +595,11 @@ void Domain_Filter::pencil_sketch(const Mat &img, Mat &sketch, Mat &color_res, f
         if(i==0)
         {
             sketch = pen_res.clone();
-
-            for(int k = 0; k < h; k++)
-                for(int j = 0; j < w; j++)
-                {
-                    Y_channel.at<float>(k,j) = color_sketch.at<float>(k,j*channel+0);
-                    U_channel.at<float>(k,j) = color_sketch.at<float>(k,j*channel+1);
-                    V_channel.at<float>(k,j) = color_sketch.at<float>(k,j*channel+2);
-                }
-
-
-            for(int k=0;k<h;k++)
-                for(int j=0;j<w;j++)
-                    Y_channel.at<float>(k,j) = pen_res.at<float>(k,j);
-
-            // cvMerge(Y_channel,U_channel,V_channel,0,color_sketch);
-            for(int k = 0; k < h; k++)
-                for(int j = 0; j < w; j++)
-                {
-                    color_sketch.at<float>(k,j*channel+0) = Y_channel.at<float>(k,j);
-                    color_sketch.at<float>(k,j*channel+1) = U_channel.at<float>(k,j);
-                    color_sketch.at<float>(k,j*channel+2) = V_channel.at<float>(k,j);
-                }
-
+            split(color_sketch,YUV_channel);
+            pen_res.copyTo(YUV_channel[0]);
+            merge(YUV_channel,color_sketch);
             cvtColor(color_sketch,color_res,COLOR_YCrCb2BGR);
-
         }
 
     }
-
 }
index c49f4a2..63f4530 100644 (file)
@@ -41,7 +41,6 @@
 
 #include "precomp.hpp"
 #include "opencv2/photo.hpp"
-#include "opencv2/highgui.hpp"
 #include <iostream>
 #include <stdlib.h>
 
@@ -100,30 +99,17 @@ void cv::seamlessClone(InputArray _src, InputArray _dst, InputArray _mask, Point
         exit(0);
     }
 
-    for(int i=minx, k=minxd;i<(minx+lenx);i++,k++)
-        for(int j=miny,l=minyd ;j<(miny+leny);j++,l++)
-        {
-            dst_mask.at<uchar>(k,l) = gray.at<uchar>(i,j);
-        }
+    Rect roi_d(minyd,minxd,leny,lenx);
+    Rect roi_s(miny,minx,leny,lenx);
 
-    int channel = 3;
+    Mat destinationROI = dst_mask(roi_d);
+    Mat sourceROI = cs_mask(roi_s);
 
-    for(int i=minx;i<(minx+lenx);i++)
-        for(int j=miny;j<(miny+leny);j++)
-        {
-            for(int c=0;c<3;c++)
-            {
-                if(gray.at<uchar>(i,j) == 255)
-                    cs_mask.at<uchar>(i,j*channel+c) = src.at<uchar>(i,j*channel+c);
-            }
-        }
+    gray(roi_s).copyTo(destinationROI);
+    src(roi_s).copyTo(sourceROI,gray(roi_s));
 
-    for(int i=minx, k=minxd;i<(minx+lenx);i++,k++)
-        for(int j=miny,l=minyd ;j<(miny+leny);j++,l++)
-        {
-            for(int c=0;c<channel;c++)
-                cd_mask.at<uchar>(k,l*channel+c) = cs_mask.at<uchar>(i,j*channel+c);
-        }
+    destinationROI = cd_mask(roi_d);
+    cs_mask(roi_s).copyTo(destinationROI);
 
     Cloning obj;
     obj.normal_clone(dest,cd_mask,dst_mask,blend,flags);
index 90bd2e4..612d764 100644 (file)
@@ -555,8 +555,6 @@ void Cloning::local_color_change(Mat &I, Mat &mask, Mat &wmask, Mat &cloned, flo
 
 void Cloning::illum_change(Mat &I, Mat &mask, Mat &wmask, Mat &cloned, float alpha, float beta)
 {
-    int channel = I.channels();
-
     initialization(I,mask,wmask);
 
     array_product(srx32,sgx,smask);
@@ -565,18 +563,20 @@ void Cloning::illum_change(Mat &I, Mat &mask, Mat &wmask, Mat &cloned, float alp
     Mat mag = Mat(I.size(),CV_32FC3);
     magnitude(srx32,sry32,mag);
 
-    for(int i=0;i < I.size().height; i++)
-        for(int j=0; j < I.size().width; j++)
-            for(int c=0;c < channel;++c)
-            {
-                if(srx32.at<float>(i,j*channel+c) != 0)
-                {
-                    srx32.at<float>(i,j*channel+c) =
-                        pow(alpha,beta)*srx32.at<float>(i,j*channel+c)*pow(mag.at<float>(i,j*channel+c),-1*beta);
-                    sry32.at<float>(i,j*channel+c) =
-                        pow(alpha,beta)*sry32.at<float>(i,j*channel+c)*pow(mag.at<float>(i,j*channel+c),-1*beta);
-                }
-            }
+    Mat multX, multY, multx_temp, multy_temp;
+
+    multiply(srx32,pow(alpha,beta),multX);
+    pow(mag,-1*beta, multx_temp);
+    multiply(multX,multx_temp,srx32);
+
+    multiply(sry32,pow(alpha,beta),multY);
+    pow(mag,-1*beta, multy_temp);
+    multiply(multY,multy_temp,sry32);
+
+    Mat zeroMask = (srx32 != 0);
+
+    srx32.copyTo(srx32, zeroMask);
+    sry32.copyTo(sry32, zeroMask);
 
     evaluate(I,wmask,cloned);
 }
@@ -584,23 +584,16 @@ void Cloning::illum_change(Mat &I, Mat &mask, Mat &wmask, Mat &cloned, float alp
 void Cloning::texture_flatten(Mat &I, Mat &mask, Mat &wmask, double low_threshold,
         double high_threshold, int kernel_size, Mat &cloned)
 {
-    int channel = mask.channels();
-
     initialization(I,mask,wmask);
 
     Mat out = Mat(mask.size(),CV_8UC1);
     Canny(mask,out,low_threshold,high_threshold,kernel_size);
 
-    for(int i=0;i<mask.size().height;i++)
-        for(int j=0;j<mask.size().width;j++)
-            for(int c=0;c<channel;c++)
-            {
-                if(out.at<uchar>(i,j) != 255)
-               {
-                    sgx.at<float>(i,j*channel+c) = 0.0;
-                    sgy.at<float>(i,j*channel+c) = 0.0;
-               }
-            }
+    Mat zeros(sgx.size(), CV_32FC3);
+    zeros.setTo(0);
+    Mat zerosMask = (out != 255);
+    zeros.copyTo(sgx, zerosMask);
+    zeros.copyTo(sgy, zerosMask);
 
     array_product(srx32,sgx,smask);
     array_product(sry32,sgy,smask);