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

Paul McIntosh paul.mcintosh at internetscooter.com
Mon Oct 29 01:47:49 EDT 2012


Ok - almost there - the following will no doubt help someone in future - I
have attached an adapted picture from
http://vis.computer.org/vis2004/dvd/tutorial/tut_5/notes_1.pdf

 

Correct me if I am wrong but the concept of Parallel Projection means that
the scale is like making the size of the camera bigger, there is no
perspective, so if you move the camera 1000 metres back, the picture will be
exactly the same (assuming the object is still in the clipping plane). If
you want to fit something that is "yay big" in the picture, you need to make
sure that the camera is "yay big" also. So if I have a sphere with radius
0.05, I set the parallel scale to 0.05 and the camera will reach that height
above and below (VTK code snippet below)... so far I think I am right J

 

//    // Description:

//    // Set/Get the scaling used for a parallel projection, i.e. the height

//    // of the viewport in world-coordinate distances. The default is 1.

//    // Note that the "scale" parameter works as an "inverse scale" ---

//    // larger numbers produce smaller images.

//    // This method has no effect in perspective projection mode.

//    void SetParallelScale(double scale);

//    vtkGetMacro(ParallelScale,double);

 

 

//    How camera is configured from:

//    D:\vtk-5.8.0\Rendering\vtkCamera.cxx

//    if ( this->ParallelProjection)

//      {

//      // set up a rectangular parallelipiped

//      double width = this->ParallelScale * aspect;

//      double height = this->ParallelScale;

//      double xmin = ( this->WindowCenter[0] - 1.0 ) * width;

//      double xmax = ( this->WindowCenter[0] + 1.0 ) * width;

//      double ymin = ( this->WindowCenter[1] - 1.0 ) * height;

//      double ymax = ( this->WindowCenter[1] + 1.0 ) * height;

//      this->ProjectionTransform->Ortho( xmin, xmax, ymin, ymax,

//                                        this->ClippingRange[0],

//                                        this->ClippingRange[1] );

//      }

 

Paul

 

From: Paul McIntosh [mailto:paul.mcintosh at internetscooter.com] 
Sent: Friday, 26 October 2012 1:05 PM
To: 'vtkusers at vtk.org'
Subject: RE: [vtkusers] Projected area calculation - i.e. frontal area

 

Ok - I have all the bits in place now load STL -> render STL -> capture
screenshot -> count black pixels. Now I just have to make the count mean
something useful ;)

 

Has anyone got some suggestions of how to make, for example, 10000 pixels =
1 metre.  I'll keep playing and will eventually nut it out but any
suggestions to speed up this process is most welcome.

 

Below is my test data dimensions and below that is the camera properties
that produce the attached image. 

 

Reading: D:\scooter\vespalabs\3D_Scan_small\201205_attempt\Horn_Cover.stl

xmin: -1.41477 xmax: 216.216

ymin: -72.8783 ymax: 62.79

zmin: 234.3 zmax: 632.477

position x:107.401  y:-5.04419 z:433.389

 

    // Camera

    vtkSmartPointer<vtkCamera> camera = vtkSmartPointer<vtkCamera>::New();

    camera->SetFocalPoint(centre);

    camera->SetPosition(centre[0],centre[1],centre[2]+500);

    camera->SetParallelProjection(1);

    camera->SetParallelScale(100);

 

 

Cheers,

 

Paul

 

From: vtkusers-bounces at vtk.org [mailto:vtkusers-bounces at vtk.org] On Behalf
Of Paul McIntosh
Sent: Thursday, 25 October 2012 8:08 AM
To: 'Goodwin Lawlor'
Cc: vtkusers at vtk.org
Subject: Re: [vtkusers] Projected area calculation - i.e. frontal area

 

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/20121029/5dbee471/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: parallelprojection.png
Type: image/png
Size: 66798 bytes
Desc: not available
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20121029/5dbee471/attachment.png>


More information about the vtkusers mailing list