C++ MultiArray::shape方法代码示例

Fields FieldAlgorithms::fieldsByLaplasian(MultiArray<2, float> & image) 
{
    MultiArray<2, float> valley(image.shape());
    Pyramid pyramid(image);
    for (int i = 3; i > 0; i--)
    {
        MultiArray<2, float> resized(pyramid.get(i));
        MultiArray<2, float> tmpArr(resized.shape());
        laplacianOfGaussianMultiArray(resized, tmpArr, 1.0);
        valley += pyramid.toOriginalSize(tmpArr);
    }
    MultiArray<2, float> peak = valley * -1;
    //remove DC component
    float thrhld = valley[argMax(valley)] * 0.3;
    threshold(valley, valley, thrhld);
    thrhld = peak[argMax(peak)] * 0.3;
    threshold(peak, peak, thrhld);
    

    //Edge as Gradients
    MultiArray<2, float> edgeField(image.shape());
    gaussianGradientMagnitude(image, edgeField, 1.0);

    //Localize
    std::vector<Shape2> valleyLocals(localizePOI(image));

    //Result as Field Object
    Fields fields(valley, valleyLocals, peak, edgeField, image);
    //A heuristic initialization of the localization, as a priori known position
    Shape2 nextToIris = Shape2(image.width() / 2, image.height() / 2);
    fields.specializedIrisValley = localizeByFollowingLocalMaxima(image, nextToIris);
    return fields;
}

void mark(
    MultiArray& m,const Sequence& pts,
    const point& p,int dx,int dy)
{
    for(int k=p.x<0||p.y<0||p.x>=m.shape()[0]||p.y>=m.shape()[1]?1:0,
            k_end=distZ2inf(p,pts.begin(),pts.end()); k<k_end; ++k) {
        int x=p.x-k*dx;
        int y=p.y-k*dy;
        if(x>=0&&y>=0&&x<m.shape()[0]&&y<m.shape()[1])m[x][y]=true;
    }
}

void test() {
    using namespace vigra;
    typedef Thresholding<3, float>::V V;

    MultiArray<3, float> data(V(24,33,40));
    FillRandom<float, float*>::fillRandom(data.data(), data.data()+data.size());
    {
        HDF5File f("test.h5", HDF5File::Open);
        f.write("test", data);
    }

    SourceHDF5<3, float> source("test.h5", "test");
    SinkHDF5<3, vigra::UInt8> sink("thresh.h5", "thresh");
    sink.setBlockShape(V(10,10,10));

    Thresholding<3, float> bs(&source, V(6,4,7));

    bs.run(0.5, 0, 1, &sink);

    HDF5File f("thresh.h5", HDF5File::Open);
    MultiArray<3, UInt8> r;
    f.readAndResize("thresh", r);

    MultiArray<3, UInt8> t(data.shape());
    transformMultiArray(srcMultiArrayRange(data), destMultiArray(t), Threshold<float, UInt8>(-std::numeric_limits<float>::max(), 0.5, 1, 0));
    shouldEqual(t.shape(), r.shape());
    shouldEqualSequence(r.begin(), r.end(), t.begin());
}

   MultiArray<2, float > FieldAlgorithms::matchGradients(MultiArray<2, TinyVector<float, 2> > &imageGradients, MultiArray<2, TinyVector<float, 2> > &mask)
{
    MultiArray<2, float >  dest(imageGradients.shape());
    dest = 0;
    //to make sure, that the mask is always fully within the image, consider the mask shape
    int diff = (mask.width() -1) / 2;
    for (int x = diff; x + diff < imageGradients.width(); x ++) 
    {
        for (int y = diff; y + diff < imageGradients.height(); y++)
        {
            //The masks center is at point x,y now
            //Every vector of the image currently 'covered' by the mask, is compared to its corresponding vector in the mask.
            float vals = 0;
            for (int xM = 0; xM < mask.width(); xM++)
            {
                for (int yM = 0; yM < mask.height(); yM++)
                {
                    TinyVector<float, 2> imageVal = imageGradients((x - diff) + xM, (y - diff) + yM);
                    TinyVector<float, 2> maskVal = mask(xM, yM);
                    vals += compareVectors(imageVal, maskVal);
                }
            }
            int res = vals / (mask.size());
            dest(x,y) = res;
        }
    }
    return dest;
};

