C++ McDArray::size方法代码示例

void HxCPDSpatialGraphWarp::resamplePairs(McDArray<McVec3f>& p1,
                                          McDArray<McVec3f>& p2) {
    int numPairs = p1.size();
    for (int i = p1.size() - 1; i > -1; i--) {
        bool resetI = false;
        for (int j = i - 1; j > -1; j--) {
            // std::cout<<"\nCompare "<<i<<" and "<<j;
            McVec3f set1Coord = p1[i];
            McVec3f set2Coord = p1[j];
            float dist = (set1Coord - set2Coord).length();
            if (dist < portSampleDist.getValue()) {
                std::cout << "\ndist between " << i << " and " << j << " is "
                          << dist;
            }
            if (dist < portSampleDist.getValue()) {
                p1.remove(j, 1);
                p2.remove(j, 1);
                resetI = true;
            }
        }
        if (resetI) {
            i = p1.size() - 1;
        }
    }
    std::cout << "\n" << p1.size() << " point from " << numPairs
              << " left after resampling.";
}

void BruteForceOptMatching::getMaxProbAssignments(
    const BP& ia, const FactorGraph& fg, const ConnectedFactorGraph& graph,
    McDArray<McVec2i>& pairs) {

    for (int i = 0; i < graph.variables.size(); i++) {

        McDArray<int> possibleAssignments;
        getAssignmentsForVariable(graph.variables[i], possibleAssignments);
        Factor belief =
            ia.belief(Var(graph.variables[i], possibleAssignments.size() + 1));
        float maxVal = -1 * FLT_MAX;
        int maxIdx = -1;

        for (int j = 0; j < possibleAssignments.size() + 1; j++) {
            if (belief.get(j) > maxVal) {
                maxVal = belief.get(j);
                maxIdx = j;
            }
        }
        int indexOfAssignmentInVertexList =
            mapVariableAssignmentToIndexInVertexList(graph.variables[i],
                                                     maxIdx);

        McVec2i pair =
            McVec2i(graph.variables[i], indexOfAssignmentInVertexList);

        pairs.append(pair);
    }
    outputSingleFactorValues(graph);

    //    std::vector<std::size_t> maxes= ia.findMaximum();
    //    vector<std::size_t>::iterator it=maxes.begin();
}

void HxIteratePointMatchingUntilConvergence::addEvidenceForNodes(
    const McDArray<int>& nodesToAssign) {
    HxSpatialGraph* graph = hxconnection_cast<HxSpatialGraph>(portData);
    const EdgeVertexAttribute* evidenceAttrib =
        dynamic_cast<const EdgeVertexAttribute*>(graph->findAttribute(
            HxSpatialGraph::VERTEX, "UserDefinedMatchings"));
    McDArray<int> pairNodes(nodesToAssign.size());
    pairNodes.fill(-1);

    for (int j = 0; j < nodesToAssign.size(); j++) {
        int nodeJLabel = evidenceAttrib->getIntDataAtIdx(nodesToAssign[j]);
        if (nodeJLabel == 0) {
            pairNodes[j] = -1;
            continue;
        }
        for (int i = 0; i < graph->getNumVertices(); i++) {
            if (i != nodesToAssign[j]) {
                int nodeILabel = evidenceAttrib->getIntDataAtIdx(i);
                if (nodeILabel == nodeJLabel) {
                    pairNodes[j] = i;
                }
            }
        }
    }
    for (int j = 0; j < nodesToAssign.size(); j++) {
        addEvidence(nodesToAssign[j], pairNodes[j]);
    }
}

static McDArray<McVec2d> asVec2dArray(const McDArray<McVec3f>& a) {
    McDArray<McVec2d> b;
    b.resize(a.size());
    for (long i = 0; i < a.size(); i++) {
        b[i].x = a[i].x;
        b[i].y = a[i].y;
    }
    return b;
}

