surfaceRegistration/preProcessing.cpp at master · Tonsty/surfaceRegistration · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#include "preProcessing.h"

void freeData(double **finalData)
{
    for (int i=0; i<3; i++)
        delete [] finalData[i];
    delete [] finalData;
}

double absoluteDouble (double number) {
    if (number < 0) {
        return -number;
    }
    else return number;
}

vtkSmartPointer<vtkPolyData> TransformScaleTranslate(vtkSmartPointer<vtkPolyData> polydata, double value)
{
    vtkSmartPointer<vtkTransform> customTransform = vtkSmartPointer<vtkTransform>::New();
    customTransform->Scale(value, value, value);

    vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
#if VTK_MAJOR_VERSION <= 5
    transformFilter->SetInput(polydata);
#else
    transformFilter->SetInputData(polydata);
#endif
    transformFilter->SetTransform(customTransform);
    transformFilter->Update();

    return transformFilter->GetOutput();
}

vtkSmartPointer<vtkOctreePointLocator> buildOctree(vtkSmartPointer<vtkPolyData> newdata, vtkSmartPointer<vtkPoints> points,
                                                   vtkSmartPointer<vtkPolyData> olddata) {
    newdata->SetPoints(points);
    newdata->SetPolys(olddata->GetPolys());

    vtkSmartPointer<vtkOctreePointLocator> octree = vtkSmartPointer<vtkOctreePointLocator>::New();
    octree->SetDataSet(newdata);
    octree->BuildLocator();

    return octree;
}

double applyPCA(int pointCountTarget, int pointCountTemplate,
                         vtkSmartPointer<vtkPolyData> dataTarget, vtkSmartPointer<vtkPolyData> dataTemplate,
                         vtkSmartPointer<vtkPolyData> polydataTarget, vtkSmartPointer<vtkPolyData> polydataTemplate, vtkSmartPointer<vtkPolyData> transformedTemplate,
                         int maxDistance) {

    double maxCoordinate = 0;

    double **finalDataTarget;
    finalDataTarget = PCAStatistics(dataTarget, pointCountTarget, 0);

    vtkSmartPointer<vtkPoints> pointsTarget = vtkSmartPointer<vtkPoints>::New();

    for(int i = 0; i < pointCountTarget; ++i)
    {
        pointsTarget->InsertNextPoint ( finalDataTarget[0][i], finalDataTarget[1][i], finalDataTarget[2][i] );
        for (int j=0; j<3; j++) {
            if (absoluteDouble(finalDataTarget[j][i]) > maxCoordinate) {
                maxCoordinate = absoluteDouble(finalDataTarget[j][i]);
            }
        }
    }

    // Building the octree

    vtkSmartPointer<vtkOctreePointLocator> octreeTarget = buildOctree(polydataTarget, pointsTarget, dataTarget);

    // Used to rotate the template after the PCA if needed

    double **finalDataTemplate;
    double distance;
    double overallDistance = 0;
    double minDistance = 1000000;
    int minOrientation = 0;

    vtkSmartPointer<vtkPoints> pointsTemplate = vtkSmartPointer<vtkPoints>::New();

    for (int orientation = 0; orientation < 8; orientation++) {

        finalDataTemplate = PCAStatistics(dataTemplate, pointCountTemplate, orientation);

        for(int i = 0; i < pointCountTemplate; ++i)
        {
            pointsTemplate->InsertNextPoint( finalDataTemplate[0][i], finalDataTemplate[1][i], finalDataTemplate[2][i] );
        }

        for(int i = 0; i < pointCountTemplate; ++i)
        {
            octreeTarget->FindClosestPoint( finalDataTemplate[0][i], finalDataTemplate[1][i], finalDataTemplate[2][i], distance );
            overallDistance += distance;
        }

        if (minDistance > overallDistance) {
            minDistance = overallDistance;
            minOrientation = orientation;
        }

        overallDistance = 0;
        pointsTemplate = vtkSmartPointer<vtkPoints>::New();
    }

    finalDataTemplate = PCAStatistics(dataTemplate, pointCountTemplate, minOrientation);

    // The initialization of the transformed template is the same as the actual template

    vtkSmartPointer<vtkPoints> transformedPointsTemplate = vtkSmartPointer<vtkPoints>::New();

    for(int i = 0; i < pointCountTemplate; ++i)
    {
        pointsTemplate->InsertPoint( i, finalDataTemplate[0][i], finalDataTemplate[1][i], finalDataTemplate[2][i] );
        transformedPointsTemplate->InsertPoint( i, finalDataTemplate[0][i], finalDataTemplate[1][i], finalDataTemplate[2][i] );
        for (int j=0; j<3; j++) {
            if (absoluteDouble(finalDataTemplate[j][i]) > maxCoordinate) {
                maxCoordinate = absoluteDouble(finalDataTemplate[j][i]);
            }
        }
    }

    freeData(finalDataTarget);
    freeData(finalDataTemplate);

    // Initialization of all polydata

    polydataTarget->SetPoints(pointsTarget);
    polydataTarget->SetPolys(dataTarget->GetPolys());

    polydataTemplate->SetPoints(pointsTemplate);
    polydataTemplate->SetPolys(dataTemplate->GetPolys());

    transformedTemplate->SetPoints(transformedPointsTemplate);
    transformedTemplate->SetPolys(dataTemplate->GetPolys());

    return maxCoordinate;
}