/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4521 by jfenwick, Mon Aug 26 11:51:30 2013 UTC revision 4526 by jfenwick, Mon Sep 2 06:34:25 2013 UTC
# Line 3919  namespace Line 3919  namespace
3919      // Calculates a guassian blur colvolution matrix for 2D      // Calculates a guassian blur colvolution matrix for 2D
3920      double* get2DGauss(unsigned radius, double sigma)      double* get2DGauss(unsigned radius, double sigma)
3921      {      {
3922          double* arr=new double[2*(radius*2+1)];          double* arr=new double[(radius*2+1)*(radius*2+1)];
3923          double common=M_1_PI*0.5*1/(sigma*sigma);          double common=M_1_PI*0.5*1/(sigma*sigma);
3924      double total=0;      double total=0;
3925      for (size_t y=-radius;y<=radius;++y)      int r=static_cast<int>(radius);
3926        for (int y=-r;y<=r;++y)
3927      {      {
3928          for (size_t x=-radius;x<=radius;++x)          for (int x=-r;x<=r;++x)
3929          {          {        
3930              arr[x+y*(radius*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));              arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));
3931              total+=arr[x+y*(radius*2+1)];  // cout << (x+y*(r*2+1)) << " " << arr[(x+r)+(y+r)*(r*2+1)] << endl;
3932                total+=arr[(x+r)+(y+r)*(r*2+1)];
3933          }          }
3934      }      }
3935      double invtotal=1/total;      double invtotal=1/total;
3936      for (size_t p=0;p<(2*(radius*2+1));++p)  //cout << "Inv total is "    << invtotal << endl;
3937        for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)
3938      {      {
3939          arr[p]*=invtotal;          arr[p]*=invtotal;
3940    //cout << p << "->" << arr[p] << endl;      
3941      }      }
3942      return arr;      return arr;
3943      }      }
# Line 3942  namespace Line 3946  namespace
3946      // (xp, yp) are the coords in the source matrix not the destination matrix      // (xp, yp) are the coords in the source matrix not the destination matrix
3947      double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)      double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)
3948      {      {
3949            size_t bx=xp-radius, by=yp-radius;
3950        size_t sbase=bx+by*width;
3951      double result=0;      double result=0;
3952      for (size_t y=-radius;y<=radius;++y)      for (int y=0;y<2*radius+1;++y)
3953      {      {    
3954          for (size_t x=-radius;x<=radius;++x)          for (int x=0;x<2*radius+1;++x)
3955          {          {
3956              result+=conv[x+y*(2*radius+1)] * source[xp + yp*width];              result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];
3957          }          }
3958      }      }
3959          return result;                return result;      
# Line 3986  escript::Data Rectangle::randomFill(long Line 3992  escript::Data Rectangle::randomFill(long
3992      {      {
3993          throw RipleyException("Sigma must be a postive floating point number.");          throw RipleyException("Sigma must be a postive floating point number.");
3994      }          }    
3995        size_t numpoints[2];
3996        numpoints[0]=m_ownNE[0]+1;
3997        numpoints[1]=m_ownNE[1]+1;
3998      size_t padding=max((unsigned)max((m_NE[0]-m_ownNE[0])/2, (m_NE[1]-m_ownNE[1])/2), radius);      size_t padding=max((unsigned)max((m_NE[0]-m_ownNE[0])/2, (m_NE[1]-m_ownNE[1])/2), radius);
3999      size_t width=(m_ownNE[0]+2*padding);          size_t width=(numpoints[0]+2*padding);      // width of one row in points
4000      size_t dsize=width * (m_ownNE[1]+2*padding);      size_t height=(numpoints[1]+2*padding); // height of one row in points
4001        size_t dsize=width * height; // size of padded source grid
4002            
4003      double* src=new double[dsize];      double* src=new double[dsize];
4004      esysUtils::randomFillArray(seed, src, dsize);        esysUtils::randomFillArray(seed, src, dsize);  
4005            
4006      // Now we need to copy the regions owned by other ranks over here        // Now we need to copy the regions owned by other ranks over here  
4007        
4008      // Lets call that done for now      // Lets call that done for now
# Line 4001  escript::Data Rectangle::randomFill(long Line 4011  escript::Data Rectangle::randomFill(long
4011      // don't need to check for exwrite because we just made it      // don't need to check for exwrite because we just made it
4012      escript::DataVector& dv=resdat.getExpandedVectorReference();      escript::DataVector& dv=resdat.getExpandedVectorReference();
4013      double* convolution=get2DGauss(radius, sigma);      double* convolution=get2DGauss(radius, sigma);
4014      for (size_t y=0;y<m_ownNE[1];++y)          for (size_t y=0;y<(m_ownNE[1]+1);++y)    
4015      {      {
4016          for (size_t x=0;x<m_ownNE[0];++x)          for (size_t x=0;x<(m_ownNE[0]+1);++x)
4017      {      {    
4018          dv[x+y*width]=Convolve2D(convolution, src, x, y, radius, width);          dv[x+y*(m_ownNE[0]+1)]=Convolve2D(convolution, src, x+radius, y+radius, radius, width);
4019      }      }
4020      }      }
4021        delete[] convolution;
4022      return resdat;      return resdat;
4023  }  }
4024    

Legend:
Removed from v.4521  
changed lines
  Added in v.4526

  ViewVC Help
Powered by ViewVC 1.1.26