void BruteForceOptMatching::computeAngleProbs(McDArray<float>& angles) {
    float sumAngles = 0;
    for (int i = 0; i < angles.size(); i++) {
        angles[i] = 180.0 - angles[i];
        angles[i] = exp((float)(-angles[i] / 10.0));
        sumAngles += angles[i];
    }
    for (int i = 0; i < angles.size(); i++)
        angles[i] /= sumAngles;
}

void BruteForceOptMatching::computeProjDistProbs(McDArray<float>& dists) {
    float sumDists = 0.0;
    for (int i = 0; i < dists.size(); i++) {
        if (dists[i] != FLT_MAX) {
            dists[i] = exp((float)(-dists[i] / 100.0));
        } else
            dists[i] = 0;
        sumDists += dists[i];
    }
    for (int i = 0; i < dists.size(); i++)
        dists[i] /= sumDists;
}

// creates all singleton factors for one connected component
// does not set the values
void BruteForceOptMatching::createSingletonFactors(
    const McDArray<int>& connectedComp, vector<Factor>& singletonFactorList) {

    for (int i = 0; i < connectedComp.size(); i++) {
        McDArray<int> assigmentsForIthVariable;
        int curVar = connectedComp[i];
        getAssignmentsForVariable(curVar, assigmentsForIthVariable);

        Var pi(curVar, assigmentsForIthVariable.size() + 1);
        Factor faci(pi);
        singletonFactorList.push_back(faci);
    }
}

// Computes all variables that form a connected component in the
// adjacenceMatrix.
// The connected component chosen is arbitrary - it takes the first it finds.
bool BruteForceOptMatching::getConnectedComponent(
    const McDMatrix<int>& adjacenceMatrix,
    McDMatrix<int>& adjacenceMatrixWithoutConnctedComponent,
    McDArray<int>& connComp) {
    adjacenceMatrixWithoutConnctedComponent.resize(adjacenceMatrix.nRows(),
                                                   adjacenceMatrix.nCols());
    memcpy(adjacenceMatrixWithoutConnctedComponent.dataPtr(),
           adjacenceMatrix.dataPtr(),
           sizeof(int) * adjacenceMatrix.nRows() * adjacenceMatrix.nCols());

    // find first a startpoint
    int start = -1;
    connComp.resize(0);
    for (int i = 0; i < adjacenceMatrix.nRows(); i++) {
        for (int j = i; j < adjacenceMatrix.nCols(); j++) {
            if (adjacenceMatrix[i][j] == 1) {
                start = i;
                break;
            }
        }
    }
    if (start == -1)
        return false;

    McDArray<int> queue;
    queue.append(start);
    connComp.clear();
    while (queue.size() > 0) {
        int cur = queue.last();
        connComp.append(cur);
        queue.pop_back();
        for (int i = 0; i < adjacenceMatrixWithoutConnctedComponent.nCols();
             i++) {
            if (adjacenceMatrixWithoutConnctedComponent[cur][i] == 1) {
                queue.push(i);
                adjacenceMatrixWithoutConnctedComponent[cur][i] = 0;
            }
        }
    }
    // remove duplicates
    connComp.sort(&mcStandardCompare);
    int cur = connComp.last();
    for (int i = connComp.size() - 2; i >= 0; i--) {
        if (cur == connComp[i])
            connComp.remove(i, 1);
        else
            cur = connComp[i];
    }
    return true;
}

void BruteForceOptMatching::compute3dDistProbs(McDArray<float>& dists) {
    float sumDists = 0.0;
    // This should be a user parameter...
    float probParam = 500.0;
    for (int i = 0; i < dists.size(); i++) {
        if (dists[i] != FLT_MAX) {
            dists[i] = exp(-dists[i] / probParam);
        } else
            dists[i] = 0;
        sumDists += dists[i];
    }
    for (int i = 0; i < dists.size(); i++)
        dists[i] /= sumDists;
}