Fields FieldAlgorithms::fieldsByErosionDilation(MultiArray<2, float> & image) 
{
    Shape2 shape = image.shape();
    MultiArray<2, float> source(image);
    MultiArray<2, float> t1(shape);
    MultiArray<2, float> t2(shape);

    double radius = 9; //Sinnvoll?

    multiGrayscaleErosion(source, t1, radius);
    multiGrayscaleDilation(source, t2, radius);

    MultiArray<2, float> psi_e = (t1- t2) * -1;
    //Opening of source
    multiGrayscaleErosion(source, t1, radius);
    multiGrayscaleDilation(t1, t2, radius);

    MultiArray<2, float> psi_p = source - t2;
    MultiArray<2, float> psi_pSmooth(psi_p.shape());
    gaussianSmoothMultiArray(psi_p, psi_pSmooth, 8.0);
    //Opening of source
    MultiArray<2, float> inverse = source * -1;
    multiGrayscaleErosion(inverse, t1, 9);
    multiGrayscaleDilation(t1, t2, 9);
    MultiArray<2, float> psi_v = t2 - inverse; 
    MultiArray<2, float> psi_vSmooth(psi_v.shape());
    gaussianSmoothMultiArray(psi_v, psi_vSmooth, 3.0);
    std::vector<Shape2> valleyLocals(0);
    Fields fields(psi_vSmooth, valleyLocals, psi_pSmooth, psi_e, image);
    return fields;
};

MultiArray<2, float > FieldAlgorithms::morphologyByGradientPattern(MultiArray<2, float> & image, MultiArray<2,float> &mask) {
    
    Shape2 shape(image.shape());
    MultiArray<2, TinyVector<float, 2>> imageGradients(shape);
    MultiArray<2, TinyVector<float, 2>> maskGradients(mask.shape());
    gaussianGradientMultiArray(image, imageGradients, 2.0);
    gaussianGradientMultiArray(mask, maskGradients, 2.0);

    //Threshold gradients with small magnitude, 
    //because we only consider directions from now on.
    MultiArray<2, float> magnitudes(shape);
    gaussianGradientMagnitude(image, magnitudes, 3.0);
    float thrhld = magnitudes[argMax(magnitudes)] * 0.3;
    thresholdGrad(magnitudes, imageGradients, thrhld);
    //The actual machting:
    return matchGradients(imageGradients, maskGradients);
   };

Fields FieldAlgorithms::fieldsByGradientPattern(MultiArray<2, float> & image) 
{
    //Get Masks for Valley and Peak
    //not efficient, but can be scaled easly
    vigra::ImageImportInfo valleyInfo("../images/valleyMask.png");
    MultiArray<2, float>  valleyArray(valleyInfo.shape());  
    importImage(valleyInfo, valleyArray);
    MultiArray<2, float > valleyMask(9,9);
    resizeImageNoInterpolation(valleyArray, valleyMask);
    vigra::ImageImportInfo peakInfo("../images/peakMask.png");
    MultiArray<2, float>  peakArray(peakInfo.shape());  
    importImage(peakInfo, peakArray);
    MultiArray<2, float > peakMask(9,9);
    resizeImageNoInterpolation(peakArray, peakMask);

    Pyramid pyramid(image);
    MultiArray<2, float> valleyField(image.shape());
    MultiArray<2, float> peakField(image.shape());
    //go 3 octaves
    for (int i = 3; i > 0; i--)
    {
        MultiArray<2, float> resized(pyramid.get(i));
        MultiArray<2, float> tmpArr = morphologyByGradientPattern(resized, valleyMask);
        valleyField += pyramid.toOriginalSize(tmpArr);
        tmpArr = morphologyByGradientPattern(resized, peakMask);
        peakField += pyramid.toOriginalSize(tmpArr);
    }

    //Edge as Gradients
    MultiArray<2, float> edgeField(image.shape());
    gaussianGradientMagnitude(image, edgeField, 1.0);

    //Localize, smooth for increased range of interaction (following local maxima)
    MultiArray<2, float> smoothed(valleyField.shape());
    gaussianSmoothMultiArray(valleyField, smoothed, 6.0);
    std::vector<Shape2> valleyLocals(localizePOI(smoothed));

    //Result as Field Object
    Fields fields(valleyField, valleyLocals, peakField, edgeField, image);
    //A heuristic initialization of the localization, as a priori known position
    Shape2 nextToIris = Shape2(image.width() / 2 + 10, image.height() / 2);
    fields.specializedIrisValley = Shape2(image.width() / 2 + 10, image.height() / 2 + 15);//localizeByFollowingLocalMaxima(valleyField, nextToIris);
    return fields;
};

