6 * $Date: 2012-07-05 16:45:20 +0200 (Do, 05. Jul 2012) $
7 ***************************************************************/
10 * \brief Implementation of class LinearQuadtreeExpansion.
12 * \author Martin Gronemann
15 * This file is part of the Open Graph Drawing Framework (OGDF).
19 * See README.txt in the root directory of the OGDF installation for details.
22 * This program is free software; you can redistribute it and/or
23 * modify it under the terms of the GNU General Public License
24 * Version 2 or 3 as published by the Free Software Foundation;
25 * see the file LICENSE.txt included in the packaging of this file
29 * This program is distributed in the hope that it will be useful,
30 * but WITHOUT ANY WARRANTY; without even the implied warranty of
31 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 * GNU General Public License for more details.
35 * You should have received a copy of the GNU General Public
36 * License along with this program; if not, write to the Free
37 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
38 * Boston, MA 02110-1301, USA.
40 * \see http://www.gnu.org/copyleft/gpl.html
41 ***************************************************************/
43 #include "LinearQuadtreeExpansion.h"
44 #include "ComplexDouble.h"
48 using namespace ogdf::sse
;
52 LinearQuadtreeExpansion::LinearQuadtreeExpansion(__uint32 precision
, const LinearQuadtree
& tree
) : m_tree(tree
), m_numCoeff(precision
), binCoef(2*m_numCoeff
)
54 m_numExp
= m_tree
.maxNumberOfNodes();
59 LinearQuadtreeExpansion::~LinearQuadtreeExpansion(void)
65 void LinearQuadtreeExpansion::allocate()
67 m_multiExp
= (double*)MALLOC_16(m_numCoeff
*sizeof(double)*2*m_numExp
);
68 m_localExp
= (double*)MALLOC_16(m_numCoeff
*sizeof(double)*2*m_numExp
);
72 void LinearQuadtreeExpansion::deallocate()
79 void LinearQuadtreeExpansion::P2M(__uint32 point
, __uint32 receiver
)
81 double* receiv_coeff
= m_multiExp
+ receiver
*(m_numCoeff
<<1);
82 const double q
= (double)m_tree
.pointSize(point
);
83 const double x
= (double)m_tree
.pointX(point
);
84 const double y
= (double)m_tree
.pointY(point
);
85 const double centerX
= (double)m_tree
.nodeX(receiver
);
86 const double centerY
= (double)m_tree
.nodeY(receiver
);
92 ComplexDouble
delta(ComplexDouble(x
, y
) - ComplexDouble(centerX
, centerY
));
94 ComplexDouble
delta_k(delta
);
95 for (__uint32 k
=1; k
<m_numCoeff
; k
++)
97 ak
.load(receiv_coeff
+(k
<<1));
98 ak
-= delta_k
*(q
/(double)k
);
99 ak
.store(receiv_coeff
+(k
<<1));
105 void LinearQuadtreeExpansion::L2P(__uint32 source
, __uint32 point
, float& fx
, float& fy
)
107 const double* source_coeff
= m_localExp
+ source
*(m_numCoeff
<< 1);
108 const double x
= (double)m_tree
.pointX(point
);
109 const double y
= (double)m_tree
.pointY(point
);
110 const double centerX
= (double)m_tree
.nodeX(source
);
111 const double centerY
= (double)m_tree
.nodeY(source
);
114 ComplexDouble
delta(ComplexDouble(x
,y
) - ComplexDouble(centerX
, centerY
));
115 ComplexDouble
delta_k(1);
116 for (__uint32 k
=1; k
<m_numCoeff
; k
++)
118 ak
.load(source_coeff
+(k
<<1));
119 res
+= ak
*delta_k
*(double)k
;
124 res
.store_unaligned(resTemp
);
126 fx
-= ((float)resTemp
[0]);
127 fy
-= ((float)resTemp
[1]);
131 void LinearQuadtreeExpansion::M2M(__uint32 source
, __uint32 receiver
)
133 double* receiv_coeff
= m_multiExp
+ receiver
*(m_numCoeff
<<1);
134 double* source_coeff
= m_multiExp
+ source
*(m_numCoeff
<<1);
136 const double center_x_source
= (double)m_tree
.nodeX(source
);
137 const double center_y_source
= (double)m_tree
.nodeY(source
);
138 const double center_x_receiver
= (double)m_tree
.nodeX(receiver
);
139 const double center_y_receiver
= (double)m_tree
.nodeY(receiver
);
141 ComplexDouble
delta(ComplexDouble(center_x_source
, center_y_source
) - ComplexDouble(center_x_receiver
, center_y_receiver
));
143 ComplexDouble
a(source_coeff
);
144 ComplexDouble
b(receiv_coeff
);
146 b
.store(receiv_coeff
);
147 for (__uint32 l
=1;l
<m_numCoeff
;l
++)
149 b
.load(receiv_coeff
+(l
<<1));
150 ComplexDouble
delta_k(1.0,0.0);
151 for (__uint32 k
=0;k
<l
;k
++)
153 a
.load(source_coeff
+((l
-k
)<<1));
154 b
+= a
*delta_k
*binCoef
.value(l
-1, k
);
157 a
.load(source_coeff
);
158 b
-= a
*delta_k
*(1/(double)l
);
159 b
.store(receiv_coeff
+(l
<<1));
164 void LinearQuadtreeExpansion::L2L(__uint32 source
, __uint32 receiver
)
166 double* receiv_coeff
= m_localExp
+ receiver
*(m_numCoeff
<<1);
167 double* source_coeff
= m_localExp
+ source
*(m_numCoeff
<<1);
169 const double center_x_source
= (double)m_tree
.nodeX(source
);
170 const double center_y_source
= (double)m_tree
.nodeY(source
);
171 const double center_x_receiver
= (double)m_tree
.nodeX(receiver
);
172 const double center_y_receiver
= (double)m_tree
.nodeY(receiver
);
174 ComplexDouble
center_receiver(center_x_receiver
, center_y_receiver
);
175 ComplexDouble
center_source(center_x_source
, center_y_source
);
176 ComplexDouble
delta(center_source
- center_receiver
);
178 for (__uint32 l
=0;l
<m_numCoeff
;l
++)
180 ComplexDouble
b(receiv_coeff
+(l
<<1));
181 ComplexDouble
delta_k(1.0,0.0);
182 for (__uint32 k
=l
;k
<m_numCoeff
;k
++)
184 ComplexDouble
a(source_coeff
+(k
<<1));
185 b
+= a
*delta_k
*binCoef
.value(k
, l
);
188 b
.store(receiv_coeff
+(l
<<1));
193 void LinearQuadtreeExpansion::M2L(__uint32 source
, __uint32 receiver
)
195 double* receiv_coeff
= m_localExp
+ receiver
*(m_numCoeff
<<1);
196 double* source_coeff
= m_multiExp
+ source
*(m_numCoeff
<<1);
198 const float center_x_source
= (float)m_tree
.nodeX(source
);
199 const float center_y_source
= (float)m_tree
.nodeY(source
);
200 const float center_x_receiver
= (float)m_tree
.nodeX(receiver
);
201 const float center_y_receiver
= (float)m_tree
.nodeY(receiver
);
203 ComplexDouble
center_receiver(center_x_receiver
, center_y_receiver
);
204 ComplexDouble
center_source(center_x_source
, center_y_source
);
205 ComplexDouble
delta0(center_source
- center_receiver
);
207 ComplexDouble delta1
= -delta0
;
208 ComplexDouble
delta1_l(delta1
);
210 ComplexDouble
a0(source_coeff
);
213 for (__uint32 l
=1;l
<m_numCoeff
;l
++)
215 b
.load(receiv_coeff
+(l
<<1));
216 sum
= a0
*(-1/(double)l
);
217 ComplexDouble
delta0_k(delta0
);
218 for (__uint32 k
=1;k
<m_numCoeff
;k
++)
220 a
.load(source_coeff
+(k
<<1));
221 sum
+= (a
*binCoef
.value(l
+k
-1, k
-1)) / delta0_k
;
226 b
.store(receiv_coeff
+(l
<<1));
231 b
.load(receiv_coeff
);
232 //std::complex<double> sdelta(m_center[(receiver<<1)] - m_center[(source<<1)], m_center[(receiver<<1)+1] - m_center[(source<<1) +1]);//, m_x[receiver] - m_x[source], m_y[receiver] - m_y[source]);
233 /*std::complex<double> sdelta(center_x_receiver - center_x_source,
234 center_y_receiver - center_y_source);//, m_x[receiver] - m_x[source], m_y[receiver] - m_y[source]);
236 if ((sdelta.real() <=0) && (sdelta.imag() == 0)) //no cont. compl. log fct exists !!!
237 sdelta = log(sdelta + 0.00000001);
239 sdelta = log(sdelta);
242 double r
= delta1
.length();//= sqrt(sdelta.real()*sdelta.real()+sdelta.imag()*sdelta.imag());
243 double phi
= atan((center_x_receiver
- center_x_source
)/(center_y_receiver
- center_y_source
));
244 // sum = a0*log(z1 - z0)
245 b
+= a0
*ComplexDouble(log(r
), phi
);
247 //b += a0*ComplexDouble(sdelta.real(), sdelta.imag());
249 for (__uint32 k
=1;k
<m_numCoeff
;k
++)
251 a
.load(source_coeff
+(k
<<1));
256 b
.store(receiv_coeff
);
259 } // end of namespace ogdf