void BruteForceOptMatching::getSingletonProbs(
    const McDMatrix<float>& angleMatrix,
    const McDMatrix<float>& projDistanceMatrix,
    const McDMatrix<float>& distanceMatrix3d,
    const McDMatrix<int>& variableAssignmentMat,
    const McDArray<int>& assignments, const int varLabel,
    McDArray<double>& probs) {
    McDArray<float> floatProbs;
    getSingletonProbs(angleMatrix, projDistanceMatrix, distanceMatrix3d,
                      variableAssignmentMat, assignments, varLabel, floatProbs);

    probs.resize(floatProbs.size());
    for (int i = 0; i < floatProbs.size(); ++i)
        probs[i] = floatProbs[i];
}

void BruteForceOptMatching::projectToPlaneApproxDirection(
    const McDArray<McVec3f>& vertices, const McDArray<McVec3f>& directions,
    const float planeZ, McDArray<McVec3f>& result) {

    for (int i = 0; i < vertices.size(); i++) {
        McPlane theZPlane(McVec3f(0, 0, 1), planeZ);
        McVec3f vertexCoord = vertices[i];
        McVec3f dir = directions[i] * -1;
        dir.normalize();
        float angle = dir.angle(McVec3f(0, 0, 1));
        if (angle > M_PI / 2.0)
            angle = M_PI - angle;

        McLine theLine(vertexCoord, vertexCoord + directions[i]);

        McVec3f intersectionPoint;
        bool intersected = theZPlane.intersect(theLine, intersectionPoint);
        // if(fabs(angle)<0.1)
        //   cout<<"\n Angle for vertex "<<i <<" too low: "<< angle;
        if (intersected && (fabs(angle) < (M_PI / 2.0 - M_PI / 8.0))) {
            result.append(intersectionPoint);
        } else {
            result.append(McVec3f(vertexCoord.x, vertexCoord.y, planeZ));
        }
    }
}

float BruteForceOptMatching::getMedianZ(const McDArray<McVec3f>& vertices) {

    if (!vertices.size())
        return -1 * FLT_MAX;
    McDArray<float> zs;
    float mean = 0.0;

    for (int i = 0; i < vertices.size(); i++) {
        zs.append(vertices[i].z);
        mean += vertices[i].z;
    }
    zs.sort(&mcStandardCompare);
    int medianIdx = zs.size() / 2.0;
    // cout <<"MeanZ: "<<mean/vertices.size();
    // return mean/vertices.size();
    return zs[medianIdx];
}

void HxCPDSpatialGraphWarp::preparePoints(McDArray<McVec3f>& p1,
                                          McDArray<McVec3f>& p2,
                                          SpatialGraphSelection& slice2,
                                          const HxSpatialGraph* spatialGraph) {
    ma::SliceSelector selectionHelper(spatialGraph, "TransformInfo");

    ma::EndPointParams params;
    params.endPointRegion = 30;
    params.projectionPlane = selectionHelper.computeMidPlane(0, 1);
    params.projectionType = ma::P_ORTHOGONAL;
    params.refSliceNum = 0;
    params.transSliceNum = 1;
    params.useAbsoluteValueForEndPointRegion = false;
    params.maxDistForAngle = 2000;
    params.angleToPlaneFilter = 0.01;
    SpatialGraphSelection slice1;
    ma::FacingPointSets pr =
        ma::projectEndPoints(spatialGraph, slice1, slice2, params);
    McDArray<McVec3f> refCoords = pr.ref.positions;
    McDArray<McVec3f> transCoords = pr.trans.positions;
    mcassert(refCoords.size() == slice1.getNumSelectedVertices());
    mcassert(transCoords.size() == slice2.getNumSelectedVertices());

    p1.resize(refCoords.size());
    for (int i = 0; i < refCoords.size(); i++) {
        McVec3f coord = refCoords[i];
        p1[i] = McVec3f(coord.x, coord.y, 0);
    }
    p2.resize(transCoords.size());
    for (int i = 0; i < transCoords.size(); i++) {
        McVec3f coord = transCoords[i];
        p2[i] = McVec3f(coord.x, coord.y, 0);
    }
    mcassert(p2.size() == slice2.getNumSelectedVertices());
}

