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

[Goodwin Lawlor]

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
> 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/20121024/ffc95954/attachment.htm>