[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