void BruteForceOptMatching::checkAmbiguities(const BP& ia,
                                             const FactorGraph& fg,
                                             const ConnectedFactorGraph& graph,
                                             McDArray<int>& ambiguities) {
    for (int h = 0; h < graph.variables.size(); h++) {
        McDArray<int> possibleAssignments;
        getAssignmentsForVariable(graph.variables[h], possibleAssignments);
        Factor belief =
            ia.belief(Var(graph.variables[h], possibleAssignments.size() + 1));

        float maxProb = belief.max();
        int countSame = 0;

        for (int k = 0; k < possibleAssignments.size() + 1; k++) {
            float curProb = belief.get(k);
            if (fabs(curProb - maxProb) < 0.1)
                countSame++;
        }

        /////
        cout << "\n Belief for var " << graph.variables[h] << "\n";
        for (int k = 0; k < possibleAssignments.size() + 1; k++) {
            float curProb = belief.get(k);
            cout << curProb << " ";
        }
        cout << "\n";

        ////

        if (countSame > 1) {
            // oh no! We found an ambiguos assignment!

            ambiguities.append(graph.variables[h]);

            // print it out:
            cout << "Found an ambiguous assignemnt to variable "
                 << graph.variables[h] << "\n";
            for (int k = 0; k < possibleAssignments.size() + 1; k++) {
                float curProb = belief.get(k);
                cout << curProb << " ";
            }
            cout << "\n";
        }
    }
}

float AssignEModulus::averageVoxels( const OBLelement3D* FEelement,
									const McDArray<McVec3f>& vertices,
									HxUniformScalarField3* field )
{
	if ( !FEelement ) { return 0; }

	McVec3f pcoords;

	// Create an instance of HxLoc3Regular instead of HxLocation3,
	// because we know we are dealing with a uniform scalar field and
	// we want access to the voxelindices
	//HxLoc3Regular* location = image->coords()->createLocation();
	HxLoc3Regular* location = (HxLoc3Regular*)field->createLocation();

	// determine bounding box of voxel indices within element
	location->set( vertices[0][0], vertices[0][1], vertices[0][2] );
	int indexbbox[6] = { location->ix, location->ix,
						location->iy, location->iy,
						location->iz, location->iz };
	for ( int el = 1; el < vertices.size(); el++ )
	{
		location->move( vertices[el][0], vertices[el][1], vertices[el][2] );
		if ( location->ix < indexbbox[0] ) { indexbbox[0] = location->ix; }
		if ( location->ix > indexbbox[1] ) { indexbbox[1] = location->ix; }
		if ( location->iy < indexbbox[2] ) { indexbbox[2] = location->iy; }
		if ( location->iy > indexbbox[3] ) { indexbbox[3] = location->iy; }
		if ( location->iz < indexbbox[4] ) { indexbbox[4] = location->iz; }
		if ( location->iz > indexbbox[5] ) { indexbbox[5] = location->iz; }
	}

	float sum = 0.0;
	int count = 0;
	for ( int k = indexbbox[4]; k <= indexbbox[5]; k++ )
	{
		for ( int j = indexbbox[2]; j <= indexbbox[3]; j++ )
		{
			for ( int i = indexbbox[0]; i <= indexbbox[1]; i++ )
			{
				float point[3];
				field->lattice.coords()->pos( i, j, k, point );
				int inside = FEelement->getIsoParam( vertices,
					McVec3f(point[0],point[1], point[2]), pcoords );
				if ( inside == 1 )
				{
					location->move( point );
					float value;
					field->eval( location, &value );
					sum += value;
					count++;
				}
			}
		}
	}
	delete location;

	return ( count > 0 ) ? ( sum/count ) : 0.0;
}