void Deformation::drawFunctions(MultiArray<2, int> &distanceToCenter)
{
    std::vector<MultiArray<2,float>> functions = getFit(distanceToCenter);
    MultiArray<2,float> resEdge(distanceToCenter.shape());
    MultiArray<2,float> resValley(distanceToCenter.shape());
    MultiArray<2,float> resBoth(distanceToCenter.shape());
    for (int i = 0; i<distanceToCenter.width() / 2;i++)
    {
        int yEdge = functions[0][i] > distanceToCenter.height() ? distanceToCenter.height() -1 : (distanceToCenter.height() - (functions[0][i] -1));
        int yValley = functions[1][i] > distanceToCenter.height() ? distanceToCenter.height() -1 : (distanceToCenter.height() - (functions[1][i] -1));
        int yBoth = functions[2][i] > distanceToCenter.height() ? distanceToCenter.height() -1 : (distanceToCenter.height() - (functions[2][i] -1));
        resEdge(i, yEdge) = 1;
        resValley(i, yValley) = 1;
        resBoth(i, yBoth) = 1;
    }
    exportImage(resEdge,"./../images/results/edgeFunction.png");
    exportImage(resValley,"./../images/results/valleyFunction.png");
    exportImage(resBoth ,"./../images/results/bothFunction.png");
    std::cout << "\n expected Radius: ";
    int rad = argMax(functions[2]);
    std::cout << rad;
    functions[2][rad] = 0;
    rad = argMax(functions[2]);
    std::cout << "\n next best Radius: ";
    std::cout << rad;
    functions[2][rad] = 0;
    rad = argMax(functions[2]);
    std::cout << "\n next best Radius: ";
    std::cout << rad;
    functions[2][rad] = 0;
    rad = argMax(functions[2]);
    std::cout << "\n next best Radius: ";
    std::cout << rad;
    functions[2][rad] = 0;
    rad = argMax(functions[2]);
    std::cout << "\n";
    
}

void fitPSF(DataParams &params, MultiArray<2, double> &ps) {
    double minval=9999999, maxval = 0;
    int size = 2*ps.shape(0)*ps.shape(1);
    std::vector<double> data2(2*ps.shape(0)*ps.shape(1));
    for(int i=0, counter = 0;i<ps.shape(0);++i) {
        for(int j=0;j<ps.shape(1);++j,++counter) {
            //std::cout<<i<<" "<<j<<" "<<counter<<std::endl;
            data2[2*counter] = std::sqrt(std::abs(std::pow(i-ps.shape(0)/2,2)+std::pow(j-ps.shape(1)/2,2)));
            data2[2*counter+1] = ps(i,j);
            if (ps(i,j)< minval){minval = ps(i,j);}
            if (ps(i,j)> maxval){maxval = ps(i,j);}
        }
    }
    double sigma = 2.0, scale = maxval - minval, offset = minval;
    //std::cout<<sigma<<" "<<scale<<" "<<offset<<std::endl;
    fitGaussian(&data2[0], size/2, sigma, scale, offset);
    params.setSigma(sigma);

}

double fitPSF2D(MultiArray<2, double> &ps, double &sigma) {
    double minval=9999999, maxval = 0;
	int w = ps.shape(0), h = ps.shape(1);
    int size = ps.shape(0)*ps.shape(1);
    std::vector<double> data2(3*ps.shape(0)*ps.shape(1));
    //std::ofstream outputRoi;
	//outputRoi.open("c:\\tmp\\outputRoi.txt");
	for(int i=0, counter = 0;i<ps.shape(0);++i) {
        for(int j=0;j<ps.shape(1);++j,++counter) {
            //std::cout<<i<<" "<<j<<" "<<counter<<std::endl;
		//	outputRoi <<i<<" "<<j<<" "<<ps(i,j)<<std::endl;
            data2[3*counter] = i-ps.shape(0)/2;
			data2[3*counter+1] = j-ps.shape(1)/2;
            data2[3*counter+2] = ps(i,j);
            if (ps(i,j)< minval){minval = ps(i,j);}
            if (ps(i,j)> maxval){maxval = ps(i,j);}
        }
    }
	//outputRoi.close();
    double scale = maxval - minval, offset = minval, x0=0, y0=0;
	sigma = 2.0;
    //std::cout<<"sigma: "<<sigma<<" scale: "<<scale<<" offset:"<<offset<<std::endl;
    fitGaussian2D(&data2[0], size, sigma, scale, offset, x0, y0);

	double SSE = 0, SSD = 0, sumY = 0, sumY2 = 0;
	for(int i= w/2-2, counter = 0;i<w/2+2;++i) {
        for(int j=h/2-2;j<h/2+2;++j,++counter) {
			double xs = sq((i-ps.shape(0)/2 - x0) / sigma)+ sq((j-ps.shape(1)/2 - y0) / sigma);
            double e = std::exp(-0.5 * xs);
            double r = ps(i,j) - (scale * e + offset);
			SSE += r*r;
			sumY+= ps(i,j);
			sumY2+= sq(ps(i,j));
        }
    }
	SSD = sumY2 - sq(sumY) / size;
	double error = 1-SSE/SSD;
	sigma = std::abs(sigma);

	//std::cout<<"sigma: "<<sigma<<" scale: "<<scale<<" offset:"<<offset<< " x0: " << x0<<" y0: "<<y0<<" Error: "<<error<<std::endl;
	//std::cin.get();
	
	return error;
}