1 /*************************************************************************
2  *
3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4  *
5  * Copyright 2000, 2010 Oracle and/or its affiliates.
6  *
7  * OpenOffice.org - a multi-platform office productivity suite
8  *
9  * This file is part of OpenOffice.org.
10  *
11  * OpenOffice.org is free software: you can redistribute it and/or modify
12  * it under the terms of the GNU Lesser General Public License version 3
13  * only, as published by the Free Software Foundation.
14  *
15  * OpenOffice.org is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18  * GNU Lesser General Public License version 3 for more details
19  * (a copy is included in the LICENSE file that accompanied this code).
20  *
21  * You should have received a copy of the GNU Lesser General Public License
22  * version 3 along with OpenOffice.org.  If not, see
23  * <http://www.openoffice.org/license.html>
24  * for a copy of the LGPLv3 License.
25  *
26  ************************************************************************/
27 
28 #ifndef _HOMMATRIX_TEMPLATE_HXX
29 #define _HOMMATRIX_TEMPLATE_HXX
30 
31 #include <sal/types.h>
32 #include <basegfx/numeric/ftools.hxx>
33 #include <math.h>
34 #include <string.h>
35 
36 namespace basegfx
37 {
38     namespace internal
39     {
40 
41         inline double implGetDefaultValue(sal_uInt16 nRow, sal_uInt16 nColumn)
42         {
43             if(nRow == nColumn)
44                 return 1.0;
45             return 0.0;
46         }
47 
48         template < unsigned int _RowSize > class ImplMatLine
49         {
50             enum { RowSize = _RowSize };
51 
52             double											mfValue[RowSize];
53 
54         public:
55             ImplMatLine()
56             {
57             }
58 
59             ImplMatLine(sal_uInt16 nRow, ImplMatLine< RowSize >* pToBeCopied = 0L)
60             {
61                 if(pToBeCopied)
62                 {
63                     memcpy(&mfValue, pToBeCopied, sizeof(double) * RowSize);
64                 }
65                 else
66                 {
67                     for(sal_uInt16 a(0); a < RowSize; a++)
68                     {
69                         mfValue[a] = implGetDefaultValue(nRow, a);
70                     }
71                 }
72             }
73 
74             double get(sal_uInt16 nColumn) const
75             {
76                 return mfValue[nColumn];
77             }
78 
79             void set(sal_uInt16 nColumn, const double& rValue)
80             {
81                 mfValue[nColumn] = rValue;
82             }
83         };
84 
85         template < unsigned int _RowSize > class ImplHomMatrixTemplate
86         {
87             enum { RowSize = _RowSize };
88 
89             ImplMatLine< RowSize >							maLine[RowSize - 1];
90             ImplMatLine< RowSize >*							mpLine;
91 
92         public:
93             // Is last line used?
94             bool isLastLineDefault() const
95             {
96                 if(!mpLine)
97                     return true;
98 
99                 for(sal_uInt16 a(0); a < RowSize; a++)
100                 {
101                     const double fDefault(implGetDefaultValue((RowSize - 1), a));
102                     const double fLineValue(mpLine->get(a));
103 
104                     if(!::basegfx::fTools::equal(fDefault, fLineValue))
105                     {
106                         return false;
107                     }
108                 }
109 
110                 // reset last line, it equals default
111                 delete ((ImplHomMatrixTemplate< RowSize >*)this)->mpLine;
112                 ((ImplHomMatrixTemplate< RowSize >*)this)->mpLine = 0L;
113 
114                 return true;
115             }
116 
117             ImplHomMatrixTemplate()
118                 :	mpLine(0L)
119             {
120                 // complete initialization with identity matrix, all lines
121                 // were initialized with a trailing 1 followed by 0's.
122                 for(sal_uInt16 a(0); a < RowSize-1; a++)
123                 {
124                     for(sal_uInt16 b(0); b < RowSize; b++)
125                         maLine[a].set(b, implGetDefaultValue(a, b) );
126                 }
127             }
128 
129             ImplHomMatrixTemplate(const ImplHomMatrixTemplate& rToBeCopied)
130                 :	mpLine(0L)
131             {
132                 // complete initialization using copy
133                 for(sal_uInt16 a(0); a < (RowSize - 1); a++)
134                 {
135                     memcpy(&maLine[a], &rToBeCopied.maLine[a], sizeof(ImplMatLine< RowSize >));
136                 }
137 
138                 if(rToBeCopied.mpLine)
139                 {
140                     mpLine = new ImplMatLine< RowSize >((RowSize - 1), rToBeCopied.mpLine);
141                 }
142             }
143 
144             ~ImplHomMatrixTemplate()
145             {
146                 if(mpLine)
147                 {
148                     delete mpLine;
149                 }
150             }
151 
152             sal_uInt16 getEdgeLength() const { return RowSize; }
153 
154             double get(sal_uInt16 nRow, sal_uInt16 nColumn) const
155             {
156                 if(nRow < (RowSize - 1))
157                 {
158                     return maLine[nRow].get(nColumn);
159                 }
160 
161                 if(mpLine)
162                 {
163                     return mpLine->get(nColumn);
164                 }
165 
166                 return implGetDefaultValue((RowSize - 1), nColumn);
167             }
168 
169             void set(sal_uInt16 nRow, sal_uInt16 nColumn, const double& rValue)
170             {
171                 if(nRow < (RowSize - 1))
172                 {
173                     maLine[nRow].set(nColumn, rValue);
174                 }
175                 else if(mpLine)
176                 {
177                     mpLine->set(nColumn, rValue);
178                 }
179                 else
180                 {
181                     const double fDefault(implGetDefaultValue((RowSize - 1), nColumn));
182 
183                     if(!::basegfx::fTools::equal(fDefault, rValue))
184                     {
185                         mpLine = new ImplMatLine< RowSize >((RowSize - 1), 0L);
186                         mpLine->set(nColumn, rValue);
187                     }
188                 }
189             }
190 
191             void testLastLine()
192             {
193                 if(mpLine)
194                 {
195                     bool bNecessary(false);
196 
197                     for(sal_uInt16 a(0);!bNecessary && a < RowSize; a++)
198                     {
199                         const double fDefault(implGetDefaultValue((RowSize - 1), a));
200                         const double fLineValue(mpLine->get(a));
201 
202                         if(!::basegfx::fTools::equal(fDefault, fLineValue))
203                         {
204                             bNecessary = true;
205                         }
206                     }
207 
208                     if(!bNecessary)
209                     {
210                         delete mpLine;
211                         mpLine = 0L;
212                     }
213                 }
214             }
215 
216             // Left-upper decompositon
217             bool ludcmp(sal_uInt16 nIndex[], sal_Int16& nParity)
218             {
219                 double fBig, fSum, fDum;
220                 double fStorage[RowSize];
221                 sal_uInt16 a, b, c;
222 
223                 // #i30874# Initialize nAMax (compiler warns)
224                 sal_uInt16 nAMax = 0;
225 
226                 nParity = 1;
227 
228                 // Calc the max of each line. If a line is empty,
229                 // stop immediately since matrix is not invertible then.
230                 for(a = 0; a < RowSize; a++)
231                 {
232                     fBig = 0.0;
233 
234                     for(b = 0; b < RowSize; b++)
235                     {
236                         double fTemp(fabs(get(a, b)));
237 
238                         if(::basegfx::fTools::more(fTemp, fBig))
239                         {
240                             fBig = fTemp;
241                         }
242                     }
243 
244                     if(::basegfx::fTools::equalZero(fBig))
245                     {
246                         return false;
247                     }
248 
249                     fStorage[a] = 1.0 / fBig;
250                 }
251 
252                 // start normalizing
253                 for(b = 0; b < RowSize; b++)
254                 {
255                     for(a = 0; a < b; a++)
256                     {
257                         fSum = get(a, b);
258 
259                         for(c = 0; c < a; c++)
260                         {
261                             fSum -= get(a, c) * get(c, b);
262                         }
263 
264                         set(a, b, fSum);
265                     }
266 
267                     fBig = 0.0;
268 
269                     for(a = b; a < RowSize; a++)
270                     {
271                         fSum = get(a, b);
272 
273                         for(c = 0; c < b; c++)
274                         {
275                             fSum -= get(a, c) * get(c, b);
276                         }
277 
278                         set(a, b, fSum);
279                         fDum = fStorage[a] * fabs(fSum);
280 
281                         if(::basegfx::fTools::moreOrEqual(fDum, fBig))
282                         {
283                             fBig = fDum;
284                             nAMax = a;
285                         }
286                     }
287 
288                     if(b != nAMax)
289                     {
290                         for(c = 0; c < RowSize; c++)
291                         {
292                             fDum = get(nAMax, c);
293                             set(nAMax, c, get(b, c));
294                             set(b, c, fDum);
295                         }
296 
297                         nParity = -nParity;
298                         fStorage[nAMax] = fStorage[b];
299                     }
300 
301                     nIndex[b] = nAMax;
302 
303                     // here the failure of precision occurs
304                     const double fValBB(fabs(get(b, b)));
305 
306                     if(::basegfx::fTools::equalZero(fValBB))
307                     {
308                         return false;
309                     }
310 
311                     if(b != (RowSize - 1))
312                     {
313                         fDum = 1.0 / get(b, b);
314 
315                         for(a = b + 1; a < RowSize; a++)
316                         {
317                             set(a, b, get(a, b) * fDum);
318                         }
319                     }
320                 }
321 
322                 return true;
323             }
324 
325             void lubksb(const sal_uInt16 nIndex[], double fRow[]) const
326             {
327                 sal_uInt16 b, ip;
328                 sal_Int16 a, a2 = -1;
329                 double fSum;
330 
331                 for(a = 0; a < RowSize; a++)
332                 {
333                     ip = nIndex[a];
334                     fSum = fRow[ip];
335                     fRow[ip] = fRow[a];
336 
337                     if(a2 >= 0)
338                     {
339                         for(b = a2; b < a; b++)
340                         {
341                             fSum -= get(a, b) * fRow[b];
342                         }
343                     }
344                     else if(!::basegfx::fTools::equalZero(fSum))
345                     {
346                         a2 = a;
347                     }
348 
349                     fRow[a] = fSum;
350                 }
351 
352                 for(a = (RowSize - 1); a >= 0; a--)
353                 {
354                     fSum = fRow[a];
355 
356                     for(b = a + 1; b < RowSize; b++)
357                     {
358                         fSum -= get(a, b) * fRow[b];
359                     }
360 
361                     const double fValueAA(get(a, a));
362 
363                     if(!::basegfx::fTools::equalZero(fValueAA))
364                     {
365                         fRow[a] = fSum / get(a, a);
366                     }
367                 }
368             }
369 
370             bool isIdentity() const
371             {
372                 // last line needs no testing if not existing
373                 const sal_uInt16 nMaxLine(
374                     sal::static_int_cast<sal_uInt16>(mpLine ? RowSize : (RowSize - 1)) );
375 
376                 for(sal_uInt16 a(0); a < nMaxLine; a++)
377                 {
378                     for(sal_uInt16 b(0); b < RowSize; b++)
379                     {
380                         const double fDefault(implGetDefaultValue(a, b));
381                         const double fValueAB(get(a, b));
382 
383                         if(!::basegfx::fTools::equal(fDefault, fValueAB))
384                         {
385                             return false;
386                         }
387                     }
388                 }
389 
390                 return true;
391             }
392 
393             bool isInvertible() const
394             {
395                 ImplHomMatrixTemplate aWork(*this);
396                 sal_uInt16 nIndex[RowSize];
397                 sal_Int16 nParity;
398 
399                 return aWork.ludcmp(nIndex, nParity);
400             }
401 
402             bool isNormalized() const
403             {
404                 if(!mpLine)
405                     return true;
406 
407                 const double fHomValue(get((RowSize - 1), (RowSize - 1)));
408 
409                 if(::basegfx::fTools::equalZero(fHomValue))
410                 {
411                     return true;
412                 }
413 
414                 const double fOne(1.0);
415 
416                 if(::basegfx::fTools::equal(fOne, fHomValue))
417                 {
418                     return true;
419                 }
420 
421                 return false;
422             }
423 
424             void doInvert(const ImplHomMatrixTemplate& rWork, const sal_uInt16 nIndex[])
425             {
426                 double fArray[RowSize];
427 
428                 for(sal_uInt16 a(0); a < RowSize; a++)
429                 {
430                     // prepare line
431 		    sal_uInt16 b;
432                     for( b = 0; b < RowSize; b++)
433                     {
434                         fArray[b] = implGetDefaultValue(a, b);
435                     }
436 
437                     // expand line
438                     rWork.lubksb(nIndex, fArray);
439 
440                     // copy line transposed to this matrix
441                     for( b = 0; b < RowSize; b++)
442                     {
443                         set(b, a, fArray[b]);
444                     }
445                 }
446 
447                 // evtl. get rid of last matrix line
448                 testLastLine();
449             }
450 
451             void doNormalize()
452             {
453                 if(mpLine)
454                 {
455                     const double fHomValue(get((RowSize - 1), (RowSize - 1)));
456 
457                     for(sal_uInt16 a(0); a < RowSize; a++)
458                     {
459                         for(sal_uInt16 b(0); b < RowSize; b++)
460                         {
461                             set(a, b, get(a, b) / fHomValue);
462                         }
463                     }
464 
465                     // evtl. get rid of last matrix line
466                     testLastLine();
467                 }
468             }
469 
470             double doDeterminant() const
471             {
472                 ImplHomMatrixTemplate aWork(*this);
473                 sal_uInt16 nIndex[RowSize];
474                 sal_Int16 nParity;
475                 double fRetval(0.0);
476 
477                 if(aWork.ludcmp(nIndex, nParity))
478                 {
479                     fRetval = (double)nParity;
480 
481                     // last line needs no multiply if not existing; default value would be 1.
482                     const sal_uInt16 nMaxLine(
483                         sal::static_int_cast<sal_uInt16>(aWork.mpLine ? RowSize : (RowSize - 1)) );
484 
485                     for(sal_uInt16 a(0); a < nMaxLine; a++)
486                     {
487                         fRetval *= aWork.get(a, a);
488                     }
489                 }
490 
491                 return fRetval;
492             }
493 
494             double doTrace() const
495             {
496                 double fTrace = (mpLine) ? 0.0 : 1.0;
497                 const sal_uInt16 nMaxLine(
498                     sal::static_int_cast<sal_uInt16>(mpLine ? RowSize : (RowSize - 1)) );
499 
500                 for(sal_uInt16 a(0); a < nMaxLine; a++)
501                 {
502                     fTrace += get(a, a);
503                 }
504 
505                 return fTrace;
506             }
507 
508             void doTranspose()
509             {
510                 for(sal_uInt16 a(0); a < (RowSize - 1); a++)
511                 {
512                     for(sal_uInt16 b(a + 1); b < RowSize; b++)
513                     {
514                         const double fTemp(get(a, b));
515                         set(a, b, get(b, a));
516                         set(b, a, fTemp);
517                     }
518                 }
519 
520                 testLastLine();
521             }
522 
523             void doAddMatrix(const ImplHomMatrixTemplate& rMat)
524             {
525                 for(sal_uInt16 a(0); a < RowSize; a++)
526                 {
527                     for(sal_uInt16 b(0); b < RowSize; b++)
528                     {
529                         set(a, b, get(a, b) + rMat.get(a, b));
530                     }
531                 }
532 
533                 testLastLine();
534             }
535 
536             void doSubMatrix(const ImplHomMatrixTemplate& rMat)
537             {
538                 for(sal_uInt16 a(0); a < RowSize; a++)
539                 {
540                     for(sal_uInt16 b(0); b < RowSize; b++)
541                     {
542                         set(a, b, get(a, b) - rMat.get(a, b));
543                     }
544                 }
545 
546                 testLastLine();
547             }
548 
549             void doMulMatrix(const double& rfValue)
550             {
551                 for(sal_uInt16 a(0); a < RowSize; a++)
552                 {
553                     for(sal_uInt16 b(0); b < RowSize; b++)
554                     {
555                         set(a, b, get(a, b) * rfValue);
556                     }
557                 }
558 
559                 testLastLine();
560             }
561 
562             void doMulMatrix(const ImplHomMatrixTemplate& rMat)
563             {
564                 // create a copy as source for the original values
565                 const ImplHomMatrixTemplate aCopy(*this);
566 
567                 // TODO: maybe optimize cases where last line is [0 0 1].
568 
569                 double fValue(0.0);
570 
571                 for(sal_uInt16 a(0); a < RowSize; ++a)
572                 {
573                     for(sal_uInt16 b(0); b < RowSize; ++b)
574                     {
575                         fValue = 0.0;
576 
577                         for(sal_uInt16 c(0); c < RowSize; ++c)
578                             fValue += aCopy.get(c, b) * rMat.get(a, c);
579 
580                         set(a, b, fValue);
581                     }
582                 }
583 
584                 testLastLine();
585             }
586 
587             bool isEqual(const ImplHomMatrixTemplate& rMat) const
588             {
589                 const sal_uInt16 nMaxLine(
590                     sal::static_int_cast<sal_uInt16>((mpLine || rMat.mpLine) ? RowSize : (RowSize - 1)) );
591 
592                 for(sal_uInt16 a(0); a < nMaxLine; a++)
593                 {
594                     for(sal_uInt16 b(0); b < RowSize; b++)
595                     {
596                         const double fValueA(get(a, b));
597                         const double fValueB(rMat.get(a, b));
598 
599                         if(!::basegfx::fTools::equal(fValueA, fValueB))
600                         {
601                             return false;
602                         }
603                     }
604                 }
605 
606                 return true;
607             }
608         };
609 
610     } // namespace internal
611 } // namespace basegfx
612 
613 #endif /* _HOMMATRIX_TEMPLATE_HXX */
614