C++ MatrixView::colsize方法代码示例

    void SBMoffat::SBMoffatImpl::fillXValue(tmv::MatrixView<double> val,
                                            double x0, double dx, double dxy,
                                            double y0, double dy, double dyx) const
    {
        dbg<<"SBMoffat fillXValue\n";
        dbg<<"x = "<<x0<<" + i * "<<dx<<" + j * "<<dxy<<std::endl;
        dbg<<"y = "<<y0<<" + i * "<<dyx<<" + j * "<<dy<<std::endl;
        assert(val.stepi() == 1);
        assert(val.canLinearize());
        const int m = val.colsize();
        const int n = val.rowsize();
        typedef tmv::VIt<double,1,tmv::NonConj> It;

        x0 *= _inv_rD;
        dx *= _inv_rD;
        dxy *= _inv_rD;
        y0 *= _inv_rD;
        dy *= _inv_rD;
        dyx *= _inv_rD;

        It valit = val.linearView().begin();
        for (int j=0;j<n;++j,x0+=dxy,y0+=dy) {
            double x = x0;
            double y = y0;
            for (int i=0;i<m;++i,x+=dx,y+=dyx) {
                double rsq = x*x + y*y;
                if (rsq > _maxRrD_sq) *valit++ = 0.;
                else *valit++ = _norm / _pow_beta(1.+rsq, _beta);
            }
        }
    }

    void SBExponential::SBExponentialImpl::fillXValue(tmv::MatrixView<double> val,
                                                      double x0, double dx, double dxy,
                                                      double y0, double dy, double dyx) const
    {
        dbg<<"SBExponential fillXValue\n";
        dbg<<"x = "<<x0<<" + i * "<<dx<<" + j * "<<dxy<<std::endl;
        dbg<<"y = "<<y0<<" + i * "<<dyx<<" + j * "<<dy<<std::endl;
        assert(val.stepi() == 1);
        assert(val.canLinearize());
        const int m = val.colsize();
        const int n = val.rowsize();
        typedef tmv::VIt<double,1,tmv::NonConj> It;

        x0 *= _inv_r0;
        dx *= _inv_r0;
        dxy *= _inv_r0;
        y0 *= _inv_r0;
        dy *= _inv_r0;
        dyx *= _inv_r0;

        It valit = val.linearView().begin();
        for (int j=0;j<n;++j,x0+=dxy,y0+=dy) {
            double x = x0;
            double y = y0;
            for (int i=0;i<m;++i,x+=dx,y+=dyx) *valit++ = _norm * std::exp(-sqrt(x*x + y*y));
        }
    }

    void SBTopHat::SBTopHatImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                            double kx0, double dkx, double dkxy,
                                            double ky0, double dky, double dkyx) const
    {
        dbg<<"SBTopHat fillKValue\n";
        dbg<<"kx = "<<kx0<<" + i * "<<dkx<<" + j * "<<dkxy<<std::endl;
        dbg<<"ky = "<<ky0<<" + i * "<<dkyx<<" + j * "<<dky<<std::endl;
        assert(val.stepi() == 1);
        assert(val.canLinearize());
        const int m = val.colsize();
        const int n = val.rowsize();
        typedef tmv::VIt<std::complex<double>,1,tmv::NonConj> It;

        kx0 *= _r0;
        dkx *= _r0;
        dkxy *= _r0;
        ky0 *= _r0;
        dky *= _r0;
        dkyx *= _r0;

        It valit = val.linearView().begin();
        for (int j=0;j<n;++j,kx0+=dkxy,ky0+=dky) {
            double kx = kx0;
            double ky = ky0;
            for (int i=0;i<m;++i,kx+=dkx,ky+=dkyx) {
                double ksq = kx*kx + ky*ky;
                *valit++ = kValue2(ksq);
            }
        }
    }

    void SBTopHat::SBTopHatImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                            double kx0, double dkx, int izero,
                                            double ky0, double dky, int jzero) const
    {
        dbg<<"SBTopHat fillKValue\n";
        dbg<<"kx = "<<kx0<<" + i * "<<dkx<<", izero = "<<izero<<std::endl;
        dbg<<"ky = "<<ky0<<" + j * "<<dky<<", jzero = "<<jzero<<std::endl;
        if (izero != 0 || jzero != 0) {
            xdbg<<"Use Quadrant\n";
            fillKValueQuadrant(val,kx0,dkx,izero,ky0,dky,jzero);
        } else {
            xdbg<<"Non-Quadrant\n";
            assert(val.stepi() == 1);
            const int m = val.colsize();
            const int n = val.rowsize();
            typedef tmv::VIt<std::complex<double>,1,tmv::NonConj> It;

            kx0 *= _r0;
            dkx *= _r0;
            ky0 *= _r0;
            dky *= _r0;

            for (int j=0;j<n;++j,ky0+=dky) {
                double kx = kx0;
                double kysq = ky0*ky0;
                It valit = val.col(j).begin();
                for (int i=0;i<m;++i,kx+=dkx) {
                    double ksq = kx*kx + kysq;
                    *valit++ = kValue2(ksq);
                }
            }
        }
    }

    void SBAiry::SBAiryImpl::fillXValue(tmv::MatrixView<double> val,
                                        double x0, double dx, int izero,
                                        double y0, double dy, int jzero) const
    {
        dbg<<"SBAiry fillXValue\n";
        dbg<<"x = "<<x0<<" + i * "<<dx<<", izero = "<<izero<<std::endl;
        dbg<<"y = "<<y0<<" + j * "<<dy<<", jzero = "<<jzero<<std::endl;
        if (izero != 0 || jzero != 0) {
            xdbg<<"Use Quadrant\n";
            fillXValueQuadrant(val,x0,dx,izero,y0,dy,jzero);
        } else {
            xdbg<<"Non-Quadrant\n";
            assert(val.stepi() == 1);
            const int m = val.colsize();
            const int n = val.rowsize();
            typedef tmv::VIt<double,1,tmv::NonConj> It;

            x0 *= _D;
            dx *= _D;
            y0 *= _D;
            dy *= _D;

            for (int j=0;j<n;++j,y0+=dy) {
                double x = x0;
                double ysq = y0*y0;
                It valit = val.col(j).begin();
                for (int i=0;i<m;++i,x+=dx)
                    *valit++ = _xnorm * _info->xValue(sqrt(x*x + ysq));
            }
        }
    }

    void SBMoffat::SBMoffatImpl::fillXValue(tmv::MatrixView<double> val,
                                            double x0, double dx, int izero,
                                            double y0, double dy, int jzero) const
    {
        dbg<<"SBMoffat fillXValue\n";
        dbg<<"x = "<<x0<<" + i * "<<dx<<", izero = "<<izero<<std::endl;
        dbg<<"y = "<<y0<<" + j * "<<dy<<", jzero = "<<jzero<<std::endl;
        if (izero != 0 || jzero != 0) {
            xdbg<<"Use Quadrant\n";
            fillXValueQuadrant(val,x0,dx,izero,y0,dy,jzero);
        } else {
            xdbg<<"Non-Quadrant\n";
            assert(val.stepi() == 1);
            const int m = val.colsize();
            const int n = val.rowsize();
            typedef tmv::VIt<double,1,tmv::NonConj> It;

            x0 *= _inv_rD;
            dx *= _inv_rD;
            y0 *= _inv_rD;
            dy *= _inv_rD;

            for (int j=0;j<n;++j,y0+=dy) {
                double x = x0;
                double ysq = y0*y0;
                It valit = val.col(j).begin();
                for (int i=0;i<m;++i,x+=dx) {
                    double rsq = x*x + ysq;
                    if (rsq > _maxRrD_sq) *valit++ = 0.;
                    else *valit++ = _norm / _pow_beta(1.+rsq, _beta);
                }
            }
        }
    }

 void SBShapelet::SBShapeletImpl::fillXValue(
     tmv::MatrixView<double> val,
     const tmv::Matrix<double>& x, const tmv::Matrix<double>& y) const
 {
     dbg<<"order = "<<_bvec.getOrder()<<", sigma = "<<_sigma<<std::endl;
     xdbg<<"fillXValue with bvec = "<<_bvec<<std::endl;
     assert(val.stepi() == 1);
     assert(val.canLinearize());
     const int m = val.colsize();
     const int n = val.rowsize();
     tmv::Matrix<double> psi(m*n,_bvec.size());
     LVector::basis(x.constLinearView(),y.constLinearView(),psi.view(),
                    _bvec.getOrder(),_sigma);
     val.linearView() = psi * _bvec.rVector();
 }

 void SBAdd::SBAddImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                   double kx0, double dkx, double dkxy,
                                   double ky0, double dky, double dkyx) const
 {
     dbg<<"SBAdd fillKValue\n";
     dbg<<"kx = "<<kx0<<" + i * "<<dkx<<" + j * "<<dkxy<<std::endl;
     dbg<<"ky = "<<ky0<<" + i * "<<dkyx<<" + j * "<<dky<<std::endl;
     ConstIter pptr = _plist.begin();
     assert(pptr != _plist.end());
     GetImpl(*pptr)->fillKValue(val,kx0,dkx,dkxy,ky0,dky,dkyx);
     if (++pptr != _plist.end()) {
         tmv::Matrix<std::complex<double> > val2(val.colsize(),val.rowsize());
         for (; pptr != _plist.end(); ++pptr) {
             GetImpl(*pptr)->fillKValue(val2.view(),kx0,dkx,dkxy,ky0,dky,dkyx);
             val += val2;
         }
     }
 }

 void SBConvolve::SBConvolveImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                             double kx0, double dkx, int izero,
                                             double ky0, double dky, int jzero) const
 {
     dbg<<"SBConvolve fillKValue\n";
     dbg<<"kx = "<<kx0<<" + i * "<<dkx<<", izero = "<<izero<<std::endl;
     dbg<<"ky = "<<ky0<<" + j * "<<dky<<", jzero = "<<jzero<<std::endl;
     ConstIter pptr = _plist.begin();
     assert(pptr != _plist.end());
     GetImpl(*pptr)->fillKValue(val,kx0,dkx,izero,ky0,dky,jzero);
     if (++pptr != _plist.end()) {
         tmv::Matrix<std::complex<double> > val2(val.colsize(),val.rowsize());
         for (; pptr != _plist.end(); ++pptr) {
             GetImpl(*pptr)->fillKValue(val2.view(),kx0,dkx,izero,ky0,dky,jzero);
             val = ElemProd(val,val2);
         }
     }
 }

 void SBConvolve::SBConvolveImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                             double x0, double dx, double dxy,
                                             double y0, double dy, double dyx) const
 {
     dbg<<"SBConvolve fillKValue\n";
     dbg<<"x = "<<x0<<" + ix * "<<dx<<" + iy * "<<dxy<<std::endl;
     dbg<<"y = "<<y0<<" + ix * "<<dyx<<" + iy * "<<dy<<std::endl;
     ConstIter pptr = _plist.begin();
     assert(pptr != _plist.end());
     GetImpl(*pptr)->fillKValue(val,x0,dx,dxy,y0,dy,dyx);
     if (++pptr != _plist.end()) {
         tmv::Matrix<std::complex<double> > val2(val.colsize(),val.rowsize());
         for (; pptr != _plist.end(); ++pptr) {
             GetImpl(*pptr)->fillKValue(val2.view(),x0,dx,dxy,y0,dy,dyx);
             val = ElemProd(val,val2);
         }
     }
 }

    void SBExponential::SBExponentialImpl::fillKValue(tmv::MatrixView<std::complex<double> > val,
                                                      double kx0, double dkx, int izero,
                                                      double ky0, double dky, int jzero) const
    {
        dbg<<"SBExponential fillKValue\n";
        dbg<<"kx = "<<kx0<<" + i * "<<dkx<<", izero = "<<izero<<std::endl;
        dbg<<"ky = "<<ky0<<" + j * "<<dky<<", jzero = "<<jzero<<std::endl;
        if (izero != 0 || jzero != 0) {
            xdbg<<"Use Quadrant\n";
            fillKValueQuadrant(val,kx0,dkx,izero,ky0,dky,jzero);
        } else {
            xdbg<<"Non-Quadrant\n";
            assert(val.stepi() == 1);
            const int m = val.colsize();
            const int n = val.rowsize();
            typedef tmv::VIt<std::complex<double>,1,tmv::NonConj> It;

            kx0 *= _r0;
            dkx *= _r0;
            ky0 *= _r0;
            dky *= _r0;

            for (int j=0;j<n;++j,ky0+=dky) {
                double kx = kx0;
                double kysq = ky0*ky0;
                It valit = val.col(j).begin();
                for (int i=0;i<m;++i,kx+=dkx) {
                    double ksq = kx*kx + kysq;
                    if (ksq > _ksq_max) {
                        *valit++ = 0.;
                    } else if (ksq < _ksq_min) {
                        *valit++ = _flux * (1. - 1.5*ksq*(1. - 1.25*ksq));
                    } else {
                        double temp = 1. + ksq;
                        *valit++ =  _flux/(temp*sqrt(temp));
                    }
                }
            }
        }
    }

    void SBBox::SBBoxImpl::fillXValue(tmv::MatrixView<double> val,
                                      double x0, double dx, int izero,
                                      double y0, double dy, int jzero) const
    {
        dbg<<"SBBox fillXValue\n";
        dbg<<"x = "<<x0<<" + i * "<<dx<<", izero = "<<izero<<std::endl;
        dbg<<"y = "<<y0<<" + j * "<<dy<<", jzero = "<<jzero<<std::endl;

        assert(val.stepi() == 1);
        const int m = val.colsize();
        const int n = val.rowsize();
        typedef tmv::VIt<double,1,tmv::NonConj> It;

        // It will be useful to do everything in units of dx,dy
        x0 /= dx;
        double wo2 = _wo2 / std::abs(dx);
        y0 /= dy;
        double ho2 = _ho2 / std::abs(dy);
        xdbg<<"x0,y0 -> "<<x0<<','<<y0<<std::endl;
        xdbg<<"width,height -> "<<wo2*2.<<','<<ho2*2.<<std::endl;

        // Start by setting everything to zero
        val.setZero();

        // Then fill the interior with _norm:
        // Fill pixels where:
        //     x0 + ix >= -width/2
        //     x0 + ix < width/2
        //     y0 + iy >= -width/2
        //     y0 + iy < width/2

        int ix1 = std::max(0, int(std::ceil(-wo2 - x0)));
        int ix2 = std::min(m, int(std::ceil(wo2 - x0)));
        int iy1 = std::max(0, int(std::ceil(-ho2 - y0)));
        int iy2 = std::min(n, int(std::ceil(ho2 - y0)));

        if (ix1 < ix2 && iy1 < iy2)
            val.subMatrix(ix1,ix2,iy1,iy2).setAllTo(_norm);

#if 0
        // We used to implement this by making the pixels that cross the edge have a
        // fractional flux value appropriate for the fraction of the box that goes through
        // each pixel.  However, this isn't actually correct.  SBProfile objects are always
        // rendered as the local surface brightness at the center of the pixel.  To get
        // the right flux, you need to convolve by a Pixel.  So if someone renders a Box
        // without convolving by a pixel, it is inconsistent to do this differently than we
        // do all the other SBProfile types.  However, since it was an involved calculation
        // and someone might actually care to resurrect it in a different guise at some point,
        // I'm leaving it here, just commented out.

        // We need to make sure the pixels where the edges of the box fall only get
        // a fraction of the flux.
        //
        // We divide up the range into 3 sections in x:
        //    left of the box where val = 0
        //    in the box where val = _norm
        //    right of the box where val = 0 again
        //
        // ... and 3 sections in y:
        //    below the box where val = 0
        //    in the box where val = _norm
        //    above the box where val = 0 again
        //
        // Furthermore, we have to calculate the correct values for the pixels on the border.

        int ix_left, ix_right, iy_bottom, iy_top;
        double x_left, x_right, y_bottom, y_top;

        // Find the x edges:
        double tmp = 0.5*width + 0.5;
        ix_left = int(-tmp-x0+1);
        ix_right = int(tmp-x0);

        // If the box goes off the image, it's ok, but it will cause us problems
        // later on if we don't change it.  Just use ix_left = 0.
        if (ix_left < 0) { ix_left = 0; x_left = 1.; }

        // If the whole box is off the image, just zero and return.
        else if (ix_left >= m) { val.setZero(); return; }

        // Normal case: calculate the fractional flux in the edge
        else x_left = tmp+x0+ix_left;

        // Now the right side.
        if (ix_right >= m) { ix_right = m-1; x_right = 1.; }
        else if (ix_right < 0) { val.setZero(); return; }
        else x_right = tmp-x0-ix_right;
        xdbg<<"ix_left = "<<ix_left<<" with partial flux "<<x_left<<std::endl;
        xdbg<<"ix_right = "<<ix_right<<" with partial flux "<<x_right<<std::endl;

        // Repeat for y values
        tmp = 0.5*height + 0.5;
        iy_bottom = int(-tmp-y0+1);
        iy_top = int(tmp-y0);

        if (iy_bottom < 0) { iy_bottom = 0; y_bottom = 1.; }
        else if (iy_bottom >= n) { val.setZero(); return; }
        else y_bottom = tmp+y0+iy_bottom;

        if (iy_top >= n) { iy_top = n-1; y_top = 1.; }
//.........这里部分代码省略.........