1/*
2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3%                                                                             %
4%                                                                             %
5%                                                                             %
6%                  M   M   AAA   TTTTT  RRRR   IIIII  X   X                   %
7%                  MM MM  A   A    T    R   R    I     X X                    %
8%                  M M M  AAAAA    T    RRRR     I      X                     %
9%                  M   M  A   A    T    R R      I     X X                    %
10%                  M   M  A   A    T    R  R   IIIII  X   X                   %
11%                                                                             %
12%                                                                             %
13%                         MagickCore Matrix Methods                           %
14%                                                                             %
15%                            Software Design                                  %
16%                                 Cristy                                      %
17%                              August 2007                                    %
18%                                                                             %
19%                                                                             %
20%  Copyright 1999-2016 ImageMagick Studio LLC, a non-profit organization      %
21%  dedicated to making software imaging solutions freely available.           %
22%                                                                             %
23%  You may not use this file except in compliance with the License.  You may  %
24%  obtain a copy of the License at                                            %
25%                                                                             %
26%    http://www.imagemagick.org/script/license.php                            %
27%                                                                             %
28%  Unless required by applicable law or agreed to in writing, software        %
29%  distributed under the License is distributed on an "AS IS" BASIS,          %
30%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
31%  See the License for the specific language governing permissions and        %
32%  limitations under the License.                                             %
33%                                                                             %
34%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35%
36%
37*/
38
39/*
40  Include declarations.
41*/
42#include "MagickCore/studio.h"
43#include "MagickCore/blob.h"
44#include "MagickCore/blob-private.h"
45#include "MagickCore/cache.h"
46#include "MagickCore/exception.h"
47#include "MagickCore/exception-private.h"
48#include "MagickCore/image-private.h"
49#include "MagickCore/matrix.h"
50#include "MagickCore/memory_.h"
51#include "MagickCore/pixel-accessor.h"
52#include "MagickCore/pixel-private.h"
53#include "MagickCore/resource_.h"
54#include "MagickCore/semaphore.h"
55#include "MagickCore/thread-private.h"
56#include "MagickCore/utility.h"
57
58/*
59  Typedef declaration.
60*/
61struct _MatrixInfo
62{
63  CacheType
64    type;
65
66  size_t
67    columns,
68    rows,
69    stride;
70
71  MagickSizeType
72    length;
73
74  MagickBooleanType
75    mapped,
76    synchronize;
77
78  char
79    path[MagickPathExtent];
80
81  int
82    file;
83
84  void
85    *elements;
86
87  SemaphoreInfo
88    *semaphore;
89
90  size_t
91    signature;
92};
93
94/*
95%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96%                                                                             %
97%                                                                             %
98%                                                                             %
99%   A c q u i r e M a t r i x I n f o                                         %
100%                                                                             %
101%                                                                             %
102%                                                                             %
103%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
104%
105%  AcquireMatrixInfo() allocates the ImageInfo structure.
106%
107%  The format of the AcquireMatrixInfo method is:
108%
109%      MatrixInfo *AcquireMatrixInfo(const size_t columns,const size_t rows,
110%        const size_t stride,ExceptionInfo *exception)
111%
112%  A description of each parameter follows:
113%
114%    o columns: the matrix columns.
115%
116%    o rows: the matrix rows.
117%
118%    o stride: the matrix stride.
119%
120%    o exception: return any errors or warnings in this structure.
121%
122*/
123
124#if defined(SIGBUS)
125static void MatrixSignalHandler(int status)
126{
127  ThrowFatalException(CacheFatalError,"UnableToExtendMatrixCache");
128}
129#endif
130
131static inline MagickOffsetType WriteMatrixElements(
132  const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset,
133  const MagickSizeType length,const unsigned char *magick_restrict buffer)
134{
135  register MagickOffsetType
136    i;
137
138  ssize_t
139    count;
140
141#if !defined(MAGICKCORE_HAVE_PWRITE)
142  LockSemaphoreInfo(matrix_info->semaphore);
143  if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
144    {
145      UnlockSemaphoreInfo(matrix_info->semaphore);
146      return((MagickOffsetType) -1);
147    }
148#endif
149  count=0;
150  for (i=0; i < (MagickOffsetType) length; i+=count)
151  {
152#if !defined(MAGICKCORE_HAVE_PWRITE)
153    count=write(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
154      (MagickSizeType) SSIZE_MAX));
155#else
156    count=pwrite(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
157      (MagickSizeType) SSIZE_MAX),(off_t) (offset+i));
158#endif
159    if (count <= 0)
160      {
161        count=0;
162        if (errno != EINTR)
163          break;
164      }
165  }
166#if !defined(MAGICKCORE_HAVE_PWRITE)
167  UnlockSemaphoreInfo(matrix_info->semaphore);
168#endif
169  return(i);
170}
171
172static MagickBooleanType SetMatrixExtent(
173  MatrixInfo *magick_restrict matrix_info,
174  MagickSizeType length)
175{
176  MagickOffsetType
177    count,
178    extent,
179    offset;
180
181  if (length != (MagickSizeType) ((MagickOffsetType) length))
182    return(MagickFalse);
183  offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END);
184  if (offset < 0)
185    return(MagickFalse);
186  if ((MagickSizeType) offset >= length)
187    return(MagickTrue);
188  extent=(MagickOffsetType) length-1;
189  count=WriteMatrixElements(matrix_info,extent,1,(const unsigned char *) "");
190#if defined(MAGICKCORE_HAVE_POSIX_FALLOCATE)
191  if (matrix_info->synchronize != MagickFalse)
192    (void) posix_fallocate(matrix_info->file,offset+1,extent-offset);
193#endif
194#if defined(SIGBUS)
195  (void) signal(SIGBUS,MatrixSignalHandler);
196#endif
197  return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue);
198}
199
200MagickExport MatrixInfo *AcquireMatrixInfo(const size_t columns,
201  const size_t rows,const size_t stride,ExceptionInfo *exception)
202{
203  char
204    *synchronize;
205
206  MagickBooleanType
207    status;
208
209  MatrixInfo
210    *matrix_info;
211
212  matrix_info=(MatrixInfo *) AcquireMagickMemory(sizeof(*matrix_info));
213  if (matrix_info == (MatrixInfo *) NULL)
214    return((MatrixInfo *) NULL);
215  (void) ResetMagickMemory(matrix_info,0,sizeof(*matrix_info));
216  matrix_info->signature=MagickCoreSignature;
217  matrix_info->columns=columns;
218  matrix_info->rows=rows;
219  matrix_info->stride=stride;
220  matrix_info->semaphore=AcquireSemaphoreInfo();
221  synchronize=GetEnvironmentValue("MAGICK_SYNCHRONIZE");
222  if (synchronize != (const char *) NULL)
223    {
224      matrix_info->synchronize=IsStringTrue(synchronize);
225      synchronize=DestroyString(synchronize);
226    }
227  matrix_info->length=(MagickSizeType) columns*rows*stride;
228  if (matrix_info->columns != (size_t) (matrix_info->length/rows/stride))
229    {
230      (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
231        "CacheResourcesExhausted","`%s'","matrix cache");
232      return(DestroyMatrixInfo(matrix_info));
233    }
234  matrix_info->type=MemoryCache;
235  status=AcquireMagickResource(AreaResource,matrix_info->length);
236  if ((status != MagickFalse) &&
237      (matrix_info->length == (MagickSizeType) ((size_t) matrix_info->length)))
238    {
239      status=AcquireMagickResource(MemoryResource,matrix_info->length);
240      if (status != MagickFalse)
241        {
242          matrix_info->mapped=MagickFalse;
243          matrix_info->elements=AcquireMagickMemory((size_t)
244            matrix_info->length);
245          if (matrix_info->elements == NULL)
246            {
247              matrix_info->mapped=MagickTrue;
248              matrix_info->elements=MapBlob(-1,IOMode,0,(size_t)
249                matrix_info->length);
250            }
251          if (matrix_info->elements == (unsigned short *) NULL)
252            RelinquishMagickResource(MemoryResource,matrix_info->length);
253        }
254    }
255  matrix_info->file=(-1);
256  if (matrix_info->elements == (unsigned short *) NULL)
257    {
258      status=AcquireMagickResource(DiskResource,matrix_info->length);
259      if (status == MagickFalse)
260        {
261          (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
262            "CacheResourcesExhausted","`%s'","matrix cache");
263          return(DestroyMatrixInfo(matrix_info));
264        }
265      matrix_info->type=DiskCache;
266      (void) AcquireMagickResource(MemoryResource,matrix_info->length);
267      matrix_info->file=AcquireUniqueFileResource(matrix_info->path);
268      if (matrix_info->file == -1)
269        return(DestroyMatrixInfo(matrix_info));
270      status=AcquireMagickResource(MapResource,matrix_info->length);
271      if (status != MagickFalse)
272        {
273          status=SetMatrixExtent(matrix_info,matrix_info->length);
274          if (status != MagickFalse)
275            {
276              matrix_info->elements=(void *) MapBlob(matrix_info->file,IOMode,0,
277                (size_t) matrix_info->length);
278              if (matrix_info->elements != NULL)
279                matrix_info->type=MapCache;
280              else
281                RelinquishMagickResource(MapResource,matrix_info->length);
282            }
283        }
284    }
285  return(matrix_info);
286}
287
288/*
289%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
290%                                                                             %
291%                                                                             %
292%                                                                             %
293%   A c q u i r e M a g i c k M a t r i x                                     %
294%                                                                             %
295%                                                                             %
296%                                                                             %
297%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
298%
299%  AcquireMagickMatrix() allocates and returns a matrix in the form of an
300%  array of pointers to an array of doubles, with all values pre-set to zero.
301%
302%  This used to generate the two dimensional matrix, and vectors required
303%  for the GaussJordanElimination() method below, solving some system of
304%  simultanious equations.
305%
306%  The format of the AcquireMagickMatrix method is:
307%
308%      double **AcquireMagickMatrix(const size_t number_rows,
309%        const size_t size)
310%
311%  A description of each parameter follows:
312%
313%    o number_rows: the number pointers for the array of pointers
314%      (first dimension).
315%
316%    o size: the size of the array of doubles each pointer points to
317%      (second dimension).
318%
319*/
320MagickExport double **AcquireMagickMatrix(const size_t number_rows,
321  const size_t size)
322{
323  double
324    **matrix;
325
326  register ssize_t
327    i,
328    j;
329
330  matrix=(double **) AcquireQuantumMemory(number_rows,sizeof(*matrix));
331  if (matrix == (double **) NULL)
332    return((double **) NULL);
333  for (i=0; i < (ssize_t) number_rows; i++)
334  {
335    matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i]));
336    if (matrix[i] == (double *) NULL)
337    {
338      for (j=0; j < i; j++)
339        matrix[j]=(double *) RelinquishMagickMemory(matrix[j]);
340      matrix=(double **) RelinquishMagickMemory(matrix);
341      return((double **) NULL);
342    }
343    for (j=0; j < (ssize_t) size; j++)
344      matrix[i][j]=0.0;
345  }
346  return(matrix);
347}
348
349/*
350%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
351%                                                                             %
352%                                                                             %
353%                                                                             %
354%   D e s t r o y M a t r i x I n f o                                         %
355%                                                                             %
356%                                                                             %
357%                                                                             %
358%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
359%
360%  DestroyMatrixInfo() dereferences a matrix, deallocating memory associated
361%  with the matrix.
362%
363%  The format of the DestroyImage method is:
364%
365%      MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info)
366%
367%  A description of each parameter follows:
368%
369%    o matrix_info: the matrix.
370%
371*/
372MagickExport MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info)
373{
374  assert(matrix_info != (MatrixInfo *) NULL);
375  assert(matrix_info->signature == MagickCoreSignature);
376  LockSemaphoreInfo(matrix_info->semaphore);
377  switch (matrix_info->type)
378  {
379    case MemoryCache:
380    {
381      if (matrix_info->mapped == MagickFalse)
382        matrix_info->elements=RelinquishMagickMemory(matrix_info->elements);
383      else
384        {
385          (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length);
386          matrix_info->elements=(unsigned short *) NULL;
387        }
388      RelinquishMagickResource(MemoryResource,matrix_info->length);
389      break;
390    }
391    case MapCache:
392    {
393      (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length);
394      matrix_info->elements=NULL;
395      RelinquishMagickResource(MapResource,matrix_info->length);
396    }
397    case DiskCache:
398    {
399      if (matrix_info->file != -1)
400        (void) close(matrix_info->file);
401      (void) RelinquishUniqueFileResource(matrix_info->path);
402      RelinquishMagickResource(DiskResource,matrix_info->length);
403      break;
404    }
405    default:
406      break;
407  }
408  UnlockSemaphoreInfo(matrix_info->semaphore);
409  RelinquishSemaphoreInfo(&matrix_info->semaphore);
410  return((MatrixInfo *) RelinquishMagickMemory(matrix_info));
411}
412
413/*
414%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
415%                                                                             %
416%                                                                             %
417%                                                                             %
418+   G a u s s J o r d a n E l i m i n a t i o n                               %
419%                                                                             %
420%                                                                             %
421%                                                                             %
422%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
423%
424%  GaussJordanElimination() returns a matrix in reduced row echelon form,
425%  while simultaneously reducing and thus solving the augumented results
426%  matrix.
427%
428%  See also  http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
429%
430%  The format of the GaussJordanElimination method is:
431%
432%      MagickBooleanType GaussJordanElimination(double **matrix,
433%        double **vectors,const size_t rank,const size_t number_vectors)
434%
435%  A description of each parameter follows:
436%
437%    o matrix: the matrix to be reduced, as an 'array of row pointers'.
438%
439%    o vectors: the additional matrix argumenting the matrix for row reduction.
440%             Producing an 'array of column vectors'.
441%
442%    o rank:  The size of the matrix (both rows and columns).
443%             Also represents the number terms that need to be solved.
444%
445%    o number_vectors: Number of vectors columns, argumenting the above matrix.
446%             Usally 1, but can be more for more complex equation solving.
447%
448%  Note that the 'matrix' is given as a 'array of row pointers' of rank size.
449%  That is values can be assigned as   matrix[row][column]   where 'row' is
450%  typically the equation, and 'column' is the term of the equation.
451%  That is the matrix is in the form of a 'row first array'.
452%
453%  However 'vectors' is a 'array of column pointers' which can have any number
454%  of columns, with each column array the same 'rank' size as 'matrix'.
455%
456%  This allows for simpler handling of the results, especially is only one
457%  column 'vector' is all that is required to produce the desired solution.
458%
459%  For example, the 'vectors' can consist of a pointer to a simple array of
460%  doubles.  when only one set of simultanious equations is to be solved from
461%  the given set of coefficient weighted terms.
462%
463%     double **matrix = AcquireMagickMatrix(8UL,8UL);
464%     double coefficents[8];
465%     ...
466%     GaussJordanElimination(matrix, &coefficents, 8UL, 1UL);
467%
468%  However by specifing more 'columns' (as an 'array of vector columns',
469%  you can use this function to solve a set of 'separable' equations.
470%
471%  For example a distortion function where    u = U(x,y)   v = V(x,y)
472%  And the functions U() and V() have separate coefficents, but are being
473%  generated from a common x,y->u,v  data set.
474%
475%  Another example is generation of a color gradient from a set of colors at
476%  specific coordients, such as a list x,y -> r,g,b,a.
477%
478%  You can also use the 'vectors' to generate an inverse of the given 'matrix'
479%  though as a 'column first array' rather than a 'row first array'. For
480%  details see http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
481%
482*/
483MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix,
484  double **vectors,const size_t rank,const size_t number_vectors)
485{
486#define GaussJordanSwap(x,y) \
487{ \
488  if ((x) != (y)) \
489    { \
490      (x)+=(y); \
491      (y)=(x)-(y); \
492      (x)=(x)-(y); \
493    } \
494}
495
496  double
497    max,
498    scale;
499
500  register ssize_t
501    i,
502    j,
503    k;
504
505  ssize_t
506    column,
507    *columns,
508    *pivots,
509    row,
510    *rows;
511
512  columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns));
513  rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows));
514  pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots));
515  if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) ||
516      (pivots == (ssize_t *) NULL))
517    {
518      if (pivots != (ssize_t *) NULL)
519        pivots=(ssize_t *) RelinquishMagickMemory(pivots);
520      if (columns != (ssize_t *) NULL)
521        columns=(ssize_t *) RelinquishMagickMemory(columns);
522      if (rows != (ssize_t *) NULL)
523        rows=(ssize_t *) RelinquishMagickMemory(rows);
524      return(MagickFalse);
525    }
526  (void) ResetMagickMemory(columns,0,rank*sizeof(*columns));
527  (void) ResetMagickMemory(rows,0,rank*sizeof(*rows));
528  (void) ResetMagickMemory(pivots,0,rank*sizeof(*pivots));
529  column=0;
530  row=0;
531  for (i=0; i < (ssize_t) rank; i++)
532  {
533    max=0.0;
534    for (j=0; j < (ssize_t) rank; j++)
535      if (pivots[j] != 1)
536        {
537          for (k=0; k < (ssize_t) rank; k++)
538            if (pivots[k] != 0)
539              {
540                if (pivots[k] > 1)
541                  return(MagickFalse);
542              }
543            else
544              if (fabs(matrix[j][k]) >= max)
545                {
546                  max=fabs(matrix[j][k]);
547                  row=j;
548                  column=k;
549                }
550        }
551    pivots[column]++;
552    if (row != column)
553      {
554        for (k=0; k < (ssize_t) rank; k++)
555          GaussJordanSwap(matrix[row][k],matrix[column][k]);
556        for (k=0; k < (ssize_t) number_vectors; k++)
557          GaussJordanSwap(vectors[k][row],vectors[k][column]);
558      }
559    rows[i]=row;
560    columns[i]=column;
561    if (matrix[column][column] == 0.0)
562      return(MagickFalse);  /* sigularity */
563    scale=PerceptibleReciprocal(matrix[column][column]);
564    matrix[column][column]=1.0;
565    for (j=0; j < (ssize_t) rank; j++)
566      matrix[column][j]*=scale;
567    for (j=0; j < (ssize_t) number_vectors; j++)
568      vectors[j][column]*=scale;
569    for (j=0; j < (ssize_t) rank; j++)
570      if (j != column)
571        {
572          scale=matrix[j][column];
573          matrix[j][column]=0.0;
574          for (k=0; k < (ssize_t) rank; k++)
575            matrix[j][k]-=scale*matrix[column][k];
576          for (k=0; k < (ssize_t) number_vectors; k++)
577            vectors[k][j]-=scale*vectors[k][column];
578        }
579  }
580  for (j=(ssize_t) rank-1; j >= 0; j--)
581    if (columns[j] != rows[j])
582      for (i=0; i < (ssize_t) rank; i++)
583        GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
584  pivots=(ssize_t *) RelinquishMagickMemory(pivots);
585  rows=(ssize_t *) RelinquishMagickMemory(rows);
586  columns=(ssize_t *) RelinquishMagickMemory(columns);
587  return(MagickTrue);
588}
589
590/*
591%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
592%                                                                             %
593%                                                                             %
594%                                                                             %
595%   G e t M a t r i x C o l u m n s                                           %
596%                                                                             %
597%                                                                             %
598%                                                                             %
599%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
600%
601%  GetMatrixColumns() returns the number of columns in the matrix.
602%
603%  The format of the GetMatrixColumns method is:
604%
605%      size_t GetMatrixColumns(const MatrixInfo *matrix_info)
606%
607%  A description of each parameter follows:
608%
609%    o matrix_info: the matrix.
610%
611*/
612MagickExport size_t GetMatrixColumns(const MatrixInfo *matrix_info)
613{
614  assert(matrix_info != (MatrixInfo *) NULL);
615  assert(matrix_info->signature == MagickCoreSignature);
616  return(matrix_info->columns);
617}
618
619/*
620%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
621%                                                                             %
622%                                                                             %
623%                                                                             %
624%   G e t M a t r i x E l e m e n t                                           %
625%                                                                             %
626%                                                                             %
627%                                                                             %
628%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
629%
630%  GetMatrixElement() returns the specifed element in the matrix.
631%
632%  The format of the GetMatrixElement method is:
633%
634%      MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
635%        const ssize_t x,const ssize_t y,void *value)
636%
637%  A description of each parameter follows:
638%
639%    o matrix_info: the matrix columns.
640%
641%    o x: the matrix x-offset.
642%
643%    o y: the matrix y-offset.
644%
645%    o value: return the matrix element in this buffer.
646%
647*/
648
649static inline ssize_t EdgeX(const ssize_t x,const size_t columns)
650{
651  if (x < 0L)
652    return(0L);
653  if (x >= (ssize_t) columns)
654    return((ssize_t) (columns-1));
655  return(x);
656}
657
658static inline ssize_t EdgeY(const ssize_t y,const size_t rows)
659{
660  if (y < 0L)
661    return(0L);
662  if (y >= (ssize_t) rows)
663    return((ssize_t) (rows-1));
664  return(y);
665}
666
667static inline MagickOffsetType ReadMatrixElements(
668  const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset,
669  const MagickSizeType length,unsigned char *magick_restrict buffer)
670{
671  register MagickOffsetType
672    i;
673
674  ssize_t
675    count;
676
677#if !defined(MAGICKCORE_HAVE_PREAD)
678  LockSemaphoreInfo(matrix_info->semaphore);
679  if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
680    {
681      UnlockSemaphoreInfo(matrix_info->semaphore);
682      return((MagickOffsetType) -1);
683    }
684#endif
685  count=0;
686  for (i=0; i < (MagickOffsetType) length; i+=count)
687  {
688#if !defined(MAGICKCORE_HAVE_PREAD)
689    count=read(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
690      (MagickSizeType) SSIZE_MAX));
691#else
692    count=pread(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
693      (MagickSizeType) SSIZE_MAX),(off_t) (offset+i));
694#endif
695    if (count <= 0)
696      {
697        count=0;
698        if (errno != EINTR)
699          break;
700      }
701  }
702#if !defined(MAGICKCORE_HAVE_PREAD)
703  UnlockSemaphoreInfo(matrix_info->semaphore);
704#endif
705  return(i);
706}
707
708MagickExport MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
709  const ssize_t x,const ssize_t y,void *value)
710{
711  MagickOffsetType
712    count,
713    i;
714
715  assert(matrix_info != (const MatrixInfo *) NULL);
716  assert(matrix_info->signature == MagickCoreSignature);
717  i=(MagickOffsetType) EdgeY(y,matrix_info->rows)*matrix_info->columns+
718    EdgeX(x,matrix_info->columns);
719  if (matrix_info->type != DiskCache)
720    {
721      (void) memcpy(value,(unsigned char *) matrix_info->elements+i*
722        matrix_info->stride,matrix_info->stride);
723      return(MagickTrue);
724    }
725  count=ReadMatrixElements(matrix_info,i*matrix_info->stride,
726    matrix_info->stride,(unsigned char *) value);
727  if (count != (MagickOffsetType) matrix_info->stride)
728    return(MagickFalse);
729  return(MagickTrue);
730}
731
732/*
733%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
734%                                                                             %
735%                                                                             %
736%                                                                             %
737%   G e t M a t r i x R o w s                                                 %
738%                                                                             %
739%                                                                             %
740%                                                                             %
741%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
742%
743%  GetMatrixRows() returns the number of rows in the matrix.
744%
745%  The format of the GetMatrixRows method is:
746%
747%      size_t GetMatrixRows(const MatrixInfo *matrix_info)
748%
749%  A description of each parameter follows:
750%
751%    o matrix_info: the matrix.
752%
753*/
754MagickExport size_t GetMatrixRows(const MatrixInfo *matrix_info)
755{
756  assert(matrix_info != (const MatrixInfo *) NULL);
757  assert(matrix_info->signature == MagickCoreSignature);
758  return(matrix_info->rows);
759}
760
761/*
762%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
763%                                                                             %
764%                                                                             %
765%                                                                             %
766+   L e a s t S q u a r e s A d d T e r m s                                   %
767%                                                                             %
768%                                                                             %
769%                                                                             %
770%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
771%
772%  LeastSquaresAddTerms() adds one set of terms and associate results to the
773%  given matrix and vectors for solving using least-squares function fitting.
774%
775%  The format of the AcquireMagickMatrix method is:
776%
777%      void LeastSquaresAddTerms(double **matrix,double **vectors,
778%        const double *terms,const double *results,const size_t rank,
779%        const size_t number_vectors);
780%
781%  A description of each parameter follows:
782%
783%    o matrix: the square matrix to add given terms/results to.
784%
785%    o vectors: the result vectors to add terms/results to.
786%
787%    o terms: the pre-calculated terms (without the unknown coefficent
788%             weights) that forms the equation being added.
789%
790%    o results: the result(s) that should be generated from the given terms
791%               weighted by the yet-to-be-solved coefficents.
792%
793%    o rank: the rank or size of the dimensions of the square matrix.
794%            Also the length of vectors, and number of terms being added.
795%
796%    o number_vectors: Number of result vectors, and number or results being
797%      added.  Also represents the number of separable systems of equations
798%      that is being solved.
799%
800%  Example of use...
801%
802%     2 dimensional Affine Equations (which are separable)
803%         c0*x + c2*y + c4*1 => u
804%         c1*x + c3*y + c5*1 => v
805%
806%     double **matrix = AcquireMagickMatrix(3UL,3UL);
807%     double **vectors = AcquireMagickMatrix(2UL,3UL);
808%     double terms[3], results[2];
809%     ...
810%     for each given x,y -> u,v
811%        terms[0] = x;
812%        terms[1] = y;
813%        terms[2] = 1;
814%        results[0] = u;
815%        results[1] = v;
816%        LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL);
817%     ...
818%     if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) {
819%       c0 = vectors[0][0];
820%       c2 = vectors[0][1];
821%       c4 = vectors[0][2];
822%       c1 = vectors[1][0];
823%       c3 = vectors[1][1];
824%       c5 = vectors[1][2];
825%     }
826%     else
827%       printf("Matrix unsolvable\n);
828%     RelinquishMagickMatrix(matrix,3UL);
829%     RelinquishMagickMatrix(vectors,2UL);
830%
831*/
832MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors,
833  const double *terms,const double *results,const size_t rank,
834  const size_t number_vectors)
835{
836  register ssize_t
837    i,
838    j;
839
840  for (j=0; j < (ssize_t) rank; j++)
841  {
842    for (i=0; i < (ssize_t) rank; i++)
843      matrix[i][j]+=terms[i]*terms[j];
844    for (i=0; i < (ssize_t) number_vectors; i++)
845      vectors[i][j]+=results[i]*terms[j];
846  }
847}
848
849/*
850%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
851%                                                                             %
852%                                                                             %
853%                                                                             %
854%   M a t r i x T o I m a g e                                                 %
855%                                                                             %
856%                                                                             %
857%                                                                             %
858%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
859%
860%  MatrixToImage() returns a matrix as an image.  The matrix elements must be
861%  of type double otherwise nonsense is returned.
862%
863%  The format of the MatrixToImage method is:
864%
865%      Image *MatrixToImage(const MatrixInfo *matrix_info,
866%        ExceptionInfo *exception)
867%
868%  A description of each parameter follows:
869%
870%    o matrix_info: the matrix.
871%
872%    o exception: return any errors or warnings in this structure.
873%
874*/
875MagickExport Image *MatrixToImage(const MatrixInfo *matrix_info,
876  ExceptionInfo *exception)
877{
878  CacheView
879    *image_view;
880
881  double
882    max_value,
883    min_value,
884    scale_factor,
885    value;
886
887  Image
888    *image;
889
890  MagickBooleanType
891    status;
892
893  ssize_t
894    y;
895
896  assert(matrix_info != (const MatrixInfo *) NULL);
897  assert(matrix_info->signature == MagickCoreSignature);
898  assert(exception != (ExceptionInfo *) NULL);
899  assert(exception->signature == MagickCoreSignature);
900  if (matrix_info->stride < sizeof(double))
901    return((Image *) NULL);
902  /*
903    Determine range of matrix.
904  */
905  (void) GetMatrixElement(matrix_info,0,0,&value);
906  min_value=value;
907  max_value=value;
908  for (y=0; y < (ssize_t) matrix_info->rows; y++)
909  {
910    register ssize_t
911      x;
912
913    for (x=0; x < (ssize_t) matrix_info->columns; x++)
914    {
915      if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
916        continue;
917      if (value < min_value)
918        min_value=value;
919      else
920        if (value > max_value)
921          max_value=value;
922    }
923  }
924  if ((min_value == 0.0) && (max_value == 0.0))
925    scale_factor=0;
926  else
927    if (min_value == max_value)
928      {
929        scale_factor=(double) QuantumRange/min_value;
930        min_value=0;
931      }
932    else
933      scale_factor=(double) QuantumRange/(max_value-min_value);
934  /*
935    Convert matrix to image.
936  */
937  image=AcquireImage((ImageInfo *) NULL,exception);
938  image->columns=matrix_info->columns;
939  image->rows=matrix_info->rows;
940  image->colorspace=GRAYColorspace;
941  status=MagickTrue;
942  image_view=AcquireAuthenticCacheView(image,exception);
943#if defined(MAGICKCORE_OPENMP_SUPPORT)
944  #pragma omp parallel for schedule(static,4) shared(status) \
945    magick_threads(image,image,image->rows,1)
946#endif
947  for (y=0; y < (ssize_t) image->rows; y++)
948  {
949    double
950      value;
951
952    register Quantum
953      *q;
954
955    register ssize_t
956      x;
957
958    if (status == MagickFalse)
959      continue;
960    q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
961    if (q == (Quantum *) NULL)
962      {
963        status=MagickFalse;
964        continue;
965      }
966    for (x=0; x < (ssize_t) image->columns; x++)
967    {
968      if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
969        continue;
970      value=scale_factor*(value-min_value);
971      *q=ClampToQuantum(value);
972      q+=GetPixelChannels(image);
973    }
974    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
975      status=MagickFalse;
976  }
977  image_view=DestroyCacheView(image_view);
978  if (status == MagickFalse)
979    image=DestroyImage(image);
980  return(image);
981}
982
983/*
984%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
985%                                                                             %
986%                                                                             %
987%                                                                             %
988%   N u l l M a t r i x                                                       %
989%                                                                             %
990%                                                                             %
991%                                                                             %
992%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
993%
994%  NullMatrix() sets all elements of the matrix to zero.
995%
996%  The format of the ResetMagickMemory method is:
997%
998%      MagickBooleanType *NullMatrix(MatrixInfo *matrix_info)
999%
1000%  A description of each parameter follows:
1001%
1002%    o matrix_info: the matrix.
1003%
1004*/
1005MagickExport MagickBooleanType NullMatrix(MatrixInfo *matrix_info)
1006{
1007  register ssize_t
1008    x;
1009
1010  ssize_t
1011    count,
1012    y;
1013
1014  unsigned char
1015    value;
1016
1017  assert(matrix_info != (const MatrixInfo *) NULL);
1018  assert(matrix_info->signature == MagickCoreSignature);
1019  if (matrix_info->type != DiskCache)
1020    {
1021      (void) ResetMagickMemory(matrix_info->elements,0,(size_t)
1022        matrix_info->length);
1023      return(MagickTrue);
1024    }
1025  value=0;
1026  (void) lseek(matrix_info->file,0,SEEK_SET);
1027  for (y=0; y < (ssize_t) matrix_info->rows; y++)
1028  {
1029    for (x=0; x < (ssize_t) matrix_info->length; x++)
1030    {
1031      count=write(matrix_info->file,&value,sizeof(value));
1032      if (count != (ssize_t) sizeof(value))
1033        break;
1034    }
1035    if (x < (ssize_t) matrix_info->length)
1036      break;
1037  }
1038  return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
1039}
1040
1041/*
1042%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1043%                                                                             %
1044%                                                                             %
1045%                                                                             %
1046%   R e l i n q u i s h M a g i c k M a t r i x                               %
1047%                                                                             %
1048%                                                                             %
1049%                                                                             %
1050%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1051%
1052%  RelinquishMagickMatrix() frees the previously acquired matrix (array of
1053%  pointers to arrays of doubles).
1054%
1055%  The format of the RelinquishMagickMatrix method is:
1056%
1057%      double **RelinquishMagickMatrix(double **matrix,
1058%        const size_t number_rows)
1059%
1060%  A description of each parameter follows:
1061%
1062%    o matrix: the matrix to relinquish
1063%
1064%    o number_rows: the first dimension of the acquired matrix (number of
1065%      pointers)
1066%
1067*/
1068MagickExport double **RelinquishMagickMatrix(double **matrix,
1069  const size_t number_rows)
1070{
1071  register ssize_t
1072    i;
1073
1074  if (matrix == (double **) NULL )
1075    return(matrix);
1076  for (i=0; i < (ssize_t) number_rows; i++)
1077     matrix[i]=(double *) RelinquishMagickMemory(matrix[i]);
1078  matrix=(double **) RelinquishMagickMemory(matrix);
1079  return(matrix);
1080}
1081
1082/*
1083%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1084%                                                                             %
1085%                                                                             %
1086%                                                                             %
1087%   S e t M a t r i x E l e m e n t                                           %
1088%                                                                             %
1089%                                                                             %
1090%                                                                             %
1091%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1092%
1093%  SetMatrixElement() sets the specifed element in the matrix.
1094%
1095%  The format of the SetMatrixElement method is:
1096%
1097%      MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
1098%        const ssize_t x,const ssize_t y,void *value)
1099%
1100%  A description of each parameter follows:
1101%
1102%    o matrix_info: the matrix columns.
1103%
1104%    o x: the matrix x-offset.
1105%
1106%    o y: the matrix y-offset.
1107%
1108%    o value: set the matrix element to this value.
1109%
1110*/
1111
1112MagickExport MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
1113  const ssize_t x,const ssize_t y,const void *value)
1114{
1115  MagickOffsetType
1116    count,
1117    i;
1118
1119  assert(matrix_info != (const MatrixInfo *) NULL);
1120  assert(matrix_info->signature == MagickCoreSignature);
1121  i=(MagickOffsetType) y*matrix_info->columns+x;
1122  if ((i < 0) ||
1123      ((MagickSizeType) (i*matrix_info->stride) >= matrix_info->length))
1124    return(MagickFalse);
1125  if (matrix_info->type != DiskCache)
1126    {
1127      (void) memcpy((unsigned char *) matrix_info->elements+i*
1128        matrix_info->stride,value,matrix_info->stride);
1129      return(MagickTrue);
1130    }
1131  count=WriteMatrixElements(matrix_info,i*matrix_info->stride,
1132    matrix_info->stride,(unsigned char *) value);
1133  if (count != (MagickOffsetType) matrix_info->stride)
1134    return(MagickFalse);
1135  return(MagickTrue);
1136}
1137