[vtkusers] Projected area calculation - i.e. frontal area

[Paul McIntosh]

HI Goodwin,

 

The meshes are all over the place and overlapping - if I could just select
triangles facing forward, then I could transform them flat and do the area
that way. Selecting only triangles facing forward and not obscured by other
triangles is the hard part. Also select visible surfaces doesn't work
because the mesh is much finer than what is being rendered.

 

Cheers,

 

Paul

 

From: Goodwin Lawlor [mailto:goodwin.lawlor.lists at gmail.com] 
Sent: Wednesday, 24 October 2012 10:45 PM
To: Paul McIntosh
Cc: vtkusers at vtk.org
Subject: Re: [vtkusers] Projected area calculation - i.e. frontal area

 

Hi Paul,

 

I may have misunderstood but could you not just project each triangle in the
stl mesh onto the frontal plane (using vtkPlane::ProjectPoint() ), then sum
the projected areas?

 

You may have to prune the mesh first of triangles not "facing" the plane...
you could use vtkVectorDot to do this. You'll also have to deal with
occluded triangles (if you have them)

 

hth,

 

Goodwin

On Wed, Oct 24, 2012 at 2:18 PM, Paul McIntosh
<paul.mcintosh at internetscooter.com> wrote:

Hi All,

I have been working on a utility for calculating frontal area of an STL mesh
for CFD. I thought I had it cracked my scanning a line through the mesh at
evenly spaced distances and then counting the number of hits (based on the
OBBTree example). However on closer inspection the areas I was getting were
always larger than they should be, with reason being the OBBTree method is
not precise and I was getting double/triple hits at edges. I have included
the code below as it would probably help someone solve some other problem ;)

Ok - back to the drawing board... the next obvious way of doing this (and
probably quicker) is to render with parallel projection and then count
pixels in the result.

Is there an easy way to take sceenshot so that pixels are a known size? Any
pointers welcome!

Cheers,

Paul
---
www.internetscooter.com


https://github.com/internetscooter/Vespa-Labs/blob/master/VespaCFD/Calculate
<https://github.com/internetscooter/Vespa-Labs/blob/master/VespaCFD/Calculat
eFrontalArea/CalculateFrontalArea.cxx> 
FrontalArea/CalculateFrontalArea.cxx

#include <vtkPolyData.h>
#include <vtkTriangle.h>
#include <vtkLineSource.h>
#include <vtkSTLReader.h>
#include <vtkSmartPointer.h>
#include <vtkPolyDataMapper.h>
#include <vtkIdTypeArray.h>
#include <vtkSelectionNode.h>
#include <vtkActor.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkSelection.h>
#include <vtkExtractSelection.h>
#include <vtkDataSetMapper.h>
#include <vtkProperty.h>
#include <vtkObjectFactory.h>
#include <vtkCellArray.h>
#include <vtkCell.h>
#include <vtkInformation.h>
#include <vtkExtractSelectedPolyDataIds.h>
#include <vtkOBBTree.h>


// C++
#include <iostream>
using namespace std;

// just to make things nicer later and save some documentation reading ;)
struct boundingBox {
    double xmin;
    double xmax;
    double ymin;
    double ymax;
    double zmin;
    double zmax;
};

// this projects a series of lines through a mesh and looks for
intersections to determine frontal area
// lines are projected based on the bounding box and the resolution is
determined by the steps
// For my application I just need to look at frontal area along the X
direction. The follow code could easily be
// adapted to take a direction as a command line option and produce the area
from that .
// So this currently "scans" yMin->yMax,zMin->ZMax, with a line from Xmin to
Xmax
// ref:
http://www.cmake.org/Wiki/VTK/Examples/Cxx/DataStructures/OBBTree_IntersectW
ithLine
int main(int argc, char *argv[])
{
        // check and get the stl input file provided
    if ( argc != 2 )
    {
        cout << "Required parameters: Filename" << endl;
        return EXIT_FAILURE;
    }
    std::string inputfile = argv[1];

    // later this should be a command line option - hard coded for now
    double resolution = 0.005; // step value while scanning - will need to
do some more calculations later if changed
    //char direction = "X";  // direction of scan

    // read STL and print out some info
    std::cout << "Reading: " << inputfile << std::endl;
    vtkSmartPointer<vtkSTLReader> stlReader =
vtkSmartPointer<vtkSTLReader>::New();
    stlReader->SetFileName(inputfile.c_str());
        vtkSmartPointer<vtkPolyData> polydata =
vtkSmartPointer<vtkPolyData>::New();
        polydata = stlReader->GetOutput();
    polydata->Update();
    //    Debug info if needed:
    cout << "Cells: " << polydata->GetNumberOfCells() << endl;
    cout << "Points: " << polydata->GetNumberOfPoints() << endl;
    cout << "Polys: " << polydata->GetNumberOfPolys() << endl;
    cout << "Verts: " << polydata->GetNumberOfVerts() << endl;
    polydata->Print(cout);

    // makes sense to only scan in an area the object exists, the bounding
box will tell us this
    double bounds[6];
    boundingBox boxBounds;
    polydata->GetBounds(bounds);
    boxBounds.xmin = bounds[0] - resolution*5;
    boxBounds.xmax = bounds[1] + resolution*5;
    boxBounds.ymin = bounds[2] - resolution*5;
    boxBounds.ymax = bounds[3] + resolution*5;
    boxBounds.zmin = bounds[4] - resolution*5;
    boxBounds.zmax = bounds[5] + resolution*5;

    //    Debug info if needed:
    std::cout  << "xmin: " << boxBounds.xmin << " "
               << "xmax: " << boxBounds.xmax << std::endl
               << "ymin: " << boxBounds.ymin << " "
               << "ymax: " << boxBounds.ymax << std::endl
               << "zmin: " << boxBounds.zmin << " "
               << "zmax: " << boxBounds.zmax << std::endl;

    // build OBB Tree locator
    vtkSmartPointer<vtkOBBTree> tree = vtkSmartPointer<vtkOBBTree>::New();
    tree->SetDataSet(polydata);
    tree->BuildLocator();

    // calculate area
    double area = 0;
    for (double zdir = boxBounds.zmin; zdir < boxBounds.zmax; zdir = zdir +
resolution)
    {
        for (double ydir = boxBounds.ymin; ydir < boxBounds.ymax; ydir =
ydir + resolution)
        {
            // Line to intersect with
            double p1[3] = {boxBounds.xmin, ydir + resolution/2, zdir +
resolution/2,};
            double p2[3] = {boxBounds.xmax, ydir + resolution/2, zdir +
resolution/2,};
            //Find intersection points
            vtkSmartPointer<vtkPoints> intersectPoints =
vtkSmartPointer<vtkPoints>::New();
            tree->IntersectWithLine(p1, p2, intersectPoints, NULL);
            // this bit will draw a rough approximation of the shape
            if (intersectPoints->GetNumberOfPoints() > 0)
            {
                cout << "#";
                area++;
            }
            else
            {
                cout << " ";
            }
        }
        cout << std::endl;
    }
    // output how many hits, which for a resolution of 1 should equal
something like 1mm^2
    cout << "area: " << area * resolution * resolution << std::endl;

    return EXIT_SUCCESS;
}



_______________________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Please keep messages on-topic and check the VTK FAQ at:
http://www.vtk.org/Wiki/VTK_FAQ

Follow this link to subscribe/unsubscribe:
http://www.vtk.org/mailman/listinfo/vtkusers

 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20121025/0f2f2a6e/attachment.htm>