Commit ccd60e52 authored by Aron Helser's avatar Aron Helser Committed by Kitware Robot

Merge topic 'moments_remote_module'

8776e904 MomentInvariants final remote repo location and hash
0aeb5bab Remote: parallel code included in single module
aa8b5698 Remote: change ParallelMoments to use branch
3249d4e8 Remote Modules, add MomentInvariants, and ParallelMI
700215a8 Remove MomentInvariants filter, prep to make remote
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Merge-request: !4301
parents cca3d092 8776e904
set(Module_SRCS
vtkMomentsTensor.h
vtkComputeMoments.cxx
vtkMomentInvariants.cxx
vtkMomentsHelper.cxx
vtkSimilarityBalls.cxx
vtkReconstructFromMoments.cxx
)
vtk_module_library(vtkFiltersMomentInvariants ${Module_SRCS})
vtk_add_test_python(
NO_VALID
patternDetectionTestSimple.py
)
"""
we cutout a random point from the given dataset,
rotate it, scale it, multiply with a constant, add a constant
and look for it.
The test is considered successful if the original position is the global maximum of similarity.
Inaccuracies come from the resolution of the dataset and the pattern as well as the coarse resolutions.
I chose the numbers to perfectly match and the rotation angle to be a multiple of 90 degree. So that the sampling error is suppressed.
"""
import vtk
try:
import numpy as np
except ImportError:
print("Numpy (http://numpy.scipy.org) not found.")
print("This test requires numpy!")
from vtk.test import Testing
Testing.skip()
import sys
import os, shutil
import math
from patternDetectionHelper import *
import random
from vtk.util.misc import vtkGetDataRoot
VTK_DATA_ROOT = vtkGetDataRoot()
numberOfIntegrationSteps = 0
order = 2
angleResolution = 8
eps = 1e-2
numberOfMoments = 5 #4*n+5
nameOfPointData = "values"
def doTest(filePath):
if not os.path.isfile(filePath):
print("File: " + filePath + " does not exist.")
sys.exit(1)
else:
(dataset, bounds, dimension) = readDataset(filePath)
npBounds = np.array(bounds)
if (dataset.GetDimensions()[0]-1) % (numberOfMoments-1) != 0:
print("Warning! The numberOfMoments will not produce a true coarse version of the input dataset.")
# produce coarse grid where the moments will be computed
if dimension == 3:
momentDataCoarse = createCoarseDataset(bounds, numberOfMoments, numberOfMoments, numberOfMoments)
else:
momentDataCoarse = createCoarseDataset(bounds, numberOfMoments, numberOfMoments, 0)
# produce pattern
if dimension == 3:
radiusPattern = np.amin(npBounds[1::2]-npBounds[0::2]) / 10
offset = max(radiusPattern, momentDataCoarse.GetSpacing()[0], momentDataCoarse.GetSpacing()[1], momentDataCoarse.GetSpacing()[2])
patternPosition = [random.uniform(bounds[2*i]+offset, bounds[2*i+1]-offset) for i in range(3)]
else:
radiusPattern = np.amin((npBounds[1::2]-npBounds[0::2])[:-1]) / 10
offset = max(radiusPattern, momentDataCoarse.GetSpacing()[0], momentDataCoarse.GetSpacing()[1])
patternPosition = [random.uniform(bounds[2*i]+offset, bounds[2*i+1]-offset) for i in range(3)]
patternPosition[2] = bounds[5]
patternPosition = momentDataCoarse.GetPoint(momentDataCoarse.FindPoint(patternPosition))
# rotate, multiply constant, change size, add constant
rotationAngle = 2 * math.pi * 0.25 * random.randint(0,3)
factorOuter = random.uniform(0.5,2)
factorInner = random.uniform(0.5,2)
summand = [random.uniform(-1,1) for i in range(9)]
# rotationAngle = 0
# factorInner= 1
# factorOuter= 1
# summand = [0 for i in range(9)]
if dimension == 2:
summand[2] = summand[5] = summand[6] = summand[7] = summand[8] = 0
pattern = scaleDataset(shiftDataset(rotateDatasetExact(cutoutPattern(dataset, dimension, patternPosition, radiusPattern), rotationAngle, nameOfPointData),summand, nameOfPointData),factorOuter, nameOfPointData)
pattern.SetSpacing(pattern.GetSpacing()[0]*factorInner, pattern.GetSpacing()[1]*factorInner, pattern.GetSpacing()[2]*factorInner)
print("patternPosition=", patternPosition, " rotationAngle=", rotationAngle, " factorInner=", factorInner, " factorOuter=", factorOuter, " summand=", summand)
# compute moments
momentsAlgo = vtk.vtkComputeMoments()
momentsAlgo.SetFieldData(dataset)
momentsAlgo.SetGridData(momentDataCoarse)
momentsAlgo.SetNameOfPointData(nameOfPointData)
momentsAlgo.SetNumberOfIntegrationSteps(numberOfIntegrationSteps)
momentsAlgo.SetOrder(order)
momentsAlgo.SetRadiiArray([radiusPattern,0,0,0,0,0,0,0,0,0])
momentsAlgo.Update()
# pattern detetction
invariantsAlgo = vtk.vtkMomentInvariants()
invariantsAlgo.SetMomentData(momentsAlgo.GetOutput())
invariantsAlgo.SetPatternData(pattern)
invariantsAlgo.SetNameOfPointData(nameOfPointData)
invariantsAlgo.SetIsScaling(1)
invariantsAlgo.SetIsTranslation(1)
invariantsAlgo.SetIsRotation(1)
invariantsAlgo.SetIsReflection(1)
invariantsAlgo.SetNumberOfIntegrationSteps(numberOfIntegrationSteps)
invariantsAlgo.SetAngleResolution(angleResolution)
invariantsAlgo.SetEps(eps)
invariantsAlgo.Update()
# print invariantsAlgo
# detetction of the local maxima and provide visualization with solid balls and hollow spheres
ballsAlgo = vtk.vtkSimilarityBalls()
ballsAlgo.SetSimilarityData(invariantsAlgo.GetOutput())
ballsAlgo.SetGridData(dataset)
# ballsAlgo.SetKindOfMaxima(2)
ballsAlgo.Update()
print("maxSimilarity =", ballsAlgo.GetOutput(0).GetPointData().GetArray("localMaxValue").GetRange()[1])
# print ballsAlgo.GetOutput().FindPoint(patternPosition)
if ballsAlgo.GetOutput(0).GetPointData().GetArray("localMaxValue").GetTuple1( ballsAlgo.GetOutput().FindPoint(patternPosition)) == ballsAlgo.GetOutput(0).GetPointData().GetArray("localMaxValue").GetRange()[1]:
return
else:
sys.exit(1)
files = ["2DScalar.vti", "2DMatrix.vti", "2DVector.vti", "3DScalar.vti"]
for f in files:
filePath = VTK_DATA_ROOT + "/Data/" + f
doTest(filePath)
sys.exit(0)
vtk_module(vtkFiltersMomentInvariants
TEST_DEPENDS
vtkIOXML
DEPENDS
vtkFiltersCore
vtkImagingCore
PRIVATE_DEPENDS
vtkeigen
vtkkissfft
)
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
set(Module_SRCS
vtkPComputeMoments.cxx
)
vtk_module_library(vtkFiltersParallelMomentInvariants ${Module_SRCS})
# Tests with data
ExternalData_Expand_Arguments(VTKData _
"DATA{../Data/Input/,RECURSE:,REGEX:.*}"
)
# We don't use `vtk_add_test_python_mpi` since we need to customize the
# arguments and number of ranks these tests are run on.
ExternalData_add_test(VTKData
NAME TestParalleComputeShort-p4
COMMAND ${VTK_MPIRUN_EXE} ${VTK_MPI_PRENUMPROC_FLAGS} ${VTK_MPI_NUMPROC_FLAG} 4 ${VTK_MPI_PREFLAGS}
$<TARGET_FILE:pvtkpython>
${CMAKE_CURRENT_SOURCE_DIR}/TestParalleComputeShort.py
DATA{../Data/Input/2DScalar.vti}
DATA{../Data/Input/2DScalar.vtm}
)
ExternalData_add_test(VTKData
NAME TestParalleComputeShort-p8
COMMAND ${VTK_MPIRUN_EXE} ${VTK_MPI_PRENUMPROC_FLAGS} ${VTK_MPI_NUMPROC_FLAG} 8 ${VTK_MPI_PREFLAGS}
$<TARGET_FILE:pvtkpython>
${CMAKE_CURRENT_SOURCE_DIR}/TestParalleComputeShort.py
DATA{../Data/Input/3DScalar.vti}
DATA{../Data/Input/3DScalar.vtm}
)
"""
we cutout a random point from the given dataset,
rotate it, scale it, multiply with a constant, add a constant
and look for it.
The test is considered successful if the original position is the global maximum of similarity.
Inaccuracies come from the resolution of the dataset and the pattern as well as the coarse resolutions.
I chose the numbers to perfectly match and the rotation angle to be a multiple of 90 degree. So that the sampling error is suppressed.
"""
import vtk
import numpy as np
import sys
import os, shutil
import math
from patternDetectionHelper import *
numberOfIntegrationSteps = 0
order = 2
numberOfMoments = 0
nameOfPointData = "values"
radius = 0.1
def compute(whole, split):
# prepare
controller = vtk.vtkMPIController()
vtk.vtkMultiProcessController.SetGlobalController(controller)
rank = controller.GetLocalProcessId()
nprocs = controller.GetNumberOfProcesses()
# print('rank={0}'.format(rank))
# print('nprocs={0}'.format(nprocs))
# whole
(wholeDataset, bounds, dimension) = readDataset(whole)
if pow(2,dimension) != nprocs:
print 'the number of procs must be 2^dimension'
sys.exit()
# coarse
coarseData = vtk.vtkImageData()
coarseData.CopyStructure(wholeDataset)
coarseData.SetSpacing(0.1,0.1,0.1)
if dimension == 2:
coarseData.SetDimensions(11,11,1)
else:
coarseData.SetDimensions(11,11,11)
# print data
# print coarseData
# compute moments
momentsFilter = vtk.vtkComputeMoments()
momentsFilter.SetFieldData(wholeDataset)
momentsFilter.SetGridData(coarseData)
momentsFilter.SetNameOfPointData(nameOfPointData)
momentsFilter.SetUseFFT(0)
momentsFilter.SetNumberOfIntegrationSteps(numberOfIntegrationSteps)
momentsFilter.SetOrder(order)
momentsFilter.SetRadiiArray([radius,0,0,0,0,0,0,0,0,0])
momentsFilter.Update()
# if rank ==0: print momentsFilter
# this is to make sure each proc only loads their piece
ids_to_read = [rank]
reqs = vtk.vtkInformation()
reqs.Set(vtk.vtkCompositeDataPipeline.UPDATE_COMPOSITE_INDICES(), ids_to_read, 1)
reader = vtk.vtkXMLMultiBlockDataReader()
reader.SetFileName(split)
reader.Update(reqs)
# print('reader={0}'.format(reader.GetOutput()))
multiblock = vtk.vtkMultiBlockDataSet.SafeDownCast(reader.GetOutput())
# print('multiblock={0}'.format(multiblock))
multipiece = vtk.vtkMultiPieceDataSet.SafeDownCast(multiblock.GetBlock(0))
# print('multipiece={0}'.format(multipiece))
data = vtk.vtkImageData.SafeDownCast(multipiece.GetPiece(rank))
# print('data={0}'.format(data))
# coarse
coarseData = vtk.vtkImageData()
coarseData.CopyStructure(data)
coarseData.SetSpacing(0.1,0.1,0.1)
if dimension == 2:
coarseData.SetDimensions(6,6,1)
else:
coarseData.SetDimensions(6,6,6)
# print data
# coarseData
# parallel
pMomentsFilter = vtk.vtkPComputeMoments()
pMomentsFilter.SetFieldData(data)
pMomentsFilter.SetGridData(coarseData)
pMomentsFilter.SetNameOfPointData(nameOfPointData)
pMomentsFilter.SetNumberOfIntegrationSteps(numberOfIntegrationSteps)
pMomentsFilter.SetRadiiArray([radius, 0, 0, 0, 0, 0, 0, 0, 0, 0 ])
pMomentsFilter.Update(reqs)
# if rank ==0: print pMomentsFilter
# check for difference except for the global boundary (the probe and parallel probe filters behave differently there. So, there is a difference if numberOfIntegrationSteps > 0 )
diff = 0
for i in xrange(pMomentsFilter.GetOutput().GetPointData().GetNumberOfArrays()):
diffArray = vtk.vtkDoubleArray()
diffArray.SetName( pMomentsFilter.GetOutput().GetPointData().GetArrayName(i) )
diffArray.SetNumberOfComponents( 1 )
diffArray.SetNumberOfTuples( pMomentsFilter.GetOutput().GetNumberOfPoints() )
diffArray.Fill( 0.0 )
momentsArray = momentsFilter.GetOutput().GetPointData().GetArray(i)
pMomentsArray = pMomentsFilter.GetOutput().GetPointData().GetArray(i)
for j in xrange(pMomentsFilter.GetOutput().GetNumberOfPoints()):
diffArray.SetTuple1(j, abs(momentsArray.GetTuple1(momentsFilter.GetOutput().FindPoint(pMomentsFilter.GetOutput().GetPoint(j))) - pMomentsArray.GetTuple1(j)))
if (dimension == 2 and pMomentsFilter.GetOutput().GetPoint(j)[0] > 0 and pMomentsFilter.GetOutput().GetPoint(j)[0] < 1 and pMomentsFilter.GetOutput().GetPoint(j)[1] > 0 and pMomentsFilter.GetOutput().GetPoint(j)[1] < 1) or (dimension == 3 and pMomentsFilter.GetOutput().GetPoint(j)[0] > 0 and pMomentsFilter.GetOutput().GetPoint(j)[0] < 1 and pMomentsFilter.GetOutput().GetPoint(j)[1] > 0 and pMomentsFilter.GetOutput().GetPoint(j)[1] < 1 and pMomentsFilter.GetOutput().GetPoint(j)[2] > 0 and pMomentsFilter.GetOutput().GetPoint(j)[2] < 1):
# if diffArray.GetTuple1(j) > 1e-10:
# print rank, i, j, pMomentsFilter.GetOutput().GetPoint(j), diffArray.GetTuple(j)
diff = max(diff, diffArray.GetTuple1(j))
if diff > 1e-10:
print "test failed, maxdiff =", diff
else:
print "test successful"
if rank == 0:
print 'all done!'
if __name__ == '__main__':
if len(sys.argv) < 2:
print 'usage: <file whole> <file split>'
else:
compute(sys.argv[1], sys.argv[2])
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment